Transcendental functions for 'long double', from Paolo Bonzini.
authorBruno Haible <bruno@clisp.org>
Tue, 18 Feb 2003 17:05:23 +0000 (17:05 +0000)
committerBruno Haible <bruno@clisp.org>
Tue, 18 Feb 2003 17:05:23 +0000 (17:05 +0000)
22 files changed:
ChangeLog
lib/ChangeLog
lib/acosl.c [new file with mode: 0644]
lib/asinl.c [new file with mode: 0644]
lib/atanl.c [new file with mode: 0644]
lib/ceill.c [new file with mode: 0644]
lib/cosl.c [new file with mode: 0644]
lib/expl.c [new file with mode: 0644]
lib/floorl.c [new file with mode: 0644]
lib/frexpl.c [new file with mode: 0644]
lib/ldexpl.c [new file with mode: 0644]
lib/logl.c [new file with mode: 0644]
lib/mathl.h [new file with mode: 0644]
lib/sincosl.c [new file with mode: 0644]
lib/sinl.c [new file with mode: 0644]
lib/sqrtl.c [new file with mode: 0644]
lib/tanl.c [new file with mode: 0644]
lib/trigl.c [new file with mode: 0644]
lib/trigl.h [new file with mode: 0644]
m4/ChangeLog
m4/mathl.m4 [new file with mode: 0644]
modules/mathl [new file with mode: 0644]

index 7129a62..95bec6f 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,7 @@
+2003-02-18  Paolo Bonzino  <bonzini@gnu.org>
+
+       * modules/mathl: New file.
+
 2003-02-17  Bruno Haible  <bruno@clisp.org>
 
        * modules/mkdtemp: New module.
index 1279061..553a8a1 100644 (file)
@@ -1,3 +1,23 @@
+2003-02-18  Paolo Bonzino  <bonzini@gnu.org>
+
+       * mathl.h: New file.
+       * acosl.c: New file.
+       * asinl.c: New file.
+       * atanl.c: New file.
+       * ceill.c: New file.
+       * cosl.c: New file.
+       * expl.c: New file.
+       * floorl.c: New file.
+       * frexpl.c: New file.
+       * ldexpl.c: New file.
+       * logl.c: New file.
+       * sincosl.c: New file.
+       * sinl.c: New file.
+       * sqrtl.c: New file.
+       * tanl.c: New file.
+       * trigl.c: New file.
+       * trigl.h: New file.
+
 2003-02-17  Bruno Haible  <bruno@clisp.org>
 
        * mkdtemp.h: New file, from GNU gettext.
diff --git a/lib/acosl.c b/lib/acosl.c
new file mode 100644 (file)
index 0000000..d430254
--- /dev/null
@@ -0,0 +1,224 @@
+/*
+ * ====================================================
+ * 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 "mathl.h"
+
+/*
+  Long double expansions contributed by
+  Stephen L. Moshier <moshier@na-net.ornl.gov>
+*/
+
+/* asin(x)
+ * Method :
+ *     Since  asin(x) = x + x^3/6 + x^5*3/40 + x^7*15/336 + ...
+ *     we approximate asin(x) on [0,0.5] by
+ *             asin(x) = x + x*x^2*R(x^2)
+ *      Between .5 and .625 the approximation is
+ *              asin(0.5625 + x) = asin(0.5625) + x rS(x) / sS(x)
+ *     For x in [0.625,1]
+ *             asin(x) = pi/2-2*asin(sqrt((1-x)/2))
+ *
+ * Special cases:
+ *     if x is NaN, return x itself;
+ *     if |x|>1, return NaN with invalid signal.
+ *
+ */
+
+
+#include <math.h>
+
+static const long double
+  one = 1.0L,
+  huge = 1.0e+4932L,
+  pi =      3.1415926535897932384626433832795028841972L,
+  pio2_hi = 1.5707963267948966192313216916397514420986L,
+  pio2_lo = 4.3359050650618905123985220130216759843812E-35L,
+  pio4_hi = 7.8539816339744830961566084581987569936977E-1L,
+
+       /* coefficient for R(x^2) */
+
+  /* asin(x) = x + x^3 pS(x^2) / qS(x^2)
+     0 <= x <= 0.5
+     peak relative error 1.9e-35  */
+  pS0 = -8.358099012470680544198472400254596543711E2L,
+  pS1 =  3.674973957689619490312782828051860366493E3L,
+  pS2 = -6.730729094812979665807581609853656623219E3L,
+  pS3 =  6.643843795209060298375552684423454077633E3L,
+  pS4 = -3.817341990928606692235481812252049415993E3L,
+  pS5 =  1.284635388402653715636722822195716476156E3L,
+  pS6 = -2.410736125231549204856567737329112037867E2L,
+  pS7 =  2.219191969382402856557594215833622156220E1L,
+  pS8 = -7.249056260830627156600112195061001036533E-1L,
+  pS9 =  1.055923570937755300061509030361395604448E-3L,
+
+  qS0 = -5.014859407482408326519083440151745519205E3L,
+  qS1 =  2.430653047950480068881028451580393430537E4L,
+  qS2 = -4.997904737193653607449250593976069726962E4L,
+  qS3 =  5.675712336110456923807959930107347511086E4L,
+  qS4 = -3.881523118339661268482937768522572588022E4L,
+  qS5 =  1.634202194895541569749717032234510811216E4L,
+  qS6 = -4.151452662440709301601820849901296953752E3L,
+  qS7 =  5.956050864057192019085175976175695342168E2L,
+  qS8 = -4.175375777334867025769346564600396877176E1L,
+  /* 1.000000000000000000000000000000000000000E0 */
+
+  /* asin(0.5625 + x) = asin(0.5625) + x rS(x) / sS(x)
+     -0.0625 <= x <= 0.0625
+     peak relative error 3.3e-35  */
+  rS0 = -5.619049346208901520945464704848780243887E0L,
+  rS1 =  4.460504162777731472539175700169871920352E1L,
+  rS2 = -1.317669505315409261479577040530751477488E2L,
+  rS3 =  1.626532582423661989632442410808596009227E2L,
+  rS4 = -3.144806644195158614904369445440583873264E1L,
+  rS5 = -9.806674443470740708765165604769099559553E1L,
+  rS6 =  5.708468492052010816555762842394927806920E1L,
+  rS7 =  1.396540499232262112248553357962639431922E1L,
+  rS8 = -1.126243289311910363001762058295832610344E1L,
+  rS9 = -4.956179821329901954211277873774472383512E-1L,
+  rS10 =  3.313227657082367169241333738391762525780E-1L,
+
+  sS0 = -4.645814742084009935700221277307007679325E0L,
+  sS1 =  3.879074822457694323970438316317961918430E1L,
+  sS2 = -1.221986588013474694623973554726201001066E2L,
+  sS3 =  1.658821150347718105012079876756201905822E2L,
+  sS4 = -4.804379630977558197953176474426239748977E1L,
+  sS5 = -1.004296417397316948114344573811562952793E2L,
+  sS6 =  7.530281592861320234941101403870010111138E1L,
+  sS7 =  1.270735595411673647119592092304357226607E1L,
+  sS8 = -1.815144839646376500705105967064792930282E1L,
+  sS9 = -7.821597334910963922204235247786840828217E-2L,
+  /*  1.000000000000000000000000000000000000000E0 */
+
+ asinr5625 =  5.9740641664535021430381036628424864397707E-1L;
+
+
+long double
+acosl (long double x)
+{
+  long double t, p, q;
+
+  if (x < 0.0L)
+    {
+      t = pi - acosl(-x);
+      if (huge + x > one) /* return with inexact */
+        return t;
+    }
+
+  if (x >= 1.0L)       /* |x|>= 1 */
+    {
+      if (x == 1.0L)
+       return 0.0L;   /* return zero */
+
+      return (x - x) / (x - x);        /* asin(|x|>1) is NaN */
+    }
+
+  else if (x < 0.5L) /* |x| < 0.5 */
+    {
+      if (x < 0.000000000000000006938893903907228377647697925567626953125L) /* |x| < 2**-57 */
+       /* acos(0)=+-pi/2 with inexact */
+       return x * pio2_hi + x * pio2_lo;
+
+      t = x * x;
+      p = (((((((((pS9 * t
+                   + pS8) * t
+                  + pS7) * t
+                 + pS6) * t
+                + pS5) * t
+               + pS4) * t
+              + pS3) * t
+             + pS2) * t
+            + pS1) * t
+           + pS0) * t;
+
+      q = (((((((( t
+                  + qS8) * t
+                 + qS7) * t
+                + qS6) * t
+               + qS5) * t
+              + qS4) * t
+             + qS3) * t
+            + qS2) * t
+           + qS1) * t
+        + qS0;
+
+      return pio2_hi - (x + x * (p / q) - pio2_lo);
+    }
+
+  else if (x < 0.625) /* 0.625 */
+    {
+      t = x - 0.5625;
+      p = ((((((((((rS10 * t
+                   + rS9) * t
+                  + rS8) * t
+                 + rS7) * t
+                + rS6) * t
+               + rS5) * t
+              + rS4) * t
+             + rS3) * t
+            + rS2) * t
+           + rS1) * t
+          + rS0) * t;
+
+      q = ((((((((( t
+                   + sS9) * t
+                 + sS8) * t
+                + sS7) * t
+               + sS6) * t
+              + sS5) * t
+             + sS4) * t
+            + sS3) * t
+           + sS2) * t
+          + sS1) * t
+       + sS0;
+
+      return (pio2_hi - asinr5625) - (p / q - pio2_lo);
+    }
+  else
+    return 2 * asinl(sqrtl((1-x)/2));
+}
+
+#if 0
+main()
+{
+  printf ("%.18Lg %.18Lg\n",
+          acosl(1.0L),
+          1.5707963267948966192313216916397514420984L - 
+          1.5707963267948966192313216916397514420984L);
+  printf ("%.18Lg %.18Lg\n",
+          acosl(0.7071067811865475244008443621048490392848L),
+          1.5707963267948966192313216916397514420984L - 
+          0.7853981633974483096156608458198757210492L);
+  printf ("%.18Lg %.18Lg\n",
+          acosl(0.5L),
+          1.5707963267948966192313216916397514420984L - 
+          0.5235987755982988730771072305465838140328L);
+  printf ("%.18Lg %.18Lg\n",
+          acosl(0.3090169943749474241022934171828190588600L),
+          1.5707963267948966192313216916397514420984L - 
+          0.3141592653589793238462643383279502884196L);
+  printf ("%.18Lg %.18Lg\n",
+          acosl(-1.0L),
+          1.5707963267948966192313216916397514420984L - 
+          -1.5707963267948966192313216916397514420984L);
+  printf ("%.18Lg %.18Lg\n",
+          acosl(-0.7071067811865475244008443621048490392848L),
+          1.5707963267948966192313216916397514420984L - 
+          -0.7853981633974483096156608458198757210492L);
+  printf ("%.18Lg %.18Lg\n",
+          acosl(-0.5L),
+          1.5707963267948966192313216916397514420984L - 
+          -0.5235987755982988730771072305465838140328L);
+  printf ("%.18Lg %.18Lg\n",
+          acosl(-0.3090169943749474241022934171828190588600L),
+          1.5707963267948966192313216916397514420984L - 
+          -0.3141592653589793238462643383279502884196L);
+}
+#endif
diff --git a/lib/asinl.c b/lib/asinl.c
new file mode 100644 (file)
index 0000000..9230240
--- /dev/null
@@ -0,0 +1,218 @@
+/*
+ * ====================================================
+ * 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 "mathl.h"
+
+/*
+  Long double expansions contributed by
+  Stephen L. Moshier <moshier@na-net.ornl.gov>
+*/
+
+/* asin(x)
+ * Method :
+ *     Since  asin(x) = x + x^3/6 + x^5*3/40 + x^7*15/336 + ...
+ *     we approximate asin(x) on [0,0.5] by
+ *             asin(x) = x + x*x^2*R(x^2)
+ *      Between .5 and .625 the approximation is
+ *              asin(0.5625 + x) = asin(0.5625) + x rS(x) / sS(x)
+ *     For x in [0.625,1]
+ *             asin(x) = pi/2-2*asin(sqrt((1-x)/2))
+ *
+ * Special cases:
+ *     if x is NaN, return x itself;
+ *     if |x|>1, return NaN with invalid signal.
+ *
+ */
+
+
+#include <math.h>
+
+static const long double
+  one = 1.0L,
+  huge = 1.0e+4932L,
+  pio2_hi = 1.5707963267948966192313216916397514420986L,
+  pio2_lo = 4.3359050650618905123985220130216759843812E-35L,
+  pio4_hi = 7.8539816339744830961566084581987569936977E-1L,
+
+       /* coefficient for R(x^2) */
+
+  /* asin(x) = x + x^3 pS(x^2) / qS(x^2)
+     0 <= x <= 0.5
+     peak relative error 1.9e-35  */
+  pS0 = -8.358099012470680544198472400254596543711E2L,
+  pS1 =  3.674973957689619490312782828051860366493E3L,
+  pS2 = -6.730729094812979665807581609853656623219E3L,
+  pS3 =  6.643843795209060298375552684423454077633E3L,
+  pS4 = -3.817341990928606692235481812252049415993E3L,
+  pS5 =  1.284635388402653715636722822195716476156E3L,
+  pS6 = -2.410736125231549204856567737329112037867E2L,
+  pS7 =  2.219191969382402856557594215833622156220E1L,
+  pS8 = -7.249056260830627156600112195061001036533E-1L,
+  pS9 =  1.055923570937755300061509030361395604448E-3L,
+
+  qS0 = -5.014859407482408326519083440151745519205E3L,
+  qS1 =  2.430653047950480068881028451580393430537E4L,
+  qS2 = -4.997904737193653607449250593976069726962E4L,
+  qS3 =  5.675712336110456923807959930107347511086E4L,
+  qS4 = -3.881523118339661268482937768522572588022E4L,
+  qS5 =  1.634202194895541569749717032234510811216E4L,
+  qS6 = -4.151452662440709301601820849901296953752E3L,
+  qS7 =  5.956050864057192019085175976175695342168E2L,
+  qS8 = -4.175375777334867025769346564600396877176E1L,
+  /* 1.000000000000000000000000000000000000000E0 */
+
+  /* asin(0.5625 + x) = asin(0.5625) + x rS(x) / sS(x)
+     -0.0625 <= x <= 0.0625
+     peak relative error 3.3e-35  */
+  rS0 = -5.619049346208901520945464704848780243887E0L,
+  rS1 =  4.460504162777731472539175700169871920352E1L,
+  rS2 = -1.317669505315409261479577040530751477488E2L,
+  rS3 =  1.626532582423661989632442410808596009227E2L,
+  rS4 = -3.144806644195158614904369445440583873264E1L,
+  rS5 = -9.806674443470740708765165604769099559553E1L,
+  rS6 =  5.708468492052010816555762842394927806920E1L,
+  rS7 =  1.396540499232262112248553357962639431922E1L,
+  rS8 = -1.126243289311910363001762058295832610344E1L,
+  rS9 = -4.956179821329901954211277873774472383512E-1L,
+  rS10 =  3.313227657082367169241333738391762525780E-1L,
+
+  sS0 = -4.645814742084009935700221277307007679325E0L,
+  sS1 =  3.879074822457694323970438316317961918430E1L,
+  sS2 = -1.221986588013474694623973554726201001066E2L,
+  sS3 =  1.658821150347718105012079876756201905822E2L,
+  sS4 = -4.804379630977558197953176474426239748977E1L,
+  sS5 = -1.004296417397316948114344573811562952793E2L,
+  sS6 =  7.530281592861320234941101403870010111138E1L,
+  sS7 =  1.270735595411673647119592092304357226607E1L,
+  sS8 = -1.815144839646376500705105967064792930282E1L,
+  sS9 = -7.821597334910963922204235247786840828217E-2L,
+  /*  1.000000000000000000000000000000000000000E0 */
+
+ asinr5625 =  5.9740641664535021430381036628424864397707E-1L;
+
+
+long double
+asinl (long double x)
+{
+  long double y, t, p, q;
+  int sign;
+
+  sign = 1;
+  y = x;
+  if (x < 0.0L)
+    {
+      sign = -1;
+      y = -x;
+    }
+
+  if (y >= 1.0L)       /* |x|>= 1 */
+    {
+      if (y == 1.0L)
+       /* asin(1)=+-pi/2 with inexact */
+       return x * pio2_hi + x * pio2_lo;
+
+      return (x - x) / (x - x);        /* asin(|x|>1) is NaN */
+    }
+  else if (y < 0.5L) /* |x| < 0.5 */
+    {
+      if (y < 0.000000000000000006938893903907228377647697925567626953125L) /* |x| < 2**-57 */
+       if (huge + y > one)
+         return y;             /* return x with inexact if x!=0 */
+
+      t = x * x;
+      p = (((((((((pS9 * t
+                   + pS8) * t
+                  + pS7) * t
+                 + pS6) * t
+                + pS5) * t
+               + pS4) * t
+              + pS3) * t
+             + pS2) * t
+            + pS1) * t
+           + pS0) * t;
+
+      q = (((((((( t
+                  + qS8) * t
+                 + qS7) * t
+                + qS6) * t
+               + qS5) * t
+              + qS4) * t
+             + qS3) * t
+            + qS2) * t
+           + qS1) * t
+        + qS0;
+
+      return x + x * (p / q);
+    }
+
+  else if (y < 0.625) /* 0.625 */
+    {
+      t = y - 0.5625;
+      p = ((((((((((rS10 * t
+                   + rS9) * t
+                  + rS8) * t
+                 + rS7) * t
+                + rS6) * t
+               + rS5) * t
+              + rS4) * t
+             + rS3) * t
+            + rS2) * t
+           + rS1) * t
+          + rS0) * t;
+
+      q = ((((((((( t
+                   + sS9) * t
+                 + sS8) * t
+                + sS7) * t
+               + sS6) * t
+              + sS5) * t
+             + sS4) * t
+            + sS3) * t
+           + sS2) * t
+          + sS1) * t
+       + sS0;
+      t = asinr5625 + p / q;
+    }
+  else
+    t = pio2_hi + pio2_lo - 2 * asinl(sqrtl((1-y)/2));
+
+  return t * sign;
+}
+
+#if 0
+main()
+{
+  printf ("%.18Lg %.18Lg\n",
+          asinl(1.0L),
+          1.5707963267948966192313216916397514420984L);
+  printf ("%.18Lg %.18Lg\n",
+          asinl(0.7071067811865475244008443621048490392848L),
+          0.7853981633974483096156608458198757210492L);
+  printf ("%.18Lg %.18Lg\n",
+          asinl(0.5L),
+          0.5235987755982988730771072305465838140328L);
+  printf ("%.18Lg %.18Lg\n",
+          asinl(0.3090169943749474241022934171828190588600L),
+          0.3141592653589793238462643383279502884196L);
+  printf ("%.18Lg %.18Lg\n",
+          asinl(-1.0L),
+          -1.5707963267948966192313216916397514420984L);
+  printf ("%.18Lg %.18Lg\n",
+          asinl(-0.7071067811865475244008443621048490392848L),
+          -0.7853981633974483096156608458198757210492L);
+  printf ("%.18Lg %.18Lg\n",
+          asinl(-0.5L),
+          -0.5235987755982988730771072305465838140328L);
+  printf ("%.18Lg %.18Lg\n",
+          asinl(-0.3090169943749474241022934171828190588600L),
+          -0.3141592653589793238462643383279502884196L);
+}
+#endif
diff --git a/lib/atanl.c b/lib/atanl.c
new file mode 100644 (file)
index 0000000..2959e0b
--- /dev/null
@@ -0,0 +1,200 @@
+/*                                                     s_atanl.c
+ *
+ *     Inverse circular tangent for 128-bit long double precision
+ *      (arctangent)
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * long double x, y, atanl();
+ *
+ * y = atanl( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns radian angle between -pi/2 and +pi/2 whose tangent is x.
+ *
+ * The function uses a rational approximation of the form
+ * t + t^3 P(t^2)/Q(t^2), optimized for |t| < 0.09375.
+ *
+ * The argument is reduced using the identity
+ *    arctan x - arctan u  =  arctan ((x-u)/(1 + ux))
+ * and an 83-entry lookup table for arctan u, with u = 0, 1/8, ..., 10.25.
+ * Use of the table improves the execution speed of the routine.
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE      -19, 19       4e5       1.7e-34     5.4e-35
+ *
+ *
+ * WARNING:
+ *
+ * This program uses integer operations on bit fields of floating-point
+ * numbers.  It does not work with data structures other than the
+ * structure assumed.
+ *
+ */
+
+/* Copyright 2001 by Stephen L. Moshier <moshier@na-net.ornl.gov>  */
+
+#include "mathl.h"
+
+#include <math.h>
+
+/* arctan(k/8), k = 0, ..., 82 */
+static const long double atantbl[84] = {
+  0.0000000000000000000000000000000000000000E0L,
+  1.2435499454676143503135484916387102557317E-1L, /* arctan(0.125)  */
+  2.4497866312686415417208248121127581091414E-1L,
+  3.5877067027057222039592006392646049977698E-1L,
+  4.6364760900080611621425623146121440202854E-1L,
+  5.5859931534356243597150821640166127034645E-1L,
+  6.4350110879328438680280922871732263804151E-1L,
+  7.1882999962162450541701415152590465395142E-1L,
+  7.8539816339744830961566084581987572104929E-1L,
+  8.4415398611317100251784414827164750652594E-1L,
+  8.9605538457134395617480071802993782702458E-1L,
+  9.4200004037946366473793717053459358607166E-1L,
+  9.8279372324732906798571061101466601449688E-1L,
+  1.0191413442663497346383429170230636487744E0L,
+  1.0516502125483736674598673120862998296302E0L,
+  1.0808390005411683108871567292171998202703E0L,
+  1.1071487177940905030170654601785370400700E0L,
+  1.1309537439791604464709335155363278047493E0L,
+  1.1525719972156675180401498626127513797495E0L,
+  1.1722738811284763866005949441337046149712E0L,
+  1.1902899496825317329277337748293183376012E0L,
+  1.2068173702852525303955115800565576303133E0L,
+  1.2220253232109896370417417439225704908830E0L,
+  1.2360594894780819419094519711090786987027E0L,
+  1.2490457723982544258299170772810901230778E0L,
+  1.2610933822524404193139408812473357720101E0L,
+  1.2722973952087173412961937498224804940684E0L,
+  1.2827408797442707473628852511364955306249E0L,
+  1.2924966677897852679030914214070816845853E0L,
+  1.3016288340091961438047858503666855921414E0L,
+  1.3101939350475556342564376891719053122733E0L,
+  1.3182420510168370498593302023271362531155E0L,
+  1.3258176636680324650592392104284756311844E0L,
+  1.3329603993374458675538498697331558093700E0L,
+  1.3397056595989995393283037525895557411039E0L,
+  1.3460851583802539310489409282517796256512E0L,
+  1.3521273809209546571891479413898128509842E0L,
+  1.3578579772154994751124898859640585287459E0L,
+  1.3633001003596939542892985278250991189943E0L,
+  1.3684746984165928776366381936948529556191E0L,
+  1.3734007669450158608612719264449611486510E0L,
+  1.3780955681325110444536609641291551522494E0L,
+  1.3825748214901258580599674177685685125566E0L,
+  1.3868528702577214543289381097042486034883E0L,
+  1.3909428270024183486427686943836432060856E0L,
+  1.3948567013423687823948122092044222644895E0L,
+  1.3986055122719575950126700816114282335732E0L,
+  1.4021993871854670105330304794336492676944E0L,
+  1.4056476493802697809521934019958079881002E0L,
+  1.4089588955564736949699075250792569287156E0L,
+  1.4121410646084952153676136718584891599630E0L,
+  1.4152014988178669079462550975833894394929E0L,
+  1.4181469983996314594038603039700989523716E0L,
+  1.4209838702219992566633046424614466661176E0L,
+  1.4237179714064941189018190466107297503086E0L,
+  1.4263547484202526397918060597281265695725E0L,
+  1.4288992721907326964184700745371983590908E0L,
+  1.4313562697035588982240194668401779312122E0L,
+  1.4337301524847089866404719096698873648610E0L,
+  1.4360250423171655234964275337155008780675E0L,
+  1.4382447944982225979614042479354815855386E0L,
+  1.4403930189057632173997301031392126865694E0L,
+  1.4424730991091018200252920599377292525125E0L,
+  1.4444882097316563655148453598508037025938E0L,
+  1.4464413322481351841999668424758804165254E0L,
+  1.4483352693775551917970437843145232637695E0L,
+  1.4501726582147939000905940595923466567576E0L,
+  1.4519559822271314199339700039142990228105E0L,
+  1.4536875822280323362423034480994649820285E0L,
+  1.4553696664279718992423082296859928222270E0L,
+  1.4570043196511885530074841089245667532358E0L,
+  1.4585935117976422128825857356750737658039E0L,
+  1.4601391056210009726721818194296893361233E0L,
+  1.4616428638860188872060496086383008594310E0L,
+  1.4631064559620759326975975316301202111560E0L,
+  1.4645314639038178118428450961503371619177E0L,
+  1.4659193880646627234129855241049975398470E0L,
+  1.4672716522843522691530527207287398276197E0L,
+  1.4685896086876430842559640450619880951144E0L,
+  1.4698745421276027686510391411132998919794E0L,
+  1.4711276743037345918528755717617308518553E0L,
+  1.4723501675822635384916444186631899205983E0L,
+  1.4735431285433308455179928682541563973416E0L, /* arctan(10.25) */
+  1.5707963267948966192313216916397514420986E0L  /* pi/2 */
+};
+
+
+/* arctan t = t + t^3 p(t^2) / q(t^2)
+   |t| <= 0.09375
+   peak relative error 5.3e-37 */
+
+static const long double
+  p0 = -4.283708356338736809269381409828726405572E1L,
+  p1 = -8.636132499244548540964557273544599863825E1L,
+  p2 = -5.713554848244551350855604111031839613216E1L,
+  p3 = -1.371405711877433266573835355036413750118E1L,
+  p4 = -8.638214309119210906997318946650189640184E-1L,
+  q0 = 1.285112506901621042780814422948906537959E2L,
+  q1 = 3.361907253914337187957855834229672347089E2L,
+  q2 = 3.180448303864130128268191635189365331680E2L,
+  q3 = 1.307244136980865800160844625025280344686E2L,
+  q4 = 2.173623741810414221251136181221172551416E1L;
+  /* q5 = 1.000000000000000000000000000000000000000E0 */
+
+
+long double
+atanl (long double x)
+{
+  int k, sign;
+  long double t, u, p, q;
+
+  sign = (x < 0) ? -1 : 1;
+
+  /* Check for zero or NaN.  */
+  if (x != x || x == 0.0)
+    return x + x;
+
+  /* Check for infinity.  */
+  if (x + x == x)
+    return sign * atantbl[83];
+
+  x *= sign;
+
+  if (x >= 10.25) /* 10.25 */
+    {
+      k = 83;
+      t = -1.0/x;
+    }
+  else
+    {
+      /* Index of nearest table element.
+        Roundoff to integer is asymmetrical to avoid cancellation when t < 0
+         (cf. fdlibm). */
+      k = 8.0 * x + 0.25;
+      u = 0.125 * k;
+      /* Small arctan argument.  */
+      t = (x - u) / (1.0 + x * u);
+    }
+
+  /* Arctan of small argument t.  */
+  u = t * t;
+  p =     ((((p4 * u) + p3) * u + p2) * u + p1) * u + p0;
+  q = ((((u + q4) * u + q3) * u + q2) * u + q1) * u + q0;
+  u = t * u * p / q  +  t;
+
+  /* arctan x = arctan u  +  arctan t */
+  return sign * (atantbl[k] + u);
+}
diff --git a/lib/ceill.c b/lib/ceill.c
new file mode 100644 (file)
index 0000000..94ffd6f
--- /dev/null
@@ -0,0 +1,45 @@
+/* Emulation for ceill.
+   Contributed by Paolo Bonzini
+
+   Copyright 2002, 2003 Free Software Foundation, Inc.
+
+   This file is part of gnulib.
+
+   gnulib is free software; you can redistribute it and/or modify it
+   under the terms of the GNU Lesser General Public License as published
+   by the Free Software Foundation; either version 2.1, or (at your option)
+   any later version.
+
+   gnulib is distributed in the hope that it will be useful, but WITHOUT
+   ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+   FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
+   License for more details.
+
+   You should have received a copy of the GNU Lesser General Public License
+   along with gnulib; see the file COPYING.LIB.  If not, write to the Free
+   Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307,
+   USA.
+ */
+
+#include <float.h>
+
+#include "mathl.h"
+
+/* To compute the integer part of X, sum a big enough
+   integer so that the precision of the floating point
+   number is exactly 1.  */
+
+long double
+ceill(long double x)
+{
+  long double y;
+  if (x < 0.0L)
+    y = -(1.0L / LDBL_EPSILON - x - 1.0 / LDBL_EPSILON);
+  else
+    y = 1.0L / LDBL_EPSILON + x - 1.0 / LDBL_EPSILON;
+
+  if (y < x)
+    return y + 1.0L;
+  else
+    return y;
+}
diff --git a/lib/cosl.c b/lib/cosl.c
new file mode 100644 (file)
index 0000000..f019eee
--- /dev/null
@@ -0,0 +1,97 @@
+/* s_cosl.c -- long double version of s_sin.c.
+ * Conversion to long double by Jakub Jelinek, jj@ultra.linux.cz.
+ */
+
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+
+/* sinl(x)
+ * Return sine function of x.
+ *
+ * kernel function:
+ *     __kernel_sinl           ... sine function on [-pi/4,pi/4]
+ *     __kernel_cosl           ... cose function on [-pi/4,pi/4]
+ *     __ieee754_rem_pio2l     ... argument reduction routine
+ *
+ * Method.
+ *      Let S,C and T denote the sin, cos and tan respectively on
+ *     [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2
+ *     in [-pi/4 , +pi/4], and let n = k mod 4.
+ *     We have
+ *
+ *          n        sin(x)      cos(x)        tan(x)
+ *     ----------------------------------------------------------
+ *         0          S           C             T
+ *         1          C          -S            -1/T
+ *         2         -S          -C             T
+ *         3         -C           S            -1/T
+ *     ----------------------------------------------------------
+ *
+ * Special cases:
+ *      Let trig be any of sin, cos, or tan.
+ *      trig(+-INF)  is NaN, with signals;
+ *      trig(NaN)    is that NaN;
+ *
+ * Accuracy:
+ *     TRIG(x) returns trig(x) nearly rounded
+ */
+
+#include <math.h>
+
+#include "mathl.h"
+
+#include "trigl.h"
+#ifdef HAVE_SINL
+#include "trigl.c"
+#include "sincosl.c"
+#endif
+
+long double cosl(long double x)
+{
+       long double y[2],z=0.0L;
+       int n;
+
+    /* |x| ~< pi/4 */
+        if(x >= -0.7853981633974483096156608458198757210492 &&
+           x <= 0.7853981633974483096156608458198757210492)
+          return kernel_cosl(x, z);
+
+    /* sinl(Inf or NaN) is NaN, sinl(0) is 0 */
+        else if ((x + x == x && x != 0.0) || x != x)
+          return x-x;           /* NaN */
+
+    /* argument reduction needed */
+       else {
+           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);
+           }
+       }
+}
+
+#if 0
+int
+main ()
+{
+  printf ("%.16Lg\n", cosl(0.7853981633974483096156608458198757210492));
+  printf ("%.16Lg\n", cosl(0.7853981633974483096156608458198757210492 *29));
+  printf ("%.16Lg\n", cosl(0.7853981633974483096156608458198757210492 *2));
+  printf ("%.16Lg\n", cosl(0.7853981633974483096156608458198757210492 *30));
+  printf ("%.16Lg\n", cosl(0.7853981633974483096156608458198757210492 *4));
+  printf ("%.16Lg\n", cosl(0.7853981633974483096156608458198757210492 *32));
+  printf ("%.16Lg\n", cosl(0.7853981633974483096156608458198757210492 *2/3));
+  printf ("%.16Lg\n", cosl(0.7853981633974483096156608458198757210492 *4/3));
+}
+#endif
diff --git a/lib/expl.c b/lib/expl.c
new file mode 100644 (file)
index 0000000..122e189
--- /dev/null
@@ -0,0 +1,136 @@
+/* Emulation for expl.
+   Contributed by Paolo Bonzini
+
+   Copyright 2002, 2003 Free Software Foundation, Inc.
+
+   This file is part of gnulib.
+
+   gnulib is free software; you can redistribute it and/or modify it
+   under the terms of the GNU Lesser General Public License as published
+   by the Free Software Foundation; either version 2.1, or (at your option)
+   any later version.
+
+   gnulib is distributed in the hope that it will be useful, but WITHOUT
+   ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+   FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
+   License for more details.
+
+   You should have received a copy of the GNU Lesser General Public License
+   along with gnulib; see the file COPYING.LIB.  If not, write to the Free
+   Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307,
+   USA.
+ */
+
+#include <float.h>
+#include <math.h>
+
+#include "mathl.h"
+
+static const long double C[] = {
+/* Chebyshev polynom coeficients for (exp(x)-1)/x */
+#define P1 C[0]
+#define P2 C[1]
+#define P3 C[2]
+#define P4 C[3]
+#define P5 C[4]
+#define P6 C[5]
+ 0.5L,
+ 1.66666666666666666666666666666666683E-01L,
+ 4.16666666666666666666654902320001674E-02L,
+ 8.33333333333333333333314659767198461E-03L,
+ 1.38888888889899438565058018857254025E-03L,
+ 1.98412698413981650382436541785404286E-04L,
+
+/* Smallest integer x for which e^x overflows.  */
+#define himark C[6]
+ 11356.523406294143949491931077970765L,
+
+/* Largest integer x for which e^x underflows.  */
+#define lomark C[7]
+-11433.4627433362978788372438434526231L,
+
+/* very small number */
+#define TINY C[8]
+ 1.0e-4900L,
+
+/* 2^16383 */
+#define TWO16383 C[9]
+ 5.94865747678615882542879663314003565E+4931L};
+
+long double
+expl (long double x)
+{
+  /* Check for usual case.  */
+  if (x < himark && x > lomark)
+    {
+      int exponent;
+      long double t, x22;
+      int k = 1;
+      long double result = 1.0;
+
+      /* Compute an integer power of e with a granularity of 0.125.  */
+      exponent = (int) floorl (x * 8.0L);
+      x -= exponent / 8.0L;
+
+      if (x > 0.0625)
+       {
+         exponent++;
+         x -= 0.125L;
+       }
+
+      if (exponent < 0)
+        {
+          t = 0.8824969025845954028648921432290507362220L; /* e^-0.25 */
+         exponent = -exponent;
+       }
+      else
+        t = 1.1331484530668263168290072278117938725655L; /* e^0.25 */
+
+      while (exponent)
+        {
+          if (exponent & k)
+            {
+              result *= t;
+              exponent ^= k;
+            }
+          t *= t;
+          k <<= 1;
+        }
+
+      /* Approximate (e^x - 1)/x, using a seventh-degree polynomial,
+        with maximum error in [-2^-16-2^-53,2^-16+2^-53]
+        less than 4.8e-39.  */
+      x22 = x + x*x*(P1+x*(P2+x*(P3+x*(P4+x*(P5+x*P6)))));
+
+      return result + result * x22;
+    }
+  /* Exceptional cases:  */
+  else if (x < himark)
+    {
+      if (x + x == x)
+       /* e^-inf == 0, with no error.  */
+       return 0;
+      else
+       /* Underflow */
+       return TINY * TINY;
+    }
+  else
+    /* Return x, if x is a NaN or Inf; or overflow, otherwise.  */
+    return TWO16383*x;
+}
+
+#if 0
+int
+main ()
+{
+  printf ("%.16Lg\n", expl(1.0L));
+  printf ("%.16Lg\n", expl(-1.0L));
+  printf ("%.16Lg\n", expl(2.0L));
+  printf ("%.16Lg\n", expl(4.0L));
+  printf ("%.16Lg\n", expl(-2.0L));
+  printf ("%.16Lg\n", expl(-4.0L));
+  printf ("%.16Lg\n", expl(0.0625L));
+  printf ("%.16Lg\n", expl(0.3L));
+  printf ("%.16Lg\n", expl(0.6L));
+}
+#endif
diff --git a/lib/floorl.c b/lib/floorl.c
new file mode 100644 (file)
index 0000000..0f1f0cb
--- /dev/null
@@ -0,0 +1,45 @@
+/* Emulation for floorl.
+   Contributed by Paolo Bonzini
+
+   Copyright 2002, 2003 Free Software Foundation, Inc.
+
+   This file is part of gnulib.
+
+   gnulib is free software; you can redistribute it and/or modify it
+   under the terms of the GNU Lesser General Public License as published
+   by the Free Software Foundation; either version 2.1, or (at your option)
+   any later version.
+
+   gnulib is distributed in the hope that it will be useful, but WITHOUT
+   ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+   FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
+   License for more details.
+
+   You should have received a copy of the GNU Lesser General Public License
+   along with gnulib; see the file COPYING.LIB.  If not, write to the Free
+   Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307,
+   USA.
+ */
+
+#include <float.h>
+
+#include "mathl.h"
+
+/* To compute the integer part of X, sum a big enough
+   integer so that the precision of the floating point
+   number is exactly 1.  */
+
+long double
+floorl(long double x)
+{
+  long double y;
+  if (x < 0.0L)
+    y = -(1.0L / LDBL_EPSILON - x - 1.0 / LDBL_EPSILON);
+  else
+    y = 1.0L / LDBL_EPSILON + x - 1.0 / LDBL_EPSILON;
+
+  if (y > x)
+    return y - 1.0L;
+  else
+    return y;
+}
diff --git a/lib/frexpl.c b/lib/frexpl.c
new file mode 100644 (file)
index 0000000..cdc3ad8
--- /dev/null
@@ -0,0 +1,103 @@
+/* Emulation for frexpl.
+   Contributed by Paolo Bonzini
+
+   Copyright 2002, 2003 Free Software Foundation, Inc.
+
+   This file is part of gnulib.
+
+   gnulib is free software; you can redistribute it and/or modify it
+   under the terms of the GNU Lesser General Public License as published
+   by the Free Software Foundation; either version 2.1, or (at your option)
+   any later version.
+
+   gnulib is distributed in the hope that it will be useful, but WITHOUT
+   ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+   FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
+   License for more details.
+
+   You should have received a copy of the GNU Lesser General Public License
+   along with gnulib; see the file COPYING.LIB.  If not, write to the Free
+   Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307,
+   USA.
+ */
+
+#include <float.h>
+#include <math.h>
+
+#include "mathl.h"
+
+/* Binary search.  Quite inefficient but portable. */
+long double
+frexpl(long double x, int *exp)
+{
+  long double exponents[20], *next;
+  int exponent, bit;
+
+  /* Check for zero, nan and infinity. */
+  if (x != x || x + x == x )
+    {
+      *exp = 0;
+      return x;
+    }
+
+  if (x < 0)
+    return -frexpl(-x, exp);
+
+  exponent = 0;
+  if (x > 1.0)
+    {
+      for (next = exponents, exponents[0] = 2.0L, bit = 1;
+          *next <= x + x;
+          bit <<= 1, next[1] = next[0] * next[0], next++);
+
+      for (; next >= exponents; bit >>= 1, next--)
+       if (x + x >= *next)
+         {
+           x /= *next;
+           exponent |= bit;
+         }
+
+    }
+
+  else if (x < 0.5)
+    {
+      for (next = exponents, exponents[0] = 0.5L, bit = 1;
+          *next > x;
+          bit <<= 1, next[1] = next[0] * next[0], next++);
+
+      for (; next >= exponents; bit >>= 1, next--)
+       if (x < *next)
+         {
+           x /= *next;
+           exponent |= bit;
+         }
+
+      exponent = -exponent;
+    }
+
+  *exp = exponent;
+  return x;
+}
+
+#if 0
+int
+main()
+{
+  long double x;
+  int y;
+  x = frexpl(0.0L / 0.0L, &y); printf ("%.6Lg %d\n", x, y);
+  x = frexpl(1.0L / 0.0L, &y); printf ("%.6Lg %d\n", x, y);
+  x = frexpl(-1.0L / 0.0L, &y); printf ("%.6Lg %d\n", x, y);
+  x = frexpl(0.5L, &y); printf ("%.6Lg %d\n", x, y);
+  x = frexpl(0.75L, &y); printf ("%.6Lg %d\n", x, y);
+  x = frexpl(3.6L, &y); printf ("%.6Lg %d\n", x, y);
+  x = frexpl(17.8L, &y); printf ("%.6Lg %d\n", x, y);
+  x = frexpl(8.0L, &y); printf ("%.6Lg %d\n", x, y);
+  x = frexpl(0.3L, &y); printf ("%.6Lg %d\n", x, y);
+  x = frexpl(0.49L, &y); printf ("%.6Lg %d\n", x, y);
+  x = frexpl(0.049L, &y); printf ("%.6Lg %d\n", x, y);
+  x = frexpl(0.0245L, &y); printf ("%.6Lg %d\n", x, y);
+  x = frexpl(0.0625L, &y); printf ("%.6Lg %d\n", x, y);
+}
+#endif
+
diff --git a/lib/ldexpl.c b/lib/ldexpl.c
new file mode 100644 (file)
index 0000000..cd9bd21
--- /dev/null
@@ -0,0 +1,63 @@
+/* Emulation for ldexpl.
+   Contributed by Paolo Bonzini
+
+   Copyright 2002, 2003 Free Software Foundation, Inc.
+
+   This file is part of gnulib.
+
+   gnulib is free software; you can redistribute it and/or modify it
+   under the terms of the GNU Lesser General Public License as published
+   by the Free Software Foundation; either version 2.1, or (at your option)
+   any later version.
+
+   gnulib is distributed in the hope that it will be useful, but WITHOUT
+   ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+   FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
+   License for more details.
+
+   You should have received a copy of the GNU Lesser General Public License
+   along with gnulib; see the file COPYING.LIB.  If not, write to the Free
+   Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307,
+   USA.
+ */
+
+#include <float.h>
+#include <math.h>
+
+#include "mathl.h"
+
+long double
+ldexpl(long double x, int exp)
+{
+  long double factor;
+  int bit;
+
+  /* Check for zero, nan and infinity. */
+  if (x != x || x + x == x )
+    return x;
+
+  if (exp < 0)
+    {
+      exp = -exp;
+      factor = 0.5L;
+    }
+  else
+    factor = 2.0L;
+
+  for (bit = 1; bit <= exp; bit <<= 1, factor *= factor)
+    if (exp & bit)
+      x *= factor;
+
+  return x;
+}
+
+#if 0
+int
+main()
+{
+  long double x;
+  int y;
+  for (y = 0; y < 29; y++)
+    printf ("%5d %.16Lg %.16Lg\n", y, ldexpl(0.8L, y), ldexpl(0.8L, -y) * ldexpl(0.8L, y));
+}
+#endif
diff --git a/lib/logl.c b/lib/logl.c
new file mode 100644 (file)
index 0000000..9048e45
--- /dev/null
@@ -0,0 +1,238 @@
+/*                                                     logll.c
+ *
+ * Natural logarithm for 128-bit long double precision.
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * long double x, y, logl();
+ *
+ * y = logl( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns the base e (2.718...) logarithm of x.
+ *
+ * The argument is separated into its exponent and fractional
+ * parts.  Use of a lookup table increases the speed of the routine.
+ * The program uses logarithms tabulated at intervals of 1/128 to
+ * cover the domain from approximately 0.7 to 1.4.
+ *
+ * On the interval [-1/128, +1/128] the logarithm of 1+x is approximated by
+ *     log(1+x) = x - 0.5 x^2 + x^3 P(x) .
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE   0.875, 1.125   100000      1.2e-34    4.1e-35
+ *    IEEE   0.125, 8       100000      1.2e-34    4.1e-35
+ *
+ *
+ * WARNING:
+ *
+ * This program uses integer operations on bit fields of floating-point
+ * numbers.  It does not work with data structures other than the
+ * structure assumed.
+ *
+ */
+
+#include <math.h>
+
+#include "mathl.h"
+
+/* Copyright 2001 by Stephen L. Moshier <moshier@na-net.ornl.gov> */
+
+/* log(1+x) = x - .5 x^2 + x^3 l(x)
+   -.0078125 <= x <= +.0078125
+   peak relative error 1.2e-37 */
+static const long double
+l3 =   3.333333333333333333333333333333336096926E-1L,
+l4 =  -2.499999999999999999999999999486853077002E-1L,
+l5 =   1.999999999999999999999999998515277861905E-1L,
+l6 =  -1.666666666666666666666798448356171665678E-1L,
+l7 =   1.428571428571428571428808945895490721564E-1L,
+l8 =  -1.249999999999999987884655626377588149000E-1L,
+l9 =   1.111111111111111093947834982832456459186E-1L,
+l10 = -1.000000000000532974938900317952530453248E-1L,
+l11 =  9.090909090915566247008015301349979892689E-2L,
+l12 = -8.333333211818065121250921925397567745734E-2L,
+l13 =  7.692307559897661630807048686258659316091E-2L,
+l14 = -7.144242754190814657241902218399056829264E-2L,
+l15 =  6.668057591071739754844678883223432347481E-2L;
+
+/* Lookup table of ln(t) - (t-1)
+    t = 0.5 + (k+26)/128)
+    k = 0, ..., 91   */
+static const long double logtbl[92] = {
+-5.5345593589352099112142921677820359632418E-2L,
+-5.2108257402767124761784665198737642086148E-2L,
+-4.8991686870576856279407775480686721935120E-2L,
+-4.5993270766361228596215288742353061431071E-2L,
+-4.3110481649613269682442058976885699556950E-2L,
+-4.0340872319076331310838085093194799765520E-2L,
+-3.7682072451780927439219005993827431503510E-2L,
+-3.5131785416234343803903228503274262719586E-2L,
+-3.2687785249045246292687241862699949178831E-2L,
+-3.0347913785027239068190798397055267411813E-2L,
+-2.8110077931525797884641940838507561326298E-2L,
+-2.5972247078357715036426583294246819637618E-2L,
+-2.3932450635346084858612873953407168217307E-2L,
+-2.1988775689981395152022535153795155900240E-2L,
+-2.0139364778244501615441044267387667496733E-2L,
+-1.8382413762093794819267536615342902718324E-2L,
+-1.6716169807550022358923589720001638093023E-2L,
+-1.5138929457710992616226033183958974965355E-2L,
+-1.3649036795397472900424896523305726435029E-2L,
+-1.2244881690473465543308397998034325468152E-2L,
+-1.0924898127200937840689817557742469105693E-2L,
+-9.6875626072830301572839422532631079809328E-3L,
+-8.5313926245226231463436209313499745894157E-3L,
+-7.4549452072765973384933565912143044991706E-3L,
+-6.4568155251217050991200599386801665681310E-3L,
+-5.5356355563671005131126851708522185605193E-3L,
+-4.6900728132525199028885749289712348829878E-3L,
+-3.9188291218610470766469347968659624282519E-3L,
+-3.2206394539524058873423550293617843896540E-3L,
+-2.5942708080877805657374888909297113032132E-3L,
+-2.0385211375711716729239156839929281289086E-3L,
+-1.5522183228760777967376942769773768850872E-3L,
+-1.1342191863606077520036253234446621373191E-3L,
+-7.8340854719967065861624024730268350459991E-4L,
+-4.9869831458030115699628274852562992756174E-4L,
+-2.7902661731604211834685052867305795169688E-4L,
+-1.2335696813916860754951146082826952093496E-4L,
+-3.0677461025892873184042490943581654591817E-5L,
+#define ZERO logtbl[38]
+ 0.0000000000000000000000000000000000000000E0L,
+-3.0359557945051052537099938863236321874198E-5L,
+-1.2081346403474584914595395755316412213151E-4L,
+-2.7044071846562177120083903771008342059094E-4L,
+-4.7834133324631162897179240322783590830326E-4L,
+-7.4363569786340080624467487620270965403695E-4L,
+-1.0654639687057968333207323853366578860679E-3L,
+-1.4429854811877171341298062134712230604279E-3L,
+-1.8753781835651574193938679595797367137975E-3L,
+-2.3618380914922506054347222273705859653658E-3L,
+-2.9015787624124743013946600163375853631299E-3L,
+-3.4938307889254087318399313316921940859043E-3L,
+-4.1378413103128673800485306215154712148146E-3L,
+-4.8328735414488877044289435125365629849599E-3L,
+-5.5782063183564351739381962360253116934243E-3L,
+-6.3731336597098858051938306767880719015261E-3L,
+-7.2169643436165454612058905294782949315193E-3L,
+-8.1090214990427641365934846191367315083867E-3L,
+-9.0486422112807274112838713105168375482480E-3L,
+-1.0035177140880864314674126398350812606841E-2L,
+-1.1067990155502102718064936259435676477423E-2L,
+-1.2146457974158024928196575103115488672416E-2L,
+-1.3269969823361415906628825374158424754308E-2L,
+-1.4437927104692837124388550722759686270765E-2L,
+-1.5649743073340777659901053944852735064621E-2L,
+-1.6904842527181702880599758489058031645317E-2L,
+-1.8202661505988007336096407340750378994209E-2L,
+-1.9542647000370545390701192438691126552961E-2L,
+-2.0924256670080119637427928803038530924742E-2L,
+-2.2346958571309108496179613803760727786257E-2L,
+-2.3810230892650362330447187267648486279460E-2L,
+-2.5313561699385640380910474255652501521033E-2L,
+-2.6856448685790244233704909690165496625399E-2L,
+-2.8438398935154170008519274953860128449036E-2L,
+-3.0058928687233090922411781058956589863039E-2L,
+-3.1717563112854831855692484086486099896614E-2L,
+-3.3413836095418743219397234253475252001090E-2L,
+-3.5147290019036555862676702093393332533702E-2L,
+-3.6917475563073933027920505457688955423688E-2L,
+-3.8723951502862058660874073462456610731178E-2L,
+-4.0566284516358241168330505467000838017425E-2L,
+-4.2444048996543693813649967076598766917965E-2L,
+-4.4356826869355401653098777649745233339196E-2L,
+-4.6304207416957323121106944474331029996141E-2L,
+-4.8285787106164123613318093945035804818364E-2L,
+-5.0301169421838218987124461766244507342648E-2L,
+-5.2349964705088137924875459464622098310997E-2L,
+-5.4431789996103111613753440311680967840214E-2L,
+-5.6546268881465384189752786409400404404794E-2L,
+-5.8693031345788023909329239565012647817664E-2L,
+-6.0871713627532018185577188079210189048340E-2L,
+-6.3081958078862169742820420185833800925568E-2L,
+-6.5323413029406789694910800219643791556918E-2L,
+-6.7595732653791419081537811574227049288168E-2L
+};
+
+/* ln(2) = ln2a + ln2b with extended precision. */
+static const long double
+  ln2a = 6.93145751953125e-1L,
+  ln2b = 1.4286068203094172321214581765680755001344E-6L;
+
+long double
+logl(long double x)
+{
+  long double z, y, w, u, t;
+  unsigned int m;
+  int k, e;
+
+  /* Check for IEEE special cases.  */
+
+  /* log(0) = -infinity. */
+  if (x == 0.0L)
+    return -0.5L / 0.0L;
+
+  /* log ( x < 0 ) = NaN */
+  if (x < 0.0L)
+    return (x - x) / (x - x);
+
+  /* log (infinity or NaN) */
+  if (x + x == x || x != x)
+    return x + x;
+
+  /* Extract exponent and reduce domain to 0.703125 <= u < 1.40625  */
+  x = frexpl(x, &e);
+  if (x < 0.703125L)
+    {
+      x += x;
+      e--;
+    }
+
+  /* On this interval the table is not used due to cancellation error.  */
+  if ((x <= 1.0078125L) && (x >= 0.9921875L))
+    {
+      k = 64;
+      t = 1.0L;
+      z = x - 1.0L;
+    }
+  else
+    {
+      k = floorl((x - 0.5L) * 128.0L);
+      t = 0.5L + k / 128.0L;
+      z = (x - t) / t;
+    }
+
+  /* Series expansion of log(1+z).  */
+  w = z * z;
+  y = ((((((((((((l15 * z
+                 + l14) * z
+                + l13) * z
+               + l12) * z
+              + l11) * z
+             + l10) * z
+            + l9) * z
+           + l8) * z
+          + l7) * z
+         + l6) * z
+        + l5) * z
+       + l4) * z
+       + l3) * z * w;
+  y -= 0.5 * w;
+  y += e * ln2b;  /* Base 2 exponent offset times ln(2).  */
+  y += z;
+  y += logtbl[k-26]; /* log(t) - (t-1) */
+  y += (t - 1.0L);
+  y += e * ln2a;
+  return y;
+}
diff --git a/lib/mathl.h b/lib/mathl.h
new file mode 100644 (file)
index 0000000..841b1d6
--- /dev/null
@@ -0,0 +1,41 @@
+/* Declaration for long double functions.
+   Contributed by Paolo Bonzini
+
+   Copyright 2002, 2003 Free Software Foundation, Inc.
+
+   This file is part of gnulib.
+
+   gnulib is free software; you can redistribute it and/or modify it
+   under the terms of the GNU Lesser General Public License as published
+   by the Free Software Foundation; either version 2.1, or (at your option)
+   any later version.
+
+   gnulib is distributed in the hope that it will be useful, but WITHOUT
+   ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+   FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
+   License for more details.
+
+   You should have received a copy of the GNU Lesser General Public License
+   along with gnulib; see the file COPYING.LIB.  If not, write to the Free
+   Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307,
+   USA.
+ */
+
+#ifndef GNULIB_MATHL_H
+#define GNULIB_MATHL_H
+
+extern long double acosl (long double x);
+extern long double asinl (long double x);
+extern long double atanl (long double x);
+extern long double ceill (long double x);
+extern long double cosl (long double x);
+extern long double expl (long double x);
+extern long double floorl (long double x);
+extern long double frexpl (long double x, int *exp);
+extern long double ldexpl (long double x, int exp);
+extern long double logl (long double x);
+extern long double sinl (long double x);
+extern long double sqrtl (long double x);
+extern long double tanl (long double x);
+
+#endif
diff --git a/lib/sincosl.c b/lib/sincosl.c
new file mode 100644 (file)
index 0000000..9fa8be9
--- /dev/null
@@ -0,0 +1,897 @@
+/* Quad-precision floating point trigonometric functions on <-pi/4,pi/4>.
+   Copyright (C) 1999 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+   Contributed by Jakub Jelinek <jj@ultra.linux.cz>
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+
+   The GNU C Library is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, write to the Free
+   Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
+   02111-1307 USA.  */
+
+#include <math.h>
+#include <float.h>
+
+#include "mathl.h"
+
+static const long double sin_c[] = {
+#define ONE sin_c[0]
+  1.00000000000000000000000000000000000E+00L,  /* 3fff0000000000000000000000000000 */
+
+/* cos x ~ ONE + x^2 ( SCOS1 + SCOS2 * x^2 + ... + SCOS4 * x^6 + SCOS5 * x^8 )
+   x in <0,1/256>  */
+#define SCOS1 sin_c[1]
+#define SCOS2 sin_c[2]
+#define SCOS3 sin_c[3]
+#define SCOS4 sin_c[4]
+#define SCOS5 sin_c[5]
+  -5.00000000000000000000000000000000000E-01L, /* bffe0000000000000000000000000000 */
+  4.16666666666666666666666666556146073E-02L,  /* 3ffa5555555555555555555555395023 */
+  -1.38888888888888888888309442601939728E-03L, /* bff56c16c16c16c16c16a566e42c0375 */
+  2.48015873015862382987049502531095061E-05L,  /* 3fefa01a01a019ee02dcf7da2d6d5444 */
+  -2.75573112601362126593516899592158083E-07L, /* bfe927e4f5dce637cb0b54908754bde0 */
+
+/* sin x ~ ONE * x + x^3 ( SIN1 + SIN2 * x^2 + ... + SIN7 * x^12 + SIN8 * x^14 )
+   x in <0,0.1484375>  */
+#define SIN1 sin_c[6]
+#define SIN2 sin_c[7]
+#define SIN3 sin_c[8]
+#define SIN4 sin_c[9]
+#define SIN5 sin_c[10]
+#define SIN6 sin_c[11]
+#define SIN7 sin_c[12]
+#define SIN8 sin_c[13]
+  -1.66666666666666666666666666666666538e-01L, /* bffc5555555555555555555555555550 */
+  8.33333333333333333333333333307532934e-03L,  /* 3ff811111111111111111111110e7340 */
+  -1.98412698412698412698412534478712057e-04L, /* bff2a01a01a01a01a01a019e7a626296 */
+  2.75573192239858906520896496653095890e-06L,  /* 3fec71de3a556c7338fa38527474b8f5 */
+  -2.50521083854417116999224301266655662e-08L, /* bfe5ae64567f544e16c7de65c2ea551f */
+  1.60590438367608957516841576404938118e-10L,  /* 3fde6124613a811480538a9a41957115 */
+  -7.64716343504264506714019494041582610e-13L, /* bfd6ae7f3d5aef30c7bc660b060ef365 */
+  2.81068754939739570236322404393398135e-15L,  /* 3fce9510115aabf87aceb2022a9a9180 */
+
+/* sin x ~ ONE * x + x^3 ( SSIN1 + SSIN2 * x^2 + ... + SSIN4 * x^6 + SSIN5 * x^8 )
+   x in <0,1/256>  */
+#define SSIN1 sin_c[14]
+#define SSIN2 sin_c[15]
+#define SSIN3 sin_c[16]
+#define SSIN4 sin_c[17]
+#define SSIN5 sin_c[18]
+  -1.66666666666666666666666666666666659E-01L, /* bffc5555555555555555555555555555 */
+  8.33333333333333333333333333146298442E-03L,  /* 3ff81111111111111111111110fe195d */
+  -1.98412698412698412697726277416810661E-04L, /* bff2a01a01a01a01a019e7121e080d88 */
+  2.75573192239848624174178393552189149E-06L,  /* 3fec71de3a556c640c6aaa51aa02ab41 */
+  -2.50521016467996193495359189395805639E-08L, /* bfe5ae644ee90c47dc71839de75b2787 */
+};
+
+static const long double cos_c[] = {
+/* cos x ~ ONE + x^2 ( COS1 + COS2 * x^2 + ... + COS7 * x^12 + COS8 * x^14 )
+   x in <0,0.1484375>  */
+#define COS1 cos_c[0]
+#define COS2 cos_c[1]
+#define COS3 cos_c[2]
+#define COS4 cos_c[3]
+#define COS5 cos_c[4]
+#define COS6 cos_c[5]
+#define COS7 cos_c[6]
+#define COS8 cos_c[7]
+  -4.99999999999999999999999999999999759E-01L, /* bffdfffffffffffffffffffffffffffb */
+  4.16666666666666666666666666651287795E-02L,  /* 3ffa5555555555555555555555516f30 */
+  -1.38888888888888888888888742314300284E-03L, /* bff56c16c16c16c16c16c16a463dfd0d */
+  2.48015873015873015867694002851118210E-05L,  /* 3fefa01a01a01a01a0195cebe6f3d3a5 */
+  -2.75573192239858811636614709689300351E-07L, /* bfe927e4fb7789f5aa8142a22044b51f */
+  2.08767569877762248667431926878073669E-09L,  /* 3fe21eed8eff881d1e9262d7adff4373 */
+  -1.14707451049343817400420280514614892E-11L, /* bfda9397496922a9601ed3d4ca48944b */
+  4.77810092804389587579843296923533297E-14L,  /* 3fd2ae5f8197cbcdcaf7c3fb4523414c */
+
+};
+
+#define SINCOSL_COS_HI 0
+#define SINCOSL_COS_LO 1
+#define SINCOSL_SIN_HI 2
+#define SINCOSL_SIN_LO 3
+static const long double sincosl_table[];
+
+long double
+kernel_sinl (long double x, long double y, int iy)
+{
+  long double h, l, z, sin_l, cos_l_m1;
+  int index, sign;
+
+  sign = 1;
+  if (x < 0)
+    {
+      x = -x;
+      sign = -1;
+    }
+
+  if (x < 0.148375L)           /* |x| < 0.1484375 */
+    {
+      /* Argument is small enough to approximate it by a Chebyshev
+         polynomial of degree 17.  */
+      if (x < 0.000000000000000006938893903907228377647697925567626953125L)    /* |x| < 2^-57 */
+       if (!((int) x))
+         return x;             /* generate inexact */
+
+      z = x * x;
+      return x + (x * (z * (SIN1 + z * (SIN2 + z * (SIN3 + z * (SIN4 +
+                                                               z * (SIN5 +
+                                                                    z *
+                                                                    (SIN6 +
+                                                                     z *
+                                                                     (SIN7 +
+                                                                      z *
+                                                                      SIN8)))))))));
+    }
+  else
+    {
+      /* So that we don't have to use too large polynomial,  we find
+         l and h such that x = l + h,  where fabsl(l) <= 1.0/256 with 83
+         possible values for h.  We look up cosl(h) and sinl(h) in
+         pre-computed tables,  compute cosl(l) and sinl(l) using a
+         Chebyshev polynomial of degree 10(11) and compute
+         sinl(h+l) = sinl(h)cosl(l) + cosl(h)sinl(l).  */
+      x -= 0.1484375L;
+      index = (int) (x * 128L + 0.5L);
+      h = index / 128.0L;
+      if (iy)
+       l = y - (h - x);
+      else
+       l = x - h;
+
+      z = l * l;
+      sin_l =
+       l * (ONE +
+            z * (SSIN1 +
+                 z * (SSIN2 + z * (SSIN3 + z * (SSIN4 + z * SSIN5)))));
+      cos_l_m1 =
+       z * (SCOS1 + z * (SCOS2 + z * (SCOS3 + z * (SCOS4 + z * SCOS5))));
+
+      index *= 4;
+      z =
+       sincosl_table[index + SINCOSL_SIN_HI] +
+       (sincosl_table[index + SINCOSL_SIN_LO] +
+        (sincosl_table[index + SINCOSL_SIN_HI] * cos_l_m1) +
+        (sincosl_table[index + SINCOSL_COS_HI] * sin_l));
+      return z * sign;
+    }
+}
+
+long double
+kernel_cosl (long double x, long double y)
+{
+  long double h, l, z, sin_l, cos_l_m1;
+  int index;
+
+  if (x < 0)
+    x = -x;
+
+  if (x < 0.1484375L)          /* |x| < 0.1484375 */
+    {
+      /* Argument is small enough to approximate it by a Chebyshev
+         polynomial of degree 16.  */
+      if (x < 0.000000000000000006938893903907228377647697925567626953125L)    /* |x| < 2^-57 */
+       if (!((int) x))
+         return ONE;           /* generate inexact */
+      z = x * x;
+      return ONE + (z * (COS1 + z * (COS2 + z * (COS3 + z * (COS4 +
+                                                            z * (COS5 +
+                                                                 z * (COS6 +
+                                                                      z *
+                                                                      (COS7 +
+                                                                       z *
+                                                                       COS8))))))));
+    }
+  else
+    {
+      /* So that we don't have to use too large polynomial,  we find
+         l and h such that x = l + h,  where fabsl(l) <= 1.0/256 with 83
+         possible values for h.  We look up cosl(h) and sinl(h) in
+         pre-computed tables,  compute cosl(l) and sinl(l) using a
+         Chebyshev polynomial of degree 10(11) and compute
+         sinl(h+l) = sinl(h)cosl(l) + cosl(h)sinl(l).  */
+      x -= 0.1484375L;
+      index = (int) (x * 128L + 0.5L);
+      h = index / 128.0L;
+      l = y - (h - x);
+      z = l * l;
+      sin_l =
+       l * (ONE +
+            z * (SSIN1 +
+                 z * (SSIN2 + z * (SSIN3 + z * (SSIN4 + z * SSIN5)))));
+      cos_l_m1 =
+       z * (SCOS1 + z * (SCOS2 + z * (SCOS3 + z * (SCOS4 + z * SCOS5))));
+
+      index *= 4;
+      z = sincosl_table [index + SINCOSL_COS_HI]
+          + (sincosl_table [index + SINCOSL_COS_LO]
+             - (sincosl_table [index + SINCOSL_SIN_HI] * sin_l)
+             - (sincosl_table [index + SINCOSL_COS_HI] * cos_l_m1));
+      return z;
+    }
+}
+
+
+/* For 0.1484375 + n/128.0, n=0..82 this table contains
+   first 113 bits of cosine, then at least 113 additional
+   bits and the same for sine.
+   0.1484375+82.0/128.0 is the smallest number among above defined numbers
+   larger than pi/4.
+   Computed using gmp.
+ */
+
+const long double sincosl_table[] = {
+
+/* x =  1.48437500000000000000000000000000000e-01L 3ffc3000000000000000000000000000 */
+/* cos(x) = 0.fd2f5320e1b790209b4dda2f98f79caaa7b873aff1014b0fbc5243766d03cb006bc837c4358 */
+  9.89003367927322909016887196069562069e-01L,  /* 3ffefa5ea641c36f2041369bb45f31ef */
+  2.15663692029265697782289400027743703e-35L,  /* 3f8bcaaa7b873aff1014b0fbc5243767 */
+/* sin(x) = 0.25dc50bc95711d0d9787d108fd438cf5959ee0bfb7a1e36e8b1a112968f356657420e9cc9ea */
+  1.47892995873409608580026675734609314e-01L,  /* 3ffc2ee285e4ab88e86cbc3e8847ea1c */
+  9.74950446464233268291647449768590886e-36L,  /* 3f8a9eb2b3dc17f6f43c6dd16342252d */
+
+/* x = 1.56250000000000000000000000000000000e-01 3ffc4000000000000000000000000000 */
+/* cos(x) = 0.fce1a053e621438b6d60c76e8c45bf0a9dc71aa16f922acc10e95144ec796a249813c9cb649 */
+  9.87817783816471944100503034363211317e-01L,  /* 3ffef9c340a7cc428716dac18edd188b */
+  4.74271307836705897892468107620526395e-35L,  /* 3f8cf854ee38d50b7c915660874a8a27 */
+/* sin(x) = 0.27d66258bacd96a3eb335b365c87d59438c5142bb56a489e9b8db9d36234ffdebb6bdc22d8e */
+  1.55614992773556041209920643203516258e-01L,  /* 3ffc3eb312c5d66cb51f599ad9b2e43f */
+  -7.83989563419287980121718050629497270e-36L, /* bf8a4d78e75d7a8952b6ec2c8e48c594 */
+
+/* x = 1.64062500000000000000000000000000000e-01 3ffc5000000000000000000000000000 */
+/* cos(x) = 0.fc8ffa01ba6807417e05962b0d9fdf1fddb0cc4c07d22e19e08019bffa50a6c7acdb40307a3 */
+  9.86571908399497588757337407495308409e-01L,  /* 3ffef91ff40374d00e82fc0b2c561b40 */
+  -2.47327949936985362476252401212720725e-35L, /* bf8c070112799d9fc16e8f30fbff3200 */
+/* sin(x) = 0.29cfd49b8be4f665276cab01cbf0426934906c3dd105473b226e410b1450f62e53ff7c6cce1 */
+  1.63327491736612850846866172454354370e-01L,  /* 3ffc4e7ea4dc5f27b3293b65580e5f82 */
+  1.81380344301155485770367902300754350e-36L,  /* 3f88349a48361ee882a39d913720858a */
+
+/* x = 1.71875000000000000000000000000000000e-01 3ffc6000000000000000000000000000 */
+/* cos(x) = 0.fc3a6170f767ac735d63d99a9d439e1db5e59d3ef153a4265d5855850ed82b536bf361b80e3 */
+  9.85265817718213816204294709759578994e-01L,  /* 3ffef874c2e1eecf58e6bac7b3353a87 */
+  2.26568029505818066141517497778527952e-35L,  /* 3f8be1db5e59d3ef153a4265d5855851 */
+/* sin(x) = 0.2bc89f9f424de5485de7ce03b2514952b9faf5648c3244d4736feb95dbb9da49f3b58a9253b */
+  1.71030022031395019281347969239834331e-01L,  /* 3ffc5e44fcfa126f2a42ef3e701d928a */
+  7.01395875187487608875416030203241317e-36L,  /* 3f8a2a573f5eac9186489a8e6dfd72bb */
+
+/* x = 1.79687500000000000000000000000000000e-01 3ffc7000000000000000000000000000 */
+/* cos(x) = 0.fbe0d7f7fef11e70aa43b8abf4f6a457cea20c8f3f676b47781f9821bbe9ce04b3c7b981c0b */
+  9.83899591489663972178309351416487245e-01L,  /* 3ffef7c1afeffde23ce154877157e9ed */
+  2.73414318948066207810486330723761265e-35L,  /* 3f8c22be75106479fb3b5a3bc0fcc10e */
+/* sin(x) = 0.2dc0bb80b49a97ffb34e8dd1f8db9df7af47ed2dcf58b12c8e7827e048cae929da02c04ecac */
+  1.78722113535153659375356241864180724e-01L,  /* 3ffc6e05dc05a4d4bffd9a746e8fc6dd */
+  -1.52906926517265103202547561260594148e-36L, /* bf8804285c09691853a769b8c3ec0fdc */
+
+/* x = 1.87500000000000000000000000000000000e-01 3ffc8000000000000000000000000000 */
+/* cos(x) = 0.fb835efcf670dd2ce6fe7924697eea13ea358867e9cdb3899b783f4f9f43aa5626e8b67b3bc */
+  9.82473313101255257487327683243622495e-01L,  /* 3ffef706bdf9ece1ba59cdfcf248d2fe */
+  -1.64924358891557584625463868014230342e-35L, /* bf8b5ec15ca779816324c766487c0b06 */
+/* sin(x) = 0.2fb8205f75e56a2b56a1c4792f856258769af396e0189ef72c05e4df59a6b00e4b44a6ea515 */
+  1.86403296762269884552379983103205261e-01L,  /* 3ffc7dc102fbaf2b515ab50e23c97c2b */
+  1.76460304806826780010586715975331753e-36L,  /* 3f882c3b4d79cb700c4f7b9602f26fad */
+
+/* x = 1.95312500000000000000000000000000000e-01 3ffc9000000000000000000000000000 */
+/* cos(x) = 0.fb21f7f5c156696b00ac1fe28ac5fd76674a92b4df80d9c8a46c684399005deccc41386257c */
+  9.80987069605669190469329896435309665e-01L,  /* 3ffef643efeb82acd2d601583fc5158c */
+  -1.90899259410096419886996331536278461e-36L, /* bf8844cc5ab6a5903f931badc9cbde34 */
+/* sin(x) = 0.31aec65df552876f82ece9a2356713246eba6799983d7011b0b3698d6e1da919c15d57c30c1 */
+  1.94073102892909791156055200214145404e-01L,  /* 3ffc8d7632efaa943b7c17674d11ab39 */
+  -9.67304741051998267208945242944928999e-36L, /* bf8a9b7228b30cccf851fdc9e992ce52 */
+
+/* x = 2.03125000000000000000000000000000000e-01 3ffca000000000000000000000000000 */
+/* cos(x) = 0.fabca467fb3cb8f1d069f01d8ea33ade5bfd68296ecd1cc9f7b7609bbcf3676e726c3301334 */
+  9.79440951715548359998530954502987493e-01L,  /* 3ffef57948cff67971e3a0d3e03b1d46 */
+  4.42878056591560757066844797290067990e-35L,  /* 3f8cd6f2dfeb414b7668e64fbdbb04de */
+/* sin(x) = 0.33a4a5a19d86246710f602c44df4fa513f4639ce938477aeeabb82e8e0a7ed583a188879fd4 */
+  2.01731063801638804725038151164000971e-01L,  /* 3ffc9d252d0cec31233887b016226fa8 */
+  -4.27513434754966978435151290617384120e-36L, /* bf896bb02e718c5b1ee21445511f45c8 */
+
+/* x = 2.10937500000000000000000000000000000e-01 3ffcb000000000000000000000000000 */
+/* cos(x) = 0.fa5365e8f1d3ca27be1db5d76ae64d983d7470a4ab0f4ccf65a2b8c67a380df949953a09bc1 */
+  9.77835053797959793331971572944454549e-01L,  /* 3ffef4a6cbd1e3a7944f7c3b6baed5cd */
+  -3.79207422905180416937210853779192702e-35L, /* bf8c933e145c7adaa7859984d2ea39cc */
+/* sin(x) = 0.3599b652f40ec999df12a0a4c8561de159c98d4e54555de518b97f48886f715d8df5f4f093e */
+  2.09376712085993643711890752724881652e-01L,  /* 3ffcaccdb297a0764ccef895052642b1 */
+  -1.59470287344329449965314638482515925e-36L, /* bf880f531b3958d5d5510d73a3405bbc */
+
+/* x = 2.18750000000000000000000000000000000e-01 3ffcc000000000000000000000000000 */
+/* cos(x) = 0.f9e63e1d9e8b6f6f2e296bae5b5ed9c11fd7fa2fe11e09fc7bde901abed24b6365e72f7db4e */
+  9.76169473868635276723989035435135534e-01L,  /* 3ffef3cc7c3b3d16dede5c52d75cb6be */
+  -2.87727974249481583047944860626985460e-35L, /* bf8c31f701402e80f70fb01c210b7f2a */
+/* sin(x) = 0.378df09db8c332ce0d2b53d865582e4526ea336c768f68c32b496c6d11c1cd241bb9f1da523 */
+  2.17009581095010156760578095826055396e-01L,  /* 3ffcbc6f84edc6199670695a9ec32ac1 */
+  1.07356488794216831812829549198201194e-35L,  /* 3f8ac8a4dd466d8ed1ed1865692d8da2 */
+
+/* x = 2.26562500000000000000000000000000000e-01 3ffcd000000000000000000000000000 */
+/* cos(x) = 0.f9752eba9fff6b98842beadab054a932fb0f8d5b875ae63d6b2288d09b148921aeb6e52f61b */
+  9.74444313585988980349711056045434344e-01L,  /* 3ffef2ea5d753ffed7310857d5b560a9 */
+  3.09947905955053419304514538592548333e-35L,  /* 3f8c4997d87c6adc3ad731eb59144685 */
+/* sin(x) = 0.39814cb10513453cb97b21bc1ca6a337b150c21a675ab85503bc09a436a10ab1473934e20c8 */
+  2.24629204957705292350428549796424820e-01L,  /* 3ffccc0a6588289a29e5cbd90de0e535 */
+  2.42061510849297469844695751870058679e-36L,  /* 3f889bd8a8610d33ad5c2a81de04d21b */
+
+/* x = 2.34375000000000000000000000000000000e-01 3ffce000000000000000000000000000 */
+/* cos(x) = 0.f90039843324f9b940416c1984b6cbed1fc733d97354d4265788a86150493ce657cae032674 */
+  9.72659678244912752670913058267565260e-01L,  /* 3ffef20073086649f3728082d833096e */
+  -3.91759231819314904966076958560252735e-35L, /* bf8ca09701c6613465595ecd43babcf5 */
+/* sin(x) = 0.3b73c2bf6b4b9f668ef9499c81f0d965087f1753fa64b086e58cb8470515c18c1412f8c2e02 */
+  2.32235118611511462413930877746235872e-01L,  /* 3ffcdb9e15fb5a5cfb3477ca4ce40f87 */
+  -4.96930483364191020075024624332928910e-36L, /* bf89a6bde03a2b0166d3de469cd1ee3f */
+
+/* x = 2.42187500000000000000000000000000000e-01 3ffcf000000000000000000000000000 */
+/* cos(x) = 0.f887604e2c39dbb20e4ec5825059a789ffc95b275ad9954078ba8a28d3fcfe9cc2c1d49697b */
+  9.70815676770349462947490545785046027e-01L,  /* 3ffef10ec09c5873b7641c9d8b04a0b3 */
+  2.97458820972393859125277682021202860e-35L,  /* 3f8c3c4ffe4ad93ad6ccaa03c5d45147 */
+/* sin(x) = 0.3d654aff15cb457a0fca854698aba33039a8a40626609204472d9d40309b626eccc6dff0ffa */
+  2.39826857830661564441369251810886574e-01L,  /* 3ffceb2a57f8ae5a2bd07e542a34c55d */
+  2.39867036569896287240938444445071448e-36L,  /* 3f88981cd45203133049022396cea018 */
+
+/* x = 2.50000000000000000000000000000000000e-01 3ffd0000000000000000000000000000 */
+/* cos(x) = 0.f80aa4fbef750ba783d33cb95f94f8a41426dbe79edc4a023ef9ec13c944551c0795b84fee1 */
+  9.68912421710644784144595449494189205e-01L,  /* 3ffef01549f7deea174f07a67972bf2a */
+  -5.53634706113461989398873287749326500e-36L, /* bf89d6faf649061848ed7f704184fb0e */
+/* sin(x) = 0.3f55dda9e62aed7513bd7b8e6a3d1635dd5676648d7db525898d7086af9330f03c7f285442a */
+  2.47403959254522929596848704849389203e-01L,  /* 3ffcfaaeed4f31576ba89debdc7351e9 */
+  -7.36487001108599532943597115275811618e-36L, /* bf8a39445531336e50495b4ece51ef2a */
+
+/* x = 2.57812500000000000000000000000000000e-01 3ffd0800000000000000000000000000 */
+/* cos(x) = 0.f78a098069792daabc9ee42591b7c5a68cb1ab822aeb446b3311b4ba5371b8970e2c1547ad7 */
+  9.66950029230677822008341623610531503e-01L,  /* 3ffeef141300d2f25b55793dc84b2370 */
+  -4.38972214432792412062088059990480514e-35L, /* bf8cd2cb9a72a3eea8a5dca667725a2d */
+/* sin(x) = 0.414572fd94556e6473d620271388dd47c0ba050cdb5270112e3e370e8c4705ae006426fb5d5 */
+  2.54965960415878467487556574864872628e-01L,  /* 3ffd0515cbf65155b991cf58809c4e23 */
+  2.20280377918534721005071688328074154e-35L,  /* 3f8bd47c0ba050cdb5270112e3e370e9 */
+
+/* x = 2.65625000000000000000000000000000000e-01 3ffd1000000000000000000000000000 */
+/* cos(x) = 0.f7058fde0788dfc805b8fe88789e4f4253e3c50afe8b22f41159620ab5940ff7df9557c0d1f */
+  9.64928619104771009581074665315748371e-01L,  /* 3ffeee0b1fbc0f11bf900b71fd10f13d */
+  -3.66685832670820775002475545602761113e-35L, /* bf8c85ed60e1d7a80ba6e85f7534efaa */
+/* sin(x) = 0.4334033bcd90d6604f5f36c1d4b84451a87150438275b77470b50e5b968fa7962b5ffb379b7 */
+  2.62512399769153281450949626395692931e-01L,  /* 3ffd0cd00cef364359813d7cdb0752e1 */
+  3.24923677072031064673177178571821843e-36L,  /* 3f89146a1c5410e09d6ddd1c2d4396e6 */
+
+/* x = 2.73437500000000000000000000000000000e-01 3ffd1800000000000000000000000000 */
+/* cos(x) = 0.f67d3a26af7d07aa4bd6d42af8c0067fefb96d5b46c031eff53627f215ea3242edc3f2e13eb */
+  9.62848314709379699899701093480214365e-01L,  /* 3ffeecfa744d5efa0f5497ada855f180 */
+  4.88986966383343450799422013051821394e-36L,  /* 3f899ffbee5b56d1b00c7bfd4d89fc85 */
+/* sin(x) = 0.452186aa5377ab20bbf2524f52e3a06a969f47166ab88cf88c111ad12c55941021ef3317a1a */
+  2.70042816718585031552755063618827102e-01L,  /* 3ffd14861aa94ddeac82efc9493d4b8f */
+  -2.37608892440611310321138680065803162e-35L, /* bf8bf956960b8e99547730773eee52ed */
+
+/* x = 2.81250000000000000000000000000000000e-01 3ffd2000000000000000000000000000 */
+/* cos(x) = 0.f5f10a7bb77d3dfa0c1da8b57842783280d01ce3c0f82bae3b9d623c168d2e7c29977994451 */
+  9.60709243015561903066659350581313472e-01L,  /* 3ffeebe214f76efa7bf4183b516af085 */
+  -5.87011558231583960712013351601221840e-36L, /* bf89f35fcbf8c70fc1f5147118a770fa */
+/* sin(x) = 0.470df5931ae1d946076fe0dcff47fe31bb2ede618ebc607821f8462b639e1f4298b5ae87fd3 */
+  2.77556751646336325922023446828128568e-01L,  /* 3ffd1c37d64c6b8765181dbf8373fd20 */
+  -1.35848595468998128214344668770082997e-36L, /* bf87ce44d1219e71439f87de07b9d49c */
+
+/* x = 2.89062500000000000000000000000000000e-01 3ffd2800000000000000000000000000 */
+/* cos(x) = 0.f561030ddd7a78960ea9f4a32c6521554995667f5547bafee9ec48b3155cdb0f7fd00509713 */
+  9.58511534581228627301969408154919822e-01L,  /* 3ffeeac2061bbaf4f12c1d53e94658ca */
+  2.50770779371636481145735089393154404e-35L,  /* 3f8c0aaa4cab33faaa3dd7f74f624599 */
+/* sin(x) = 0.48f948446abcd6b0f7fccb100e7a1b26eccad880b0d24b59948c7cdd49514d44b933e6985c2 */
+  2.85053745940547424587763033323252561e-01L,  /* 3ffd23e52111aaf35ac3dff32c4039e8 */
+  2.04269325885902918802700123680403749e-35L,  /* 3f8bb26eccad880b0d24b59948c7cdd5 */
+
+/* x = 2.96875000000000000000000000000000000e-01 3ffd3000000000000000000000000000 */
+/* cos(x) = 0.f4cd261d3e6c15bb369c8758630d2ac00b7ace2a51c0631bfeb39ed158ba924cc91e259c195 */
+  9.56255323543175296975599942263028361e-01L,  /* 3ffee99a4c3a7cd82b766d390eb0c61a */
+  3.21616572190865997051103645135837207e-35L,  /* 3f8c56005bd671528e0318dff59cf68b */
+/* sin(x) = 0.4ae37710fad27c8aa9c4cf96c03519b9ce07dc08a1471775499f05c29f86190aaebaeb9716e */
+  2.92533342023327543624702326493913423e-01L,  /* 3ffd2b8ddc43eb49f22aa7133e5b00d4 */
+  1.93539408668704450308003687950685128e-35L,  /* 3f8b9b9ce07dc08a1471775499f05c2a */
+
+/* x = 3.04687500000000000000000000000000000e-01 3ffd3800000000000000000000000000 */
+/* cos(x) = 0.f43575f94d4f6b272f5fb76b14d2a64ab52df1ee8ddf7c651034e5b2889305a9ea9015d758a */
+  9.53940747608894733981324795987611623e-01L,  /* 3ffee86aebf29a9ed64e5ebf6ed629a5 */
+  2.88075689052478602008395972924657164e-35L,  /* 3f8c3255a96f8f746efbe32881a72d94 */
+/* sin(x) = 0.4ccc7a50127e1de0cb6b40c302c651f7bded4f9e7702b0471ae0288d091a37391950907202f */
+  2.99995083378683051163248282011699944e-01L,  /* 3ffd3331e94049f877832dad030c0b19 */
+  1.35174265535697850139283361475571050e-35L,  /* 3f8b1f7bded4f9e7702b0471ae0288d1 */
+
+/* x = 3.12500000000000000000000000000000000e-01 3ffd4000000000000000000000000000 */
+/* cos(x) = 0.f399f500c9e9fd37ae9957263dab8877102beb569f101ee4495350868e5847d181d50d3cca2 */
+  9.51567948048172202145488217364270962e-01L,  /* 3ffee733ea0193d3fa6f5d32ae4c7b57 */
+  6.36842628598115658308749288799884606e-36L,  /* 3f8a0ee2057d6ad3e203dc892a6a10d2 */
+/* sin(x) = 0.4eb44a5da74f600207aaa090f0734e288603ffadb3eb2542a46977b105f8547128036dcf7f0 */
+  3.07438514580380850670502958201982091e-01L,  /* 3ffd3ad129769d3d80081eaa8243c1cd */
+  1.06515172423204645839241099453417152e-35L,  /* 3f8ac510c07ff5b67d64a8548d2ef621 */
+
+/* x = 3.20312500000000000000000000000000000e-01 3ffd4800000000000000000000000000 */
+/* cos(x) = 0.f2faa5a1b74e82fd61fa05f9177380e8e69b7b15a945e8e5ae1124bf3d12b0617e03af4fab5 */
+  9.49137069684463027665847421762105623e-01L,  /* 3ffee5f54b436e9d05fac3f40bf22ee7 */
+  6.84433965991637152250309190468859701e-37L,  /* 3f86d1cd36f62b528bd1cb5c22497e7a */
+/* sin(x) = 0.509adf9a7b9a5a0f638a8fa3a60a199418859f18b37169a644fdb986c21ecb00133853bc35b */
+  3.14863181319745250865036315126939016e-01L,  /* 3ffd426b7e69ee69683d8e2a3e8e9828 */
+  1.92431240212432926993057705062834160e-35L,  /* 3f8b99418859f18b37169a644fdb986c */
+
+/* x = 3.28125000000000000000000000000000000e-01 3ffd5000000000000000000000000000 */
+/* cos(x) = 0.f2578a595224dd2e6bfa2eb2f99cc674f5ea6f479eae2eb580186897ae3f893df1113ca06b8 */
+  9.46648260886053321846099507295532976e-01L,  /* 3ffee4af14b2a449ba5cd7f45d65f33a */
+  -4.32906339663000890941529420498824645e-35L, /* bf8ccc5850ac85c30a8e8a53ff3cbb43 */
+/* sin(x) = 0.5280326c3cf481823ba6bb08eac82c2093f2bce3c4eb4ee3dec7df41c92c8a4226098616075 */
+  3.22268630433386625687745919893188031e-01L,  /* 3ffd4a00c9b0f3d20608ee9aec23ab21 */
+  -1.49505897804759263483853908335500228e-35L, /* bf8b3df6c0d431c3b14b11c213820be3 */
+
+/* x = 3.35937500000000000000000000000000000e-01 3ffd5800000000000000000000000000 */
+/* cos(x) = 0.f1b0a5b406b526d886c55feadc8d0dcc8eb9ae2ac707051771b48e05b25b000009660bdb3e3 */
+  9.44101673557004345630017691253124860e-01L,  /* 3ffee3614b680d6a4db10d8abfd5b91a */
+  1.03812535240120229609822461172145584e-35L,  /* 3f8ab991d735c558e0e0a2ee3691c0b6 */
+/* sin(x) = 0.54643b3da29de9b357155eef0f332fb3e66c83bf4dddd9491c5eb8e103ccd92d6175220ed51 */
+  3.29654409930860171914317725126463176e-01L,  /* 3ffd5190ecf68a77a6cd5c557bbc3ccd */
+  -1.22606996784743214973082192294232854e-35L, /* bf8b04c19937c40b22226b6e3a1471f0 */
+
+/* x = 3.43750000000000000000000000000000000e-01 3ffd6000000000000000000000000000 */
+/* cos(x) = 0.f105fa4d66b607a67d44e042725204435142ac8ad54dfb0907a4f6b56b06d98ee60f19e557a */
+  9.41497463127881068644511236053670815e-01L,  /* 3ffee20bf49acd6c0f4cfa89c084e4a4 */
+  3.20709366603165602071590241054884900e-36L,  /* 3f8910d450ab22b5537ec241e93dad5b */
+/* sin(x) = 0.5646f27e8bd65cbe3a5d61ff06572290ee826d9674a00246b05ae26753cdfc90d9ce81a7d02 */
+  3.37020069022253076261281754173810024e-01L,  /* 3ffd591bc9fa2f5972f8e97587fc195d */
+  -2.21435756148839473677777545049890664e-35L, /* bf8bd6f117d92698b5ffdb94fa51d98b */
+
+/* x = 3.51562500000000000000000000000000000e-01 3ffd6800000000000000000000000000 */
+/* cos(x) = 0.f0578ad01ede707fa39c09dc6b984afef74f3dc8d0efb0f4c5a6b13771145b3e0446fe33887 */
+  9.38835788546265488632578305984712554e-01L,  /* 3ffee0af15a03dbce0ff473813b8d731 */
+  -3.98758068773974031348585072752245458e-35L, /* bf8ca808458611b978827859d2ca7644 */
+/* sin(x) = 0.582850a41e1dd46c7f602ea244cdbbbfcdfa8f3189be794dda427ce090b5f85164f1f80ac13 */
+  3.44365158145698408207172046472223747e-01L,  /* 3ffd60a14290787751b1fd80ba891337 */
+  -3.19791885005480924937758467594051927e-36L, /* bf89100c815c339d9061ac896f60c7dc */
+
+/* x = 3.59375000000000000000000000000000000e-01 3ffd7000000000000000000000000000 */
+/* cos(x) = 0.efa559f5ec3aec3a4eb03319278a2d41fcf9189462261125fe6147b078f1daa0b06750a1654 */
+  9.36116812267055290294237411019508588e-01L,  /* 3ffedf4ab3ebd875d8749d6066324f14 */
+  3.40481591236710658435409862439032162e-35L,  /* 3f8c6a0fe7c8c4a31130892ff30a3d84 */
+/* sin(x) = 0.5a084e28e35fda2776dfdbbb5531d74ced2b5d17c0b1afc4647529d50c295e36d8ceec126c1 */
+  3.51689228994814059222584896955547016e-01L,  /* 3ffd682138a38d7f689ddb7f6eed54c7 */
+  1.75293433418270210567525412802083294e-35L,  /* 3f8b74ced2b5d17c0b1afc4647529d51 */
+
+/* x = 3.67187500000000000000000000000000000e-01 3ffd7800000000000000000000000000 */
+/* cos(x) = 0.eeef6a879146af0bf9b95ea2ea0ac0d3e2e4d7e15d93f48cbd41bf8e4fded40bef69e19eafa */
+  9.33340700242548435655299229469995527e-01L,  /* 3ffeddded50f228d5e17f372bd45d416 */
+  -4.75255707251679831124800898831382223e-35L, /* bf8cf960e8d940f513605b9a15f2038e */
+/* sin(x) = 0.5be6e38ce8095542bc14ee9da0d36483e6734bcab2e07624188af5653f114eeb46738fa899d */
+  3.58991834546065053677710299152868941e-01L,  /* 3ffd6f9b8e33a025550af053ba76834e */
+  -2.06772389262723368139416970257112089e-35L, /* bf8bb7c198cb4354d1f89dbe7750a9ac */
+
+/* x = 3.75000000000000000000000000000000000e-01 3ffd8000000000000000000000000000 */
+/* cos(x) = 0.ee35bf5ccac89052cd91ddb734d3a47e262e3b609db604e217053803be0091e76daf28a89b7 */
+  9.30507621912314291149476792229555481e-01L,  /* 3ffedc6b7eb9959120a59b23bb6e69a7 */
+  2.74541088551732982573335285685416092e-35L,  /* 3f8c23f13171db04edb02710b829c01e */
+/* sin(x) = 0.5dc40955d9084f48a94675a2498de5d851320ff5528a6afb3f2e24de240fce6cbed1ba0ccd6 */
+  3.66272529086047561372909351716264177e-01L,  /* 3ffd7710255764213d22a519d6892638 */
+  -1.96768433534936592675897818253108989e-35L, /* bf8ba27aecdf00aad759504c0d1db21e */
+
+/* x = 3.82812500000000000000000000000000000e-01 3ffd8800000000000000000000000000 */
+/* cos(x) = 0.ed785b5c44741b4493c56bcb9d338a151c6f6b85d8f8aca658b28572c162b199680eb9304da */
+  9.27617750192851909628030798799961350e-01L,  /* 3ffedaf0b6b888e83689278ad7973a67 */
+  7.58520371916345756281201167126854712e-36L,  /* 3f8a42a38ded70bb1f1594cb1650ae58 */
+/* sin(x) = 0.5f9fb80f21b53649c432540a50e22c53057ff42ae0fdf1307760dc0093f99c8efeb2fbd7073 */
+  3.73530868238692946416839752660848112e-01L,  /* 3ffd7e7ee03c86d4d92710c950294389 */
+  -1.48023494778986556048879113411517128e-35L, /* bf8b3acfa800bd51f020ecf889f23ff7 */
+
+/* x = 3.90625000000000000000000000000000000e-01 3ffd9000000000000000000000000000 */
+/* cos(x) = 0.ecb7417b8d4ee3fec37aba4073aa48f1f14666006fb431d9671303c8100d10190ec8179c41d */
+  9.24671261467036098502113014560138771e-01L,  /* 3ffed96e82f71a9dc7fd86f57480e755 */
+  -4.14187124860031825108649347251175815e-35L, /* bf8cb87075cccffc825e7134c767e1bf */
+/* sin(x) = 0.6179e84a09a5258a40e9b5face03e525f8b5753cd0105d93fe6298010c3458e84d75fe420e9 */
+  3.80766408992390192057200703388896675e-01L,  /* 3ffd85e7a1282694962903a6d7eb3810 */
+  -2.02009541175208636336924533372496107e-35L, /* bf8bada074a8ac32fefa26c019d67fef */
+
+/* x = 3.98437500000000000000000000000000000e-01 3ffd9800000000000000000000000000 */
+/* cos(x) = 0.ebf274bf0bda4f62447e56a093626798d3013b5942b1abfd155aacc9dc5c6d0806a20d6b9c1 */
+  9.21668335573351918175411368202712714e-01L,  /* 3ffed7e4e97e17b49ec488fcad4126c5 */
+  -1.83587995433957622948710263541479322e-35L, /* bf8b8672cfec4a6bd4e5402eaa553362 */
+/* sin(x) = 0.6352929dd264bd44a02ea766325d8aa8bd9695fc8def3caefba5b94c9a3c873f7b2d3776ead */
+  3.87978709727025046051079690813741960e-01L,  /* 3ffd8d4a4a774992f51280ba9d98c976 */
+  8.01904783870935075844443278617586301e-36L,  /* 3f8a5517b2d2bf91bde795df74b72993 */
+
+/* x = 4.06250000000000000000000000000000000e-01 3ffda000000000000000000000000000 */
+/* cos(x) = 0.eb29f839f201fd13b93796827916a78f15c85230a4e8ea4b21558265a14367e1abb4c30695a */
+  9.18609155794918267837824977718549863e-01L,  /* 3ffed653f073e403fa27726f2d04f22d */
+  2.97608282778274433460057745798409849e-35L,  /* 3f8c3c78ae429185274752590aac132d */
+/* sin(x) = 0.6529afa7d51b129631ec197c0a840a11d7dc5368b0a47956feb285caa8371c4637ef17ef01b */
+  3.95167330240934236244832640419653657e-01L,  /* 3ffd94a6be9f546c4a58c7b065f02a10 */
+  7.57560031388312550940040194042627704e-36L,  /* 3f8a423afb8a6d16148f2adfd650b955 */
+
+/* x = 4.14062500000000000000000000000000000e-01 3ffda800000000000000000000000000 */
+/* cos(x) = 0.ea5dcf0e30cf03e6976ef0b1ec26515fba47383855c3b4055a99b5e86824b2cd1a691fdca7b */
+  9.15493908848301228563917732180221882e-01L,  /* 3ffed4bb9e1c619e07cd2edde163d84d */
+  -3.50775517955306954815090901168305659e-35L, /* bf8c75022dc63e3d51e25fd52b3250bd */
+/* sin(x) = 0.66ff380ba0144109e39a320b0a3fa5fd65ea0585bcbf9b1a769a9b0334576c658139e1a1cbe */
+  4.02331831777773111217105598880982387e-01L,  /* 3ffd9bfce02e805104278e68c82c28ff */
+  -1.95678722882848174723569916504871563e-35L, /* bf8ba029a15fa7a434064e5896564fcd */
+
+/* x = 4.21875000000000000000000000000000000e-01 3ffdb000000000000000000000000000 */
+/* cos(x) = 0.e98dfc6c6be031e60dd3089cbdd18a75b1f6b2c1e97f79225202f03dbea45b07a5ec4efc062 */
+  9.12322784872117846492029542047341734e-01L,  /* 3ffed31bf8d8d7c063cc1ba611397ba3 */
+  7.86903886556373674267948132178845568e-36L,  /* 3f8a4eb63ed6583d2fef244a405e07b8 */
+/* sin(x) = 0.68d32473143327973bc712bcc4ccddc47630d755850c0655243b205934dc49ffed8eb76adcb */
+  4.09471777053295066122694027011452236e-01L,  /* 3ffda34c91cc50cc9e5cef1c4af31333 */
+  2.23945241468457597921655785729821354e-35L,  /* 3f8bdc47630d755850c0655243b20593 */
+
+/* x = 4.29687500000000000000000000000000000e-01 3ffdb800000000000000000000000000 */
+/* cos(x) = 0.e8ba8393eca7821aa563d83491b6101189b3b101c3677f73d7bad7c10f9ee02b7ab4009739a */
+  9.09095977415431051650381735684476417e-01L,  /* 3ffed1750727d94f04354ac7b069236c */
+  1.20886014028444155733776025085677953e-35L,  /* 3f8b01189b3b101c3677f73d7bad7c11 */
+/* sin(x) = 0.6aa56d8e8249db4eb60a761fe3f9e559be456b9e13349ca99b0bfb787f22b95db3b70179615 */
+  4.16586730282041119259112448831069657e-01L,  /* 3ffdaa95b63a09276d3ad829d87f8fe8 */
+  -2.00488106831998813675438269796963612e-35L, /* bf8baa641ba9461eccb635664f404878 */
+
+/* x = 4.37500000000000000000000000000000000e-01 3ffdc000000000000000000000000000 */
+/* cos(x) = 0.e7e367d2956cfb16b6aa11e5419cd0057f5c132a6455bf064297e6a76fe2b72bb630d6d50ff */
+  9.05813683425936420744516660652700258e-01L,  /* 3ffecfc6cfa52ad9f62d6d5423ca833a */
+  -3.60950307605941169775676563004467163e-35L, /* bf8c7fd4051f66acdd5207cdeb40cac5 */
+/* sin(x) = 0.6c760c14c8585a51dbd34660ae6c52ac7036a0b40887a0b63724f8b4414348c3063a637f457 */
+  4.23676257203938010361683988031102480e-01L,  /* 3ffdb1d83053216169476f4d1982b9b1 */
+  1.40484456388654470329473096579312595e-35L,  /* 3f8b2ac7036a0b40887a0b63724f8b44 */
+
+/* x = 4.45312500000000000000000000000000000e-01 3ffdc800000000000000000000000000 */
+/* cos(x) = 0.e708ac84d4172a3e2737662213429e14021074d7e702e77d72a8f1101a7e70410df8273e9aa */
+  9.02476103237941504925183272675895999e-01L,  /* 3ffece115909a82e547c4e6ecc442685 */
+  2.26282899501344419018306295680210602e-35L,  /* 3f8be14021074d7e702e77d72a8f1102 */
+/* sin(x) = 0.6e44f8c36eb10a1c752d093c00f4d47ba446ac4c215d26b0316442f168459e677d06e7249e3 */
+  4.30739925110803197216321517850849190e-01L,  /* 3ffdb913e30dbac42871d4b424f003d3 */
+  1.54096780001629398850891218396761548e-35L,  /* 3f8b47ba446ac4c215d26b0316442f17 */
+
+/* x = 4.53125000000000000000000000000000000e-01 3ffdd000000000000000000000000000 */
+/* cos(x) = 0.e62a551594b970a770b15d41d4c0e483e47aca550111df6966f9e7ac3a94ae49e6a71eb031e */
+  8.99083440560138456216544929209379307e-01L,  /* 3ffecc54aa2b2972e14ee162ba83a982 */
+  -2.06772615490904370666670275154751976e-35L, /* bf8bb7c1b8535aafeee209699061853c */
+/* sin(x) = 0.70122c5ec5028c8cff33abf4fd340ccc382e038379b09cf04f9a52692b10b72586060cbb001 */
+  4.37777302872755132861618974702796680e-01L,  /* 3ffdc048b17b140a3233fcceafd3f4d0 */
+  9.62794364503442612477117426033922467e-36L,  /* 3f8a998705c0706f36139e09f34a4d25 */
+
+/* x = 4.60937500000000000000000000000000000e-01 3ffdd800000000000000000000000000 */
+/* cos(x) = 0.e54864fe33e8575cabf5bd0e5cf1b1a8bc7c0d5f61702450fa6b6539735820dd2603ae355d5 */
+  8.95635902463170698900570000446256350e-01L,  /* 3ffeca90c9fc67d0aeb957eb7a1cb9e3 */
+  3.73593741659866883088620495542311808e-35L,  /* 3f8c8d45e3e06afb0b812287d35b29cc */
+/* sin(x) = 0.71dd9fb1ff4677853acb970a9f6729c6e3aac247b1c57cea66c77413f1f98e8b9e98e49d851 */
+  4.44787960964527211433056012529525211e-01L,  /* 3ffdc7767ec7fd19de14eb2e5c2a7d9d */
+  -1.67187936511493678007508371613954899e-35L, /* bf8b6391c553db84e3a831599388bec1 */
+
+/* x = 4.68750000000000000000000000000000000e-01 3ffde000000000000000000000000000 */
+/* cos(x) = 0.e462dfc670d421ab3d1a15901228f146a0547011202bf5ab01f914431859aef577966bc4fa4 */
+  8.92133699366994404723900253723788575e-01L,  /* 3ffec8c5bf8ce1a843567a342b202452 */
+  -1.10771937602567314732693079264692504e-35L, /* bf8ad72bf571fddbfa814a9fc0dd779d */
+/* sin(x) = 0.73a74b8f52947b681baf6928eb3fb021769bf4779bad0e3aa9b1cdb75ec60aad9fc63ff19d5 */
+  4.51771471491683776581688750134062870e-01L,  /* 3ffdce9d2e3d4a51eda06ebda4a3acff */
+  -1.19387223016472295893794387275284505e-35L, /* bf8afbd12c81710c8a5e38aac9c64914 */
+
+/* x = 4.76562500000000000000000000000000000e-01 3ffde800000000000000000000000000 */
+/* cos(x) = 0.e379c9045f29d517c4808aa497c2057b2b3d109e76c0dc302d4d0698b36e3f0bdbf33d8e952 */
+  8.88577045028035543317609023116020980e-01L,  /* 3ffec6f39208be53aa2f890115492f84 */
+  4.12354278954664731443813655177022170e-36L,  /* 3f895ecacf44279db0370c0b5341a62d */
+/* sin(x) = 0.756f28d011d98528a44a75fc29c779bd734ecdfb582fdb74b68a4c4c4be54cfd0b2d3ad292f */
+  4.58727408216736592377295028972874773e-01L,  /* 3ffdd5bca340476614a29129d7f0a71e */
+  -4.70946994194182908929251719575431779e-36L, /* bf8990a32c4c8129f40922d25d6ceced */
+
+/* x = 4.84375000000000000000000000000000000e-01 3ffdf000000000000000000000000000 */
+/* cos(x) = 0.e28d245c58baef72225e232abc003c4366acd9eb4fc2808c2ab7fe7676cf512ac7f945ae5fb */
+  8.84966156526143291697296536966647926e-01L,  /* 3ffec51a48b8b175dee444bc46557800 */
+  4.53370570288325630442037826313462165e-35L,  /* 3f8ce21b3566cf5a7e14046155bff3b4 */
+/* sin(x) = 0.77353054ca72690d4c6e171fd99e6b39fa8e1ede5f052fd2964534c75340970a3a9cd3c5c32 */
+  4.65655346585160182681199512507546779e-01L,  /* 3ffddcd4c15329c9a43531b85c7f667a */
+  -1.56282598978971872478619772155305961e-35L, /* bf8b4c60571e121a0fad02d69bacb38b */
+
+/* x = 4.92187500000000000000000000000000000e-01 3ffdf800000000000000000000000000 */
+/* cos(x) = 0.e19cf580eeec046aa1422fa74807ecefb2a1911c94e7b5f20a00f70022d940193691e5bd790 */
+  8.81301254251340599140161908298100173e-01L,  /* 3ffec339eb01ddd808d542845f4e9010 */
+  -1.43419192312116687783945619009629445e-35L, /* bf8b3104d5e6ee36b184a0df5ff08ffe */
+/* sin(x) = 0.78f95b0560a9a3bd6df7bd981dc38c61224d08bc20631ea932e605e53b579e9e0767dfcbbcb */
+  4.72554863751304451146551317808516942e-01L,  /* 3ffde3e56c1582a68ef5b7def660770e */
+  9.31324774957768018850224267625371204e-36L,  /* 3f8a8c2449a117840c63d5265cc0bca7 */
+
+/* x = 5.00000000000000000000000000000000000e-01 3ffe0000000000000000000000000000 */
+/* cos(x) = 0.e0a94032dbea7cedbddd9da2fafad98556566b3a89f43eabd72350af3e8b19e801204d8fe2e */
+  8.77582561890372716116281582603829681e-01L,  /* 3ffec1528065b7d4f9db7bbb3b45f5f6 */
+  -2.89484960181363924855192538540698851e-35L, /* bf8c33d54d4ca62bb05e0aa146e57a86 */
+/* sin(x) = 0.7abba1d12c17bfa1d92f0d93f60ded9992f45b4fcaf13cd58b303693d2a0db47db35ae8a3a9 */
+  4.79425538604203000273287935215571402e-01L,  /* 3ffdeaee8744b05efe8764bc364fd838 */
+  -1.38426977616718318950175848639381926e-35L, /* bf8b2666d0ba4b0350ec32a74cfc96c3 */
+
+/* x = 5.07812500000000000000000000000000000e-01 3ffe0400000000000000000000000000 */
+/* cos(x) = 0.dfb20840f3a9b36f7ae2c515342890b5ec583b8366cc2b55029e95094d31112383f2553498b */
+  8.73810306413054508282556837071377159e-01L,  /* 3ffebf641081e75366def5c58a2a6851 */
+  1.25716864497849302237218128599994785e-35L,  /* 3f8b0b5ec583b8366cc2b55029e95095 */
+/* sin(x) = 0.7c7bfdaf13e5ed17212f8a7525bfb113aba6c0741b5362bb8d59282a850b63716bca0c910f0 */
+  4.86266951793275574311011306895834993e-01L,  /* 3ffdf1eff6bc4f97b45c84be29d496ff */
+  -1.12269393250914752644352376448094271e-35L, /* bf8add8a8b27f17c9593a88e54dafaaf */
+
+/* x = 5.15625000000000000000000000000000000e-01 3ffe0800000000000000000000000000 */
+/* cos(x) = 0.deb7518814a7a931bbcc88c109cd41c50bf8bb48f20ae8c36628d1d3d57574f7dc58f27d91c */
+  8.69984718058417388828915599901466243e-01L,  /* 3ffebd6ea310294f526377991182139b */
+  -4.68168638300575626782741319792183837e-35L, /* bf8cf1d7a03a25b86fa8b9e4ceb97161 */
+/* sin(x) = 0.7e3a679daaf25c676542bcb4028d0964172961c921823a4ef0c3a9070d886dbd073f6283699 */
+  4.93078685753923057265136552753487121e-01L,  /* 3ffdf8e99e76abc9719d950af2d00a34 */
+  7.06498693112535056352301101088624950e-36L,  /* 3f8a2c82e52c3924304749de187520e2 */
+
+/* x = 5.23437500000000000000000000000000000e-01 3ffe0c00000000000000000000000000 */
+/* cos(x) = 0.ddb91ff318799172bd2452d0a3889f5169c64a0094bcf0b8aa7dcf0d7640a2eba68955a80be */
+  8.66106030320656714696616831654267220e-01L,  /* 3ffebb723fe630f322e57a48a5a14711 */
+  2.35610597588322493119667003904687628e-35L,  /* 3f8bf5169c64a0094bcf0b8aa7dcf0d7 */
+/* sin(x) = 0.7ff6d8a34bd5e8fa54c97482db5159df1f24e8038419c0b448b9eea8939b5d4dfcf40900257 */
+  4.99860324733013463819556536946425724e-01L,  /* 3ffdffdb628d2f57a3e95325d20b6d45 */
+  1.94636052312235297538564591686645139e-35L,  /* 3f8b9df1f24e8038419c0b448b9eea89 */
+
+/* x = 5.31250000000000000000000000000000000e-01 3ffe1000000000000000000000000000 */
+/* cos(x) = 0.dcb7777ac420705168f31e3eb780ce9c939ecada62843b54522f5407eb7f21e556059fcd734 */
+  8.62174479934880504367162510253324274e-01L,  /* 3ffeb96eeef58840e0a2d1e63c7d6f02 */
+  -3.71556818317533582234562471835771823e-35L, /* bf8c8b1b6309a92cebde255d6e855fc1 */
+/* sin(x) = 0.81b149ce34caa5a4e650f8d09fd4d6aa74206c32ca951a93074c83b2d294d25dbb0f7fdfad2 */
+  5.06611454814257367642296000893867192e-01L,  /* 3ffe0362939c69954b49cca1f1a13faa */
+  -3.10963699824274155702706043065967062e-35L, /* bf8c4aac5efc9e69ab572b67c59be269 */
+
+/* x = 5.39062500000000000000000000000000000e-01 3ffe1400000000000000000000000000 */
+/* cos(x) = 0.dbb25c25b8260c14f6e7bc98ec991b70c65335198b0ab628bad20cc7b229d4dd62183cfa055 */
+  8.58190306862660347046629564970494649e-01L,  /* 3ffeb764b84b704c1829edcf7931d932 */
+  2.06439574601190798155563653000684861e-35L,  /* 3f8bb70c65335198b0ab628bad20cc7b */
+/* sin(x) = 0.8369b434a372da7eb5c8a71fe36ce1e0b2b493f6f5cb2e38bcaec2a556b3678c401940d1c3c */
+  5.13331663943471218288801270215706878e-01L,  /* 3ffe06d3686946e5b4fd6b914e3fc6da */
+  -2.26614796466671970772244932848067224e-35L, /* bf8be1f4d4b6c090a34d1c743513d5ab */
+
+/* x = 5.46875000000000000000000000000000000e-01 3ffe1800000000000000000000000000 */
+/* cos(x) = 0.daa9d20860827063fde51c09e855e9932e1b17143e7244fd267a899d41ae1f3bc6a0ec42e27 */
+  8.54153754277385385143451785105103176e-01L,  /* 3ffeb553a410c104e0c7fbca3813d0ac */
+  -1.68707534013095152873222061722573172e-35L, /* bf8b66cd1e4e8ebc18dbb02d9857662c */
+/* sin(x) = 0.852010f4f0800521378bd8dd614753d080c2e9e0775ffc609947b9132f5357404f464f06a58 */
+  5.20020541953727004760213699874674730e-01L,  /* 3ffe0a4021e9e1000a426f17b1bac28f */
+  -3.32415021330884924833711842866896734e-35L, /* bf8c617bf9e8b0fc45001cfb35c23767 */
+
+/* x = 5.54687500000000000000000000000000000e-01 3ffe1c00000000000000000000000000 */
+/* cos(x) = 0.d99ddd44e44a43d4d4a3a3ed95204106fd54d78e8c7684545c0da0b7c2c72be7a89b7c182ad */
+  8.50065068549420263957072899177793617e-01L,  /* 3ffeb33bba89c89487a9a94747db2a41 */
+  -4.73753917078785974356016104842568442e-35L, /* bf8cf7c81559438b9c4bdd5d1f92fa42 */
+/* sin(x) = 0.86d45935ab396cb4e421e822dee54f3562dfcefeaa782184c23401d231f5ad981a1cc195b18 */
+  5.26677680590386730710789410624833901e-01L,  /* 3ffe0da8b26b5672d969c843d045bdcb */
+  -3.67066148195515214077582496518566735e-35L, /* bf8c8654e901880aac3ef3d9ee5ff16e */
+
+/* x = 5.62500000000000000000000000000000000e-01 3ffe2000000000000000000000000000 */
+/* cos(x) = 0.d88e820b1526311dd561efbc0c1a9a5375eb26f65d246c5744b13ca26a7e0fd42556da843c8 */
+  8.45924499231067954459723078597493262e-01L,  /* 3ffeb11d04162a4c623baac3df781835 */
+  1.98054947141989878179164342925274053e-35L,  /* 3f8ba5375eb26f65d246c5744b13ca27 */
+/* sin(x) = 0.88868625b4e1dbb2313310133022527200c143a5cb16637cb7daf8ade82459ff2e98511f40f */
+  5.33302673536020173329131103308161529e-01L,  /* 3ffe110d0c4b69c3b764626620266045 */
+  -3.42715291319551615996993795226755157e-35L, /* bf8c6c6ff9f5e2d1a74ce41a41283a91 */
+
+/* x = 5.70312500000000000000000000000000000e-01 3ffe2400000000000000000000000000 */
+/* cos(x) = 0.d77bc4985e93a607c9d868b906bbc6bbe3a04258814acb0358468b826fc91bd4d814827f65e */
+  8.41732299041338366963111794309701085e-01L,  /* 3ffeaef78930bd274c0f93b0d1720d78 */
+  -4.30821936750410026005408345400225948e-35L, /* bf8cca20e2fded3bf5a9a7e53dcba3ed */
+/* sin(x) = 0.8a3690fc5bfc11bf9535e2739a8512f448a41251514bbed7fc18d530f9b4650fcbb2861b0aa */
+  5.39895116435204405041660709903993340e-01L,  /* 3ffe146d21f8b7f8237f2a6bc4e7350a */
+  1.42595803521626714477253741404712093e-35L,  /* 3f8b2f448a41251514bbed7fc18d5310 */
+
+/* x = 5.78125000000000000000000000000000000e-01 3ffe2800000000000000000000000000 */
+/* cos(x) = 0.d665a937b4ef2b1f6d51bad6d988a4419c1d7051faf31a9efa151d7631117efac03713f950a */
+  8.37488723850523685315353348917240617e-01L,  /* 3ffeaccb526f69de563edaa375adb311 */
+  2.72761997872084533045777718677326179e-35L,  /* 3f8c220ce0eb828fd798d4f7d0a8ebb2 */
+/* sin(x) = 0.8be472f9776d809af2b88171243d63d66dfceeeb739cc894e023fbc165a0e3f26ff729c5d57 */
+  5.46454606919203564403349553749411001e-01L,  /* 3ffe17c8e5f2eedb0135e57102e2487b */
+  -2.11870230730160315420936523771864858e-35L, /* bf8bc29920311148c63376b1fdc043ea */
+
+/* x = 5.85937500000000000000000000000000000e-01 3ffe2c00000000000000000000000000 */
+/* cos(x) = 0.d54c3441844897fc8f853f0655f1ba695eba9fbfd7439dbb1171d862d9d9146ca5136f825ac */
+  8.33194032664581363070224042208032321e-01L,  /* 3ffeaa98688308912ff91f0a7e0cabe3 */
+  4.39440050052045486567668031751259899e-35L,  /* 3f8cd34af5d4fdfeba1cedd88b8ec317 */
+/* sin(x) = 0.8d902565817ee7839bce3cd128060119492cd36d42d82ada30d7f8bde91324808377ddbf5d4 */
+  5.52980744630527369849695082681623667e-01L,  /* 3ffe1b204acb02fdcf07379c79a2500c */
+  8.26624790417342895897164123189984127e-37L,  /* 3f8719492cd36d42d82ada30d7f8bde9 */
+
+/* x = 5.93750000000000000000000000000000000e-01 3ffe3000000000000000000000000000 */
+/* cos(x) = 0.d42f6a1b9f0168cdf031c2f63c8d9304d86f8d34cb1d5fccb68ca0f2241427fc18d1fd5bbdf */
+  8.28848487609325734810171790119116638e-01L,  /* 3ffea85ed4373e02d19be06385ec791b */
+  1.43082508100496581719048175506239770e-35L,  /* 3f8b304d86f8d34cb1d5fccb68ca0f22 */
+/* sin(x) = 0.8f39a191b2ba6122a3fa4f41d5a3ffd421417d46f19a22230a14f7fcc8fce5c75b4b28b29d1 */
+  5.59473131247366877384844006003116688e-01L,  /* 3ffe1e7343236574c24547f49e83ab48 */
+  -1.28922620524163922306886952100992796e-37L, /* bf845ef5f415c8732eeee7af584019b8 */
+
+/* x = 6.01562500000000000000000000000000000e-01 3ffe3400000000000000000000000000 */
+/* cos(x) = 0.d30f4f392c357ab0661c5fa8a7d9b26627846fef214b1d19a22379ff9eddba087cf410eb097 */
+  8.24452353914429207485643598212356053e-01L,  /* 3ffea61e9e72586af560cc38bf514fb3 */
+  3.79160239225080026987031418939026741e-35L,  /* 3f8c93313c237f790a58e8cd111bcffd */
+/* sin(x) = 0.90e0e0d81ca678796cc92c8ea8c2815bc72ca78abe571bfa8576aacc571e096a33237e0e830 */
+  5.65931370507905990773159095689276114e-01L,  /* 3ffe21c1c1b0394cf0f2d992591d5185 */
+  1.02202775968053982310991962521535027e-36L,  /* 3f875bc72ca78abe571bfa8576aacc57 */
+
+/* x = 6.09375000000000000000000000000000000e-01 3ffe3800000000000000000000000000 */
+/* cos(x) = 0.d1ebe81a95ee752e48a26bcd32d6e922d7eb44b8ad2232f6930795e84b56317269b9dd1dfa6 */
+  8.20005899897234008255550633876556043e-01L,  /* 3ffea3d7d0352bdcea5c9144d79a65ae */
+  -1.72008811955230823416724332297991247e-35L, /* bf8b6dd2814bb4752ddcd096cf86a17b */
+/* sin(x) = 0.9285dc9bc45dd9ea3d02457bcce59c4175aab6ff7929a8d287195525fdace200dba032874fb */
+  5.72355068234507240384953706824503608e-01L,  /* 3ffe250bb93788bbb3d47a048af799cb */
+  2.12572273479933123944580199464514529e-35L,  /* 3f8bc4175aab6ff7929a8d2871955260 */
+
+/* x = 6.17187500000000000000000000000000000e-01 3ffe3c00000000000000000000000000 */
+/* cos(x) = 0.d0c5394d772228195e25736c03574707de0af1ca344b13bd3914bfe27518e9e426f5deff1e1 */
+  8.15509396946375476876345384201386217e-01L,  /* 3ffea18a729aee445032bc4ae6d806af */
+  -4.28589138410712954051679139949341961e-35L, /* bf8cc7c10fa871ae5da76216375a00ec */
+/* sin(x) = 0.94288e48bd0335fc41c4cbd2920497a8f5d1d8185c99fa0081f90c27e2a53ffdd208a0dbe69 */
+  5.78743832357770354521111378581385347e-01L,  /* 3ffe28511c917a066bf8838997a52409 */
+  1.77998063432551282609698670002456093e-35L,  /* 3f8b7a8f5d1d8185c99fa0081f90c27e */
+
+/* x = 6.25000000000000000000000000000000000e-01 3ffe4000000000000000000000000000 */
+/* cos(x) = 0.cf9b476c897c25c5bfe750dd3f308eaf7bcc1ed00179a256870f4200445043dcdb1974b5878 */
+  8.10963119505217902189534803941080724e-01L,  /* 3ffe9f368ed912f84b8b7fcea1ba7e61 */
+  1.10481292856794436426051402418804358e-35L,  /* 3f8ad5ef7983da002f344ad0e1e84009 */
+/* sin(x) = 0.95c8ef544210ec0b91c49bd2aa09e8515fa61a156ebb10f5f8c232a6445b61ebf3c2ec268f9 */
+  5.85097272940462154805399314150080459e-01L,  /* 3ffe2b91dea88421d817238937a55414 */
+  -1.78164576278056195136525335403380464e-35L, /* bf8b7aea059e5ea9144ef0a073dcd59c */
+
+/* x = 6.32812500000000000000000000000000000e-01 3ffe4400000000000000000000000000 */
+/* cos(x) = 0.ce6e171f92f2e27f32225327ec440ddaefae248413efc0e58ceee1ae369aabe73f88c87ed1a */
+  8.06367345055103913698795406077297399e-01L,  /* 3ffe9cdc2e3f25e5c4fe6444a64fd888 */
+  1.04235088143133625463876245029180850e-35L,  /* 3f8abb5df5c490827df81cb19ddc35c7 */
+/* sin(x) = 0.9766f93cd18413a6aafc1cfc6fc28abb6817bf94ce349901ae3f48c3215d3eb60acc5f78903 */
+  5.91415002201316315087000225758031236e-01L,  /* 3ffe2ecdf279a308274d55f839f8df85 */
+  8.07390238063560077355762466502569603e-36L,  /* 3f8a576d02f7f299c6932035c7e91864 */
+
+/* x = 6.40625000000000000000000000000000000e-01 3ffe4800000000000000000000000000 */
+/* cos(x) = 0.cd3dad1b5328a2e459f993f4f5108819faccbc4eeba9604e81c7adad51cc8a2561631a06826 */
+  8.01722354098418450607492605652964208e-01L,  /* 3ffe9a7b5a36a65145c8b3f327e9ea21 */
+  6.09487851305233089325627939458963741e-36L,  /* 3f8a033f599789dd752c09d038f5b5aa */
+/* sin(x) = 0.9902a58a45e27bed68412b426b675ed503f54d14c8172e0d373f42cadf04daf67319a7f94be */
+  5.97696634538701531238647618967334337e-01L,  /* 3ffe32054b148bc4f7dad0825684d6cf */
+  -2.49527608940873714527427941350461554e-35L, /* bf8c0957e0559759bf468f964605e9a9 */
+
+/* x = 6.48437500000000000000000000000000000e-01 3ffe4c00000000000000000000000000 */
+/* cos(x) = 0.cc0a0e21709883a3ff00911e11a07ee3bd7ea2b04e081be99be0264791170761ae64b8b744a */
+  7.97028430141468342004642741431945296e-01L,  /* 3ffe98141c42e1310747fe01223c2341 */
+  -8.35364432831812599727083251866305534e-37L, /* bf871c42815d4fb1f7e416641fd9b86f */
+/* sin(x) = 0.9a9bedcdf01b38d993f3d7820781de292033ead73b89e28f39313dbe3a6e463f845b5fa8490 */
+  6.03941786554156657267270287527367726e-01L,  /* 3ffe3537db9be03671b327e7af040f04 */
+  -2.54578992328947177770363936132309779e-35L, /* bf8c0eb6fe60a94623b0eb863676120e */
+
+/* x = 6.56250000000000000000000000000000000e-01 3ffe5000000000000000000000000000 */
+/* cos(x) = 0.cad33f00658fe5e8204bbc0f3a66a0e6a773f87987a780b243d7be83b3db1448ca0e0e62787 */
+  7.92285859677178543141501323781709399e-01L,  /* 3ffe95a67e00cb1fcbd04097781e74cd */
+  2.47519558228473167879248891673807645e-35L,  /* 3f8c07353b9fc3cc3d3c05921ebdf41e */
+/* sin(x) = 0.9c32cba2b14156ef05256c4f857991ca6a547cd7ceb1ac8a8e62a282bd7b9183648a462bd04 */
+  6.10150077075791371273742393566183220e-01L,  /* 3ffe386597456282adde0a4ad89f0af3 */
+  1.33842237929938963780969418369150532e-35L,  /* 3f8b1ca6a547cd7ceb1ac8a8e62a282c */
+
+/* x = 6.64062500000000000000000000000000000e-01 3ffe5400000000000000000000000000 */
+/* cos(x) = 0.c99944936cf48c8911ff93fe64b3ddb7981e414bdaf6aae1203577de44878c62bc3bc9cf7b9 */
+  7.87494932167606083931328295965533034e-01L,  /* 3ffe93328926d9e9191223ff27fcc968 */
+  -2.57915385618070637156514241185180920e-35L, /* bf8c12433f0df5a1284aa8f6fe54410e */
+/* sin(x) = 0.9dc738ad14204e689ac582d0f85826590feece34886cfefe2e08cf2bb8488d55424dc9d3525 */
+  6.16321127181550943005700433761731837e-01L,  /* 3ffe3b8e715a28409cd1358b05a1f0b0 */
+  2.88497530050197716298085892460478666e-35L,  /* 3f8c32c87f7671a44367f7f17046795e */
+
+/* x = 6.71875000000000000000000000000000000e-01 3ffe5800000000000000000000000000 */
+/* cos(x) = 0.c85c23c26ed7b6f014ef546c47929682122876bfbf157de0aff3c4247d820c746e32cd4174f */
+  7.82655940026272796930787447428139026e-01L,  /* 3ffe90b84784ddaf6de029dea8d88f25 */
+  1.69332045679237919427807771288506254e-35L,  /* 3f8b682122876bfbf157de0aff3c4248 */
+/* sin(x) = 0.9f592e9b66a9cf906a3c7aa3c10199849040c45ec3f0a747597311038101780c5f266059dbf */
+  6.22454560222343683041926705090443330e-01L,  /* 3ffe3eb25d36cd539f20d478f5478203 */
+  1.91974786921147072717621236192269859e-35L,  /* 3f8b9849040c45ec3f0a747597311038 */
+
+/* x = 6.79687500000000000000000000000000000e-01 3ffe5c00000000000000000000000000 */
+/* cos(x) = 0.c71be181ecd6875ce2da5615a03cca207d9adcb9dfb0a1d6c40a4f0056437f1a59ccddd06ee */
+  7.77769178600317903122203513685412863e-01L,  /* 3ffe8e37c303d9ad0eb9c5b4ac2b407a */
+  -4.05296033424632846931240580239929672e-35L, /* bf8caefc13291a31027af149dfad87fd */
+/* sin(x) = 0.a0e8a725d33c828c11fa50fd9e9a15ffecfad43f3e534358076b9b0f6865694842b1e8c67dc */
+  6.28550001845029662028004327939032867e-01L,  /* 3ffe41d14e4ba679051823f4a1fb3d34 */
+  1.65507421184028099672784511397428852e-35L,  /* 3f8b5ffecfad43f3e534358076b9b0f7 */
+
+/* x = 6.87500000000000000000000000000000000e-01 3ffe6000000000000000000000000000 */
+/* cos(x) = 0.c5d882d2ee48030c7c07d28e981e34804f82ed4cf93655d2365389b716de6ad44676a1cc5da */
+  7.72834946152471544810851845913425178e-01L,  /* 3ffe8bb105a5dc900618f80fa51d303c */
+  3.94975229341211664237241534741146939e-35L,  /* 3f8ca4027c176a67c9b2ae91b29c4db9 */
+/* sin(x) = 0.a2759c0e79c35582527c32b55f5405c182c66160cb1d9eb7bb0b7cdf4ad66f317bda4332914 */
+  6.34607080015269296850309914203671436e-01L,  /* 3ffe44eb381cf386ab04a4f8656abea8 */
+  4.33025916939968369326060156455927002e-36L,  /* 3f897060b1985832c767adeec2df37d3 */
+
+/* x = 6.95312500000000000000000000000000000e-01 3ffe6400000000000000000000000000 */
+/* cos(x) = 0.c4920cc2ec38fb891b38827db08884fc66371ac4c2052ca8885b981bbcfd3bb7b093ee31515 */
+  7.67853543842850365879920759114193964e-01L,  /* 3ffe89241985d871f712367104fb6111 */
+  3.75100035267325597157244776081706979e-36L,  /* 3f893f198dc6b130814b2a2216e606ef */
+/* sin(x) = 0.a400072188acf49cd6b173825e038346f105e1301afe642bcc364cea455e21e506e3e927ed8 */
+  6.40625425040230409188409779413961021e-01L,  /* 3ffe48000e431159e939ad62e704bc07 */
+  2.46542747294664049615806500747173281e-36L,  /* 3f88a37882f0980d7f3215e61b267523 */
+
+/* x = 7.03125000000000000000000000000000000e-01 3ffe6800000000000000000000000000 */
+/* cos(x) = 0.c348846bbd3631338ffe2bfe9dd1381a35b4e9c0c51b4c13fe376bad1bf5caacc4542be0aa9 */
+  7.62825275710576250507098753625429792e-01L,  /* 3ffe869108d77a6c62671ffc57fd3ba2 */
+  4.22067411888601505004748939382325080e-35L,  /* 3f8cc0d1ada74e0628da609ff1bb5d69 */
+/* sin(x) = 0.a587e23555bb08086d02b9c662cdd29316c3e9bd08d93793634a21b1810cce73bdb97a99b9e */
+  6.46604669591152370524042159882800763e-01L,  /* 3ffe4b0fc46aab761010da05738cc59c */
+  -3.41742981816219412415674365946079826e-35L, /* bf8c6b6749e0b217b9364364e5aef274 */
+
+/* x = 7.10937500000000000000000000000000000e-01 3ffe6c00000000000000000000000000 */
+/* cos(x) = 0.c1fbeef380e4ffdd5a613ec8722f643ffe814ec2343e53adb549627224fdc9f2a7b77d3d69f */
+  7.57750448655219342240234832230493361e-01L,  /* 3ffe83f7dde701c9ffbab4c27d90e45f */
+  -2.08767968311222650582659938787920125e-35L, /* bf8bbc0017eb13dcbc1ac524ab69d8de */
+/* sin(x) = 0.a70d272a76a8d4b6da0ec90712bb748b96dabf88c3079246f3db7eea6e58ead4ed0e2843303 */
+  6.52544448725765956407573982284767763e-01L,  /* 3ffe4e1a4e54ed51a96db41d920e2577 */
+  -8.61758060284379660697102362141557170e-36L, /* bf8a6e8d24a80ee79f0db721849022b2 */
+
+/* x = 7.18750000000000000000000000000000000e-01 3ffe7000000000000000000000000000 */
+/* cos(x) = 0.c0ac518c8b6ae710ba37a3eeb90cb15aebcb8bed4356fb507a48a6e97de9aa6d9660116b436 */
+  7.52629372418066476054541324847143116e-01L,  /* 3ffe8158a31916d5ce21746f47dd7219 */
+  3.71306958657663189665450864311104571e-35L,  /* 3f8c8ad75e5c5f6a1ab7da83d245374c */
+/* sin(x) = 0.a88fcfebd9a8dd47e2f3c76ef9e2439920f7e7fbe735f8bcc985491ec6f12a2d4214f8cfa99 */
+  6.58444399910567541589583954884041989e-01L,  /* 3ffe511f9fd7b351ba8fc5e78eddf3c5 */
+  -4.54412944084300330523721391865787219e-35L, /* bf8ce336f840c020c6503a19b3d5b70a */
+
+/* x = 7.26562500000000000000000000000000000e-01 3ffe7400000000000000000000000000 */
+/* cos(x) = 0.bf59b17550a4406875969296567cf3e3b4e483061877c02811c6cae85fad5a6c3da58f49292 */
+  7.47462359563216166669700384714767552e-01L,  /* 3ffe7eb362eaa14880d0eb2d252cacfa */
+  -9.11094340926220027288083639048016945e-36L, /* bf8a8389636f9f3cf107fafdc726a2f4 */
+/* sin(x) = 0.aa0fd66eddb921232c28520d3911b8a03193b47f187f1471ac216fbcd5bb81029294d3a73f1 */
+  6.64304163042946276515506587432846246e-01L,  /* 3ffe541facddbb7242465850a41a7223 */
+  4.26004843895378210155889028714676019e-35L,  /* 3f8cc5018c9da3f8c3f8a38d610b7de7 */
+
+/* x = 7.34375000000000000000000000000000000e-01 3ffe7800000000000000000000000000 */
+/* cos(x) = 0.be0413f84f2a771c614946a88cbf4da1d75a5560243de8f2283fefa0ea4a48468a52d51d8b3 */
+  7.42249725458501306991347253449610537e-01L,  /* 3ffe7c0827f09e54ee38c2928d51197f */
+  -3.78925270049800913539923473871287550e-35L, /* bf8c92f1452d54fede10b86ebe0082f9 */
+/* sin(x) = 0.ab8d34b36acd987210ed343ec65d7e3adc2e7109fce43d55c8d57dfdf55b9e01d2cc1f1b9ec */
+  6.70123380473162894654531583500648495e-01L,  /* 3ffe571a6966d59b30e421da687d8cbb */
+  -1.33165852952743729897634069393684656e-36L, /* bf87c523d18ef6031bc2aa372a82020b */
+
+/* x = 7.42187500000000000000000000000000000e-01 3ffe7c00000000000000000000000000 */
+/* cos(x) = 0.bcab7e6bfb2a14a9b122c574a376bec98ab14808c64a4e731b34047e217611013ac99c0f25d */
+  7.36991788256240741057089385586450844e-01L,  /* 3ffe7956fcd7f654295362458ae946ed */
+  4.72358938637974850573747497460125519e-35L,  /* 3f8cf64c558a404632527398d9a023f1 */
+/* sin(x) = 0.ad07e4c409d08c4fa3a9057bb0ac24b8636e74e76f51e09bd6b2319707cbd9f5e254643897a */
+  6.75901697026178809189642203142423973e-01L,  /* 3ffe5a0fc98813a1189f47520af76158 */
+  2.76252586616364878801928456702948857e-35L,  /* 3f8c25c31b73a73b7a8f04deb5918cb8 */
+
+/* x = 7.50000000000000000000000000000000000e-01 3ffe8000000000000000000000000000 */
+/* cos(x) = 0.bb4ff632a908f73ec151839cb9d993b4e0bfb8f20e7e44e6e4aee845e35575c3106dbe6fd06 */
+  7.31688868873820886311838753000084529e-01L,  /* 3ffe769fec655211ee7d82a3073973b3 */
+  1.48255637548931697184991710293198620e-35L,  /* 3f8b3b4e0bfb8f20e7e44e6e4aee845e */
+/* sin(x) = 0.ae7fe0b5fc786b2d966e1d6af140a488476747c2646425fc7533f532cd044cb10a971a49a6a */
+  6.81638760023334166733241952779893908e-01L,  /* 3ffe5cffc16bf8f0d65b2cdc3ad5e281 */
+  2.74838775935027549024224114338667371e-35L,  /* 3f8c24423b3a3e1323212fe3a99fa996 */
+
+/* x = 7.57812500000000000000000000000000000e-01 3ffe8400000000000000000000000000 */
+/* cos(x) = 0.b9f180ba77dd0751628e135a9508299012230f14becacdd14c3f8862d122de5b56d55b53360 */
+  7.26341290974108590410147630237598973e-01L,  /* 3ffe73e30174efba0ea2c51c26b52a10 */
+  3.12683579338351123545814364980658990e-35L,  /* 3f8c4c80911878a5f6566e8a61fc4317 */
+/* sin(x) = 0.aff522a954f2ba16d9defdc416e33f5e9a5dfd5a6c228e0abc4d521327ff6e2517a7b3851dd */
+  6.87334219303873534951703613035647220e-01L,  /* 3ffe5fea4552a9e5742db3bdfb882dc6 */
+  4.76739454455410744997012795035529128e-35L,  /* 3f8cfaf4d2efead361147055e26a9099 */
+
+/* x = 7.65625000000000000000000000000000000e-01 3ffe8800000000000000000000000000 */
+/* cos(x) = 0.b890237d3bb3c284b614a0539016bfa1053730bbdf940fa895e185f8e58884d3dda15e63371 */
+  7.20949380945696418043812784148447688e-01L,  /* 3ffe712046fa776785096c2940a7202d */
+  4.78691285733673379499536326050811832e-35L,  /* 3f8cfd0829b985defca07d44af0c2fc7 */
+/* sin(x) = 0.b167a4c90d63c4244cf5493b7cc23bd3c3c1225e078baa0c53d6d400b926281f537a1a260e6 */
+  6.92987727246317910281815490823048210e-01L,  /* 3ffe62cf49921ac7884899ea9276f984 */
+  4.50089871077663557180849219529189918e-35L,  /* 3f8cde9e1e0912f03c5d50629eb6a006 */
+
+/* x = 7.73437500000000000000000000000000000e-01 3ffe8c00000000000000000000000000 */
+/* cos(x) = 0.b72be40067aaf2c050dbdb7a14c3d7d4f203f6b3f0224a4afe55d6ec8e92b508fd5c5984b3b */
+  7.15513467882981573520620561289896903e-01L,  /* 3ffe6e57c800cf55e580a1b7b6f42988 */
+  -3.02191815581445336509438104625489192e-35L, /* bf8c41586fe04a607eedada80d51489c */
+/* sin(x) = 0.b2d7614b1f3aaa24df2d6e20a77e1ca3e6d838c03e29c1bcb026e6733324815fadc9eb89674 */
+  6.98598938789681741301929277107891591e-01L,  /* 3ffe65aec2963e755449be5adc414efc */
+  2.15465226809256290914423429408722521e-35L,  /* 3f8bca3e6d838c03e29c1bcb026e6733 */
+
+/* x = 7.81250000000000000000000000000000000e-01 3ffe9000000000000000000000000000 */
+/* cos(x) = 0.b5c4c7d4f7dae915ac786ccf4b1a498d3e73b6e5e74fe7519d9c53ee6d6b90e881bddfc33e1 */
+  7.10033883566079674974121643959490219e-01L,  /* 3ffe6b898fa9efb5d22b58f0d99e9635 */
+  -4.09623224763692443220896752907902465e-35L, /* bf8cb3960c6248d0c580c573131d608d */
+/* sin(x) = 0.b44452709a59752905913765434a59d111f0433eb2b133f7d103207e2aeb4aae111ddc385b3 */
+  7.04167511454533672780059509973942844e-01L,  /* 3ffe6888a4e134b2ea520b226eca8695 */
+  -2.87259372740393348676633610275598640e-35L, /* bf8c3177707de60a6a76604177e6fc0f */
+
+/* x = 7.89062500000000000000000000000000000e-01 3ffe9400000000000000000000000000 */
+/* cos(x) = 0.b45ad4975b1294cadca4cf40ec8f22a68cd14b175835239a37e63acb85e8e9505215df18140 */
+  7.04510962440574606164129481545916976e-01L,  /* 3ffe68b5a92eb6252995b9499e81d91e */
+  2.60682037357042658395360726992048803e-35L,  /* 3f8c1534668a58bac1a91cd1bf31d65c */
+/* sin(x) = 0.b5ae7285bc10cf515753847e8f8b7a30e0a580d929d770103509880680f7b8b0e8ad23b65d8 */
+  7.09693105363899724959669028139035515e-01L,  /* 3ffe6b5ce50b78219ea2aea708fd1f17 */
+  -4.37026016974122945368562319136420097e-36L, /* bf8973c7d69fc9b58a23fbf2bd9dfe60 */
+};
diff --git a/lib/sinl.c b/lib/sinl.c
new file mode 100644 (file)
index 0000000..3bec0bb
--- /dev/null
@@ -0,0 +1,101 @@
+/* s_sinl.c -- long double version of s_sin.c.
+ * Conversion to long double by Jakub Jelinek, jj@ultra.linux.cz.
+ */
+
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+
+/* sinl(x)
+ * Return sine function of x.
+ *
+ * kernel function:
+ *     __kernel_sinl           ... sine function on [-pi/4,pi/4]
+ *     __kernel_cosl           ... cose function on [-pi/4,pi/4]
+ *     __ieee754_rem_pio2l     ... argument reduction routine
+ *
+ * Method.
+ *      Let S,C and T denote the sin, cos and tan respectively on
+ *     [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2
+ *     in [-pi/4 , +pi/4], and let n = k mod 4.
+ *     We have
+ *
+ *          n        sin(x)      cos(x)        tan(x)
+ *     ----------------------------------------------------------
+ *         0          S           C             T
+ *         1          C          -S            -1/T
+ *         2         -S          -C             T
+ *         3         -C           S            -1/T
+ *     ----------------------------------------------------------
+ *
+ * Special cases:
+ *      Let trig be any of sin, cos, or tan.
+ *      trig(+-INF)  is NaN, with signals;
+ *      trig(NaN)    is that NaN;
+ *
+ * Accuracy:
+ *     TRIG(x) returns trig(x) nearly rounded
+ */
+
+#include <math.h>
+
+#include "mathl.h"
+
+#include "trigl.h"
+#include "trigl.c"
+#include "sincosl.c"
+
+long double
+sinl (long double x)
+{
+  long double y[2], z = 0.0L;
+  int n;
+
+  /* |x| ~< pi/4 */
+  if (x >= -0.7853981633974483096156608458198757210492 &&
+      x <= 0.7853981633974483096156608458198757210492)
+    return kernel_sinl (x, z, 0);
+
+    /* sinl(Inf or NaN) is NaN, sinl(0) is 0 */
+  else if (x + x == x || x != x)
+    return x - x;              /* NaN */
+
+  /* argument reduction needed */
+  else
+    {
+      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]);
+       }
+    }
+}
+
+#if 0
+int
+main ()
+{
+  printf ("%.16Lg\n", sinl(0.7853981633974483096156608458198757210492));
+  printf ("%.16Lg\n", sinl(0.7853981633974483096156608458198757210492 *29));
+  printf ("%.16Lg\n", sinl(0.7853981633974483096156608458198757210492 *2));
+  printf ("%.16Lg\n", sinl(0.7853981633974483096156608458198757210492 *30));
+  printf ("%.16Lg\n", sinl(0.7853981633974483096156608458198757210492 *4));
+  printf ("%.16Lg\n", sinl(0.7853981633974483096156608458198757210492 *32));
+  printf ("%.16Lg\n", sinl(0.7853981633974483096156608458198757210492 *2/3));
+  printf ("%.16Lg\n", sinl(0.7853981633974483096156608458198757210492 *4/3));
+}
+#endif
diff --git a/lib/sqrtl.c b/lib/sqrtl.c
new file mode 100644 (file)
index 0000000..ed914a2
--- /dev/null
@@ -0,0 +1,56 @@
+/* Emulation for sqrtl.
+   Contributed by Paolo Bonzini
+
+   Copyright 2002, 2003 Free Software Foundation, Inc.
+
+   This file is part of gnulib.
+
+   gnulib is free software; you can redistribute it and/or modify it
+   under the terms of the GNU Lesser General Public License as published
+   by the Free Software Foundation; either version 2.1, or (at your option)
+   any later version.
+
+   gnulib is distributed in the hope that it will be useful, but WITHOUT
+   ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+   FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
+   License for more details.
+
+   You should have received a copy of the GNU Lesser General Public License
+   along with gnulib; see the file COPYING.LIB.  If not, write to the Free
+   Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307,
+   USA.
+ */
+
+#include <float.h>
+#include <math.h>
+
+#include "mathl.h"
+
+/* A simple Newton-Raphson method. */
+long double
+sqrtl(long double x)
+{
+  long double delta, y;
+  int exponent;
+
+  /* Check for negative numbers */
+  if (x < 0.0L)
+    return (long double) sqrt(-1);
+
+  /* Check for zero, NANs and infinites */
+  if (x + x == x || x != x)
+    return x;
+
+  frexpl (x, &exponent);
+  y = ldexpl (x, -exponent / 2);
+
+  do
+    {
+      delta = y;
+      y = (y + x / y) * 0.5L;
+      delta -= y;
+    }
+  while (delta != 0.0L);
+
+  return y;
+}
diff --git a/lib/tanl.c b/lib/tanl.c
new file mode 100644 (file)
index 0000000..e4b974e
--- /dev/null
@@ -0,0 +1,223 @@
+/* s_tanl.c -- long double version of s_tan.c.
+ * Conversion to IEEE quad long double by Jakub Jelinek, jj@ultra.linux.cz.
+ */
+
+/* @(#)s_tan.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+
+/* tanl(x)
+ * Return tangent function of x.
+ *
+ * kernel function:
+ *     __kernel_tanl           ... tangent function on [-pi/4,pi/4]
+ *     __ieee754_rem_pio2l     ... argument reduction routine
+ *
+ * Method.
+ *      Let S,C and T denote the sin, cos and tan respectively on
+ *     [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2
+ *     in [-pi/4 , +pi/4], and let n = k mod 4.
+ *     We have
+ *
+ *          n        sin(x)      cos(x)        tan(x)
+ *     ----------------------------------------------------------
+ *         0          S           C             T
+ *         1          C          -S            -1/T
+ *         2         -S          -C             T
+ *         3         -C           S            -1/T
+ *     ----------------------------------------------------------
+ *
+ * Special cases:
+ *      Let trig be any of sin, cos, or tan.
+ *      trig(+-INF)  is NaN, with signals;
+ *      trig(NaN)    is that NaN;
+ *
+ * Accuracy:
+ *     TRIG(x) returns trig(x) nearly rounded
+ */
+
+#include <math.h>
+
+#include "mathl.h"
+
+#include "trigl.h"
+#ifdef HAVE_SINL
+#ifdef HAVE_COSL
+#include "trigl.c"
+#endif
+#endif
+
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+
+/*
+  Long double expansions contributed by
+  Stephen L. Moshier <moshier@na-net.ornl.gov>
+*/
+
+/* __kernel_tanl( x, y, k )
+ * kernel tan function on [-pi/4, pi/4], pi/4 ~ 0.7854
+ * Input x is assumed to be bounded by ~pi/4 in magnitude.
+ * Input y is the tail of x.
+ * Input k indicates whether tan (if k=1) or
+ * -1/tan (if k= -1) is returned.
+ *
+ * Algorithm
+ *     1. Since tan(-x) = -tan(x), we need only to consider positive x.
+ *     2. if x < 2^-57, return x with inexact if x!=0.
+ *     3. tan(x) is approximated by a rational form x + x^3 / 3 + x^5 R(x^2)
+ *          on [0,0.67433].
+ *
+ *        Note: tan(x+y) = tan(x) + tan'(x)*y
+ *                       ~ tan(x) + (1+x*x)*y
+ *        Therefore, for better accuracy in computing tan(x+y), let
+ *             r = x^3 * R(x^2)
+ *        then
+ *             tan(x+y) = x + (x^3 / 3 + (x^2 *(r+y)+y))
+ *
+ *      4. For x in [0.67433,pi/4],  let y = pi/4 - x, then
+ *             tan(x) = tan(pi/4-y) = (1-tan(y))/(1+tan(y))
+ *                    = 1 - 2*(tan(y) - (tan(y)^2)/(1+tan(y)))
+ */
+
+
+static const long double
+  pio4hi = 7.8539816339744830961566084581987569936977E-1L,
+  pio4lo = 2.1679525325309452561992610065108379921906E-35L,
+
+  /* tan x = x + x^3 / 3 + x^5 T(x^2)/U(x^2)
+     0 <= x <= 0.6743316650390625
+     Peak relative error 8.0e-36  */
+ TH =  3.333333333333333333333333333333333333333E-1L,
+ T0 = -1.813014711743583437742363284336855889393E7L,
+ T1 =  1.320767960008972224312740075083259247618E6L,
+ T2 = -2.626775478255838182468651821863299023956E4L,
+ T3 =  1.764573356488504935415411383687150199315E2L,
+ T4 = -3.333267763822178690794678978979803526092E-1L,
+
+ U0 = -1.359761033807687578306772463253710042010E8L,
+ U1 =  6.494370630656893175666729313065113194784E7L,
+ U2 = -4.180787672237927475505536849168729386782E6L,
+ U3 =  8.031643765106170040139966622980914621521E4L,
+ U4 = -5.323131271912475695157127875560667378597E2L;
+  /* 1.000000000000000000000000000000000000000E0 */
+
+
+long double
+kernel_tanl (long double x, long double y, int iy)
+{
+  long double z, r, v, w, s, u, u1;
+  int flag, sign;
+
+       sign = 1;
+      if (x < 0)
+       {
+         x = -x;
+         y = -y;
+         sign = -1;
+       }
+
+  if (x < 0.000000000000000006938893903907228377647697925567626953125L) /* x < 2**-57 */
+    {
+      if ((int) x == 0)
+       {                       /* generate inexact */
+         if (iy == -1 && x == 0.0)
+           return 1.0L / fabs (x);
+         else
+           return (iy == 1) ? x : -1.0L / x;
+       }
+    }
+  if (x >= 0.6743316650390625) /* |x| >= 0.6743316650390625 */
+    {
+      flag = 1;
+
+      z = pio4hi - x;
+      w = pio4lo - y;
+      x = z + w;
+      y = 0.0;
+    }
+  z = x * x;
+  r = T0 + z * (T1 + z * (T2 + z * (T3 + z * T4)));
+  v = U0 + z * (U1 + z * (U2 + z * (U3 + z * (U4 + z))));
+  r = r / v;
+
+  s = z * x;
+  r = y + z * (s * r + y);
+  r += TH * s;
+  w = x + r;
+  if (flag)
+    {
+      v = (long double) iy;
+      w = (v - 2.0 * (x - (w * w / (w + v) - r)));
+      if (sign < 0)
+       w = -w;
+      return w;
+    }
+  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 */
+      u1 = (double) w;
+      v = r - (u1 - x);
+      z = -1.0 / w;
+      u = (double) z;
+      s = 1.0 + u * u1;
+      return u + z * (s + u * v);
+    }
+}
+
+long double
+tanl (long double x)
+{
+  long double y[2], z = 0.0L;
+  int n;
+
+  /* |x| ~< pi/4 */
+  if (x >= -0.7853981633974483096156608458198757210492 &&
+      x <= 0.7853981633974483096156608458198757210492)
+    return kernel_tanl (x, z, 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 */
+  else
+    {
+      n = ieee754_rem_pio2l (x, y);
+      /* 1 -- n even, -1 -- n odd */
+      return kernel_tanl (y[0], y[1], 1 - ((n & 1) << 1));
+    }
+}
+
+#if 0
+int
+main ()
+{
+  printf ("%.16Lg\n", tanl(0.7853981633974483096156608458198757210492));
+  printf ("%.16Lg\n", tanl(-0.7853981633974483096156608458198757210492));
+  printf ("%.16Lg\n", tanl(0.7853981633974483096156608458198757210492 *3));
+  printf ("%.16Lg\n", tanl(-0.7853981633974483096156608458198757210492 *31));
+  printf ("%.16Lg\n", tanl(0.7853981633974483096156608458198757210492 / 2));
+  printf ("%.16Lg\n", tanl(0.7853981633974483096156608458198757210492 * 3/2));
+  printf ("%.16Lg\n", tanl(0.7853981633974483096156608458198757210492 * 5/2));
+}
+#endif
diff --git a/lib/trigl.c b/lib/trigl.c
new file mode 100644 (file)
index 0000000..4e33f78
--- /dev/null
@@ -0,0 +1,636 @@
+/* Quad-precision floating point argument reduction.
+   Copyright (C) 1999 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+   Contributed by Jakub Jelinek <jj@ultra.linux.cz>
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+
+   The GNU C Library is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, write to the Free
+   Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
+   02111-1307 USA.  */
+
+#include <math.h>
+#include <float.h>
+
+#include "mathl.h"
+
+/* Table of constants for 2/pi, 5628 hexadecimal digits of 2/pi */
+static const int 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,
+  0x47c419, 0xc367cd, 0xdce809, 0x2a8359, 0xc4768b, 0x961ca6,
+  0xddaf44, 0xd15719, 0x053ea5, 0xff0705, 0x3f7e33, 0xe832c2,
+  0xde4f98, 0x327dbb, 0xc33d26, 0xef6b1e, 0x5ef89f, 0x3a1f35,
+  0xcaf27f, 0x1d87f1, 0x21907c, 0x7c246a, 0xfa6ed5, 0x772d30,
+  0x433b15, 0xc614b5, 0x9d19c3, 0xc2c4ad, 0x414d2c, 0x5d000c,
+  0x467d86, 0x2d71e3, 0x9ac69b, 0x006233, 0x7cd2b4, 0x97a7b4,
+  0xd55537, 0xf63ed7, 0x1810a3, 0xfc764d, 0x2a9d64, 0xabd770,
+  0xf87c63, 0x57b07a, 0xe71517, 0x5649c0, 0xd9d63b, 0x3884a7,
+  0xcb2324, 0x778ad6, 0x23545a, 0xb91f00, 0x1b0af1, 0xdfce19,
+  0xff319f, 0x6a1e66, 0x615799, 0x47fbac, 0xd87f7e, 0xb76522,
+  0x89e832, 0x60bfe6, 0xcdc4ef, 0x09366c, 0xd43f5d, 0xd7de16,
+  0xde3b58, 0x929bde, 0x2822d2, 0xe88628, 0x4d58e2, 0x32cac6,
+  0x16e308, 0xcb7de0, 0x50c017, 0xa71df3, 0x5be018, 0x34132e,
+  0x621283, 0x014883, 0x5b8ef5, 0x7fb0ad, 0xf2e91e, 0x434a48,
+  0xd36710, 0xd8ddaa, 0x425fae, 0xce616a, 0xa4280a, 0xb499d3,
+  0xf2a606, 0x7f775c, 0x83c2a3, 0x883c61, 0x78738a, 0x5a8caf,
+  0xbdd76f, 0x63a62d, 0xcbbff4, 0xef818d, 0x67c126, 0x45ca55,
+  0x36d9ca, 0xd2a828, 0x8d61c2, 0x77c912, 0x142604, 0x9b4612,
+  0xc459c4, 0x44c5c8, 0x91b24d, 0xf31700, 0xad43d4, 0xe54929,
+  0x10d5fd, 0xfcbe00, 0xcc941e, 0xeece70, 0xf53e13, 0x80f1ec,
+  0xc3e7b3, 0x28f8c7, 0x940593, 0x3e71c1, 0xb3092e, 0xf3450b,
+  0x9c1288, 0x7b20ab, 0x9fb52e, 0xc29247, 0x2f327b, 0x6d550c,
+  0x90a772, 0x1fe76b, 0x96cb31, 0x4a1679, 0xe27941, 0x89dff4,
+  0x9794e8, 0x84e6e2, 0x973199, 0x6bed88, 0x365f5f, 0x0efdbb,
+  0xb49a48, 0x6ca467, 0x427271, 0x325d8d, 0xb8159f, 0x09e5bc,
+  0x25318d, 0x3974f7, 0x1c0530, 0x010c0d, 0x68084b, 0x58ee2c,
+  0x90aa47, 0x02e774, 0x24d6bd, 0xa67df7, 0x72486e, 0xef169f,
+  0xa6948e, 0xf691b4, 0x5153d1, 0xf20acf, 0x339820, 0x7e4bf5,
+  0x6863b2, 0x5f3edd, 0x035d40, 0x7f8985, 0x295255, 0xc06437,
+  0x10d86d, 0x324832, 0x754c5b, 0xd4714e, 0x6e5445, 0xc1090b,
+  0x69f52a, 0xd56614, 0x9d0727, 0x50045d, 0xdb3bb4, 0xc576ea,
+  0x17f987, 0x7d6b49, 0xba271d, 0x296996, 0xacccc6, 0x5414ad,
+  0x6ae290, 0x89d988, 0x50722c, 0xbea404, 0x940777, 0x7030f3,
+  0x27fc00, 0xa871ea, 0x49c266, 0x3de064, 0x83dd97, 0x973fa3,
+  0xfd9443, 0x8c860d, 0xde4131, 0x9d3992, 0x8c70dd, 0xe7b717,
+  0x3bdf08, 0x2b3715, 0xa0805c, 0x93805a, 0x921110, 0xd8e80f,
+  0xaf806c, 0x4bffdb, 0x0f9038, 0x761859, 0x15a562, 0xbbcb61,
+  0xb989c7, 0xbd4010, 0x04f2d2, 0x277549, 0xf6b6eb, 0xbb22db,
+  0xaa140a, 0x2f2689, 0x768364, 0x333b09, 0x1a940e, 0xaa3a51,
+  0xc2a31d, 0xaeedaf, 0x12265c, 0x4dc26d, 0x9c7a2d, 0x9756c0,
+  0x833f03, 0xf6f009, 0x8c402b, 0x99316d, 0x07b439, 0x15200c,
+  0x5bc3d8, 0xc492f5, 0x4badc6, 0xa5ca4e, 0xcd37a7, 0x36a9e6,
+  0x9492ab, 0x6842dd, 0xde6319, 0xef8c76, 0x528b68, 0x37dbfc,
+  0xaba1ae, 0x3115df, 0xa1ae00, 0xdafb0c, 0x664d64, 0xb705ed,
+  0x306529, 0xbf5657, 0x3aff47, 0xb9f96a, 0xf3be75, 0xdf9328,
+  0x3080ab, 0xf68c66, 0x15cb04, 0x0622fa, 0x1de4d9, 0xa4b33d,
+  0x8f1b57, 0x09cd36, 0xe9424e, 0xa4be13, 0xb52333, 0x1aaaf0,
+  0xa8654f, 0xa5c1d2, 0x0f3f0b, 0xcd785b, 0x76f923, 0x048b7b,
+  0x721789, 0x53a6c6, 0xe26e6f, 0x00ebef, 0x584a9b, 0xb7dac4,
+  0xba66aa, 0xcfcf76, 0x1d02d1, 0x2df1b1, 0xc1998c, 0x77adc3,
+  0xda4886, 0xa05df7, 0xf480c6, 0x2ff0ac, 0x9aecdd, 0xbc5c3f,
+  0x6dded0, 0x1fc790, 0xb6db2a, 0x3a25a3, 0x9aaf00, 0x9353ad,
+  0x0457b6, 0xb42d29, 0x7e804b, 0xa707da, 0x0eaa76, 0xa1597b,
+  0x2a1216, 0x2db7dc, 0xfde5fa, 0xfedb89, 0xfdbe89, 0x6c76e4,
+  0xfca906, 0x70803e, 0x156e85, 0xff87fd, 0x073e28, 0x336761,
+  0x86182a, 0xeabd4d, 0xafe7b3, 0x6e6d8f, 0x396795, 0x5bbf31,
+  0x48d784, 0x16df30, 0x432dc7, 0x356125, 0xce70c9, 0xb8cb30,
+  0xfd6cbf, 0xa200a4, 0xe46c05, 0xa0dd5a, 0x476f21, 0xd21262,
+  0x845cb9, 0x496170, 0xe0566b, 0x015299, 0x375550, 0xb7d51e,
+  0xc4f133, 0x5f6e13, 0xe4305d, 0xa92e85, 0xc3b21d, 0x3632a1,
+  0xa4b708, 0xd4b1ea, 0x21f716, 0xe4698f, 0x77ff27, 0x80030c,
+  0x2d408d, 0xa0cd4f, 0x99a520, 0xd3a2b3, 0x0a5d2f, 0x42f9b4,
+  0xcbda11, 0xd0be7d, 0xc1db9b, 0xbd17ab, 0x81a2ca, 0x5c6a08,
+  0x17552e, 0x550027, 0xf0147f, 0x8607e1, 0x640b14, 0x8d4196,
+  0xdebe87, 0x2afdda, 0xb6256b, 0x34897b, 0xfef305, 0x9ebfb9,
+  0x4f6a68, 0xa82a4a, 0x5ac44f, 0xbcf82d, 0x985ad7, 0x95c7f4,
+  0x8d4d0d, 0xa63a20, 0x5f57a4, 0xb13f14, 0x953880, 0x0120cc,
+  0x86dd71, 0xb6dec9, 0xf560bf, 0x11654d, 0x6b0701, 0xacb08c,
+  0xd0c0b2, 0x485551, 0x0efb1e, 0xc37295, 0x3b06a3, 0x3540c0,
+  0x7bdc06, 0xcc45e0, 0xfa294e, 0xc8cad6, 0x41f3e8, 0xde647c,
+  0xd8649b, 0x31bed9, 0xc397a4, 0xd45877, 0xc5e369, 0x13daf0,
+  0x3c3aba, 0x461846, 0x5f7555, 0xf5bdd2, 0xc6926e, 0x5d2eac,
+  0xed440e, 0x423e1c, 0x87c461, 0xe9fd29, 0xf3d6e7, 0xca7c22,
+  0x35916f, 0xc5e008, 0x8dd7ff, 0xe26a6e, 0xc6fdb0, 0xc10893,
+  0x745d7c, 0xb2ad6b, 0x9d6ecd, 0x7b723e, 0x6a11c6, 0xa9cff7,
+  0xdf7329, 0xbac9b5, 0x5100b7, 0x0db2e2, 0x24ba74, 0x607de5,
+  0x8ad874, 0x2c150d, 0x0c1881, 0x94667e, 0x162901, 0x767a9f,
+  0xbefdfd, 0xef4556, 0x367ed9, 0x13d9ec, 0xb9ba8b, 0xfc97c4,
+  0x27a831, 0xc36ef1, 0x36c594, 0x56a8d8, 0xb5a8b4, 0x0ecccf,
+  0x2d8912, 0x34576f, 0x89562c, 0xe3ce99, 0xb920d6, 0xaa5e6b,
+  0x9c2a3e, 0xcc5f11, 0x4a0bfd, 0xfbf4e1, 0x6d3b8e, 0x2c86e2,
+  0x84d4e9, 0xa9b4fc, 0xd1eeef, 0xc9352e, 0x61392f, 0x442138,
+  0xc8d91b, 0x0afc81, 0x6a4afb, 0xd81c2f, 0x84b453, 0x8c994e,
+  0xcc2254, 0xdc552a, 0xd6c6c0, 0x96190b, 0xb8701a, 0x649569,
+  0x605a26, 0xee523f, 0x0f117f, 0x11b5f4, 0xf5cbfc, 0x2dbc34,
+  0xeebc34, 0xcc5de8, 0x605edd, 0x9b8e67, 0xef3392, 0xb817c9,
+  0x9b5861, 0xbc57e1, 0xc68351, 0x103ed8, 0x4871dd, 0xdd1c2d,
+  0xa118af, 0x462c21, 0xd7f359, 0x987ad9, 0xc0549e, 0xfa864f,
+  0xfc0656, 0xae79e5, 0x362289, 0x22ad38, 0xdc9367, 0xaae855,
+  0x382682, 0x9be7ca, 0xa40d51, 0xb13399, 0x0ed7a9, 0x480569,
+  0xf0b265, 0xa7887f, 0x974c88, 0x36d1f9, 0xb39221, 0x4a827b,
+  0x21cf98, 0xdc9f40, 0x5547dc, 0x3a74e1, 0x42eb67, 0xdf9dfe,
+  0x5fd45e, 0xa4677b, 0x7aacba, 0xa2f655, 0x23882b, 0x55ba41,
+  0x086e59, 0x862a21, 0x834739, 0xe6e389, 0xd49ee5, 0x40fb49,
+  0xe956ff, 0xca0f1c, 0x8a59c5, 0x2bfa94, 0xc5c1d3, 0xcfc50f,
+  0xae5adb, 0x86c547, 0x624385, 0x3b8621, 0x94792c, 0x876110,
+  0x7b4c2a, 0x1a2c80, 0x12bf43, 0x902688, 0x893c78, 0xe4c4a8,
+  0x7bdbe5, 0xc23ac4, 0xeaf426, 0x8a67f7, 0xbf920d, 0x2ba365,
+  0xb1933d, 0x0b7cbd, 0xdc51a4, 0x63dd27, 0xdde169, 0x19949a,
+  0x9529a8, 0x28ce68, 0xb4ed09, 0x209f44, 0xca984e, 0x638270,
+  0x237c7e, 0x32b90f, 0x8ef5a7, 0xe75614, 0x08f121, 0x2a9db5,
+  0x4d7e6f, 0x5119a5, 0xabf9b5, 0xd6df82, 0x61dd96, 0x023616,
+  0x9f3ac4, 0xa1a283, 0x6ded72, 0x7a8d39, 0xa9b882, 0x5c326b,
+  0x5b2746, 0xed3400, 0x7700d2, 0x55f4fc, 0x4d5901, 0x8071e0,
+  0xe13f89, 0xb295f3, 0x64a8f1, 0xaea74b, 0x38fc4c, 0xeab2bb,
+  0x47270b, 0xabc3a7, 0x34ba60, 0x52dd34, 0xf8563a, 0xeb7e8a,
+  0x31bb36, 0x5895b7, 0x47f7a9, 0x94c3aa, 0xd39225, 0x1e7f3e,
+  0xd8974e, 0xbba94f, 0xd8ae01, 0xe661b4, 0x393d8e, 0xa523aa,
+  0x33068e, 0x1633b5, 0x3bb188, 0x1d3a9d, 0x4013d0, 0xcc1be5,
+  0xf862e7, 0x3bf28f, 0x39b5bf, 0x0bc235, 0x22747e, 0xa247c0,
+  0xd52d1f, 0x19add3, 0x9094df, 0x9311d0, 0xb42b25, 0x496db2,
+  0xe264b2, 0x5ef135, 0x3bc6a4, 0x1a4ad0, 0xaac92e, 0x64e886,
+  0x573091, 0x982cfb, 0x311b1a, 0x08728b, 0xbdcee1, 0x60e142,
+  0xeb641d, 0xd0bba3, 0xe559d4, 0x597b8c, 0x2a4483, 0xf332ba,
+  0xf84867, 0x2c8d1b, 0x2fa9b0, 0x50f3dd, 0xf9f573, 0xdb61b4,
+  0xfe233e, 0x6c41a6, 0xeea318, 0x775a26, 0xbc5e5c, 0xcea708,
+  0x94dc57, 0xe20196, 0xf1e839, 0xbe4851, 0x5d2d2f, 0x4e9555,
+  0xd96ec2, 0xe7d755, 0x6304e0, 0xc02e0e, 0xfc40a0, 0xbbf9b3,
+  0x7125a7, 0x222dfb, 0xf619d8, 0x838c1c, 0x6619e6, 0xb20d55,
+  0xbb5137, 0x79e809, 0xaf9149, 0x0d73de, 0x0b0da5, 0xce7f58,
+  0xac1934, 0x724667, 0x7a1a13, 0x9e26bc, 0x4555e7, 0x585cb5,
+  0x711d14, 0x486991, 0x480d60, 0x56adab, 0xd62f64, 0x96ee0c,
+  0x212ff3, 0x5d6d88, 0xa67684, 0x95651e, 0xab9e0a, 0x4ddefe,
+  0x571010, 0x836a39, 0xf8ea31, 0x9e381d, 0xeac8b1, 0xcac96b,
+  0x37f21e, 0xd505e9, 0x984743, 0x9fc56c, 0x0331b7, 0x3b8bf8,
+  0x86e56a, 0x8dc343, 0x6230e7, 0x93cfd5, 0x6a8f2d, 0x733005,
+  0x1af021, 0xa09fcb, 0x7415a1, 0xd56b23, 0x6ff725, 0x2f4bc7,
+  0xb8a591, 0x7fac59, 0x5c55de, 0x212c38, 0xb13296, 0x5cff50,
+  0x366262, 0xfa7b16, 0xf4d9a6, 0x2acfe7, 0xf07403, 0xd4d604,
+  0x6fd916, 0x31b1bf, 0xcbb450, 0x5bd7c8, 0x0ce194, 0x6bd643,
+  0x4fd91c, 0xdf4543, 0x5f3453, 0xe2b5aa, 0xc9aec8, 0x131485,
+  0xf9d2bf, 0xbadb9e, 0x76f5b9, 0xaf15cf, 0xca3182, 0x14b56d,
+  0xe9fe4d, 0x50fc35, 0xf5aed5, 0xa2d0c1, 0xc96057, 0x192eb6,
+  0xe91d92, 0x07d144, 0xaea3c6, 0x343566, 0x26d5b4, 0x3161e2,
+  0x37f1a2, 0x209eff, 0x958e23, 0x493798, 0x35f4a6, 0x4bdc02,
+  0xc2be13, 0xbe80a0, 0x0b72a3, 0x115c5f, 0x1e1bd1, 0x0db4d3,
+  0x869e85, 0x96976b, 0x2ac91f, 0x8a26c2, 0x3070f0, 0x041412,
+  0xfc9fa5, 0xf72a38, 0x9c6878, 0xe2aa76, 0x50cfe1, 0x559274,
+  0x934e38, 0x0a92f7, 0x5533f0, 0xa63db4, 0x399971, 0xe2b755,
+  0xa98a7c, 0x008f19, 0xac54d2, 0x2ea0b4, 0xf5f3e0, 0x60c849,
+  0xffd269, 0xae52ce, 0x7a5fdd, 0xe9ce06, 0xfb0ae8, 0xa50cce,
+  0xea9d3e, 0x3766dd, 0xb834f5, 0x0da090, 0x846f88, 0x4ae3d5,
+  0x099a03, 0x2eae2d, 0xfcb40a, 0xfb9b33, 0xe281dd, 0x1b16ba,
+  0xd8c0af, 0xd96b97, 0xb52dc9, 0x9c277f, 0x5951d5, 0x21ccd6,
+  0xb6496b, 0x584562, 0xb3baf2, 0xa1a5c4, 0x7ca2cf, 0xa9b93d,
+  0x7b7b89, 0x483d38,
+};
+
+static const long double c[] = {
+/* 93 bits of pi/2 */
+#define PI_2_1 c[0]
+  1.57079632679489661923132169155131424e+00L,  /* 3fff921fb54442d18469898cc5100000 */
+
+/* pi/2 - PI_2_1 */
+#define PI_2_1t c[1]
+  8.84372056613570112025531863263659260e-29L,  /* 3fa1c06e0e68948127044533e63a0106 */
+};
+
+static int kernel_rem_pio2 (double *x, double *y, int e0, int nx, int prec,
+                           const int *ipio2);
+
+int
+ieee754_rem_pio2l (long double x, long double *y)
+{
+  long double z, w, t;
+  double tx[8];
+  int exp, n;
+
+  if (x >= -0.78539816339744830961566084581987572104929234984377
+      && x < 0.78539816339744830961566084581987572104929234984377)
+    /* x in <-pi/4, pi/4> */
+    {
+      y[0] = x;
+      y[1] = 0;
+      return 0;
+    }
+
+  if (x >= 2.35619449019234492884698253745962716314787704953131
+      && x < 2.35619449019234492884698253745962716314787704953131)
+    if (x > 0)
+      {
+       /* 113 + 93 bit PI is ok */
+       z = x - PI_2_1;
+       y[0] = z - PI_2_1t;
+       y[1] = (z - y[0]) - PI_2_1t;
+       return 1;
+      }
+    else
+      {
+       /* 113 + 93 bit PI is ok */
+       z = x + PI_2_1;
+       y[0] = z + PI_2_1t;
+       y[1] = (z - y[0]) + PI_2_1t;
+       return -1;
+      }
+
+  if (x + x == x || x != x)    /* x is +=oo or NaN */
+    {
+      y[0] = x - x;
+      y[1] = y[0];
+      return 0;
+    }
+
+  /* Handle large arguments.
+     We split the 113 bits of the mantissa into 5 24bit integers
+     stored in a double array.  */
+  z = frexp (x, &exp);
+  tx[0] = floorl (z *= 16777216.0);
+  z -= tx[0];
+  tx[1] = floorl (z *= 16777216.0);
+  z -= tx[1];
+  tx[2] = floorl (z *= 16777216.0);
+  z -= tx[2];
+  tx[3] = floorl (z *= 16777216.0);
+  z -= tx[3];
+  tx[4] = floorl (z *= 16777216.0);
+
+  n = kernel_rem_pio2 (tx, tx + 5, exp - 24, tx[4] ? 5 : 4, 3, two_over_pi);
+
+  /* The result is now stored in 3 double values, we need to convert it into
+     two long double values.  */
+  t = (long double) tx[6] + (long double) tx[7];
+  w = (long double) tx[5];
+
+  if (x > 0)
+    {
+      y[0] = w + t;
+      y[1] = t - (y[0] - w);
+      return n;
+    }
+  else
+    {
+      y[0] = -(w + t);
+      y[1] = -t - (y[0] + w);
+      return -n;
+    }
+}
+
+/* @(#)k_rem_pio2.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+
+#if defined(LIBM_SCCS) && !defined(lint)
+static char rcsid[] =
+  "$NetBSD: k_rem_pio2.c,v 1.7 1995/05/10 20:46:25 jtc Exp $";
+#endif
+
+/*
+ * kernel_rem_pio2(x,y,e0,nx,prec,ipio2)
+ * double x[],y[]; int e0,nx,prec; int ipio2[];
+ *
+ * kernel_rem_pio2 return the last three digits of N with
+ *             y = x - N*pi/2
+ * so that |y| < pi/2.
+ *
+ * The method is to compute the integer (mod 8) and fraction parts of
+ * (2/pi)*x without doing the full multiplication. In general we
+ * skip the part of the product that are known to be a huge integer (
+ * more accurately, = 0 mod 8 ). Thus the number of operations are
+ * independent of the exponent of the input.
+ *
+ * (2/pi) is represented by an array of 24-bit integers in ipio2[].
+ *
+ * Input parameters:
+ *     x[]     The input value (must be positive) is broken into nx
+ *             pieces of 24-bit integers in double precision format.
+ *             x[i] will be the i-th 24 bit of x. The scaled exponent
+ *             of x[0] is given in input parameter e0 (i.e., x[0]*2^e0
+ *             match x's up to 24 bits.
+ *
+ *             Example of breaking a double positive z into x[0]+x[1]+x[2]:
+ *                     e0 = ilogb(z)-23
+ *                     z  = scalbn(z,-e0)
+ *             for i = 0,1,2
+ *                     x[i] = floor(z)
+ *                     z    = (z-x[i])*2**24
+ *
+ *
+ *     y[]     ouput result in an array of double precision numbers.
+ *             The dimension of y[] is:
+ *                     24-bit  precision       1
+ *                     53-bit  precision       2
+ *                     64-bit  precision       2
+ *                     113-bit precision       3
+ *             The actual value is the sum of them. Thus for 113-bit
+ *             precision, one may have to do something like:
+ *
+ *             long double t,w,r_head, r_tail;
+ *             t = (long double)y[2] + (long double)y[1];
+ *             w = (long double)y[0];
+ *             r_head = t+w;
+ *             r_tail = w - (r_head - t);
+ *
+ *     e0      The exponent of x[0]
+ *
+ *     nx      dimension of x[]
+ *
+ *     prec    an integer indicating the precision:
+ *                     0       24  bits (single)
+ *                     1       53  bits (double)
+ *                     2       64  bits (extended)
+ *                     3       113 bits (quad)
+ *
+ *     ipio2[]
+ *             integer array, contains the (24*i)-th to (24*i+23)-th
+ *             bit of 2/pi after binary point. The corresponding
+ *             floating value is
+ *
+ *                     ipio2[i] * 2^(-24(i+1)).
+ *
+ * External function:
+ *     double scalbn(), floor();
+ *
+ *
+ * Here is the description of some local variables:
+ *
+ *     jk      jk+1 is the initial number of terms of ipio2[] needed
+ *             in the computation. The recommended value is 2,3,4,
+ *             6 for single, double, extended,and quad.
+ *
+ *     jz      local integer variable indicating the number of
+ *             terms of ipio2[] used.
+ *
+ *     jx      nx - 1
+ *
+ *     jv      index for pointing to the suitable ipio2[] for the
+ *             computation. In general, we want
+ *                     ( 2^e0*x[0] * ipio2[jv-1]*2^(-24jv) )/8
+ *             is an integer. Thus
+ *                     e0-3-24*jv >= 0 or (e0-3)/24 >= jv
+ *             Hence jv = max(0,(e0-3)/24).
+ *
+ *     jp      jp+1 is the number of terms in PIo2[] needed, jp = jk.
+ *
+ *     q[]     double array with integral value, representing the
+ *             24-bits chunk of the product of x and 2/pi.
+ *
+ *     q0      the corresponding exponent of q[0]. Note that the
+ *             exponent for q[i] would be q0-24*i.
+ *
+ *     PIo2[]  double precision array, obtained by cutting pi/2
+ *             into 24 bits chunks.
+ *
+ *     f[]     ipio2[] in floating point
+ *
+ *     iq[]    integer array by breaking up q[] in 24-bits chunk.
+ *
+ *     fq[]    final product of x*(2/pi) in fq[0],..,fq[jk]
+ *
+ *     ih      integer. If >0 it indicates q[] is >= 0.5, hence
+ *             it also indicates the *sign* of the result.
+ *
+ */
+
+
+/*
+ * Constants:
+ * The hexadecimal values are the intended ones for the following
+ * constants. The decimal values may be used, provided that the
+ * compiler will convert from decimal to binary accurately enough
+ * to produce the hexadecimal values shown.
+ */
+
+static const int init_jk[] = { 2, 3, 4, 6 };   /* initial value for jk */
+static const double PIo2[] = {
+  1.57079625129699707031e+00,  /* 0x3FF921FB, 0x40000000 */
+  7.54978941586159635335e-08,  /* 0x3E74442D, 0x00000000 */
+  5.39030252995776476554e-15,  /* 0x3CF84698, 0x80000000 */
+  3.28200341580791294123e-22,  /* 0x3B78CC51, 0x60000000 */
+  1.27065575308067607349e-29,  /* 0x39F01B83, 0x80000000 */
+  1.22933308981111328932e-36,  /* 0x387A2520, 0x40000000 */
+  2.73370053816464559624e-44,  /* 0x36E38222, 0x80000000 */
+  2.16741683877804819444e-51,  /* 0x3569F31D, 0x00000000 */
+};
+
+static const double zero = 0.0, one = 1.0, two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
+  twon24 = 5.96046447753906250000e-08; /* 0x3E700000, 0x00000000 */
+
+int
+kernel_rem_pio2 (double *x, double *y, int e0, int nx, int prec,
+                const int *ipio2)
+{
+  int jz, jx, jv, jp, jk, carry, n, iq[20], i, j, k, m, q0, ih;
+  double z, fw, f[20], fq[20], q[20];
+
+  /* initialize jk */
+  jk = init_jk[prec];
+  jp = jk;
+
+  /* determine jx,jv,q0, note that 3>q0 */
+  jx = nx - 1;
+  jv = (e0 - 3) / 24;
+  if (jv < 0)
+    jv = 0;
+  q0 = e0 - 24 * (jv + 1);
+
+  /* set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk] */
+  j = jv - jx;
+  m = jx + jk;
+  for (i = 0; i <= m; i++, j++)
+    f[i] = (j < 0) ? zero : (double) ipio2[j];
+
+  /* compute q[0],q[1],...q[jk] */
+  for (i = 0; i <= jk; i++)
+    {
+      for (j = 0, fw = 0.0; j <= jx; j++)
+       fw += x[j] * f[jx + i - j];
+      q[i] = fw;
+    }
+
+  jz = jk;
+recompute:
+  /* distill q[] into iq[] reversingly */
+  for (i = 0, j = jz, z = q[jz]; j > 0; i++, j--)
+    {
+      fw = (double) ((int) (twon24 * z));
+      iq[i] = (int) (z - two24 * fw);
+      z = q[j - 1] + fw;
+    }
+
+  /* compute n */
+  z = ldexp (z, q0);           /* actual value of z */
+  z -= 8.0 * floor (z * 0.125);        /* trim off integer >= 8 */
+  n = (int) z;
+  z -= (double) n;
+  ih = 0;
+  if (q0 > 0)
+    {                          /* need iq[jz-1] to determine n */
+      i = (iq[jz - 1] >> (24 - q0));
+      n += i;
+      iq[jz - 1] -= i << (24 - q0);
+      ih = iq[jz - 1] >> (23 - q0);
+    }
+  else if (q0 == 0)
+    ih = iq[jz - 1] >> 23;
+  else if (z >= 0.5)
+    ih = 2;
+
+  if (ih > 0)
+    {                          /* q > 0.5 */
+      n += 1;
+      carry = 0;
+      for (i = 0; i < jz; i++)
+       {                       /* compute 1-q */
+         j = iq[i];
+         if (carry == 0)
+           {
+             if (j != 0)
+               {
+                 carry = 1;
+                 iq[i] = 0x1000000 - j;
+               }
+           }
+         else
+           iq[i] = 0xffffff - j;
+       }
+      if (q0 > 0)
+       {                       /* rare case: chance is 1 in 12 */
+         switch (q0)
+           {
+           case 1:
+             iq[jz - 1] &= 0x7fffff;
+             break;
+           case 2:
+             iq[jz - 1] &= 0x3fffff;
+             break;
+           }
+       }
+      if (ih == 2)
+       {
+         z = one - z;
+         if (carry != 0)
+           z -= ldexp (one, q0);
+       }
+    }
+
+  /* check if recomputation is needed */
+  if (z == zero)
+    {
+      j = 0;
+      for (i = jz - 1; i >= jk; i--)
+       j |= iq[i];
+      if (j == 0)
+       {                       /* need recomputation */
+         for (k = 1; iq[jk - k] == 0; k++);    /* k = no. of terms needed */
+
+         for (i = jz + 1; i <= jz + k; i++)
+           {                   /* add q[jz+1] to q[jz+k] */
+             f[jx + i] = (double) ipio2[jv + i];
+             for (j = 0, fw = 0.0; j <= jx; j++)
+               fw += x[j] * f[jx + i - j];
+             q[i] = fw;
+           }
+         jz += k;
+         goto recompute;
+       }
+    }
+
+  /* chop off zero terms */
+  if (z == 0.0)
+    {
+      jz -= 1;
+      q0 -= 24;
+      while (iq[jz] == 0)
+       {
+         jz--;
+         q0 -= 24;
+       }
+    }
+  else
+    {                          /* break z into 24-bit if necessary */
+      z = ldexp (z, -q0);
+      if (z >= two24)
+       {
+         fw = (double) ((int) (twon24 * z));
+         iq[jz] = (int) (z - two24 * fw);
+         jz += 1;
+         q0 += 24;
+         iq[jz] = (int) fw;
+       }
+      else
+       iq[jz] = (int) z;
+    }
+
+  /* convert integer "bit" chunk to floating-point value */
+  fw = ldexp (one, q0);
+  for (i = jz; i >= 0; i--)
+    {
+      q[i] = fw * (double) iq[i];
+      fw *= twon24;
+    }
+
+  /* compute PIo2[0,...,jp]*q[jz,...,0] */
+  for (i = jz; i >= 0; i--)
+    {
+      for (fw = 0.0, k = 0; k <= jp && k <= jz - i; k++)
+       fw += PIo2[k] * q[i + k];
+      fq[jz - i] = fw;
+    }
+
+  /* compress fq[] into y[] */
+  switch (prec)
+    {
+    case 0:
+      fw = 0.0;
+      for (i = jz; i >= 0; i--)
+       fw += fq[i];
+      y[0] = (ih == 0) ? fw : -fw;
+      break;
+    case 1:
+    case 2:
+      fw = 0.0;
+      for (i = jz; i >= 0; i--)
+       fw += fq[i];
+      y[0] = (ih == 0) ? fw : -fw;
+      fw = fq[0] - fw;
+      for (i = 1; i <= jz; i++)
+       fw += fq[i];
+      y[1] = (ih == 0) ? fw : -fw;
+      break;
+    case 3:                    /* painful */
+      for (i = jz; i > 0; i--)
+       {
+         fw = fq[i - 1] + fq[i];
+         fq[i] += fq[i - 1] - fw;
+         fq[i - 1] = fw;
+       }
+      for (i = jz; i > 1; i--)
+       {
+         fw = fq[i - 1] + fq[i];
+         fq[i] += fq[i - 1] - fw;
+         fq[i - 1] = fw;
+       }
+      for (fw = 0.0, i = jz; i >= 2; i--)
+       fw += fq[i];
+      if (ih == 0)
+       {
+         y[0] = fq[0];
+         y[1] = fq[1];
+         y[2] = fw;
+       }
+      else
+       {
+         y[0] = -fq[0];
+         y[1] = -fq[1];
+         y[2] = -fw;
+       }
+    }
+  return n & 7;
+}
diff --git a/lib/trigl.h b/lib/trigl.h
new file mode 100644 (file)
index 0000000..b4a7a15
--- /dev/null
@@ -0,0 +1,27 @@
+/* Declarations for sinl, cosl, tanl internal functions
+   Contributed by Paolo Bonzini
+
+   Copyright 2002, 2003 Free Software Foundation, Inc.
+
+   This file is part of gnulib.
+
+   gnulib is free software; you can redistribute it and/or modify it
+   under the terms of the GNU Lesser General Public License as published
+   by the Free Software Foundation; either version 2.1, or (at your option)
+   any later version.
+
+   gnulib is distributed in the hope that it will be useful, but WITHOUT
+   ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+   FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
+   License for more details.
+
+   You should have received a copy of the GNU Lesser General Public License
+   along with gnulib; see the file COPYING.LIB.  If not, write to the Free
+   Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307,
+   USA.
+ */
+
+extern int ieee754_rem_pio2l (long double x, long double *y);
+extern long double kernel_sinl (long double x, long double y, int iy);
+extern long double kernel_cosl (long double x, long double y);
+
index d714710..3c93617 100644 (file)
@@ -1,3 +1,7 @@
+2003-02-18  Paolo Bonzino  <bonzini@gnu.org>
+
+       * mathl.m4: New file.
+
 2003-02-17  Bruno Haible  <bruno@clisp.org>
 
        * mkdtemp.m4: New file, from GNU gettext with modifications.
diff --git a/m4/mathl.m4 b/m4/mathl.m4
new file mode 100644 (file)
index 0000000..a7d0002
--- /dev/null
@@ -0,0 +1,15 @@
+# mathl.m4 serial 1
+dnl Copyright (c) 2003 Free Software Foundation, Inc.
+dnl This file is free software, distributed under the terms of the GNU
+dnl General Public License.  As a special exception to the GNU General
+dnl Public License, this file may be distributed as part of a program
+dnl that contains a configuration script generated by Autoconf, under
+dnl the same distribution terms as the rest of the program.
+
+AC_DEFUN([gl_FUNC_LONG_DOUBLE_MATH], [
+
+AC_CHECK_LIB(m, atan)
+AC_REPLACE_FUNCS(floorl ceill sqrtl frexpl ldexpl asinl acosl atanl \
+        logl expl tanl sinl cosl)
+
+])
diff --git a/modules/mathl b/modules/mathl
new file mode 100644 (file)
index 0000000..363a240
--- /dev/null
@@ -0,0 +1,37 @@
+Description:
+C99 functions for transcendental functions with long double arguments.
+
+Files:
+lib/mathl.h
+lib/acosl.c
+lib/asinl.c
+lib/atanl.c
+lib/ceill.c
+lib/cosl.c
+lib/expl.c
+lib/floorl.c
+lib/frexpl.c
+lib/ldexpl.c
+lib/logl.c
+lib/sincosl.c
+lib/sinl.c
+lib/sqrtl.c
+lib/tanl.c
+lib/trigl.c
+lib/trigl.h
+m4/mathl.m4
+
+Depends-on:
+
+configure.ac:
+gl_FUNC_LONG_DOUBLE_MATH
+
+Makefile.am:
+noinst_HEADERS += mathl.h trigl.h trigl.c sincosl.c
+
+Include:
+"mathl.h"
+
+Maintainer:
+Paolo Bonzini  <bonzini@gnu.org>
+