Re: cosl test failure

[Paolo Bonzini]

On Sun, Mar 14, 2010 at 23:18, Bruno Haible <address@hidden> wrote:
> Hi Paolo,
>
> On 2010-01-18, in
> <http://lists.gnu.org/archive/html/bug-gnulib/2010-01/msg00256.html>,
> I observed that gnulib's cosl replacement function does not have the
> necessary accuracy. This was due to a wrong formula: The term
>  sincosl_table [index + SINCOSL_COS_HI] * cos_l_m1
> was being subtracted, when it should have been added.

I totally cannot recall, sorry.  Thanks for fixing it, it was on my
todo list but you beat me.

Paolo

> Also, in this formula and the other one for sinl(), you omitted one of
> the five summands; I cannot see a reason why.
>
> This fixes it.
>
>
> 2010-03-14  Bruno Haible  <address@hidden>
>
>        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.
>
> *** lib/trigl.h.orig    Sun Mar 14 23:11:58 2010
> --- lib/trigl.h Sun Mar 14 22:14:54 2010
> ***************
> *** 18,24 ****
>     You should have received a copy of the GNU General Public License
>     along with this program.  If not, see <http://www.gnu.org/licenses/>.  */
>
>  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);
>
> --- 18,35 ----
>     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);
> *** lib/sincosl.c.orig  Sun Mar 14 23:11:58 2010
> --- lib/sincosl.c       Sun Mar 14 23:04:53 2010
> ***************
> *** 136,146 ****
>    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;
> --- 136,147 ----
>    else
>      {
>        /* So that we don't have to use too large polynomial,  we find
> !          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(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,168 ****
>          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;
>      }
>  }
> --- 159,172 ----
>          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_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,205 ****
>    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;
> --- 199,210 ----
>    else
>      {
>        /* So that we don't have to use too large polynomial,  we find
> !          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
> !          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,222 ****
>          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;
>      }
>  }
> --- 218,231 ----
>          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_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;
>      }
>  }
>
>