From ea9bb95a5b36c0a3d2ed8fb03808745b406c2633 Mon Sep 17 00:00:00 2001 From: Szabolcs Nagy Date: Tue, 3 Sep 2013 18:50:58 +0000 Subject: [PATCH] math: long double trigonometric cleanup (cosl, sinl, sincosl, tanl) ld128 support was added to internal kernel functions (__cosl, __sinl, __tanl, __rem_pio2l) from freebsd (not tested, but should be a good start for when ld128 arch arrives) __rem_pio2l had some code cleanup, the freebsd ld128 code seems to gather the results of a large reduction with precision loss (fixed the bug but a todo comment was added for later investigation) the old copyright was removed from the non-kernel wrapper functions (cosl, sinl, sincosl, tanl) since these are trivial and the interesting parts and comments had been already rewritten. --- src/math/__cosl.c | 37 +++++++++++-- src/math/__rem_pio2l.c | 117 ++++++++++++++++++++++------------------- src/math/__sinl.c | 36 +++++++++++-- src/math/__tanl.c | 66 ++++++++++++++++++++--- src/math/cosl.c | 70 +++++------------------- src/math/sincosl.c | 19 +++---- src/math/sinl.c | 66 +++++------------------ src/math/tanl.c | 53 +++---------------- 8 files changed, 228 insertions(+), 236 deletions(-) diff --git a/src/math/__cosl.c b/src/math/__cosl.c index 9d325768..fa522ddd 100644 --- a/src/math/__cosl.c +++ b/src/math/__cosl.c @@ -1,4 +1,5 @@ /* origin: FreeBSD /usr/src/lib/msun/ld80/k_cosl.c */ +/* origin: FreeBSD /usr/src/lib/msun/ld128/k_cosl.c */ /* * ==================================================== * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. @@ -14,7 +15,8 @@ #include "libm.h" -#if LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 +#if (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 +#if LDBL_MANT_DIG == 64 /* * ld80 version of __cos.c. See __cos.c for most comments. */ @@ -43,7 +45,6 @@ */ static const long double C1 = 0.0416666666666666666136L; /* 0xaaaaaaaaaaaaaa9b.0p-68 */ - static const double C2 = -0.0013888888888888874, /* -0x16c16c16c16c10.0p-62 */ C3 = 0.000024801587301571716, /* 0x1a01a01a018e22.0p-68 */ @@ -51,13 +52,43 @@ C4 = -0.00000027557319215507120, /* -0x127e4fb7602f22.0p-74 */ C5 = 0.0000000020876754400407278, /* 0x11eed8caaeccf1.0p-81 */ C6 = -1.1470297442401303e-11, /* -0x19393412bd1529.0p-89 */ C7 = 4.7383039476436467e-14; /* 0x1aac9d9af5c43e.0p-97 */ +#define POLY(z) (z*(C1+z*(C2+z*(C3+z*(C4+z*(C5+z*(C6+z*C7))))))) +#elif LDBL_MANT_DIG == 113 +/* + * ld128 version of __cos.c. See __cos.c for most comments. + */ +/* + * Domain [-0.7854, 0.7854], range ~[-1.80e-37, 1.79e-37]: + * |cos(x) - c(x))| < 2**-122.0 + * + * 113-bit precision requires more care than 64-bit precision, since + * simple methods give a minimax polynomial with coefficient for x^2 + * that is 1 ulp below 0.5, but we want it to be precisely 0.5. See + * above for more details. + */ +static const long double +C1 = 0.04166666666666666666666666666666658424671L, +C2 = -0.001388888888888888888888888888863490893732L, +C3 = 0.00002480158730158730158730158600795304914210L, +C4 = -0.2755731922398589065255474947078934284324e-6L, +C5 = 0.2087675698786809897659225313136400793948e-8L, +C6 = -0.1147074559772972315817149986812031204775e-10L, +C7 = 0.4779477332386808976875457937252120293400e-13L; +static const double +C8 = -0.1561920696721507929516718307820958119868e-15, +C9 = 0.4110317413744594971475941557607804508039e-18, +C10 = -0.8896592467191938803288521958313920156409e-21, +C11 = 0.1601061435794535138244346256065192782581e-23; +#define POLY(z) (z*(C1+z*(C2+z*(C3+z*(C4+z*(C5+z*(C6+z*(C7+ \ + z*(C8+z*(C9+z*(C10+z*C11))))))))))) +#endif long double __cosl(long double x, long double y) { long double hz,z,r,w; z = x*x; - r = z*(C1+z*(C2+z*(C3+z*(C4+z*(C5+z*(C6+z*C7)))))); + r = POLY(z); hz = 0.5*z; w = 1.0-hz; return w + (((1.0-w)-hz) + (z*r-x*y)); diff --git a/src/math/__rem_pio2l.c b/src/math/__rem_pio2l.c index f0cb99ac..8b15b7b2 100644 --- a/src/math/__rem_pio2l.c +++ b/src/math/__rem_pio2l.c @@ -13,15 +13,22 @@ * Optimized by Bruce D. Evans. */ #include "libm.h" -#if LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 -/* ld80 version of __rem_pio2(x,y) +#if (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 +/* ld80 and ld128 version of __rem_pio2(x,y) * * return the remainder of x rem pi/2 in y[0]+y[1] * use __rem_pio2_large() for large x */ -#define BIAS (LDBL_MAX_EXP - 1) - +#if LDBL_MANT_DIG == 64 +/* u ~< 0x1p25*pi/2 */ +#define SMALL(u) (((u.i.se & 0x7fffU)<<16 | u.i.m>>48) < ((0x3fff + 25)<<16 | 0x921f>>1 | 0x8000)) +#define TOINT 0x1.8p63 +#define QUOBITS(x) ((uint32_t)(int32_t)x & 0x7fffffff) +#define ROUND1 22 +#define ROUND2 61 +#define NX 3 +#define NY 2 /* * invpio2: 64 bits of 2/pi * pio2_1: first 39 bits of pi/2 @@ -32,60 +39,61 @@ * pio2_3t: pi/2 - (pio2_1+pio2_2+pio2_3) */ static const double -two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */ pio2_1 = 1.57079632679597125389e+00, /* 0x3FF921FB, 0x54444000 */ pio2_2 = -1.07463465549783099519e-12, /* -0x12e7b967674000.0p-92 */ pio2_3 = 6.36831716351370313614e-25; /* 0x18a2e037074000.0p-133 */ - static const long double invpio2 = 6.36619772367581343076e-01L, /* 0xa2f9836e4e44152a.0p-64 */ pio2_1t = -1.07463465549719416346e-12L, /* -0x973dcb3b399d747f.0p-103 */ pio2_2t = 6.36831716351095013979e-25L, /* 0xc51701b839a25205.0p-144 */ pio2_3t = -2.75299651904407171810e-37L; /* -0xbb5bf6c7ddd660ce.0p-185 */ +#elif LDBL_MANT_DIG == 113 +/* u ~< 0x1p45*pi/2 */ +#define SMALL(u) (((u.i.se & 0x7fffU)<<16 | u.i.top) < ((0x3fff + 45)<<16 | 0x921f)) +#define TOINT 0x1.8p112 +#define QUOBITS(x) ((uint32_t)(int64_t)x & 0x7fffffff) +#define ROUND1 51 +#define ROUND2 119 +#define NX 5 +#define NY 3 +static const long double +invpio2 = 6.3661977236758134307553505349005747e-01L, /* 0x145f306dc9c882a53f84eafa3ea6a.0p-113 */ +pio2_1 = 1.5707963267948966192292994253909555e+00L, /* 0x1921fb54442d18469800000000000.0p-112 */ +pio2_1t = 2.0222662487959507323996846200947577e-21L, /* 0x13198a2e03707344a4093822299f3.0p-181 */ +pio2_2 = 2.0222662487959507323994779168837751e-21L, /* 0x13198a2e03707344a400000000000.0p-181 */ +pio2_2t = 2.0670321098263988236496903051604844e-43L, /* 0x127044533e63a0105df531d89cd91.0p-254 */ +pio2_3 = 2.0670321098263988236499468110329591e-43L, /* 0x127044533e63a0105e00000000000.0p-254 */ +pio2_3t = -2.5650587247459238361625433492959285e-65L; /* -0x159c4ec64ddaeb5f78671cbfb2210.0p-327 */ +#endif int __rem_pio2l(long double x, long double *y) { - union IEEEl2bits u,u1; + union ldshape u,uz; long double z,w,t,r,fn; - double tx[3],ty[2]; - int e0,ex,i,j,nx,n; - int16_t expsign; - - u.e = x; - expsign = u.xbits.expsign; - ex = expsign & 0x7fff; - if (ex < BIAS + 25 || (ex == BIAS + 25 && u.bits.manh < 0xc90fdaa2)) { - union IEEEl2bits u2; - int ex1; + double tx[NX],ty[NY]; + int ex,ey,n,i; - /* |x| ~< 2^25*(pi/2), medium size */ - /* Use a specialized rint() to get fn. Assume round-to-nearest. */ - fn = x*invpio2 + 0x1.8p63; - fn = fn - 0x1.8p63; -// FIXME -//#ifdef HAVE_EFFICIENT_IRINT -// n = irint(fn); -//#else - n = fn; -//#endif + u.f = x; + ex = u.i.se & 0x7fff; + if (SMALL(u)) { + /* rint(x/(pi/2)), Assume round-to-nearest. */ + fn = x*invpio2 + TOINT - TOINT; + n = QUOBITS(fn); r = x-fn*pio2_1; - w = fn*pio2_1t; /* 1st round good to 102 bit */ - j = ex; + w = fn*pio2_1t; /* 1st round good to 102/180 bits (ld80/ld128) */ y[0] = r-w; - u2.e = y[0]; - ex1 = u2.xbits.expsign & 0x7fff; - i = j-ex1; - if (i > 22) { /* 2nd iteration needed, good to 141 */ + u.f = y[0]; + ey = u.i.se & 0x7fff; + if (ex - ey > ROUND1) { /* 2nd iteration needed, good to 141/248 (ld80/ld128) */ t = r; w = fn*pio2_2; r = t-w; w = fn*pio2_2t-((t-r)-w); y[0] = r-w; - u2.e = y[0]; - ex1 = u2.xbits.expsign & 0x7fff; - i = j-ex1; - if (i > 61) { /* 3rd iteration need, 180 bits acc */ - t = r; /* will cover all possible cases */ + u.f = y[0]; + ey = u.i.se & 0x7fff; + if (ex - ey > ROUND2) { /* 3rd iteration, good to 180/316 bits */ + t = r; /* will cover all possible cases (not verified for ld128) */ w = fn*pio2_3; r = t-w; w = fn*pio2_3t-((t-r)-w); @@ -102,23 +110,26 @@ int __rem_pio2l(long double x, long double *y) y[0] = y[1] = x - x; return 0; } - /* set z = scalbn(|x|,ilogb(x)-23) */ - u1.e = x; - e0 = ex - BIAS - 23; /* e0 = ilogb(|x|)-23; */ - u1.xbits.expsign = ex - e0; - z = u1.e; - for (i=0; i<2; i++) { + /* set z = scalbn(|x|,-ilogb(x)+23) */ + uz.f = x; + uz.i.se = 0x3fff + 23; + z = uz.f; + for (i=0; i < NX - 1; i++) { tx[i] = (double)(int32_t)z; - z = (z-tx[i])*two24; + z = (z-tx[i])*0x1p24; } - tx[2] = z; - nx = 3; - while (tx[nx-1] == 0.0) - nx--; /* skip zero term */ - n = __rem_pio2_large(tx,ty,e0,nx,2); - r = (long double)ty[0] + ty[1]; - w = ty[1] - (r - ty[0]); - if (expsign < 0) { + tx[i] = z; + while (tx[i] == 0) + i--; + n = __rem_pio2_large(tx, ty, ex-0x3fff-23, i+1, NY); + w = ty[1]; + if (NY == 3) + w += ty[2]; + r = ty[0] + w; + /* TODO: for ld128 this does not follow the recommendation of the + comments of __rem_pio2_large which seem wrong if |ty[0]| > |ty[1]+ty[2]| */ + w -= r - ty[0]; + if (u.i.se >> 15) { y[0] = -r; y[1] = -w; return -n; diff --git a/src/math/__sinl.c b/src/math/__sinl.c index 068adffb..2525bbe8 100644 --- a/src/math/__sinl.c +++ b/src/math/__sinl.c @@ -1,4 +1,5 @@ /* origin: FreeBSD /usr/src/lib/msun/ld80/k_sinl.c */ +/* origin: FreeBSD /usr/src/lib/msun/ld128/k_sinl.c */ /* * ==================================================== * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. @@ -13,7 +14,8 @@ #include "libm.h" -#if LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 +#if (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 +#if LDBL_MANT_DIG == 64 /* * ld80 version of __sin.c. See __sin.c for most comments. */ @@ -23,10 +25,8 @@ * * See __cosl.c for more details about the polynomial. */ - static const long double S1 = -0.166666666666666666671L; /* -0xaaaaaaaaaaaaaaab.0p-66 */ - static const double S2 = 0.0083333333333333332, /* 0x11111111111111.0p-59 */ S3 = -0.00019841269841269427, /* -0x1a01a01a019f81.0p-65 */ @@ -35,6 +35,34 @@ S5 = -0.000000025052108218074604, /* -0x1ae64564f16cad.0p-78 */ S6 = 1.6059006598854211e-10, /* 0x161242b90243b5.0p-85 */ S7 = -7.6429779983024564e-13, /* -0x1ae42ebd1b2e00.0p-93 */ S8 = 2.6174587166648325e-15; /* 0x179372ea0b3f64.0p-101 */ +#define POLY(z) (S2+z*(S3+z*(S4+z*(S5+z*(S6+z*(S7+z*S8)))))) +#elif LDBL_MANT_DIG == 113 +/* + * ld128 version of __sin.c. See __sin.c for most comments. + */ +/* + * Domain [-0.7854, 0.7854], range ~[-1.53e-37, 1.659e-37] + * |sin(x)/x - s(x)| < 2**-122.1 + * + * See __cosl.c for more details about the polynomial. + */ +static const long double +S1 = -0.16666666666666666666666666666666666606732416116558L, +S2 = 0.0083333333333333333333333333333331135404851288270047L, +S3 = -0.00019841269841269841269841269839935785325638310428717L, +S4 = 0.27557319223985890652557316053039946268333231205686e-5L, +S5 = -0.25052108385441718775048214826384312253862930064745e-7L, +S6 = 0.16059043836821614596571832194524392581082444805729e-9L, +S7 = -0.76471637318198151807063387954939213287488216303768e-12L, +S8 = 0.28114572543451292625024967174638477283187397621303e-14L; +static const double +S9 = -0.82206352458348947812512122163446202498005154296863e-17, +S10 = 0.19572940011906109418080609928334380560135358385256e-19, +S11 = -0.38680813379701966970673724299207480965452616911420e-22, +S12 = 0.64038150078671872796678569586315881020659912139412e-25; +#define POLY(z) (S2+z*(S3+z*(S4+z*(S5+z*(S6+z*(S7+z*(S8+ \ + z*(S9+z*(S10+z*(S11+z*S12)))))))))) +#endif long double __sinl(long double x, long double y, int iy) { @@ -42,7 +70,7 @@ long double __sinl(long double x, long double y, int iy) z = x*x; v = z*x; - r = S2+z*(S3+z*(S4+z*(S5+z*(S6+z*(S7+z*S8))))); + r = POLY(z); if (iy == 0) return x+v*(S1+z*r); return x-((z*(0.5*y-v*r)-y)-v*S1); diff --git a/src/math/__tanl.c b/src/math/__tanl.c index 4b36e616..54abc3da 100644 --- a/src/math/__tanl.c +++ b/src/math/__tanl.c @@ -1,4 +1,5 @@ /* origin: FreeBSD /usr/src/lib/msun/ld80/k_tanl.c */ +/* origin: FreeBSD /usr/src/lib/msun/ld128/k_tanl.c */ /* * ==================================================== * Copyright 2004 Sun Microsystems, Inc. All Rights Reserved. @@ -12,7 +13,8 @@ #include "libm.h" -#if LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 +#if (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 +#if LDBL_MANT_DIG == 64 /* * ld80 version of __tan.c. See __tan.c for most comments. */ @@ -22,14 +24,12 @@ * * See __cosl.c for more details about the polynomial. */ - static const long double T3 = 0.333333333333333333180L, /* 0xaaaaaaaaaaaaaaa5.0p-65 */ T5 = 0.133333333333333372290L, /* 0x88888888888893c3.0p-66 */ T7 = 0.0539682539682504975744L, /* 0xdd0dd0dd0dc13ba2.0p-68 */ pio4 = 0.785398163397448309628L, /* 0xc90fdaa22168c235.0p-64 */ pio4lo = -1.25413940316708300586e-20L; /* -0xece675d1fc8f8cbb.0p-130 */ - static const double T9 = 0.021869488536312216, /* 0x1664f4882cc1c2.0p-58 */ T11 = 0.0088632355256619590, /* 0x1226e355c17612.0p-59 */ @@ -44,6 +44,59 @@ T27 = 0.0000024196006108814377, /* 0x144c0d80cc6896.0p-71 */ T29 = 0.0000078293456938132840, /* 0x106b59141a6cb3.0p-69 */ T31 = -0.0000032609076735050182, /* -0x1b5abef3ba4b59.0p-71 */ T33 = 0.0000023261313142559411; /* 0x13835436c0c87f.0p-71 */ +#define RPOLY(w) (T5 + w * (T9 + w * (T13 + w * (T17 + w * (T21 + \ + w * (T25 + w * (T29 + w * T33))))))) +#define VPOLY(w) (T7 + w * (T11 + w * (T15 + w * (T19 + w * (T23 + \ + w * (T27 + w * T31)))))) +#elif LDBL_MANT_DIG == 113 +/* + * ld128 version of __tan.c. See __tan.c for most comments. + */ +/* + * Domain [-0.67434, 0.67434], range ~[-3.37e-36, 1.982e-37] + * |tan(x)/x - t(x)| < 2**-117.8 (XXX should be ~1e-37) + * + * See __cosl.c for more details about the polynomial. + */ +static const long double +T3 = 0x1.5555555555555555555555555553p-2L, +T5 = 0x1.1111111111111111111111111eb5p-3L, +T7 = 0x1.ba1ba1ba1ba1ba1ba1ba1b694cd6p-5L, +T9 = 0x1.664f4882c10f9f32d6bbe09d8bcdp-6L, +T11 = 0x1.226e355e6c23c8f5b4f5762322eep-7L, +T13 = 0x1.d6d3d0e157ddfb5fed8e84e27b37p-9L, +T15 = 0x1.7da36452b75e2b5fce9ee7c2c92ep-10L, +T17 = 0x1.355824803674477dfcf726649efep-11L, +T19 = 0x1.f57d7734d1656e0aceb716f614c2p-13L, +T21 = 0x1.967e18afcb180ed942dfdc518d6cp-14L, +T23 = 0x1.497d8eea21e95bc7e2aa79b9f2cdp-15L, +T25 = 0x1.0b132d39f055c81be49eff7afd50p-16L, +T27 = 0x1.b0f72d33eff7bfa2fbc1059d90b6p-18L, +T29 = 0x1.5ef2daf21d1113df38d0fbc00267p-19L, +T31 = 0x1.1c77d6eac0234988cdaa04c96626p-20L, +T33 = 0x1.cd2a5a292b180e0bdd701057dfe3p-22L, +T35 = 0x1.75c7357d0298c01a31d0a6f7d518p-23L, +T37 = 0x1.2f3190f4718a9a520f98f50081fcp-24L, +pio4 = 0x1.921fb54442d18469898cc51701b8p-1L, +pio4lo = 0x1.cd129024e088a67cc74020bbea60p-116L; +static const double +T39 = 0.000000028443389121318352, /* 0x1e8a7592977938.0p-78 */ +T41 = 0.000000011981013102001973, /* 0x19baa1b1223219.0p-79 */ +T43 = 0.0000000038303578044958070, /* 0x107385dfb24529.0p-80 */ +T45 = 0.0000000034664378216909893, /* 0x1dc6c702a05262.0p-81 */ +T47 = -0.0000000015090641701997785, /* -0x19ecef3569ebb6.0p-82 */ +T49 = 0.0000000029449552300483952, /* 0x194c0668da786a.0p-81 */ +T51 = -0.0000000022006995706097711, /* -0x12e763b8845268.0p-81 */ +T53 = 0.0000000015468200913196612, /* 0x1a92fc98c29554.0p-82 */ +T55 = -0.00000000061311613386849674, /* -0x151106cbc779a9.0p-83 */ +T57 = 1.4912469681508012e-10; /* 0x147edbdba6f43a.0p-85 */ +#define RPOLY(w) (T5 + w * (T9 + w * (T13 + w * (T17 + w * (T21 + \ + w * (T25 + w * (T29 + w * (T33 + w * (T37 + w * (T41 + \ + w * (T45 + w * (T49 + w * (T53 + w * T57))))))))))))) +#define VPOLY(w) (T7 + w * (T11 + w * (T15 + w * (T19 + w * (T23 + \ + w * (T27 + w * (T31 + w * (T35 + w * (T39 + w * (T43 + \ + w * (T47 + w * (T51 + w * T55)))))))))))) +#endif long double __tanl(long double x, long double y, int odd) { long double z, r, v, w, s, a, t; @@ -62,10 +115,8 @@ long double __tanl(long double x, long double y, int odd) { } z = x * x; w = z * z; - r = T5 + w * (T9 + w * (T13 + w * (T17 + w * (T21 + - w * (T25 + w * (T29 + w * T33)))))); - v = z * (T7 + w * (T11 + w * (T15 + w * (T19 + w * (T23 + - w * (T27 + w * T31)))))); + r = RPOLY(w); + v = z * VPOLY(w); s = z * x; r = y + z * (s * (r + v) + y) + T3 * s; w = x + r; @@ -76,7 +127,6 @@ long double __tanl(long double x, long double y, int odd) { } if (!odd) return w; - /* * if allow error up to 2 ulp, simply return * -1.0 / (x+r) here diff --git a/src/math/cosl.c b/src/math/cosl.c index 0794d284..79c41c77 100644 --- a/src/math/cosl.c +++ b/src/math/cosl.c @@ -1,34 +1,3 @@ -/* origin: FreeBSD /usr/src/lib/msun/src/s_cosl.c */ -/*- - * Copyright (c) 2007 Steven G. Kargl - * All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice unmodified, this list of conditions, and the following - * disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * - * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR - * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES - * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. - * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, - * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT - * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, - * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY - * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT - * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF - * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. - */ -/* - * Limited testing on pseudorandom numbers drawn within [-2e8:4e8] shows - * an accuracy of <= 0.7412 ULP. - */ - #include "libm.h" #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 @@ -38,44 +7,33 @@ long double cosl(long double x) { #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 long double cosl(long double x) { - union IEEEl2bits z; + union ldshape u = {x}; unsigned n; - long double y[2]; - long double hi, lo; - - z.e = x; - z.bits.sign = 0; + long double y[2], hi, lo; - /* If x = NaN or Inf, then cos(x) = NaN. */ - if (z.bits.exp == 0x7fff) - return (x - x) / (x - x); - - /* |x| < (double)pi/4 */ - if (z.e < M_PI_4) { - /* |x| < 0x1p-64 */ - if (z.bits.exp < 0x3fff - 64) + u.i.se &= 0x7fff; + if (u.i.se == 0x7fff) + return x - x; + x = u.f; + if (x < M_PI_4) { + if (u.i.se < 0x3fff - LDBL_MANT_DIG) /* raise inexact if x!=0 */ return 1.0 + x; - return __cosl(z.e, 0); + return __cosl(x, 0); } - n = __rem_pio2l(x, y); hi = y[0]; lo = y[1]; switch (n & 3) { case 0: - hi = __cosl(hi, lo); - break; + return __cosl(hi, lo); case 1: - hi = -__sinl(hi, lo, 1); - break; + return -__sinl(hi, lo, 1); case 2: - hi = -__cosl(hi, lo); - break; + return -__cosl(hi, lo); case 3: - hi = __sinl(hi, lo, 1); - break; + default: + return __sinl(hi, lo, 1); } - return hi; } #endif diff --git a/src/math/sincosl.c b/src/math/sincosl.c index 5db69bd6..2c600801 100644 --- a/src/math/sincosl.c +++ b/src/math/sincosl.c @@ -9,25 +9,19 @@ void sincosl(long double x, long double *sin, long double *cos) #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 void sincosl(long double x, long double *sin, long double *cos) { - union IEEEl2bits u; + union ldshape u = {x}; unsigned n; long double y[2], s, c; - u.e = x; - u.bits.sign = 0; - - /* x = nan or inf */ - if (u.bits.exp == 0x7fff) { + u.i.se &= 0x7fff; + if (u.i.se == 0x7fff) { *sin = *cos = x - x; return; } - - /* |x| < (double)pi/4 */ - if (u.e < M_PI_4) { - /* |x| < 0x1p-64 */ - if (u.bits.exp < 0x3fff - 64) { + if (u.f < M_PI_4) { + if (u.i.se < 0x3fff - LDBL_MANT_DIG) { /* raise underflow if subnormal */ - if (u.bits.exp == 0) FORCE_EVAL(x*0x1p-120f); + if (u.i.se == 0) FORCE_EVAL(x*0x1p-120f); *sin = x; /* raise inexact if x!=0 */ *cos = 1.0 + x; @@ -37,7 +31,6 @@ void sincosl(long double x, long double *sin, long double *cos) *cos = __cosl(x, 0); return; } - n = __rem_pio2l(x, y); s = __sinl(y[0], y[1], 1); c = __cosl(y[0], y[1]); diff --git a/src/math/sinl.c b/src/math/sinl.c index 6ca99986..9c0b16ee 100644 --- a/src/math/sinl.c +++ b/src/math/sinl.c @@ -1,31 +1,3 @@ -/* origin: FreeBSD /usr/src/lib/msun/src/s_sinl.c */ -/*- - * Copyright (c) 2007 Steven G. Kargl - * All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice unmodified, this list of conditions, and the following - * disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * - * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR - * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES - * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. - * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, - * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT - * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, - * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY - * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT - * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF - * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. - */ - - #include "libm.h" #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 @@ -36,46 +8,34 @@ long double sinl(long double x) #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 long double sinl(long double x) { - union IEEEl2bits z; + union ldshape u = {x}; unsigned n; - long double y[2]; - long double hi, lo; - - z.e = x; - z.bits.sign = 0; + long double y[2], hi, lo; - /* If x = NaN or Inf, then sin(x) = NaN. */ - if (z.bits.exp == 0x7fff) - return (x - x) / (x - x); - - /* |x| < (double)pi/4 */ - if (z.e < M_PI_4) { - /* |x| < 0x1p-64 */ - if (z.bits.exp < 0x3fff - 64) { + u.i.se &= 0x7fff; + if (u.i.se == 0x7fff) + return x - x; + if (u.f < M_PI_4) { + if (u.i.se < 0x3fff - LDBL_MANT_DIG/2) { /* raise inexact if x!=0 and underflow if subnormal */ - FORCE_EVAL(z.bits.exp == 0 ? x/0x1p120f : x+0x1p120f); + FORCE_EVAL(u.i.se == 0 ? x*0x1p-120f : x+0x1p120f); return x; } return __sinl(x, 0.0, 0); } - n = __rem_pio2l(x, y); hi = y[0]; lo = y[1]; switch (n & 3) { case 0: - hi = __sinl(hi, lo, 1); - break; + return __sinl(hi, lo, 1); case 1: - hi = __cosl(hi, lo); - break; + return __cosl(hi, lo); case 2: - hi = -__sinl(hi, lo, 1); - break; + return -__sinl(hi, lo, 1); case 3: - hi = -__cosl(hi, lo); - break; + default: + return -__cosl(hi, lo); } - return hi; } #endif diff --git a/src/math/tanl.c b/src/math/tanl.c index 546c7a02..6af06712 100644 --- a/src/math/tanl.c +++ b/src/math/tanl.c @@ -1,35 +1,3 @@ -/* origin: FreeBSD /usr/src/lib/msun/src/s_tanl.c */ -/*- - * Copyright (c) 2007 Steven G. Kargl - * All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice unmodified, this list of conditions, and the following - * disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * - * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR - * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES - * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. - * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, - * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT - * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, - * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY - * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT - * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF - * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. - */ -/* - * Limited testing on pseudorandom numbers drawn within [0:4e8] shows - * an accuracy of <= 1.5 ULP where 247024 values of x out of 40 million - * possibles resulted in tan(x) that exceeded 0.5 ULP (ie., 0.6%). - */ - #include "libm.h" #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 @@ -40,28 +8,21 @@ long double tanl(long double x) #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 long double tanl(long double x) { - union IEEEl2bits z; + union ldshape u = {x}; long double y[2]; unsigned n; - z.e = x; - z.bits.sign = 0; - - /* If x = NaN or Inf, then tan(x) = NaN. */ - if (z.bits.exp == 0x7fff) - return (x - x) / (x - x); - - /* |x| < (double)pi/4 */ - if (z.e < M_PI_4) { - /* |x| < 0x1p-64 */ - if (z.bits.exp < 0x3fff - 64) { + u.i.se &= 0x7fff; + if (u.i.se == 0x7fff) + return x - x; + if (u.f < M_PI_4) { + if (u.i.se < 0x3fff - LDBL_MANT_DIG/2) { /* raise inexact if x!=0 and underflow if subnormal */ - FORCE_EVAL(z.bits.exp == 0 ? x/0x1p120f : x+0x1p120f); + FORCE_EVAL(u.i.se == 0 ? x*0x1p-120f : x+0x1p120f); return x; } return __tanl(x, 0, 0); } - n = __rem_pio2l(x, y); return __tanl(y[0], y[1], n&1); } -- 2.25.1