This is the mail archive of the
newlib@sourceware.org
mailing list for the newlib project.
RFA: Add ! LDBL_EQ_DBL implementations of cosl, sinl and tanl
- From: Nick Clifton <nickc at redhat dot com>
- To: newlib at sourceware dot org
- Date: Mon, 02 Mar 2015 17:52:53 +0000
- Subject: RFA: Add ! LDBL_EQ_DBL implementations of cosl, sinl and tanl
- Authentication-results: sourceware.org; auth=none
Hi Jeff, Hi Corinna,
Attached is a patch that implements cosl(), sinl() and tanl() when
long double is not the same as double. In the course of implementing
these functions I found that I also needed to implement the fabls()
and floorl() functions as well as several internal library support
routines.
Tested in the usual way with an x86_64-linux-gnu toolchain with
direct linking to the math functions.
OK to apply ?
Cheers
Nick
newlib/ChangeLog
2015-03-02 Nick Clifton <nickc@redhat.com>
* libm/common/cosl.c (cosl): Add implementation for when long
double is not the same as double.
* libm/common/fabls.c (fabls): Likewise.
* libm/common/floorl.c (floorl): Likewise.
* libm/common/sinl.c (sinl): Likewise.
* libm/common/tanl.c (tanl): Likewise.
* libm/common/fdlibmh.h: Add prototypes for __ieee754_rem_pio2l,
__kernel_cosl, __kernel_cosl, __kernel_sinl and __kernel_tanl.
* libm/math/e_rem_pio2l.c: New file. Implements
__ieee754_rem_pio2l.
* libm/math/k_cosl.c: New file. Implements __kernel_cosl.
* libm/math/k_sinl.c: New file. Implements __kernel_sinl.
* libm/math/k_tanl.c: New file. Implements __kernel_tanl.
* libm/math/Makefile.am (src): Add k_cosl.c, k_sinl.c, k_tanl.c
and e_rem_pio2l.c.
* libm/math/Makefile.in: Regenerate.
*** /dev/null 2015-03-02 08:03:54.706972409 +0000
--- newlib/libm/math/e_rem_pio2l.c 2015-03-02 17:25:47.115144573 +0000
***************
*** 0 ****
--- 1,149 ----
+ /*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ * Copyright (c) 2008 Steven G. Kargl, David Schultz, Bruce D. Evans.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ *
+ * Optimized by Bruce D. Evans.
+ */
+
+ #include <math.h>
+ #include "ieeefp.h"
+ #include "fdlibm.h"
+
+ #define BIAS (LDBL_MAX_EXP - 1)
+
+ static const __int32_t two_over_pi[] =
+ {
+ 0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62,
+ 0x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A,
+ 0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129,
+ 0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41,
+ 0x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8,
+ 0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF,
+ 0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5,
+ 0xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08,
+ 0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3,
+ 0x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880,
+ 0x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B,
+ };
+
+ static const double
+ zero = 0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
+ 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 */
+
+ int
+ __ieee754_rem_pio2l (long double x, long double *y)
+ {
+ union ieee_ext_u u,u1;
+ long double z,w,t,r,fn;
+ double tx[3],ty[2];
+ int e0,ex,i,j,nx,n;
+ int expsign;
+
+ u.extu_ld = x;
+ expsign = u.extu_ext.ext_sign;
+ ex = expsign & 0x7fff;
+
+ if (ex < BIAS + 25 || (ex == BIAS + 25 && u.extu_ext.ext_frach < 0xc90fdaa2))
+ {
+ union ieee_ext_u u2;
+ int ex1;
+
+ /* |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;
+ n = fn;
+ r = x - fn * pio2_1;
+ w = fn * pio2_1t; /* 1st round good to 102 bit. */
+
+ j = ex;
+ y[0] = r - w;
+ u2.extu_ld = y[0];
+ ex1 = u2.extu_ext.ext_sign & 0x7fff;
+ i = j - ex1;
+
+ if (i > 22)
+ {
+ /* 2nd iteration needed, good to 141. */
+ t = r;
+ w = fn * pio2_2;
+ r = t - w;
+ w = fn * pio2_2t - ((t - r) - w);
+ y[0] = r - w;
+ u2.extu_ld = y[0];
+ ex1 = u2.extu_ext.ext_sign & 0x7fff;
+ i = j - ex1;
+
+ if (i > 61)
+ {
+ /* 3rd iteration need, 180 bits acc. */
+ t = r; /* will cover all possible cases. */
+ w = fn * pio2_3;
+ r = t - w;
+ w = fn * pio2_3t - ((t - r) - w);
+ y[0] = r - w;
+ }
+ }
+
+ y[1] = (r - y[0]) - w;
+ return n;
+ }
+
+ /* All other (large) arguments. */
+ if (ex == 0x7fff)
+ {
+ /* x is inf or NaN. */
+ y[0] = y[1] = x - x;
+ return 0;
+ }
+
+ /* Set z = scalbn(|x|,ilogb(x)-23). */
+ u1.extu_ld = x;
+ e0 = ex - BIAS - 23; /* e0 = ilogb(|x|)-23; */
+ u1.extu_ext.ext_sign = ex - e0;
+ z = u1.extu_ld;
+
+ for (i = 0; i < 2; i++)
+ {
+ tx[i] = (double) ((__int32_t)(z));
+ z = (z - tx[i]) * two24;
+ }
+
+ tx[2] = z;
+ nx = 3;
+
+ while (tx[nx - 1] == zero)
+ nx--; /* Skip zero term. */
+
+ n = __kernel_rem_pio2 (tx, ty, e0, nx, 2, two_over_pi);
+
+ r = (long double) ty[0] + ty[1];
+ w = ty[1] - (r - ty[0]);
+
+ if (expsign < 0)
+ {
+ y[0] = -r;
+ y[1] = -w;
+ return -n;
+ }
+
+ y[0] = r;
+ y[1] = w;
+ return n;
+ }
*** /dev/null 2015-03-02 08:03:54.706972409 +0000
--- newlib/libm/math/k_cosl.c 2015-03-02 10:58:48.373846632 +0000
***************
*** 0 ****
--- 1,54 ----
+ /*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ * Copyright (c) 2008 Steven G. Kargl, David Schultz, Bruce D. Evans.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+ #include "fdlibm.h"
+
+ /*
+ * 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
+ * ../ld80/k_cosl.c for more details.
+ */
+
+ static const double one = 1.0;
+ 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;
+
+ long double
+ __kernel_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 + z * (C8 + z * (C9 + z * (C10 + z * C11))))))))));
+ hz = 0.5 * z;
+
+ w = one - hz;
+
+ return w + (((one - w) - hz) + (z * r - x * y));
+ }
*** /dev/null 2015-03-02 08:03:54.706972409 +0000
--- newlib/libm/math/k_sinl.c 2015-03-02 11:04:59.885792456 +0000
***************
*** 0 ****
--- 1,51 ----
+ /*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ * Copyright (c) 2008 Steven G. Kargl, David Schultz, Bruce D. Evans.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+ #include "fdlibm.h"
+
+ /*
+ * Domain [-0.7854, 0.7854], range ~[-1.53e-37, 1.659e-37]
+ * |sin(x)/x - s(x)| < 2**-122.1
+ */
+
+ static const double half = 0.5;
+
+ 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;
+
+ long double
+ __kernel_sinl (long double x, long double y, int iy)
+ {
+ long double z, r, v;
+
+ z = x * x;
+ v = z * x;
+ r = S2 + z * (S3 + z * (S4 + z * (S5 + z * (S6 + z * (S7 + z * (S8
+ + z * (S9 + z * (S10 + z * (S11 + z * S12)))))))));
+
+ if (iy == 0)
+ return x + v * (S1 + z * r);
+
+ return x -((z * (half * y - v * r) - y) - v * S1);
+ }
*** /dev/null 2015-03-02 08:03:54.706972409 +0000
--- newlib/libm/math/k_tanl.c 2015-03-02 11:10:00.318365714 +0000
***************
*** 0 ****
--- 1,120 ----
+ /*
+ * ====================================================
+ * Copyright 2004 Sun Microsystems, Inc. All Rights Reserved.
+ * Copyright (c) 2008 Steven G. Kargl, David Schultz, Bruce D. Evans.
+ *
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+ #include "fdlibm.h"
+
+ /*
+ * 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)
+ */
+
+ 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 */
+
+ long double
+ __kernel_tanl (long double x, long double y, int iy)
+ {
+ long double z, r, v, w, s;
+ long double osign;
+ int i;
+
+ iy = (iy == 1 ? -1 : 1); /* XXX recover original interface. */
+ osign = (x >= 0 ? 1.0 : -1.0); /* XXX slow, probably wrong for -0. */
+
+ if (fabsl (x) >= 0.67434)
+ {
+ if (x < 0)
+ {
+ x = -x;
+ y = -y;
+ }
+ z = pio4 - x;
+ w = pio4lo - y;
+ x = z + w;
+ y = 0.0;
+ i = 1;
+ }
+ else
+ i = 0;
+
+ z = x * x;
+ w = z * z;
+ r = 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))))))))))));
+
+ v = z * (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))))))))))));
+
+ s = z * x;
+ r = y + z * (s * (r + v) + y);
+ r += T3 * s;
+ w = x + r;
+
+ if (i == 1)
+ {
+ v = (long double) iy;
+
+ return osign *
+ (v - 2.0 * (x - (w * w / (w + v) - r)));
+ }
+ if (iy == 1)
+ return w;
+ else
+ {
+ /*
+ * If allow error up to 2 ulp, simply return -1.0 / (x+r) here.
+ */
+ /* compute -1.0 / (x+r) accurately. */
+ long double a, t;
+
+ z = w;
+ z = z + 0x1p32 - 0x1p32;
+ v = r - (z - x); /* z+v = r+x */
+ t = a = -1.0 / w; /* a = -1.0/w */
+ t = t + 0x1p32 - 0x1p32;
+ s = 1.0 + t * z;
+ return t + a * (s + t * v);
+ }
+ }
Index: newlib/libm/common/cosl.c
===================================================================
RCS file: /cvs/src/src/newlib/libm/common/cosl.c,v
retrieving revision 1.2
diff -u -3 -p -r1.2 cosl.c
--- newlib/libm/common/cosl.c 17 Apr 2009 22:15:42 -0000 1.2
+++ newlib/libm/common/cosl.c 2 Mar 2015 17:40:35 -0000
@@ -38,5 +38,44 @@ cosl (long double x)
{
return cos(x);
}
-#endif
+#else
+
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+#include "fdlibm.h"
+
+long double
+cosl (long double x)
+{
+ long double y[2];
+ int n;
+ /* |x| ~< pi/4 */
+ if (x >= -0.7853981633974483096156608458198757210492
+ && x <= 0.7853981633974483096156608458198757210492)
+ return __kernel_cosl (x, 0.0L);
+
+ else if ((x + x == x && x != 0.0) || x != x)
+ return x - x; /* NaN */
+
+ /* Argument reduction needed. */
+ n = __ieee754_rem_pio2l (x, y);
+
+ switch (n & 3)
+ {
+ case 0: return __kernel_cosl (y[0], y[1]);
+ case 1: return - __kernel_sinl (y[0], y[1], 1);
+ case 2: return - __kernel_cosl (y[0], y[1]);
+ default: return __kernel_sinl (y[0], y[1], 1);
+ }
+}
+#endif
Index: newlib/libm/common/fabsl.c
===================================================================
RCS file: /cvs/src/src/newlib/libm/common/fabsl.c,v
retrieving revision 1.2
diff -u -3 -p -r1.2 fabsl.c
--- newlib/libm/common/fabsl.c 17 Apr 2009 22:15:42 -0000 1.2
+++ newlib/libm/common/fabsl.c 2 Mar 2015 17:40:35 -0000
@@ -29,6 +29,7 @@ POSSIBILITY OF SUCH DAMAGE.
*/
#include <math.h>
+#include "ieeefp.h"
#include "local.h"
/* On platforms where long double is as wide as double. */
@@ -38,5 +39,13 @@ fabsl (long double x)
{
return fabs(x);
}
-#endif
+#else
+long double
+fabsl (long double x)
+{
+ union ieee_ext_u u = { x };
+ u.extu_ext.ext_sign = 0;
+ return u.extu_ld;
+}
+#endif
Index: newlib/libm/common/fdlibm.h
===================================================================
RCS file: /cvs/src/src/newlib/libm/common/fdlibm.h,v
retrieving revision 1.10
diff -u -3 -p -r1.10 fdlibm.h
--- newlib/libm/common/fdlibm.h 6 Feb 2015 16:14:03 -0000 1.10
+++ newlib/libm/common/fdlibm.h 2 Mar 2015 17:40:35 -0000
@@ -186,6 +186,12 @@ extern double __kernel_cos __P((double,d
extern double __kernel_tan __P((double,double,int));
extern int __kernel_rem_pio2 __P((double*,double*,int,int,int,const __int32_t*));
+extern long double __kernel_sinl __P((long double, long double, int));
+extern long double __kernel_cosl __P((long double, long double));
+extern long double __kernel_tanl __P((long double, long double, int));
+extern int __kernel_rem_pio2l __P((long double *, long double *, int, int, int, const __int32_t *));
+extern __int32_t __ieee754_rem_pio2l __P((long double, long double *));
+
/* Undocumented float functions. */
#ifdef _SCALB_INT
extern float scalbf __P((float, int));
Index: newlib/libm/common/floorl.c
===================================================================
RCS file: /cvs/src/src/newlib/libm/common/floorl.c,v
retrieving revision 1.2
diff -u -3 -p -r1.2 floorl.c
--- newlib/libm/common/floorl.c 17 Apr 2009 22:15:42 -0000 1.2
+++ newlib/libm/common/floorl.c 2 Mar 2015 17:40:35 -0000
@@ -38,5 +38,121 @@ floorl (long double x)
{
return floor(x);
}
+#else
+
+#include "ieeefp.h"
+
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+#ifdef LDBL_IMPLICIT_NBIT
+
+#define MANH_SIZE (EXT_FRACHBITS + 1)
+#define INC_MANH(u, c) \
+ do \
+ { \
+ __ieee_ext_field_type o = u.extu_ext.ext_frach; \
+ u.extu_ext.ext_frach += (c); \
+ if (u.extu_ext.ext_frach < o) \
+ u.extu_ext.ext_exp ++; \
+ } \
+ while (0)
+
+#else
+
+#define MANH_SIZE EXT_FRACHBITS
+#define INC_MANH(u, c) \
+ do \
+ { \
+ __ieee_ext_field_type o = u.extu_ext.ext_frach; \
+ u.extu_ext.ext_frach += (c); \
+ if (u.extu_ext.ext_frach < o) \
+ { \
+ u.extu_ext.ext_exp++; \
+ u.extu_ext.ext_frach |= 1LLU << (MANH_SIZE - 1); \
+ } \
+ } \
+ while (0)
+
#endif
+static const long double huge = 1.0e300;
+
+long double
+floorl (long double x)
+{
+ union ieee_ext_u u = { x };
+ int e = u.extu_ext.ext_exp - LDBL_MAX_EXP + 1;
+
+ if (e < MANH_SIZE - 1)
+ {
+ if (e < 0)
+ {
+ /* Raise inexact if x != 0. */
+ if (huge + x > 0.0)
+ if (u.extu_ext.ext_exp > 0 || (u.extu_ext.ext_frach | u.extu_ext.ext_fracl) != 0)
+ u.extu_ld = u.extu_ext.ext_sign ? -1.0L : 0.0L;
+ }
+ else
+ {
+ unsigned long long m = ((1LLU << MANH_SIZE) - 1) >> (e + 1);
+
+ if (((u.extu_ext.ext_frach & m) | u.extu_ext.ext_fracl) == 0)
+ return x; /* x is integral. */
+
+ if (u.extu_ext.ext_sign)
+ {
+#ifdef LDBL_IMPLICIT_NBIT
+ if (e == 0)
+ u.extu_ext.ext_exp++;
+ else
+#endif
+ INC_MANH (u, 1LLU << (MANH_SIZE - e - 1));
+ }
+
+ if (huge + x > 0.0)
+ {
+ /* Raise inexact flag. */
+ u.extu_ext.ext_frach &= ~m;
+ u.extu_ext.ext_fracl = 0;
+ }
+ }
+ }
+ else if (e < LDBL_MANT_DIG - 1)
+ {
+ __ieee_ext_field_type m = (__ieee_ext_field_type)-1 >> (64 - LDBL_MANT_DIG + e + 1);
+
+ if ((u.extu_ext.ext_fracl & m) == 0)
+ return x; /* x is integral. */
+
+ if (u.extu_ext.ext_sign)
+ {
+ if (e == MANH_SIZE - 1)
+ INC_MANH (u, 1);
+ else
+ {
+ __ieee_ext_field_type o = u.extu_ext.ext_fracl;
+
+ u.extu_ext.ext_fracl += 1llu << (LDBL_MANT_DIG - e - 1);
+ if (u.extu_ext.ext_fracl < o)
+ /* Got a carry. */
+ INC_MANH (u, 1);
+ }
+ }
+
+ if (huge + x > 0.0L)
+ /* Raise inexact flag. */
+ u.extu_ext.ext_fracl &= ~m;
+ }
+
+ return u.extu_ld;
+}
+#endif
Index: newlib/libm/common/sinl.c
===================================================================
RCS file: /cvs/src/src/newlib/libm/common/sinl.c,v
retrieving revision 1.2
diff -u -3 -p -r1.2 sinl.c
--- newlib/libm/common/sinl.c 17 Apr 2009 22:15:43 -0000 1.2
+++ newlib/libm/common/sinl.c 2 Mar 2015 17:40:35 -0000
@@ -38,5 +38,45 @@ sinl (long double x)
{
return sin(x);
}
+#else
+
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+#include "fdlibm.h"
+
+long double
+sinl (long double x)
+{
+ long double y[2];
+ int n;
+
+ /* |x| ~< pi/4 */
+ if (x >= -0.7853981633974483096156608458198757210492
+ && x <= 0.7853981633974483096156608458198757210492)
+ return __kernel_sinl (x, 0.0L, 0);
+
+ else if (x + x == x || x != x)
+ return x - x; /* NaN */
+
+ /* Argument reduction needed. */
+ n = __ieee754_rem_pio2l (x, y);
+
+ switch (n & 3)
+ {
+ case 0: return __kernel_sinl (y[0], y[1], 1);
+ case 1: return __kernel_cosl (y[0], y[1]);
+ case 2: return - __kernel_sinl (y[0], y[1], 1);
+ default:return - __kernel_cosl (y[0], y[1]);
+ }
+}
#endif
Index: newlib/libm/common/tanl.c
===================================================================
RCS file: /cvs/src/src/newlib/libm/common/tanl.c,v
retrieving revision 1.2
diff -u -3 -p -r1.2 tanl.c
--- newlib/libm/common/tanl.c 17 Apr 2009 22:15:43 -0000 1.2
+++ newlib/libm/common/tanl.c 2 Mar 2015 17:40:35 -0000
@@ -38,5 +38,41 @@ tanl (long double x)
{
return tan(x);
}
+#else
+
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+#include "fdlibm.h"
+
+long double
+tanl (long double x)
+{
+ long double y[2];
+ int n;
+
+ /* |x| ~< pi/4 */
+ if (x >= -0.7853981633974483096156608458198757210492
+ && x <= 0.7853981633974483096156608458198757210492)
+ return __kernel_tanl (x, 0.0L, 1);
+
+ /* tanl(Inf or NaN) is NaN, tanl(0) is 0. */
+ else if (x + x == x || x != x)
+ return x - x; /* NaN */
+
+ /* Argument reduction needed. */
+ n = __ieee754_rem_pio2l (x, y);
+
+ return __kernel_tanl (y[0], y[1], 1 - ((n & 1) << 1));
+}
+
#endif
Index: newlib/libm/math/Makefile.am
===================================================================
RCS file: /cvs/src/src/newlib/libm/math/Makefile.am,v
retrieving revision 1.10
diff -u -3 -p -r1.10 Makefile.am
--- newlib/libm/math/Makefile.am 6 Feb 2015 16:14:03 -0000 1.10
+++ newlib/libm/math/Makefile.am 2 Mar 2015 17:40:35 -0000
@@ -6,11 +6,13 @@ INCLUDES = -I$(srcdir)/../common $(NEWLI
src = k_standard.c k_rem_pio2.c \
k_cos.c k_sin.c k_tan.c \
+ k_cosl.c k_sinl.c k_tanl.c \
e_acos.c e_acosh.c e_asin.c e_atan2.c \
e_atanh.c e_cosh.c e_exp.c e_fmod.c \
er_gamma.c e_hypot.c e_j0.c \
e_j1.c e_jn.c er_lgamma.c \
- e_log.c e_log10.c e_pow.c e_rem_pio2.c e_remainder.c \
+ e_log.c e_log10.c e_pow.c e_rem_pio2.c \
+ e_rem_pio2l.c e_remainder.c \
e_scalb.c e_sinh.c e_sqrt.c \
w_acos.c w_acosh.c w_asin.c w_atan2.c \
w_atanh.c w_cosh.c w_exp.c w_fmod.c \