--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,11 @@
+2010-03-14  Bruno Haible  <bruno@clisp.org>
+
+	Fix values returned by sinl, cosl.
+	* lib/trigl.h: Add specification comments.
+	* lib/sincosl.c (kernel_sinl, kernel_cosl): Fix comments and formula
+	that combines the values from the precomputed table with the values of
+	the Chebyshev polynomials.
+
 2010-03-14  Bruno Haible  <bruno@clisp.org>
 
 	Fix compilation error when modules 'posix_spawn[p]' are not used.

--- a/lib/sincosl.c
+++ b/lib/sincosl.c
@@ -136,11 +136,12 @@
   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
+         k and l such that x = k + l,  where fabsl(l) <= 1.0/256 with 83
+         possible values for k.  We look up cosl(k) and sinl(k) 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).  */
+         sinl(k+l) = sinl(k)cosl(l) + cosl(k)sinl(l).
+         Furthermore write k = 0.1484375 + h.  */
       x -= 0.1484375L;
       index = (int) (x * 128L + 0.5L);
       h = index / 128.0L;
@@ -158,11 +159,14 @@
         z * (SCOS1 + z * (SCOS2 + z * (SCOS3 + z * (SCOS4 + z * SCOS5))));
 
       index *= 4;
+      /* We rely on this expression not being "contracted" by the compiler
+         (cf. ISO C 99 section 6.5 paragraph 8).  */
       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));
+        sincosl_table[index + SINCOSL_SIN_HI]
+        + (sincosl_table[index + SINCOSL_COS_HI] * sin_l
+           + (sincosl_table[index + SINCOSL_SIN_HI] * cos_l_m1
+              + (sincosl_table[index + SINCOSL_SIN_LO] * (1 + cos_l_m1)
+                 + sincosl_table[index + SINCOSL_COS_LO] * sin_l)));
       return z * sign;
     }
 }
@@ -195,11 +199,12 @@
   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
+         k and l such that x = k + l,  where fabsl(l) <= 1.0/256 with 83
+         possible values for k.  We look up cosl(k) and sinl(k) 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).  */
+         cosl(k+l) = cosl(k)cosl(l) - sinl(k)sinl(l).
+         Furthermore write k = 0.1484375 + h.  */
       x -= 0.1484375L;
       index = (int) (x * 128L + 0.5L);
       h = index / 128.0L;
@@ -213,10 +218,14 @@
         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));
+      /* We rely on this expression not being "contracted" by the compiler
+         (cf. ISO C 99 section 6.5 paragraph 8).  */
+      z =
+        sincosl_table [index + SINCOSL_COS_HI]
+        - (sincosl_table [index + SINCOSL_SIN_HI] * sin_l
+           - (sincosl_table [index + SINCOSL_COS_HI] * cos_l_m1
+              + (sincosl_table [index + SINCOSL_COS_LO] * (1 + cos_l_m1)
+                 - sincosl_table [index + SINCOSL_SIN_LO] * sin_l)));
       return z;
     }
 }

--- a/lib/trigl.h
+++ b/lib/trigl.h
@@ -18,7 +18,18 @@
    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.  */
 
+/* Decompose x into x = k * π/2 + r
+   where k is an integer and abs(r) <= π/4.
+   Store r in y[0] and y[1] (main part in y[0], small additional part in
+   y[1], r = y[0] + y[1]).
+   Return k.  */
 extern int ieee754_rem_pio2l (long double x, long double *y);
+
+/* Compute and return sinl (x + y), where x is the main part and y is the
+   small additional part of a floating-point number.
+   iy is 0 when y is known to be 0.0, otherwise iy is 1.  */
 extern long double kernel_sinl (long double x, long double y, int iy);
+
+/* Compute and return cosl (x + y), where x is the main part and y is the
+   small additional part of a floating-point number.  */
 extern long double kernel_cosl (long double x, long double y);
-