AuroraMiddleware/glibc
1
0
Fork 0
mirror of https://sourceware.org/git/glibc.git synced 2024-11-08 14:20:07 +00:00
Code Issues Packages Projects Releases Wiki Activity

Remove unused math files

Browse Source 
Remove empty files due to the sin/cos improvements: k_sinf.c, k_cosf.c,
k_cos.c, k_sin.c.  After the tanf change s_rem_pio2f.c and k_rem_pio2f.c
(and the ia64, m68k and powerpc equivalents) are no longer used,
so remove them.  All e_rem_pio2.c files were already empty or commented
out, so remove them too.  Passes build-many-glibcs.

	* math/Makefile: Remove empty files k_sin(f).c, k_cos(f).c.
	Remove unused files e_rem_pio2(f).c, k_rem_pio2f.c.
	* sysdeps/i386/fpu/e_rem_pio2.c: Delete file.
	* sysdeps/ia64/fpu/e_rem_pio2.c: Likewise.
	* sysdeps/ia64/fpu/e_rem_pio2f.c: Likewise.
	* sysdeps/ia64/fpu/k_rem_pio2f.c: Likewise.
	* sysdeps/ieee754/dbl-64/e_rem_pio2.c: Likewise.
	* sysdeps/ieee754/dbl-64/k_cos.c: Likewise.
	* sysdeps/ieee754/dbl-64/k_sin.c: Likewise.
	* sysdeps/ieee754/flt-32/e_rem_pio2f.c: Likewise.
	* sysdeps/ieee754/flt-32/k_cosf.c: Likewise.
	* sysdeps/ieee754/flt-32/k_rem_pio2f.c: Likewise.
	* sysdeps/ieee754/flt-32/k_sinf.c: Likewise.
	* sysdeps/m68k/m680x0/fpu/e_rem_pio2.c: Likewise
	* sysdeps/m68k/m680x0/fpu/e_rem_pio2f.c: Likewise
	* sysdeps/m68k/m680x0/fpu/k_rem_pio2f.c: Likewise
	* sysdeps/powerpc/fpu/e_rem_pio2f.c: Likewise.
	* sysdeps/powerpc/fpu/k_rem_pio2f.c: Likewise.
This commit is contained in:
Wilco Dijkstra 2018-08-24 15:31:25 +01:00
parent 60bcac09c0
commit ca3aac57ef
18 changed files with 30 additions and 1078 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

21
ChangeLog
View File

@ -1,3 +1,24 @@
2018-08-24 Wilco Dijkstra <wdijkstr@arm.com>
* math/Makefile: Remove empty files k_sin(f).c, k_cos(f).c.
Remove unused files e_rem_pio2(f).c, k_rem_pio2f.c.
* sysdeps/i386/fpu/e_rem_pio2.c: Delete file.
* sysdeps/ia64/fpu/e_rem_pio2.c: Likewise.
* sysdeps/ia64/fpu/e_rem_pio2f.c: Likewise.
* sysdeps/ia64/fpu/k_rem_pio2f.c: Likewise.
* sysdeps/ieee754/dbl-64/e_rem_pio2.c: Likewise.
* sysdeps/ieee754/dbl-64/k_cos.c: Likewise.
* sysdeps/ieee754/dbl-64/k_sin.c: Likewise.
* sysdeps/ieee754/flt-32/e_rem_pio2f.c: Likewise.
* sysdeps/ieee754/flt-32/k_cosf.c: Likewise.
* sysdeps/ieee754/flt-32/k_rem_pio2f.c: Likewise.
* sysdeps/ieee754/flt-32/k_sinf.c: Likewise.
* sysdeps/m68k/m680x0/fpu/e_rem_pio2.c: Likewise
* sysdeps/m68k/m680x0/fpu/e_rem_pio2f.c: Likewise
* sysdeps/m68k/m680x0/fpu/k_rem_pio2f.c: Likewise
* sysdeps/powerpc/fpu/e_rem_pio2f.c: Likewise.
* sysdeps/powerpc/fpu/k_rem_pio2f.c: Likewise.
2018-08-23 Joseph Myers <joseph@codesourcery.com>
* sysdeps/generic/math-tests-exceptions.h: New file.

16
math/Makefile
View File

@ -64,9 +64,9 @@ gen-libm-calls = cargF conjF cimagF crealF cabsF s_cacosF \
libm-calls = \
e_acosF e_acoshF e_asinF e_atan2F e_atanhF e_coshF e_expF e_fmodF \
e_hypotF e_j0F e_j1F e_jnF e_lgammaF_r e_logF e_log10F e_powF \
e_rem_pio2F e_remainderF e_scalbF e_sinhF e_sqrtF e_gammaF_r \
e_remainderF e_scalbF e_sinhF e_sqrtF e_gammaF_r \
e_ilogbF \
k_cosF k_sinF k_tanF s_asinhF s_atanF s_cbrtF \
k_tanF s_asinhF s_atanF s_cbrtF \
s_ceilF s_cosF s_erfF s_expm1F s_fabsF \
s_floorF s_log1pF s_logbF \
s_nextafterF s_nexttowardF s_rintF s_scalblnF \
@ -119,23 +119,25 @@ test-types-basic = ldouble double float
# long double support
type-ldouble-suffix := l
type-ldouble-routines := t_sincosl k_sincosl s_iscanonicall
type-ldouble-routines := t_sincosl k_sinl k_cosl k_sincosl s_iscanonicall \
e_rem_pio2l
type-ldouble-yes := ldouble
# double support
type-double-suffix :=
type-double-routines := branred doasin dosincos mpa mpatan2 \
mpatan mpsqrt mptan sincos32 \
sincostab k_rem_pio2
k_rem_pio2 mpatan mpsqrt mptan sincos32 \
sincostab
# float support
type-float-suffix := f
type-float-routines := k_rem_pio2f math_errf e_exp2f_data e_logf_data \
type-float-routines := math_errf e_exp2f_data e_logf_data \
e_log2f_data e_powf_log2_data s_sincosf_data
# _Float128 support
type-float128-suffix := f128
type-float128-routines := t_sincosf128 k_sincosf128
type-float128-routines := t_sincosf128 k_sinf128 k_cosf128 k_sincosf128 \
e_rem_pio2f128
type-float128-yes := float128
# _Float64x may be supported, only as an alias type.

3
sysdeps/i386/fpu/e_rem_pio2.c
View File

@ -1,3 +0,0 @@
/* Empty. This file is only meant to avoid compiling the file with the
same name in the libm-ieee754 directory. The code is not used since
there is an assembler version for all users of this file. */

1
sysdeps/ia64/fpu/e_rem_pio2.c
View File

@ -1 +0,0 @@
/* Not needed. */

1
sysdeps/ia64/fpu/e_rem_pio2f.c
View File

@ -1 +0,0 @@
/* Not needed. */

1
sysdeps/ia64/fpu/k_rem_pio2f.c
View File

@ -1 +0,0 @@
/* Not needed. */

193
sysdeps/ieee754/dbl-64/e_rem_pio2.c
View File

@ -1,193 +0,0 @@
#ifdef NOT_NEEDED_ANYMORE
/*
* ====================================================
* 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.
* ====================================================
*/
/* __ieee754_rem_pio2(x,y)
*
* return the remainder of x rem pi/2 in y[0]+y[1]
* use __kernel_rem_pio2()
*/
#include <math.h>
#include <math_private.h>
/*
* Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi
*/
static const int32_t two_over_pi[] = {
0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62,
0x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A,
0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129,
0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41,
0x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8,
0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF,
0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5,
0xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08,
0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3,
0x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880,
0x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B,
};
static const int32_t npio2_hw[] = {
0x3FF921FB, 0x400921FB, 0x4012D97C, 0x401921FB, 0x401F6A7A, 0x4022D97C,
0x4025FDBB, 0x402921FB, 0x402C463A, 0x402F6A7A, 0x4031475C, 0x4032D97C,
0x40346B9C, 0x4035FDBB, 0x40378FDB, 0x403921FB, 0x403AB41B, 0x403C463A,
0x403DD85A, 0x403F6A7A, 0x40407E4C, 0x4041475C, 0x4042106C, 0x4042D97C,
0x4043A28C, 0x40446B9C, 0x404534AC, 0x4045FDBB, 0x4046C6CB, 0x40478FDB,
0x404858EB, 0x404921FB,
};
/*
* invpio2: 53 bits of 2/pi
* pio2_1: first 33 bit of pi/2
* pio2_1t: pi/2 - pio2_1
* pio2_2: second 33 bit of pi/2
* pio2_2t: pi/2 - (pio2_1+pio2_2)
* pio2_3: third 33 bit of pi/2
* pio2_3t: pi/2 - (pio2_1+pio2_2+pio2_3)
*/
static const double
zero = 0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
half = 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
invpio2 = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
pio2_1 = 1.57079632673412561417e+00, /* 0x3FF921FB, 0x54400000 */
pio2_1t = 6.07710050650619224932e-11, /* 0x3DD0B461, 0x1A626331 */
pio2_2 = 6.07710050630396597660e-11, /* 0x3DD0B461, 0x1A600000 */
pio2_2t = 2.02226624879595063154e-21, /* 0x3BA3198A, 0x2E037073 */
pio2_3 = 2.02226624871116645580e-21, /* 0x3BA3198A, 0x2E000000 */
pio2_3t = 8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */
int32_t
__ieee754_rem_pio2 (double x, double *y)
{
double z, w, t, r, fn;
double tx[3];
int32_t e0, i, j, nx, n, ix, hx;
uint32_t low;
GET_HIGH_WORD (hx, x); /* high word of x */
ix = hx & 0x7fffffff;
if (ix <= 0x3fe921fb) /* |x| ~<= pi/4 , no need for reduction */
{
y[0] = x; y[1] = 0; return 0;
}
if (ix < 0x4002d97c) /* |x| < 3pi/4, special case with n=+-1 */
{
if (hx > 0)
{
z = x - pio2_1;
if (ix != 0x3ff921fb) /* 33+53 bit pi is good enough */
{
y[0] = z - pio2_1t;
y[1] = (z - y[0]) - pio2_1t;
}
else /* near pi/2, use 33+33+53 bit pi */
{
z -= pio2_2;
y[0] = z - pio2_2t;
y[1] = (z - y[0]) - pio2_2t;
}
return 1;
}
else /* negative x */
{
z = x + pio2_1;
if (ix != 0x3ff921fb) /* 33+53 bit pi is good enough */
{
y[0] = z + pio2_1t;
y[1] = (z - y[0]) + pio2_1t;
}
else /* near pi/2, use 33+33+53 bit pi */
{
z += pio2_2;
y[0] = z + pio2_2t;
y[1] = (z - y[0]) + pio2_2t;
}
return -1;
}
}
if (ix <= 0x413921fb) /* |x| ~<= 2^19*(pi/2), medium size */
{
t = fabs (x);
n = (int32_t) (t * invpio2 + half);
fn = (double) n;
r = t - fn * pio2_1;
w = fn * pio2_1t; /* 1st round good to 85 bit */
if (n < 32 && ix != npio2_hw[n - 1])
{
y[0] = r - w; /* quick check no cancellation */
}
else
{
uint32_t high;
j = ix >> 20;
y[0] = r - w;
GET_HIGH_WORD (high, y[0]);
i = j - ((high >> 20) & 0x7ff);
if (i > 16) /* 2nd iteration needed, good to 118 */
{
t = r;
w = fn * pio2_2;
r = t - w;
w = fn * pio2_2t - ((t - r) - w);
y[0] = r - w;
GET_HIGH_WORD (high, y[0]);
i = j - ((high >> 20) & 0x7ff);
if (i > 49) /* 3rd iteration need, 151 bits acc */
{
t = r; /* will cover all possible cases */
w = fn * pio2_3;
r = t - w;
w = fn * pio2_3t - ((t - r) - w);
y[0] = r - w;
}
}
}
y[1] = (r - y[0]) - w;
if (hx < 0)
{
y[0] = -y[0]; y[1] = -y[1]; return -n;
}
else
return n;
}
/*
* all other (large) arguments
*/
if (ix >= 0x7ff00000) /* x is inf or NaN */
{
y[0] = y[1] = x - x; return 0;
}
/* set z = scalbn(|x|,ilogb(x)-23) */
GET_LOW_WORD (low, x);
SET_LOW_WORD (z, low);
e0 = (ix >> 20) - 1046; /* e0 = ilogb(z)-23; */
SET_HIGH_WORD (z, ix - ((int32_t) (e0 << 20)));
for (i = 0; i < 2; i++)
{
tx[i] = (double) ((int32_t) (z));
z = (z - tx[i]) * two24;
}
tx[2] = z;
nx = 3;
while (tx[nx - 1] == zero)
nx--; /* skip zero term */
n = __kernel_rem_pio2 (tx, y, e0, nx, 2, two_over_pi);
if (hx < 0)
{
y[0] = -y[0]; y[1] = -y[1]; return -n;
}
return n;
}
#endif

1
sysdeps/ieee754/dbl-64/k_cos.c
View File

@ -1 +0,0 @@
/* Not needed anymore. */

1
sysdeps/ieee754/dbl-64/k_sin.c
View File

@ -1 +0,0 @@
/* Not needed anymore. */

179
sysdeps/ieee754/flt-32/e_rem_pio2f.c
View File

@ -1,179 +0,0 @@
/* e_rem_pio2f.c -- float version of e_rem_pio2.c
* Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
*/
/*
* ====================================================
* 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: e_rem_pio2f.c,v 1.5 1995/05/10 20:46:03 jtc Exp $";
#endif
/* __ieee754_rem_pio2f(x,y)
*
* return the remainder of x rem pi/2 in y[0]+y[1]
* use __kernel_rem_pio2f()
*/
#include <math.h>
#include <math_private.h>
/*
* Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi
*/
static const int32_t two_over_pi[] = {
0xA2, 0xF9, 0x83, 0x6E, 0x4E, 0x44, 0x15, 0x29, 0xFC,
0x27, 0x57, 0xD1, 0xF5, 0x34, 0xDD, 0xC0, 0xDB, 0x62,
0x95, 0x99, 0x3C, 0x43, 0x90, 0x41, 0xFE, 0x51, 0x63,
0xAB, 0xDE, 0xBB, 0xC5, 0x61, 0xB7, 0x24, 0x6E, 0x3A,
0x42, 0x4D, 0xD2, 0xE0, 0x06, 0x49, 0x2E, 0xEA, 0x09,
0xD1, 0x92, 0x1C, 0xFE, 0x1D, 0xEB, 0x1C, 0xB1, 0x29,
0xA7, 0x3E, 0xE8, 0x82, 0x35, 0xF5, 0x2E, 0xBB, 0x44,
0x84, 0xE9, 0x9C, 0x70, 0x26, 0xB4, 0x5F, 0x7E, 0x41,
0x39, 0x91, 0xD6, 0x39, 0x83, 0x53, 0x39, 0xF4, 0x9C,
0x84, 0x5F, 0x8B, 0xBD, 0xF9, 0x28, 0x3B, 0x1F, 0xF8,
0x97, 0xFF, 0xDE, 0x05, 0x98, 0x0F, 0xEF, 0x2F, 0x11,
0x8B, 0x5A, 0x0A, 0x6D, 0x1F, 0x6D, 0x36, 0x7E, 0xCF,
0x27, 0xCB, 0x09, 0xB7, 0x4F, 0x46, 0x3F, 0x66, 0x9E,
0x5F, 0xEA, 0x2D, 0x75, 0x27, 0xBA, 0xC7, 0xEB, 0xE5,
0xF1, 0x7B, 0x3D, 0x07, 0x39, 0xF7, 0x8A, 0x52, 0x92,
0xEA, 0x6B, 0xFB, 0x5F, 0xB1, 0x1F, 0x8D, 0x5D, 0x08,
0x56, 0x03, 0x30, 0x46, 0xFC, 0x7B, 0x6B, 0xAB, 0xF0,
0xCF, 0xBC, 0x20, 0x9A, 0xF4, 0x36, 0x1D, 0xA9, 0xE3,
0x91, 0x61, 0x5E, 0xE6, 0x1B, 0x08, 0x65, 0x99, 0x85,
0x5F, 0x14, 0xA0, 0x68, 0x40, 0x8D, 0xFF, 0xD8, 0x80,
0x4D, 0x73, 0x27, 0x31, 0x06, 0x06, 0x15, 0x56, 0xCA,
0x73, 0xA8, 0xC9, 0x60, 0xE2, 0x7B, 0xC0, 0x8C, 0x6B,
};
/* This array is like the one in e_rem_pio2.c, but the numbers are
single precision and the last 8 bits are forced to 0. */
static const int32_t npio2_hw[] = {
0x3fc90f00, 0x40490f00, 0x4096cb00, 0x40c90f00, 0x40fb5300, 0x4116cb00,
0x412fed00, 0x41490f00, 0x41623100, 0x417b5300, 0x418a3a00, 0x4196cb00,
0x41a35c00, 0x41afed00, 0x41bc7e00, 0x41c90f00, 0x41d5a000, 0x41e23100,
0x41eec200, 0x41fb5300, 0x4203f200, 0x420a3a00, 0x42108300, 0x4216cb00,
0x421d1400, 0x42235c00, 0x4229a500, 0x422fed00, 0x42363600, 0x423c7e00,
0x4242c700, 0x42490f00
};
/*
* invpio2: 24 bits of 2/pi
* pio2_1: first 17 bit of pi/2
* pio2_1t: pi/2 - pio2_1
* pio2_2: second 17 bit of pi/2
* pio2_2t: pi/2 - (pio2_1+pio2_2)
* pio2_3: third 17 bit of pi/2
* pio2_3t: pi/2 - (pio2_1+pio2_2+pio2_3)
*/
static const float
zero = 0.0000000000e+00, /* 0x00000000 */
half = 5.0000000000e-01, /* 0x3f000000 */
two8 = 2.5600000000e+02, /* 0x43800000 */
invpio2 = 6.3661980629e-01, /* 0x3f22f984 */
pio2_1 = 1.5707855225e+00, /* 0x3fc90f80 */
pio2_1t = 1.0804334124e-05, /* 0x37354443 */
pio2_2 = 1.0804273188e-05, /* 0x37354400 */
pio2_2t = 6.0770999344e-11, /* 0x2e85a308 */
pio2_3 = 6.0770943833e-11, /* 0x2e85a300 */
pio2_3t = 6.1232342629e-17; /* 0x248d3132 */
int32_t __ieee754_rem_pio2f(float x, float *y)
{
float z,w,t,r,fn;
float tx[3];
int32_t e0,i,j,nx,n,ix,hx;
GET_FLOAT_WORD(hx,x);
ix = hx&0x7fffffff;
if(ix<=0x3f490fd8) /* |x| ~<= pi/4 , no need for reduction */
{y[0] = x; y[1] = 0; return 0;}
if(ix<0x4016cbe4) { /* |x| < 3pi/4, special case with n=+-1 */
if(hx>0) {
z = x - pio2_1;
if((ix&0xffffffc0)!=0x3fc90fc0) { /* 24+24 bit pi OK */
y[0] = z - pio2_1t;
y[1] = (z-y[0])-pio2_1t;
} else { /* near pi/2, use 24+24+24 bit pi */
z -= pio2_2;
y[0] = z - pio2_2t;
y[1] = (z-y[0])-pio2_2t;
}
return 1;
} else { /* negative x */
z = x + pio2_1;
if((ix&0xffffffc0)!=0x3fc90fc0) { /* 24+24 bit pi OK */
y[0] = z + pio2_1t;
y[1] = (z-y[0])+pio2_1t;
} else { /* near pi/2, use 24+24+24 bit pi */
z += pio2_2;
y[0] = z + pio2_2t;
y[1] = (z-y[0])+pio2_2t;
}
return -1;
}
}
if(ix<=0x43490f80) { /* |x| ~<= 2^7*(pi/2), medium size */
t = fabsf(x);
n = (int32_t) (t*invpio2+half);
fn = (float)n;
r = t-fn*pio2_1;
w = fn*pio2_1t; /* 1st round good to 40 bit */
if(n<32&&(int32_t)(ix&0xffffff00)!=npio2_hw[n-1]) {
y[0] = r-w; /* quick check no cancellation */
} else {
uint32_t high;
j = ix>>23;
y[0] = r-w;
GET_FLOAT_WORD(high,y[0]);
i = j-((high>>23)&0xff);
if(i>8) { /* 2nd iteration needed, good to 57 */
t = r;
w = fn*pio2_2;
r = t-w;
w = fn*pio2_2t-((t-r)-w);
y[0] = r-w;
GET_FLOAT_WORD(high,y[0]);
i = j-((high>>23)&0xff);
if(i>25) { /* 3rd iteration need, 74 bits acc */
t = r; /* will cover all possible cases */
w = fn*pio2_3;
r = t-w;
w = fn*pio2_3t-((t-r)-w);
y[0] = r-w;
}
}
}
y[1] = (r-y[0])-w;
if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;}
else return n;
}
/*
* all other (large) arguments
*/
if(ix>=0x7f800000) { /* x is inf or NaN */
y[0]=y[1]=x-x; return 0;
}
/* set z = scalbn(|x|,ilogb(x)-7) */
e0 = (ix>>23)-134; /* e0 = ilogb(z)-7; */
SET_FLOAT_WORD(z, ix - ((int32_t)(e0<<23)));
for(i=0;i<2;i++) {
tx[i] = (float)((int32_t)(z));
z = (z-tx[i])*two8;
}
tx[2] = z;
nx = 3;
while(tx[nx-1]==zero) nx--; /* skip zero term */
n = __kernel_rem_pio2f(tx,y,e0,nx,2,two_over_pi);
if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;}
return n;
}

1
sysdeps/ieee754/flt-32/k_cosf.c
View File

@ -1 +0,0 @@
/* Not needed. */

219
sysdeps/ieee754/flt-32/k_rem_pio2f.c
View File

@ -1,219 +0,0 @@
/* k_rem_pio2f.c -- float version of k_rem_pio2.c
* Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
*/
/*
* ====================================================
* 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_pio2f.c,v 1.4 1995/05/10 20:46:28 jtc Exp $";
#endif
#include <math.h>
#include <math-narrow-eval.h>
#include <math_private.h>
#include <libc-diag.h>
/* In the float version, the input parameter x contains 8 bit
integers, not 24 bit integers. 113 bit precision is not supported. */
static const int init_jk[] = {4,7,9}; /* initial value for jk */
static const float PIo2[] = {
1.5703125000e+00, /* 0x3fc90000 */
4.5776367188e-04, /* 0x39f00000 */
2.5987625122e-05, /* 0x37da0000 */
7.5437128544e-08, /* 0x33a20000 */
6.0026650317e-11, /* 0x2e840000 */
7.3896444519e-13, /* 0x2b500000 */
5.3845816694e-15, /* 0x27c20000 */
5.6378512969e-18, /* 0x22d00000 */
8.3009228831e-20, /* 0x1fc40000 */
3.2756352257e-22, /* 0x1bc60000 */
6.3331015649e-25, /* 0x17440000 */
};
static const float
zero = 0.0,
one = 1.0,
two8 = 2.5600000000e+02, /* 0x43800000 */
twon8 = 3.9062500000e-03; /* 0x3b800000 */
int __kernel_rem_pio2f(float *x, float *y, int e0, int nx, int prec, const int32_t *ipio2)
{
int32_t jz,jx,jv,jp,jk,carry,n,iq[20],i,j,k,m,q0,ih;
float 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)/8; if(jv<0) jv=0;
q0 = e0-8*(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 : (float) 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 = (float)((int32_t)(twon8* z));
iq[i] = (int32_t)(z-two8*fw);
z = q[j-1]+fw;
}
/* compute n */
z = __scalbnf(z,q0); /* actual value of z */
z -= (float)8.0*__floorf(z*(float)0.125); /* trim off integer >= 8 */
n = (int32_t) z;
z -= (float)n;
ih = 0;
if(q0>0) { /* need iq[jz-1] to determine n */
i = (iq[jz-1]>>(8-q0)); n += i;
iq[jz-1] -= i<<(8-q0);
ih = iq[jz-1]>>(7-q0);
}
else if(q0==0) ih = iq[jz-1]>>7;
else if(z>=(float)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] = 0x100- j;
}
} else iq[i] = 0xff - j;
}
if(q0>0) { /* rare case: chance is 1 in 12 */
switch(q0) {
case 1:
iq[jz-1] &= 0x7f; break;
case 2:
iq[jz-1] &= 0x3f; break;
}
}
if(ih==2) {
z = one - z;
if(carry!=0) z -= __scalbnf(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 */
/* On s390x gcc 6.1 -O3 produces the warning "array subscript is
below array bounds [-Werror=array-bounds]". Only
__ieee754_rem_pio2f calls __kernel_rem_pio2f for normal
numbers and |x| ~> 2^7*(pi/2). Thus x can't be zero and
ipio2 is not zero, too. Thus not all iq[] values can't be
zero. */
DIAG_PUSH_NEEDS_COMMENT;
DIAG_IGNORE_NEEDS_COMMENT (6.1, "-Warray-bounds");
for(k=1;iq[jk-k]==0;k++); /* k = no. of terms needed */
DIAG_POP_NEEDS_COMMENT;
for(i=jz+1;i<=jz+k;i++) { /* add q[jz+1] to q[jz+k] */
f[jx+i] = (float) 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==(float)0.0) {
jz -= 1; q0 -= 8;
while(iq[jz]==0) { jz--; q0-=8;}
} else { /* break z into 8-bit if necessary */
z = __scalbnf(z,-q0);
if(z>=two8) {
fw = (float)((int32_t)(twon8*z));
iq[jz] = (int32_t)(z-two8*fw);
jz += 1; q0 += 8;
iq[jz] = (int32_t) fw;
} else iq[jz] = (int32_t) z ;
}
/* convert integer "bit" chunk to floating-point value */
fw = __scalbnf(one,q0);
for(i=jz;i>=0;i--) {
q[i] = fw*(float)iq[i]; fw*=twon8;
}
/* 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:;
float fv = 0.0;
for (i=jz;i>=0;i--) fv = math_narrow_eval (fv + fq[i]);
y[0] = (ih==0)? fv: -fv;
/* GCC mainline (to be GCC 9), as of 2018-05-22 on
i686, warns that fq[0] may be used uninitialized.
This is not possible because jz is always
nonnegative when the above loop initializing fq is
executed, because the result is never zero to full
precision (this function is not called for zero
arguments). */
DIAG_PUSH_NEEDS_COMMENT;
DIAG_IGNORE_NEEDS_COMMENT (9, "-Wmaybe-uninitialized");
fv = math_narrow_eval (fq[0]-fv);
DIAG_POP_NEEDS_COMMENT;
for (i=1;i<=jz;i++) fv = math_narrow_eval (fv + fq[i]);
y[1] = (ih==0)? fv: -fv;
break;
case 3: /* painful */
for (i=jz;i>0;i--) {
float fv = math_narrow_eval (fq[i-1]+fq[i]);
fq[i] += fq[i-1]-fv;
fq[i-1] = fv;
}
for (i=jz;i>1;i--) {
float fv = math_narrow_eval (fq[i-1]+fq[i]);
fq[i] += fq[i-1]-fv;
fq[i-1] = fv;
}
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;
}

1
sysdeps/ieee754/flt-32/k_sinf.c
View File

@ -1 +0,0 @@
/* Not needed. */

3
sysdeps/m68k/m680x0/fpu/e_rem_pio2.c
View File

@ -1,3 +0,0 @@
/* Empty. This file is only meant to avoid compiling the file with the
same name in the libm-ieee754 directory. The code is not used since
there is an assembler version for all users of this file. */

3
sysdeps/m68k/m680x0/fpu/e_rem_pio2f.c
View File

@ -1,3 +0,0 @@
/* Empty. This file is only meant to avoid compiling the file with the
same name in the libm-ieee754 directory. The code is not used since
there is an assembler version for all users of this file. */

3
sysdeps/m68k/m680x0/fpu/k_rem_pio2f.c
View File

@ -1,3 +0,0 @@
/* Empty. This file is only meant to avoid compiling the file with the
same name in the libm-ieee754 directory. The code is not used since
there is an assembler version for all users of this file. */

188
sysdeps/powerpc/fpu/e_rem_pio2f.c
View File

@ -1,188 +0,0 @@
/* e_rem_pio2f.c -- float version of e_rem_pio2.c
Copyright (C) 2011-2018 Free Software Foundation, Inc.
This file is part of the GNU C Library.
Contributed by Adhemerval Zanella <azanella@br.ibm.com>, 2011
The GNU C Library is free software; you can redistribute it and/or
modify it under the terms of the GNU Library General Public License as
published by the Free Software Foundation; either version 2 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
Library General Public License for more details.
You should have received a copy of the GNU Library General Public
License along with the GNU C Library; see the file COPYING.LIB. If
not, see <http://www.gnu.org/licenses/>. */
#include <math.h>
#include <math_private.h>
#include "s_float_bitwise.h"
/* defined in sysdeps/powerpc/fpu/k_rem_pio2f.c */
int __fp_kernel_rem_pio2f (float *x, float *y, float e0, int32_t nx);
/* __ieee754_rem_pio2f(x,y)
*
* return the remainder of x rem pi/2 in y[0]+y[1]
*/
static const float npio2_hw[] = {
1.57077026e+00, 3.14154053e+00, 4.71228027e+00, 6.28308105e+00,
7.85388184e+00, 9.42456055e+00, 1.09953613e+01, 1.25661621e+01,
1.41369629e+01, 1.57077637e+01, 1.72783203e+01, 1.88491211e+01,
2.04199219e+01, 2.19907227e+01, 2.35615234e+01, 2.51323242e+01,
2.67031250e+01, 2.82739258e+01, 2.98447266e+01, 3.14155273e+01,
3.29863281e+01, 3.45566406e+01, 3.61279297e+01, 3.76982422e+01,
3.92695312e+01, 4.08398438e+01, 4.24111328e+01, 4.39814453e+01,
4.55527344e+01, 4.71230469e+01, 4.86943359e+01, 5.02646484e+01
};
static const float zero = 0.0000000000e+00;
static const float two8 = 2.5600000000e+02;
static const float half = 5.0000000000e-01;
static const float invpio2 = 6.3661980629e-01;
static const float pio2_1 = 1.5707855225e+00;
static const float pio2_1t = 1.0804334124e-05;
static const float pio2_2 = 1.0804273188e-05;
static const float pio2_2t = 6.0770999344e-11;
static const float pio2_3 = 6.0770943833e-11;
static const float pio2_3t = 6.1232342629e-17;
static const float pio4 = 7.8539801e-01;
static const float pio3_4 = 2.3561945e+00;
static const float pio2_24b = 1.5707951e+00;
static const float pio2_2e7 = 2.0106054e+02;
int32_t
__ieee754_rem_pio2f (float x, float *y)
{
float ax, z, n, r, w, t, e0;
float tx[3];
int32_t i, nx;
ax = __builtin_fabsf (x);
if (ax <= pio4)
{
y[0] = x;
y[1] = 0;
return 0;
}
if (ax < pio3_4)
{
if (x > 0)
{
z = x - pio2_1;
if (!__float_and_test28 (ax, pio2_24b))
{
y[0] = z - pio2_1t;
y[1] = (z - y[0]) - pio2_1t;
}
else
{
z -= pio2_2;
y[0] = z - pio2_2t;
y[1] = (z - y[0]) - pio2_2t;
}
return 1;
}
else
{
z = x + pio2_1;
if (!__float_and_test28 (ax, pio2_24b))
{
y[0] = z + pio2_1t;
y[1] = (z - y[0]) + pio2_1t;
}
else
{
z += pio2_2;
y[0] = z + pio2_2t;
y[1] = (z - y[0]) + pio2_2t;
}
return -1;
}
}
if (ax <= pio2_2e7)
{
n = __floorf (ax * invpio2 + half);
i = (int32_t) n;
r = ax - n * pio2_1;
w = n * pio2_1t; /* 1st round good to 40 bit */
if (i < 32 && !__float_and_test24 (ax, npio2_hw[i - 1]))
{
y[0] = r - w;
}
else
{
float i, j;
j = __float_and8 (ax);
y[0] = r - w;
i = __float_and8 (y[0]);
if (j / i > 256.0 || j / i < 3.9062500e-3)
{ /* 2nd iterations needed, good to 57 */
t = r;
w = n * pio2_2;
r = t - w;
w = n * pio2_2t - ((t - r) - w);
y[0] = r - w;
i = __float_and8 (y[0]);
if (j / i > 33554432 || j / i < 2.9802322e-8)
{ /* 3rd iteration needed, 74 bits acc */
t = r;
w = n * pio2_3;
r = t - w;
w = n * pio2_3t - ((t - r) - w);
y[0] = r - w;
}
}
}
y[1] = (r - y[0]) - w;
if (x < 0)
{
y[0] = -y[0];
y[1] = -y[1];
return -i;
}
else
{
return i;
}
}
/* all other (large) arguments */
if (isnanf (x) || isinff (x))
{
y[0] = y[1] = x - x;
return 0;
}
/* set z = scalbn(|x|,ilogb(x)-7) */
e0 = __float_and8 (ax / 128.0);
z = ax / e0;
tx[0] = __floorf (z);
z = (z - tx[0]) * two8;
tx[1] = __floorf (z);
z = (z - tx[1]) * two8;
tx[2] = __floorf (z);
nx = 3;
while (tx[nx - 1] == zero)
nx--;
i = __fp_kernel_rem_pio2f (tx, y, e0, nx);
if (x < 0)
{
y[0] = -y[0];
y[1] = -y[1];
return -i;
}
return i;
}

273
sysdeps/powerpc/fpu/k_rem_pio2f.c
View File

@ -1,273 +0,0 @@
/* k_rem_pio2f.c -- float version of e_rem_pio2.c
Copyright (C) 2011-2018 Free Software Foundation, Inc.
This file is part of the GNU C Library.
Contributed by Adhemerval Zanella <azanella@br.ibm.com>, 2011
The GNU C Library is free software; you can redistribute it and/or
modify it under the terms of the GNU Library General Public License as
published by the Free Software Foundation; either version 2 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
Library General Public License for more details.
You should have received a copy of the GNU Library General Public
License along with the GNU C Library; see the file COPYING.LIB. If
not, see <http://www.gnu.org/licenses/>. */
#include <math.h>
#include <math_private.h>
#include "s_float_bitwise.h"
static const float two_over_pi[] = {
1.62000000e+02, 2.49000000e+02, 1.31000000e+02, 1.10000000e+02,
7.80000000e+01, 6.80000000e+01, 2.10000000e+01, 4.10000000e+01,
2.52000000e+02, 3.90000000e+01, 8.70000000e+01, 2.09000000e+02,
2.45000000e+02, 5.20000000e+01, 2.21000000e+02, 1.92000000e+02,
2.19000000e+02, 9.80000000e+01, 1.49000000e+02, 1.53000000e+02,
6.00000000e+01, 6.70000000e+01, 1.44000000e+02, 6.50000000e+01,
2.54000000e+02, 8.10000000e+01, 9.90000000e+01, 1.71000000e+02,
2.22000000e+02, 1.87000000e+02, 1.97000000e+02, 9.70000000e+01,
1.83000000e+02, 3.60000000e+01, 1.10000000e+02, 5.80000000e+01,
6.60000000e+01, 7.70000000e+01, 2.10000000e+02, 2.24000000e+02,
6.00000000e+00, 7.30000000e+01, 4.60000000e+01, 2.34000000e+02,
9.00000000e+00, 2.09000000e+02, 1.46000000e+02, 2.80000000e+01,
2.54000000e+02, 2.90000000e+01, 2.35000000e+02, 2.80000000e+01,
1.77000000e+02, 4.10000000e+01, 1.67000000e+02, 6.20000000e+01,
2.32000000e+02, 1.30000000e+02, 5.30000000e+01, 2.45000000e+02,
4.60000000e+01, 1.87000000e+02, 6.80000000e+01, 1.32000000e+02,
2.33000000e+02, 1.56000000e+02, 1.12000000e+02, 3.80000000e+01,
1.80000000e+02, 9.50000000e+01, 1.26000000e+02, 6.50000000e+01,
5.70000000e+01, 1.45000000e+02, 2.14000000e+02, 5.70000000e+01,
1.31000000e+02, 8.30000000e+01, 5.70000000e+01, 2.44000000e+02,
1.56000000e+02, 1.32000000e+02, 9.50000000e+01, 1.39000000e+02,
1.89000000e+02, 2.49000000e+02, 4.00000000e+01, 5.90000000e+01,
3.10000000e+01, 2.48000000e+02, 1.51000000e+02, 2.55000000e+02,
2.22000000e+02, 5.00000000e+00, 1.52000000e+02, 1.50000000e+01,
2.39000000e+02, 4.70000000e+01, 1.70000000e+01, 1.39000000e+02,
9.00000000e+01, 1.00000000e+01, 1.09000000e+02, 3.10000000e+01,
1.09000000e+02, 5.40000000e+01, 1.26000000e+02, 2.07000000e+02,
3.90000000e+01, 2.03000000e+02, 9.00000000e+00, 1.83000000e+02,
7.90000000e+01, 7.00000000e+01, 6.30000000e+01, 1.02000000e+02,
1.58000000e+02, 9.50000000e+01, 2.34000000e+02, 4.50000000e+01,
1.17000000e+02, 3.90000000e+01, 1.86000000e+02, 1.99000000e+02,
2.35000000e+02, 2.29000000e+02, 2.41000000e+02, 1.23000000e+02,
6.10000000e+01, 7.00000000e+00, 5.70000000e+01, 2.47000000e+02,
1.38000000e+02, 8.20000000e+01, 1.46000000e+02, 2.34000000e+02,
1.07000000e+02, 2.51000000e+02, 9.50000000e+01, 1.77000000e+02,
3.10000000e+01, 1.41000000e+02, 9.30000000e+01, 8.00000000e+00,
8.60000000e+01, 3.00000000e+00, 4.80000000e+01, 7.00000000e+01,
2.52000000e+02, 1.23000000e+02, 1.07000000e+02, 1.71000000e+02,
2.40000000e+02, 2.07000000e+02, 1.88000000e+02, 3.20000000e+01,
1.54000000e+02, 2.44000000e+02, 5.40000000e+01, 2.90000000e+01,
1.69000000e+02, 2.27000000e+02, 1.45000000e+02, 9.70000000e+01,
9.40000000e+01, 2.30000000e+02, 2.70000000e+01, 8.00000000e+00,
1.01000000e+02, 1.53000000e+02, 1.33000000e+02, 9.50000000e+01,
2.00000000e+01, 1.60000000e+02, 1.04000000e+02, 6.40000000e+01,
1.41000000e+02, 2.55000000e+02, 2.16000000e+02, 1.28000000e+02,
7.70000000e+01, 1.15000000e+02, 3.90000000e+01, 4.90000000e+01,
6.00000000e+00, 6.00000000e+00, 2.10000000e+01, 8.60000000e+01,
2.02000000e+02, 1.15000000e+02, 1.68000000e+02, 2.01000000e+02,
9.60000000e+01, 2.26000000e+02, 1.23000000e+02, 1.92000000e+02,
1.40000000e+02, 1.07000000e+02
};
static const float PIo2[] = {
1.5703125000e+00, /* 0x3fc90000 */
4.5776367188e-04, /* 0x39f00000 */
2.5987625122e-05, /* 0x37da0000 */
7.5437128544e-08, /* 0x33a20000 */
6.0026650317e-11, /* 0x2e840000 */
7.3896444519e-13, /* 0x2b500000 */
5.3845816694e-15, /* 0x27c20000 */
5.6378512969e-18, /* 0x22d00000 */
8.3009228831e-20, /* 0x1fc40000 */
3.2756352257e-22, /* 0x1bc60000 */
6.3331015649e-25, /* 0x17440000 */
};
static const float zero = 0.0000000000e+00;
static const float one = 1.0000000000;
static const float twon8 = 3.9062500000e-03;
static const float two8 = 2.5600000000e+02;
int32_t
__fp_kernel_rem_pio2f (float *x, float *y, float e0, int32_t nx)
{
int32_t jz, jx, jv, jp, jk, carry, n, iq[20], i, j, k, m, q0, ih, exp;
float z, fw, f[20], fq[20], q[20];
/* initialize jk */
jp = jk = 9;
/* determine jx,jv,q0, note that 3>q0 */
jx = nx - 1;
exp = __float_get_exp (e0) - 127;
jv = (exp - 3) / 8;
if (jv < 0)
jv = 0;
q0 = exp - 8 * (jv + 1);
/* set up f[0] to f[jx+jk] where f[jx+jk] = two_over_pi[jv+jk] */
j = jv - jx;
m = jx + jk;
for (i = 0; i <= m; i++, j++)
f[i] = (j < 0) ? zero : two_over_pi[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 = __truncf (twon8 * z);
iq[i] = (int32_t) (z - two8 * fw);
z = q[j - 1] + fw;
}
/* compute n */
z = __scalbnf (z, q0); /* actual value of z */
z -= 8.0 * __floorf (z * 0.125); /* trim off integer >= 8 */
n = (int32_t) z;
z -= __truncf (z);
ih = 0;
if (q0 > 0)
{ /* need iq[jz-1] to determine n */
i = (iq[jz - 1] >> (8 - q0));
n += i;
iq[jz - 1] -= i << (8 - q0);
ih = iq[jz - 1] >> (7 - q0);
}
else if (q0 == 0)
ih = iq[jz - 1] >> 7;
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] = 0x100 - j;
}
}
else
iq[i] = 0xff - j;
}
if (q0 > 0)
{ /* rare case: chance is 1 in 12 */
switch (q0)
{
case 1:
iq[jz - 1] &= 0x7f;
break;
case 2:
iq[jz - 1] &= 0x3f;
break;
}
}
if (ih == 2)
{
z = one - z;
if (carry != 0)
z -= __scalbnf (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] = two_over_pi[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 -= 8;
while (iq[jz] == 0)
{
jz--;
q0 -= 8;
}
}
else
{ /* break z into 8-bit if necessary */
z = __scalbnf (z, -q0);
if (z >= two8)
{
fw = __truncf (twon8 * z);
iq[jz] = (int32_t) (z - two8 * fw);
jz += 1;
q0 += 8;
iq[jz] = (int32_t) fw;
}
else
iq[jz] = (int32_t) z;
}
/* convert integer "bit" chunk to floating-point value */
fw = __scalbnf (one, q0);
for (i = jz; i >= 0; i--)
{
q[i] = fw * (float) iq[i];
fw *= twon8;
}
/* 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[] */
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;
return n & 7;
}