125c28e83SPiotr Jasiukajtis /*
225c28e83SPiotr Jasiukajtis * CDDL HEADER START
325c28e83SPiotr Jasiukajtis *
425c28e83SPiotr Jasiukajtis * The contents of this file are subject to the terms of the
525c28e83SPiotr Jasiukajtis * Common Development and Distribution License (the "License").
625c28e83SPiotr Jasiukajtis * You may not use this file except in compliance with the License.
725c28e83SPiotr Jasiukajtis *
825c28e83SPiotr Jasiukajtis * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
925c28e83SPiotr Jasiukajtis * or http://www.opensolaris.org/os/licensing.
1025c28e83SPiotr Jasiukajtis * See the License for the specific language governing permissions
1125c28e83SPiotr Jasiukajtis * and limitations under the License.
1225c28e83SPiotr Jasiukajtis *
1325c28e83SPiotr Jasiukajtis * When distributing Covered Code, include this CDDL HEADER in each
1425c28e83SPiotr Jasiukajtis * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
1525c28e83SPiotr Jasiukajtis * If applicable, add the following below this CDDL HEADER, with the
1625c28e83SPiotr Jasiukajtis * fields enclosed by brackets "[]" replaced with your own identifying
1725c28e83SPiotr Jasiukajtis * information: Portions Copyright [yyyy] [name of copyright owner]
1825c28e83SPiotr Jasiukajtis *
1925c28e83SPiotr Jasiukajtis * CDDL HEADER END
2025c28e83SPiotr Jasiukajtis */
2125c28e83SPiotr Jasiukajtis /*
2225c28e83SPiotr Jasiukajtis * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
2325c28e83SPiotr Jasiukajtis */
2425c28e83SPiotr Jasiukajtis /*
2525c28e83SPiotr Jasiukajtis * Copyright 2005 Sun Microsystems, Inc. All rights reserved.
2625c28e83SPiotr Jasiukajtis * Use is subject to license terms.
2725c28e83SPiotr Jasiukajtis */
2825c28e83SPiotr Jasiukajtis
2925c28e83SPiotr Jasiukajtis /*
3025c28e83SPiotr Jasiukajtis * int __rem_pio2m(x,y,e0,nx,prec,ipio2)
3125c28e83SPiotr Jasiukajtis * double x[],y[]; int e0,nx,prec; const int ipio2[];
3225c28e83SPiotr Jasiukajtis *
3325c28e83SPiotr Jasiukajtis * __rem_pio2m return the last three digits of N with
3425c28e83SPiotr Jasiukajtis * y = x - N*pi/2
3525c28e83SPiotr Jasiukajtis * so that |y| < pi/4.
3625c28e83SPiotr Jasiukajtis *
3725c28e83SPiotr Jasiukajtis * The method is to compute the integer (mod 8) and fraction parts of
3825c28e83SPiotr Jasiukajtis * (2/pi)*x without doing the full multiplication. In general we
3925c28e83SPiotr Jasiukajtis * skip the part of the product that are known to be a huge integer (
4025c28e83SPiotr Jasiukajtis * more accurately, = 0 mod 8 ). Thus the number of operations are
4125c28e83SPiotr Jasiukajtis * independent of the exponent of the input.
4225c28e83SPiotr Jasiukajtis *
4325c28e83SPiotr Jasiukajtis * (2/PI) is represented by an array of 24-bit integers in ipio2[].
4425c28e83SPiotr Jasiukajtis * Here PI could as well be a machine value pi.
4525c28e83SPiotr Jasiukajtis *
4625c28e83SPiotr Jasiukajtis * Input parameters:
47*7c94ff60SToomas Soome * x[] The input value (must be positive) is broken into nx
4825c28e83SPiotr Jasiukajtis * pieces of 24-bit integers in double precision format.
4925c28e83SPiotr Jasiukajtis * x[i] will be the i-th 24 bit of x. The scaled exponent
5025c28e83SPiotr Jasiukajtis * of x[0] is given in input parameter e0 (i.e., x[0]*2^e0
5125c28e83SPiotr Jasiukajtis * match x's up to 24 bits.
5225c28e83SPiotr Jasiukajtis *
5325c28e83SPiotr Jasiukajtis * Example of breaking a double z into x[0]+x[1]+x[2]:
5425c28e83SPiotr Jasiukajtis * e0 = ilogb(z)-23
5525c28e83SPiotr Jasiukajtis * z = scalbn(z,-e0)
5625c28e83SPiotr Jasiukajtis * for i = 0,1,2
5725c28e83SPiotr Jasiukajtis * x[i] = floor(z)
5825c28e83SPiotr Jasiukajtis * z = (z-x[i])*2**24
5925c28e83SPiotr Jasiukajtis *
6025c28e83SPiotr Jasiukajtis *
6125c28e83SPiotr Jasiukajtis * y[] ouput result in an array of double precision numbers.
6225c28e83SPiotr Jasiukajtis * The dimension of y[] is:
6325c28e83SPiotr Jasiukajtis * 24-bit precision 1
6425c28e83SPiotr Jasiukajtis * 53-bit precision 2
6525c28e83SPiotr Jasiukajtis * 64-bit precision 2
6625c28e83SPiotr Jasiukajtis * 113-bit precision 3
6725c28e83SPiotr Jasiukajtis * The actual value is the sum of them. Thus for 113-bit
6825c28e83SPiotr Jasiukajtis * precsion, one may have to do something like:
6925c28e83SPiotr Jasiukajtis *
7025c28e83SPiotr Jasiukajtis * long double t,w,r_head, r_tail;
7125c28e83SPiotr Jasiukajtis * t = (long double)y[2] + (long double)y[1];
7225c28e83SPiotr Jasiukajtis * w = (long double)y[0];
7325c28e83SPiotr Jasiukajtis * r_head = t+w;
7425c28e83SPiotr Jasiukajtis * r_tail = w - (r_head - t);
7525c28e83SPiotr Jasiukajtis *
7625c28e83SPiotr Jasiukajtis * e0 The exponent of x[0]
7725c28e83SPiotr Jasiukajtis *
7825c28e83SPiotr Jasiukajtis * nx dimension of x[]
7925c28e83SPiotr Jasiukajtis *
80*7c94ff60SToomas Soome * prec an interger indicating the precision:
8125c28e83SPiotr Jasiukajtis * 0 24 bits (single)
8225c28e83SPiotr Jasiukajtis * 1 53 bits (double)
8325c28e83SPiotr Jasiukajtis * 2 64 bits (extended)
8425c28e83SPiotr Jasiukajtis * 3 113 bits (quad)
8525c28e83SPiotr Jasiukajtis *
8625c28e83SPiotr Jasiukajtis * ipio2[]
8725c28e83SPiotr Jasiukajtis * integer array, contains the (24*i)-th to (24*i+23)-th
8825c28e83SPiotr Jasiukajtis * bit of 2/pi or 2/PI after binary point. The corresponding
8925c28e83SPiotr Jasiukajtis * floating value is
9025c28e83SPiotr Jasiukajtis *
9125c28e83SPiotr Jasiukajtis * ipio2[i] * 2^(-24(i+1)).
9225c28e83SPiotr Jasiukajtis *
9325c28e83SPiotr Jasiukajtis * External function:
9425c28e83SPiotr Jasiukajtis * double scalbn( ), floor( );
9525c28e83SPiotr Jasiukajtis *
9625c28e83SPiotr Jasiukajtis *
9725c28e83SPiotr Jasiukajtis * Here is the description of some local variables:
9825c28e83SPiotr Jasiukajtis *
99*7c94ff60SToomas Soome * jk jk+1 is the initial number of terms of ipio2[] needed
10025c28e83SPiotr Jasiukajtis * in the computation. The recommended value is 3,4,4,
10125c28e83SPiotr Jasiukajtis * 6 for single, double, extended,and quad.
10225c28e83SPiotr Jasiukajtis *
103*7c94ff60SToomas Soome * jz local integer variable indicating the number of
10425c28e83SPiotr Jasiukajtis * terms of ipio2[] used.
10525c28e83SPiotr Jasiukajtis *
10625c28e83SPiotr Jasiukajtis * jx nx - 1
10725c28e83SPiotr Jasiukajtis *
10825c28e83SPiotr Jasiukajtis * jv index for pointing to the suitable ipio2[] for the
10925c28e83SPiotr Jasiukajtis * computation. In general, we want
11025c28e83SPiotr Jasiukajtis * ( 2^e0*x[0] * ipio2[jv-1]*2^(-24jv) )/8
11125c28e83SPiotr Jasiukajtis * is an integer. Thus
11225c28e83SPiotr Jasiukajtis * e0-3-24*jv >= 0 or (e0-3)/24 >= jv
11325c28e83SPiotr Jasiukajtis * Hence jv = max(0,(e0-3)/24).
11425c28e83SPiotr Jasiukajtis *
11525c28e83SPiotr Jasiukajtis * jp jp+1 is the number of terms in pio2[] needed, jp = jk.
11625c28e83SPiotr Jasiukajtis *
117*7c94ff60SToomas Soome * q[] double array with integral value, representing the
11825c28e83SPiotr Jasiukajtis * 24-bits chunk of the product of x and 2/pi.
11925c28e83SPiotr Jasiukajtis *
12025c28e83SPiotr Jasiukajtis * q0 the corresponding exponent of q[0]. Note that the
12125c28e83SPiotr Jasiukajtis * exponent for q[i] would be q0-24*i.
12225c28e83SPiotr Jasiukajtis *
12325c28e83SPiotr Jasiukajtis * pio2[] double precision array, obtained by cutting pi/2
12425c28e83SPiotr Jasiukajtis * into 24 bits chunks.
12525c28e83SPiotr Jasiukajtis *
12625c28e83SPiotr Jasiukajtis * f[] ipio2[] in floating point
12725c28e83SPiotr Jasiukajtis *
12825c28e83SPiotr Jasiukajtis * iq[] integer array by breaking up q[] in 24-bits chunk.
12925c28e83SPiotr Jasiukajtis *
13025c28e83SPiotr Jasiukajtis * fq[] final product of x*(2/pi) in fq[0],..,fq[jk]
13125c28e83SPiotr Jasiukajtis *
13225c28e83SPiotr Jasiukajtis * ih integer. If >0 it indicats q[] is >= 0.5, hence
13325c28e83SPiotr Jasiukajtis * it also indicates the *sign* of the result.
13425c28e83SPiotr Jasiukajtis *
13525c28e83SPiotr Jasiukajtis */
13625c28e83SPiotr Jasiukajtis
137*7c94ff60SToomas Soome #include <assert.h>
13825c28e83SPiotr Jasiukajtis #include "libm.h"
13925c28e83SPiotr Jasiukajtis
14025c28e83SPiotr Jasiukajtis #if defined(__i386) && !defined(__amd64)
14125c28e83SPiotr Jasiukajtis extern int __swapRP(int);
14225c28e83SPiotr Jasiukajtis #endif
14325c28e83SPiotr Jasiukajtis
14425c28e83SPiotr Jasiukajtis static const int init_jk[] = { 3, 4, 4, 6 }; /* initial value for jk */
14525c28e83SPiotr Jasiukajtis
14625c28e83SPiotr Jasiukajtis static const double pio2[] = {
14725c28e83SPiotr Jasiukajtis 1.57079625129699707031e+00,
14825c28e83SPiotr Jasiukajtis 7.54978941586159635335e-08,
14925c28e83SPiotr Jasiukajtis 5.39030252995776476554e-15,
15025c28e83SPiotr Jasiukajtis 3.28200341580791294123e-22,
15125c28e83SPiotr Jasiukajtis 1.27065575308067607349e-29,
15225c28e83SPiotr Jasiukajtis 1.22933308981111328932e-36,
15325c28e83SPiotr Jasiukajtis 2.73370053816464559624e-44,
15425c28e83SPiotr Jasiukajtis 2.16741683877804819444e-51,
15525c28e83SPiotr Jasiukajtis };
15625c28e83SPiotr Jasiukajtis
15725c28e83SPiotr Jasiukajtis static const double
15825c28e83SPiotr Jasiukajtis zero = 0.0,
15925c28e83SPiotr Jasiukajtis one = 1.0,
16025c28e83SPiotr Jasiukajtis half = 0.5,
16125c28e83SPiotr Jasiukajtis eight = 8.0,
16225c28e83SPiotr Jasiukajtis eighth = 0.125,
16325c28e83SPiotr Jasiukajtis two24 = 16777216.0,
16425c28e83SPiotr Jasiukajtis twon24 = 5.960464477539062500E-8;
16525c28e83SPiotr Jasiukajtis
16625c28e83SPiotr Jasiukajtis int
__rem_pio2m(double * x,double * y,int e0,int nx,int prec,const int * ipio2)16725c28e83SPiotr Jasiukajtis __rem_pio2m(double *x, double *y, int e0, int nx, int prec, const int *ipio2)
16825c28e83SPiotr Jasiukajtis {
16925c28e83SPiotr Jasiukajtis int jz, jx, jv, jp, jk, carry, n, iq[20];
17025c28e83SPiotr Jasiukajtis int i, j, k, m, q0, ih;
17125c28e83SPiotr Jasiukajtis double z, fw, f[20], fq[20], q[20];
17225c28e83SPiotr Jasiukajtis #if defined(__i386) && !defined(__amd64)
17325c28e83SPiotr Jasiukajtis int rp;
17425c28e83SPiotr Jasiukajtis
17525c28e83SPiotr Jasiukajtis rp = __swapRP(fp_extended);
17625c28e83SPiotr Jasiukajtis #endif
17725c28e83SPiotr Jasiukajtis
178*7c94ff60SToomas Soome fq[0] = NAN; /* Make gcc happy */
17925c28e83SPiotr Jasiukajtis /* initialize jk */
18025c28e83SPiotr Jasiukajtis jp = jk = init_jk[prec];
18125c28e83SPiotr Jasiukajtis
18225c28e83SPiotr Jasiukajtis /* determine jx,jv,q0, note that 3>q0 */
18325c28e83SPiotr Jasiukajtis jx = nx - 1;
18425c28e83SPiotr Jasiukajtis jv = (e0 - 3) / 24;
18525c28e83SPiotr Jasiukajtis if (jv < 0)
18625c28e83SPiotr Jasiukajtis jv = 0;
18725c28e83SPiotr Jasiukajtis q0 = e0 - 24 * (jv + 1);
18825c28e83SPiotr Jasiukajtis
18925c28e83SPiotr Jasiukajtis /* set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk] */
19025c28e83SPiotr Jasiukajtis j = jv - jx;
19125c28e83SPiotr Jasiukajtis m = jx + jk;
19225c28e83SPiotr Jasiukajtis for (i = 0; i <= m; i++, j++)
19325c28e83SPiotr Jasiukajtis f[i] = (j < 0)? zero : (double)ipio2[j];
19425c28e83SPiotr Jasiukajtis
19525c28e83SPiotr Jasiukajtis /* compute q[0],q[1],...q[jk] */
19625c28e83SPiotr Jasiukajtis for (i = 0; i <= jk; i++) {
19725c28e83SPiotr Jasiukajtis for (j = 0, fw = zero; j <= jx; j++)
19825c28e83SPiotr Jasiukajtis fw += x[j] * f[jx+i-j];
19925c28e83SPiotr Jasiukajtis q[i] = fw;
20025c28e83SPiotr Jasiukajtis }
20125c28e83SPiotr Jasiukajtis
20225c28e83SPiotr Jasiukajtis jz = jk;
20325c28e83SPiotr Jasiukajtis recompute:
20425c28e83SPiotr Jasiukajtis /* distill q[] into iq[] reversingly */
20525c28e83SPiotr Jasiukajtis for (i = 0, j = jz, z = q[jz]; j > 0; i++, j--) {
20625c28e83SPiotr Jasiukajtis fw = (double)((int)(twon24 * z));
20725c28e83SPiotr Jasiukajtis iq[i] = (int)(z - two24 * fw);
20825c28e83SPiotr Jasiukajtis z = q[j-1] + fw;
20925c28e83SPiotr Jasiukajtis }
21025c28e83SPiotr Jasiukajtis
21125c28e83SPiotr Jasiukajtis /* compute n */
21225c28e83SPiotr Jasiukajtis z = scalbn(z, q0); /* actual value of z */
21325c28e83SPiotr Jasiukajtis z -= eight * floor(z * eighth); /* trim off integer >= 8 */
21425c28e83SPiotr Jasiukajtis n = (int)z;
21525c28e83SPiotr Jasiukajtis z -= (double)n;
21625c28e83SPiotr Jasiukajtis ih = 0;
21725c28e83SPiotr Jasiukajtis if (q0 > 0) { /* need iq[jz-1] to determine n */
21825c28e83SPiotr Jasiukajtis i = (iq[jz-1] >> (24 - q0));
21925c28e83SPiotr Jasiukajtis n += i;
22025c28e83SPiotr Jasiukajtis iq[jz-1] -= i << (24 - q0);
22125c28e83SPiotr Jasiukajtis ih = iq[jz-1] >> (23 - q0);
22225c28e83SPiotr Jasiukajtis } else if (q0 == 0) {
22325c28e83SPiotr Jasiukajtis ih = iq[jz-1] >> 23;
22425c28e83SPiotr Jasiukajtis } else if (z >= half) {
22525c28e83SPiotr Jasiukajtis ih = 2;
22625c28e83SPiotr Jasiukajtis }
22725c28e83SPiotr Jasiukajtis
22825c28e83SPiotr Jasiukajtis if (ih > 0) { /* q > 0.5 */
22925c28e83SPiotr Jasiukajtis n += 1;
23025c28e83SPiotr Jasiukajtis carry = 0;
23125c28e83SPiotr Jasiukajtis for (i = 0; i < jz; i++) { /* compute 1-q */
23225c28e83SPiotr Jasiukajtis j = iq[i];
23325c28e83SPiotr Jasiukajtis if (carry == 0) {
23425c28e83SPiotr Jasiukajtis if (j != 0) {
23525c28e83SPiotr Jasiukajtis carry = 1;
23625c28e83SPiotr Jasiukajtis iq[i] = 0x1000000 - j;
23725c28e83SPiotr Jasiukajtis }
23825c28e83SPiotr Jasiukajtis } else {
23925c28e83SPiotr Jasiukajtis iq[i] = 0xffffff - j;
24025c28e83SPiotr Jasiukajtis }
24125c28e83SPiotr Jasiukajtis }
24225c28e83SPiotr Jasiukajtis if (q0 > 0) { /* rare case: chance is 1 in 12 */
24325c28e83SPiotr Jasiukajtis switch (q0) {
24425c28e83SPiotr Jasiukajtis case 1:
24525c28e83SPiotr Jasiukajtis iq[jz-1] &= 0x7fffff;
24625c28e83SPiotr Jasiukajtis break;
24725c28e83SPiotr Jasiukajtis case 2:
24825c28e83SPiotr Jasiukajtis iq[jz-1] &= 0x3fffff;
24925c28e83SPiotr Jasiukajtis break;
25025c28e83SPiotr Jasiukajtis }
25125c28e83SPiotr Jasiukajtis }
25225c28e83SPiotr Jasiukajtis if (ih == 2) {
25325c28e83SPiotr Jasiukajtis z = one - z;
25425c28e83SPiotr Jasiukajtis if (carry != 0)
25525c28e83SPiotr Jasiukajtis z -= scalbn(one, q0);
25625c28e83SPiotr Jasiukajtis }
25725c28e83SPiotr Jasiukajtis }
25825c28e83SPiotr Jasiukajtis
25925c28e83SPiotr Jasiukajtis /* check if recomputation is needed */
26025c28e83SPiotr Jasiukajtis if (z == zero) {
26125c28e83SPiotr Jasiukajtis j = 0;
26225c28e83SPiotr Jasiukajtis for (i = jz - 1; i >= jk; i--)
26325c28e83SPiotr Jasiukajtis j |= iq[i];
26425c28e83SPiotr Jasiukajtis if (j == 0) { /* need recomputation */
26525c28e83SPiotr Jasiukajtis /* set k to no. of terms needed */
26625c28e83SPiotr Jasiukajtis for (k = 1; iq[jk-k] == 0; k++)
26725c28e83SPiotr Jasiukajtis ;
26825c28e83SPiotr Jasiukajtis
26925c28e83SPiotr Jasiukajtis /* add q[jz+1] to q[jz+k] */
27025c28e83SPiotr Jasiukajtis for (i = jz + 1; i <= jz + k; i++) {
27125c28e83SPiotr Jasiukajtis f[jx+i] = (double)ipio2[jv+i];
27225c28e83SPiotr Jasiukajtis for (j = 0, fw = zero; j <= jx; j++)
27325c28e83SPiotr Jasiukajtis fw += x[j] * f[jx+i-j];
27425c28e83SPiotr Jasiukajtis q[i] = fw;
27525c28e83SPiotr Jasiukajtis }
27625c28e83SPiotr Jasiukajtis jz += k;
27725c28e83SPiotr Jasiukajtis goto recompute;
27825c28e83SPiotr Jasiukajtis }
27925c28e83SPiotr Jasiukajtis }
28025c28e83SPiotr Jasiukajtis
28125c28e83SPiotr Jasiukajtis /* cut out zero terms */
28225c28e83SPiotr Jasiukajtis if (z == zero) {
28325c28e83SPiotr Jasiukajtis jz -= 1;
28425c28e83SPiotr Jasiukajtis q0 -= 24;
28525c28e83SPiotr Jasiukajtis while (iq[jz] == 0) {
28625c28e83SPiotr Jasiukajtis jz--;
28725c28e83SPiotr Jasiukajtis q0 -= 24;
28825c28e83SPiotr Jasiukajtis }
28925c28e83SPiotr Jasiukajtis } else { /* break z into 24-bit if neccessary */
29025c28e83SPiotr Jasiukajtis z = scalbn(z, -q0);
29125c28e83SPiotr Jasiukajtis if (z >= two24) {
29225c28e83SPiotr Jasiukajtis fw = (double)((int)(twon24 * z));
29325c28e83SPiotr Jasiukajtis iq[jz] = (int)(z - two24 * fw);
29425c28e83SPiotr Jasiukajtis jz += 1;
29525c28e83SPiotr Jasiukajtis q0 += 24;
29625c28e83SPiotr Jasiukajtis iq[jz] = (int)fw;
29725c28e83SPiotr Jasiukajtis } else {
29825c28e83SPiotr Jasiukajtis iq[jz] = (int)z;
29925c28e83SPiotr Jasiukajtis }
30025c28e83SPiotr Jasiukajtis }
30125c28e83SPiotr Jasiukajtis
30225c28e83SPiotr Jasiukajtis /* convert integer "bit" chunk to floating-point value */
30325c28e83SPiotr Jasiukajtis fw = scalbn(one, q0);
30425c28e83SPiotr Jasiukajtis for (i = jz; i >= 0; i--) {
30525c28e83SPiotr Jasiukajtis q[i] = fw * (double)iq[i];
30625c28e83SPiotr Jasiukajtis fw *= twon24;
30725c28e83SPiotr Jasiukajtis }
30825c28e83SPiotr Jasiukajtis
30925c28e83SPiotr Jasiukajtis /* compute pio2[0,...,jp]*q[jz,...,0] */
31025c28e83SPiotr Jasiukajtis for (i = jz; i >= 0; i--) {
31125c28e83SPiotr Jasiukajtis for (fw = zero, k = 0; k <= jp && k <= jz - i; k++)
31225c28e83SPiotr Jasiukajtis fw += pio2[k] * q[i+k];
31325c28e83SPiotr Jasiukajtis fq[jz-i] = fw;
31425c28e83SPiotr Jasiukajtis }
31525c28e83SPiotr Jasiukajtis
31625c28e83SPiotr Jasiukajtis /* compress fq[] into y[] */
31725c28e83SPiotr Jasiukajtis switch (prec) {
31825c28e83SPiotr Jasiukajtis case 0:
31925c28e83SPiotr Jasiukajtis fw = zero;
32025c28e83SPiotr Jasiukajtis for (i = jz; i >= 0; i--)
32125c28e83SPiotr Jasiukajtis fw += fq[i];
32225c28e83SPiotr Jasiukajtis y[0] = (ih == 0)? fw : -fw;
32325c28e83SPiotr Jasiukajtis break;
32425c28e83SPiotr Jasiukajtis
32525c28e83SPiotr Jasiukajtis case 1:
32625c28e83SPiotr Jasiukajtis case 2:
32725c28e83SPiotr Jasiukajtis fw = zero;
32825c28e83SPiotr Jasiukajtis for (i = jz; i >= 0; i--)
32925c28e83SPiotr Jasiukajtis fw += fq[i];
33025c28e83SPiotr Jasiukajtis y[0] = (ih == 0)? fw : -fw;
331*7c94ff60SToomas Soome
332*7c94ff60SToomas Soome assert(!isnan(fq[0]));
33325c28e83SPiotr Jasiukajtis fw = fq[0] - fw;
33425c28e83SPiotr Jasiukajtis for (i = 1; i <= jz; i++)
33525c28e83SPiotr Jasiukajtis fw += fq[i];
33625c28e83SPiotr Jasiukajtis y[1] = (ih == 0)? fw : -fw;
33725c28e83SPiotr Jasiukajtis break;
33825c28e83SPiotr Jasiukajtis
33925c28e83SPiotr Jasiukajtis default:
34025c28e83SPiotr Jasiukajtis for (i = jz; i > 0; i--) {
34125c28e83SPiotr Jasiukajtis fw = fq[i-1] + fq[i];
34225c28e83SPiotr Jasiukajtis fq[i] += fq[i-1] - fw;
34325c28e83SPiotr Jasiukajtis fq[i-1] = fw;
34425c28e83SPiotr Jasiukajtis }
34525c28e83SPiotr Jasiukajtis for (i = jz; i > 1; i--) {
34625c28e83SPiotr Jasiukajtis fw = fq[i-1] + fq[i];
34725c28e83SPiotr Jasiukajtis fq[i] += fq[i-1] - fw;
34825c28e83SPiotr Jasiukajtis fq[i-1] = fw;
34925c28e83SPiotr Jasiukajtis }
35025c28e83SPiotr Jasiukajtis for (fw = zero, i = jz; i >= 2; i--)
35125c28e83SPiotr Jasiukajtis fw += fq[i];
35225c28e83SPiotr Jasiukajtis if (ih == 0) {
35325c28e83SPiotr Jasiukajtis y[0] = fq[0];
35425c28e83SPiotr Jasiukajtis y[1] = fq[1];
35525c28e83SPiotr Jasiukajtis y[2] = fw;
35625c28e83SPiotr Jasiukajtis } else {
35725c28e83SPiotr Jasiukajtis y[0] = -fq[0];
35825c28e83SPiotr Jasiukajtis y[1] = -fq[1];
35925c28e83SPiotr Jasiukajtis y[2] = -fw;
36025c28e83SPiotr Jasiukajtis }
36125c28e83SPiotr Jasiukajtis }
36225c28e83SPiotr Jasiukajtis
36325c28e83SPiotr Jasiukajtis #if defined(__i386) && !defined(__amd64)
36425c28e83SPiotr Jasiukajtis (void) __swapRP(rp);
36525c28e83SPiotr Jasiukajtis #endif
36625c28e83SPiotr Jasiukajtis return (n & 7);
36725c28e83SPiotr Jasiukajtis }
368