X-Git-Url: http://erislabs.net/gitweb/?a=blobdiff_plain;f=lib%2Ftrigl.c;h=418ded96295e9f8b19ec29dba77c146727f69e33;hb=5191b3546cfb6c163228c23f214e325ddf60d46f;hp=ed9fd51d0aeeee4bfa821a8c27ab006bd109bb6f;hpb=e1123c2f4fedae90435426a82a05cf2e3233e97d;p=gnulib.git diff --git a/lib/trigl.c b/lib/trigl.c index ed9fd51d0..418ded962 100644 --- a/lib/trigl.c +++ b/lib/trigl.c @@ -1,28 +1,31 @@ /* Quad-precision floating point argument reduction. - Copyright (C) 1999, 2007 Free Software Foundation, Inc. + Copyright (C) 1999, 2007, 2009-2013 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Jakub Jelinek - This program is free software; you can redistribute it and/or modify + This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by - the Free Software Foundation; either version 2, or (at your option) - any later version. + the Free Software Foundation; either version 3 of the License, or + (at your option) any later version. This program 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 General Public License for more details. - You should have received a copy of the GNU General Public License along - with this program; if not, write to the Free Software Foundation, - Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. */ + You should have received a copy of the GNU General Public License + along with this program. If not, see . */ #include /* Specification. */ -#include +#include "trigl.h" #include +#include + +/* Code based on glibc/sysdeps/ieee754/ldbl-128/e_rem_pio2l.c + and glibc/sysdeps/ieee754/dbl-64/k_rem_pio2.c. */ /* Table of constants for 2/pi, 5628 hexadecimal digits of 2/pi */ static const int two_over_pi[] = { @@ -188,15 +191,15 @@ static const int two_over_pi[] = { static const long double c[] = { /* 93 bits of pi/2 */ #define PI_2_1 c[0] - 1.57079632679489661923132169155131424e+00L, /* 3fff921fb54442d18469898cc5100000 */ + 1.57079632679489661923132169155131424e+00L, /* 3fff921fb54442d18469898cc5100000 */ /* pi/2 - PI_2_1 */ #define PI_2_1t c[1] - 8.84372056613570112025531863263659260e-29L, /* 3fa1c06e0e68948127044533e63a0106 */ + 8.84372056613570112025531863263659260e-29L, /* 3fa1c06e0e68948127044533e63a0106 */ }; static int kernel_rem_pio2 (double *x, double *y, int e0, int nx, int prec, - const int *ipio2); + const int *ipio2); int ieee754_rem_pio2l (long double x, long double *y) @@ -206,7 +209,7 @@ ieee754_rem_pio2l (long double x, long double *y) int exp, n; if (x >= -0.78539816339744830961566084581987572104929234984377 - && x < 0.78539816339744830961566084581987572104929234984377) + && x <= 0.78539816339744830961566084581987572104929234984377) /* x in <-pi/4, pi/4> */ { y[0] = x; @@ -214,26 +217,25 @@ ieee754_rem_pio2l (long double x, long double *y) return 0; } - if (x >= 2.35619449019234492884698253745962716314787704953131 - && x < 2.35619449019234492884698253745962716314787704953131) - if (x > 0) + if (x > 0 && x < 2.35619449019234492884698253745962716314787704953131) { - /* 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; + /* 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 + + if (x < 0 && x > -2.35619449019234492884698253745962716314787704953131) { - /* 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; + /* 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 */ + if (x + x == x) /* x is ±oo */ { y[0] = x - x; y[1] = y[0]; @@ -297,7 +299,7 @@ static char rcsid[] = * 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 + * y = x - N*pi/2 * so that |y| < pi/2. * * The method is to compute the integer (mod 8) and fraction parts of @@ -309,93 +311,93 @@ static char rcsid[] = * (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. + * 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 + * 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: + * 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); + * 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] + * e0 The exponent of x[0] * - * nx dimension of x[] + * 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) + * 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[] + * 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)). + * ipio2[i] * 2^(-24(i+1)). * * External function: - * double scalbn(), floor(); + * 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. + * 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. + * jz local integer variable indicating the number of + * terms of ipio2[] used. * - * jx nx - 1 + * 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). + * 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. + * 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. + * 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. + * 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. + * PIo2[] double precision array, obtained by cutting pi/2 + * into 24 bits chunks. * - * f[] ipio2[] in floating point + * f[] ipio2[] in floating point * - * iq[] integer array by breaking up q[] in 24-bits chunk. + * iq[] integer array by breaking up q[] in 24-bits chunk. * - * fq[] final product of x*(2/pi) in fq[0],..,fq[jk] + * 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. + * ih integer. If >0 it indicates q[] is >= 0.5, hence + * it also indicates the *sign* of the result. * */ @@ -408,24 +410,24 @@ static char rcsid[] = * to produce the hexadecimal values shown. */ -static const int init_jk[] = { 2, 3, 4, 6 }; /* initial value for jk */ +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 */ + 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 */ +static const double zero = 0.0, one = 1.0, two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */ + twon24 = 5.96046447753906250000e-08; /* 0x3E700000, 0x00000000 */ -int +static int kernel_rem_pio2 (double *x, double *y, int e0, int nx, int prec, - const int *ipio2) + 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]; @@ -451,13 +453,13 @@ kernel_rem_pio2 (double *x, double *y, int e0, int nx, int prec, for (i = 0; i <= jk; i++) { for (j = 0, fw = 0.0; j <= jx; j++) - fw += x[j] * f[jx + i - j]; + fw += x[j] * f[jx + i - j]; q[i] = fw; } jz = jk; recompute: - /* distill q[] into iq[] reversingly */ + /* distill q[] into iq[] in reverse order */ for (i = 0, j = jz, z = q[jz]; j > 0; i++, j--) { fw = (double) ((int) (twon24 * z)); @@ -466,13 +468,13 @@ recompute: } /* compute n */ - z = ldexp (z, q0); /* actual value of z */ - z -= 8.0 * floor (z * 0.125); /* trim off integer >= 8 */ + 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 */ + { /* need iq[jz-1] to determine n */ i = (iq[jz - 1] >> (24 - q0)); n += i; iq[jz - 1] -= i << (24 - q0); @@ -484,41 +486,41 @@ recompute: ih = 2; if (ih > 0) - { /* q > 0.5 */ + { /* 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; - } + { /* 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; - } - } + { /* 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); - } + { + z = one - z; + if (carry != 0) + z -= ldexp (one, q0); + } } /* check if recomputation is needed */ @@ -526,21 +528,21 @@ recompute: { j = 0; for (i = jz - 1; i >= jk; i--) - j |= iq[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; - } + { /* 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 */ @@ -549,24 +551,24 @@ recompute: jz -= 1; q0 -= 24; while (iq[jz] == 0) - { - jz--; - q0 -= 24; - } + { + jz--; + q0 -= 24; + } } else - { /* break z into 24-bit if necessary */ + { /* 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; - } + { + fw = (double) ((int) (twon24 * z)); + iq[jz] = (int) (z - two24 * fw); + jz += 1; + q0 += 24; + iq[jz] = (int) fw; + } else - iq[jz] = (int) z; + iq[jz] = (int) z; } /* convert integer "bit" chunk to floating-point value */ @@ -581,7 +583,7 @@ recompute: for (i = jz; i >= 0; i--) { for (fw = 0.0, k = 0; k <= jp && k <= jz - i; k++) - fw += PIo2[k] * q[i + k]; + fw += PIo2[k] * q[i + k]; fq[jz - i] = fw; } @@ -591,47 +593,47 @@ recompute: case 0: fw = 0.0; for (i = jz; i >= 0; i--) - fw += fq[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]; + fw += fq[i]; y[0] = (ih == 0) ? fw : -fw; fw = fq[0] - fw; for (i = 1; i <= jz; i++) - fw += fq[i]; + fw += fq[i]; y[1] = (ih == 0) ? fw : -fw; break; - case 3: /* painful */ + 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; - } + { + 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; - } + { + 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]; + fw += fq[i]; if (ih == 0) - { - y[0] = fq[0]; - y[1] = fq[1]; - y[2] = fw; - } + { + y[0] = fq[0]; + y[1] = fq[1]; + y[2] = fw; + } else - { - y[0] = -fq[0]; - y[1] = -fq[1]; - y[2] = -fw; - } + { + y[0] = -fq[0]; + y[1] = -fq[1]; + y[2] = -fw; + } } return n & 7; }