diff options
Diffstat (limited to 'libquadmath/math/j1q.c')
| -rw-r--r-- | libquadmath/math/j1q.c | 50 |
1 files changed, 25 insertions, 25 deletions
diff --git a/libquadmath/math/j1q.c b/libquadmath/math/j1q.c index 5eb705084e2..b18e881276b 100644 --- a/libquadmath/math/j1q.c +++ b/libquadmath/math/j1q.c @@ -6,9 +6,9 @@ * * SYNOPSIS: * - * __float128 x, y, j1q(); + * long double x, y, j1l(); * - * y = j1q( x ); + * y = j1l( x ); * * * @@ -52,9 +52,9 @@ * * SYNOPSIS: * - * __float128, y, y1q(); + * double x, y, y1l(); * - * y = y1q( x ); + * y = y1l( x ); * * * @@ -92,17 +92,16 @@ Lesser General Public License for more details. You should have received a copy of the GNU Lesser General Public - License along with this library; if not, write to the Free Software - Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ + License along with this library; if not, see + <http://www.gnu.org/licenses/>. */ -#include <errno.h> #include "quadmath-imp.h" /* 1 / sqrt(pi) */ static const __float128 ONEOSQPI = 5.6418958354775628694807945156077258584405E-1Q; /* 2 / pi */ static const __float128 TWOOPI = 6.3661977236758134307553505349005744813784E-1Q; -static const __float128 zero = 0.0Q; +static const __float128 zero = 0; /* J1(x) = .5x + x x^2 R(x^2) Peak relative error 1.9e-35 @@ -126,7 +125,7 @@ static const __float128 J0_2D[NJ0_2D + 1] = { 5.673775894803172808323058205986256928794E8Q, 1.080329960080981204840966206372671147224E6Q, 1.411951256636576283942477881535283304912E3Q, - /* 1.000000000000000000000000000000000000000E0Q */ + /* 1.000000000000000000000000000000000000000E0L */ }; /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2), @@ -690,9 +689,9 @@ j1q (__float128 x) if (x != x) return x + x; else - return 0.0Q; + return 0; } - if (x == 0.0Q) + if (x == 0) return x; xx = fabsq (x); if (xx <= 0x1p-58Q) @@ -703,7 +702,7 @@ j1q (__float128 x) errno = ERANGE; return ret; } - if (xx <= 2.0Q) + if (xx <= 2) { /* 0 <= x <= 2 */ z = xx * xx; @@ -723,7 +722,7 @@ j1q (__float128 x) sincosq (xx, &s, &c); ss = -s - c; cc = s - c; - if (xx <= FLT128_MAX / 2.0Q) + if (xx <= FLT128_MAX / 2) { z = cosq (xx + xx); if ((s * c) > 0) @@ -740,7 +739,7 @@ j1q (__float128 x) return z; } - xinv = 1.0Q / xx; + xinv = 1 / xx; z = xinv * xinv; if (xinv <= 0.25) { @@ -798,7 +797,7 @@ j1q (__float128 x) q = neval (z, Q2_2r3N, NQ2_2r3N) / deval (z, Q2_2r3D, NQ2_2r3D); } } - p = 1.0Q + z * p; + p = 1 + z * p; q = z * q; q = q * xinv + 0.375Q * xinv; z = ONEOSQPI * (p * cc - q * ss) / sqrtq (xx); @@ -808,11 +807,12 @@ j1q (__float128 x) } + /* Y1(x) = 2/pi * (log(x) * J1(x) - 1/x) + x R(x^2) Peak relative error 6.2e-38 0 <= x <= 2 */ #define NY0_2N 7 -static __float128 Y0_2N[NY0_2N + 1] = { +static const __float128 Y0_2N[NY0_2N + 1] = { -6.804415404830253804408698161694720833249E19Q, 1.805450517967019908027153056150465849237E19Q, -8.065747497063694098810419456383006737312E17Q, @@ -823,7 +823,7 @@ static __float128 Y0_2N[NY0_2N + 1] = { 9.541172044989995856117187515882879304461E5Q, }; #define NY0_2D 7 -static __float128 Y0_2D[NY0_2D + 1] = { +static const __float128 Y0_2D[NY0_2D + 1] = { 3.470629591820267059538637461549677594549E20Q, 4.120796439009916326855848107545425217219E18Q, 2.477653371652018249749350657387030814542E16Q, @@ -845,9 +845,9 @@ y1q (__float128 x) if (! finiteq (x)) return 1 / (x + x * x); - if (x <= 0.0Q) + if (x <= 0) { - if (x < 0.0Q) + if (x < 0) return (zero / (zero * x)); return -1 / zero; /* -inf and divide by zero exception. */ } @@ -859,10 +859,10 @@ y1q (__float128 x) errno = ERANGE; return z; } - if (xx <= 2.0Q) - { + if (xx <= 2) + { /* 0 <= x <= 2 */ - /* FIXME: SET_RESTORE_ROUNDL (FE_TONEAREST); */ + SET_RESTORE_ROUNDF128 (FE_TONEAREST); z = xx * xx; p = xx * neval (z, Y0_2N, NY0_2N) / deval (z, Y0_2D, NY0_2D); p = -TWOOPI / xx + p; @@ -879,7 +879,7 @@ y1q (__float128 x) sincosq (xx, &s, &c); ss = -s - c; cc = s - c; - if (xx <= FLT128_MAX / 2.0Q) + if (xx <= FLT128_MAX / 2) { z = cosq (xx + xx); if ((s * c) > 0) @@ -891,7 +891,7 @@ y1q (__float128 x) if (xx > 0x1p256Q) return ONEOSQPI * ss / sqrtq (xx); - xinv = 1.0Q / xx; + xinv = 1 / xx; z = xinv * xinv; if (xinv <= 0.25) { @@ -949,7 +949,7 @@ y1q (__float128 x) q = neval (z, Q2_2r3N, NQ2_2r3N) / deval (z, Q2_2r3D, NQ2_2r3D); } } - p = 1.0Q + z * p; + p = 1 + z * p; q = z * q; q = q * xinv + 0.375Q * xinv; z = ONEOSQPI * (p * ss + q * cc) / sqrtq (xx); |