diff options
Diffstat (limited to 'libquadmath/math/j1q.c')
| -rw-r--r-- | libquadmath/math/j1q.c | 103 |
1 files changed, 66 insertions, 37 deletions
diff --git a/libquadmath/math/j1q.c b/libquadmath/math/j1q.c index eb599c949a9..5eb705084e2 100644 --- a/libquadmath/math/j1q.c +++ b/libquadmath/math/j1q.c @@ -95,6 +95,7 @@ License along with this library; if not, write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ +#include <errno.h> #include "quadmath-imp.h" /* 1 / sqrt(pi) */ @@ -687,13 +688,21 @@ j1q (__float128 x) if (! finiteq (x)) { if (x != x) - return x; + return x + x; else return 0.0Q; } if (x == 0.0Q) return x; xx = fabsq (x); + if (xx <= 0x1p-58Q) + { + __float128 ret = x * 0.5Q; + math_check_force_underflow (ret); + if (ret == 0) + errno = ERANGE; + return ret; + } if (xx <= 2.0Q) { /* 0 <= x <= 2 */ @@ -705,6 +714,32 @@ j1q (__float128 x) return p; } + /* X = x - 3 pi/4 + cos(X) = cos(x) cos(3 pi/4) + sin(x) sin(3 pi/4) + = 1/sqrt(2) * (-cos(x) + sin(x)) + sin(X) = sin(x) cos(3 pi/4) - cos(x) sin(3 pi/4) + = -1/sqrt(2) * (sin(x) + cos(x)) + cf. Fdlibm. */ + sincosq (xx, &s, &c); + ss = -s - c; + cc = s - c; + if (xx <= FLT128_MAX / 2.0Q) + { + z = cosq (xx + xx); + if ((s * c) > 0) + cc = z / ss; + else + ss = z / cc; + } + + if (xx > 0x1p256Q) + { + z = ONEOSQPI * cc / sqrtq (xx); + if (x < 0) + z = -z; + return z; + } + xinv = 1.0Q / xx; z = xinv * xinv; if (xinv <= 0.25) @@ -766,20 +801,6 @@ j1q (__float128 x) p = 1.0Q + z * p; q = z * q; q = q * xinv + 0.375Q * xinv; - /* X = x - 3 pi/4 - cos(X) = cos(x) cos(3 pi/4) + sin(x) sin(3 pi/4) - = 1/sqrt(2) * (-cos(x) + sin(x)) - sin(X) = sin(x) cos(3 pi/4) - cos(x) sin(3 pi/4) - = -1/sqrt(2) * (sin(x) + cos(x)) - cf. Fdlibm. */ - sincosq (xx, &s, &c); - ss = -s - c; - cc = s - c; - z = cosq (xx + xx); - if ((s * c) > 0) - cc = z / ss; - else - ss = z / cc; z = ONEOSQPI * (p * cc - q * ss) / sqrtq (xx); if (x < 0) z = -z; @@ -823,24 +844,25 @@ y1q (__float128 x) __float128 xx, xinv, z, p, q, c, s, cc, ss; if (! finiteq (x)) - { - if (x != x) - return x; - else - return 0.0Q; - } + return 1 / (x + x * x); if (x <= 0.0Q) { if (x < 0.0Q) return (zero / (zero * x)); - return -HUGE_VALQ + x; + return -1 / zero; /* -inf and divide by zero exception. */ } xx = fabsq (x); if (xx <= 0x1p-114) - return -TWOOPI / x; + { + z = -TWOOPI / x; + if (isinfq (z)) + errno = ERANGE; + return z; + } if (xx <= 2.0Q) { /* 0 <= x <= 2 */ + /* FIXME: SET_RESTORE_ROUNDL (FE_TONEAREST); */ z = xx * xx; p = xx * neval (z, Y0_2N, NY0_2N) / deval (z, Y0_2D, NY0_2D); p = -TWOOPI / xx + p; @@ -848,6 +870,27 @@ y1q (__float128 x) return p; } + /* X = x - 3 pi/4 + cos(X) = cos(x) cos(3 pi/4) + sin(x) sin(3 pi/4) + = 1/sqrt(2) * (-cos(x) + sin(x)) + sin(X) = sin(x) cos(3 pi/4) - cos(x) sin(3 pi/4) + = -1/sqrt(2) * (sin(x) + cos(x)) + cf. Fdlibm. */ + sincosq (xx, &s, &c); + ss = -s - c; + cc = s - c; + if (xx <= FLT128_MAX / 2.0Q) + { + z = cosq (xx + xx); + if ((s * c) > 0) + cc = z / ss; + else + ss = z / cc; + } + + if (xx > 0x1p256Q) + return ONEOSQPI * ss / sqrtq (xx); + xinv = 1.0Q / xx; z = xinv * xinv; if (xinv <= 0.25) @@ -909,20 +952,6 @@ y1q (__float128 x) p = 1.0Q + z * p; q = z * q; q = q * xinv + 0.375Q * xinv; - /* X = x - 3 pi/4 - cos(X) = cos(x) cos(3 pi/4) + sin(x) sin(3 pi/4) - = 1/sqrt(2) * (-cos(x) + sin(x)) - sin(X) = sin(x) cos(3 pi/4) - cos(x) sin(3 pi/4) - = -1/sqrt(2) * (sin(x) + cos(x)) - cf. Fdlibm. */ - sincosq (xx, &s, &c); - ss = -s - c; - cc = s - c; - z = cosq (xx + xx); - if ((s * c) > 0) - cc = z / ss; - else - ss = z / cc; z = ONEOSQPI * (p * ss + q * cc) / sqrtq (xx); return z; } |