src/libm/k_rem_pio2.c
author Ryan C. Gordon <icculus@icculus.org>
Mon, 21 May 2018 12:00:21 -0400
changeset 11994 8e094f91ddab
parent 11838 5ef6e4e70103
child 12068 ce88faaf8bd2
permissions -rw-r--r--
thread: fixed compiler warnings on non-Linux systems that use pthread.

(static function rtkit_setpriority was unused, moved it in with rest of
__LINUX__ section.)
     1 /*
     2  * ====================================================
     3  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
     4  *
     5  * Developed at SunPro, a Sun Microsystems, Inc. business.
     6  * Permission to use, copy, modify, and distribute this
     7  * software is freely granted, provided that this notice
     8  * is preserved.
     9  * ====================================================
    10  */
    11 
    12 /*
    13  * __kernel_rem_pio2(x,y,e0,nx,prec,ipio2)
    14  * double x[],y[]; int e0,nx,prec; int ipio2[];
    15  *
    16  * __kernel_rem_pio2 return the last three digits of N with
    17  *		y = x - N*pi/2
    18  * so that |y| < pi/2.
    19  *
    20  * The method is to compute the integer (mod 8) and fraction parts of
    21  * (2/pi)*x without doing the full multiplication. In general we
    22  * skip the part of the product that are known to be a huge integer (
    23  * more accurately, = 0 mod 8 ). Thus the number of operations are
    24  * independent of the exponent of the input.
    25  *
    26  * (2/pi) is represented by an array of 24-bit integers in ipio2[].
    27  *
    28  * Input parameters:
    29  * 	x[]	The input value (must be positive) is broken into nx
    30  *		pieces of 24-bit integers in double precision format.
    31  *		x[i] will be the i-th 24 bit of x. The scaled exponent
    32  *		of x[0] is given in input parameter e0 (i.e., x[0]*2^e0
    33  *		match x's up to 24 bits.
    34  *
    35  *		Example of breaking a double positive z into x[0]+x[1]+x[2]:
    36  *			e0 = ilogb(z)-23
    37  *			z  = scalbn(z,-e0)
    38  *		for i = 0,1,2
    39  *			x[i] = floor(z)
    40  *			z    = (z-x[i])*2**24
    41  *
    42  *
    43  *	y[]	ouput result in an array of double precision numbers.
    44  *		The dimension of y[] is:
    45  *			24-bit  precision	1
    46  *			53-bit  precision	2
    47  *			64-bit  precision	2
    48  *			113-bit precision	3
    49  *		The actual value is the sum of them. Thus for 113-bit
    50  *		precison, one may have to do something like:
    51  *
    52  *		long double t,w,r_head, r_tail;
    53  *		t = (long double)y[2] + (long double)y[1];
    54  *		w = (long double)y[0];
    55  *		r_head = t+w;
    56  *		r_tail = w - (r_head - t);
    57  *
    58  *	e0	The exponent of x[0]
    59  *
    60  *	nx	dimension of x[]
    61  *
    62  *  	prec	an integer indicating the precision:
    63  *			0	24  bits (single)
    64  *			1	53  bits (double)
    65  *			2	64  bits (extended)
    66  *			3	113 bits (quad)
    67  *
    68  *	ipio2[]
    69  *		integer array, contains the (24*i)-th to (24*i+23)-th
    70  *		bit of 2/pi after binary point. The corresponding
    71  *		floating value is
    72  *
    73  *			ipio2[i] * 2^(-24(i+1)).
    74  *
    75  * External function:
    76  *	double scalbn(), floor();
    77  *
    78  *
    79  * Here is the description of some local variables:
    80  *
    81  * 	jk	jk+1 is the initial number of terms of ipio2[] needed
    82  *		in the computation. The recommended value is 2,3,4,
    83  *		6 for single, double, extended,and quad.
    84  *
    85  * 	jz	local integer variable indicating the number of
    86  *		terms of ipio2[] used.
    87  *
    88  *	jx	nx - 1
    89  *
    90  *	jv	index for pointing to the suitable ipio2[] for the
    91  *		computation. In general, we want
    92  *			( 2^e0*x[0] * ipio2[jv-1]*2^(-24jv) )/8
    93  *		is an integer. Thus
    94  *			e0-3-24*jv >= 0 or (e0-3)/24 >= jv
    95  *		Hence jv = max(0,(e0-3)/24).
    96  *
    97  *	jp	jp+1 is the number of terms in PIo2[] needed, jp = jk.
    98  *
    99  * 	q[]	double array with integral value, representing the
   100  *		24-bits chunk of the product of x and 2/pi.
   101  *
   102  *	q0	the corresponding exponent of q[0]. Note that the
   103  *		exponent for q[i] would be q0-24*i.
   104  *
   105  *	PIo2[]	double precision array, obtained by cutting pi/2
   106  *		into 24 bits chunks.
   107  *
   108  *	f[]	ipio2[] in floating point
   109  *
   110  *	iq[]	integer array by breaking up q[] in 24-bits chunk.
   111  *
   112  *	fq[]	final product of x*(2/pi) in fq[0],..,fq[jk]
   113  *
   114  *	ih	integer. If >0 it indicates q[] is >= 0.5, hence
   115  *		it also indicates the *sign* of the result.
   116  *
   117  */
   118 
   119 
   120 /*
   121  * Constants:
   122  * The hexadecimal values are the intended ones for the following
   123  * constants. The decimal values may be used, provided that the
   124  * compiler will convert from decimal to binary accurately enough
   125  * to produce the hexadecimal values shown.
   126  */
   127 
   128 #include "math_libm.h"
   129 #include "math_private.h"
   130 
   131 static const int init_jk[] = {2,3,4,6}; /* initial value for jk */
   132 
   133 static const double PIo2[] = {
   134   1.57079625129699707031e+00, /* 0x3FF921FB, 0x40000000 */
   135   7.54978941586159635335e-08, /* 0x3E74442D, 0x00000000 */
   136   5.39030252995776476554e-15, /* 0x3CF84698, 0x80000000 */
   137   3.28200341580791294123e-22, /* 0x3B78CC51, 0x60000000 */
   138   1.27065575308067607349e-29, /* 0x39F01B83, 0x80000000 */
   139   1.22933308981111328932e-36, /* 0x387A2520, 0x40000000 */
   140   2.73370053816464559624e-44, /* 0x36E38222, 0x80000000 */
   141   2.16741683877804819444e-51, /* 0x3569F31D, 0x00000000 */
   142 };
   143 
   144 static const double
   145 zero   = 0.0,
   146 one    = 1.0,
   147 two24   =  1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
   148 twon24  =  5.96046447753906250000e-08; /* 0x3E700000, 0x00000000 */
   149 
   150 int attribute_hidden __kernel_rem_pio2(double *x, double *y, int e0, int nx, int prec, const int32_t *ipio2)
   151 {
   152 	int32_t jz,jx,jv,jp,jk,carry,n,iq[20],i,j,k,m,q0,ih;
   153 	double z,fw,f[20],fq[20],q[20];
   154 
   155     /* initialize jk*/
   156 	jk = init_jk[prec];
   157 	jp = jk;
   158 
   159     /* determine jx,jv,q0, note that 3>q0 */
   160 	jx =  nx-1;
   161 	jv = (e0-3)/24; if(jv<0) jv=0;
   162 	q0 =  e0-24*(jv+1);
   163 
   164     /* set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk] */
   165 	j = jv-jx; m = jx+jk;
   166 	for(i=0;i<=m;i++,j++) f[i] = (j<0)? zero : (double) ipio2[j];
   167 
   168     /* compute q[0],q[1],...q[jk] */
   169 	for (i=0;i<=jk;i++) {
   170 	    for(j=0,fw=0.0;j<=jx;j++) fw += x[j]*f[jx+i-j];
   171 	    q[i] = fw;
   172 	}
   173 
   174 	jz = jk;
   175 recompute:
   176     /* distill q[] into iq[] reversingly */
   177 	for(i=0,j=jz,z=q[jz];j>0;i++,j--) {
   178 	    fw    =  (double)((int32_t)(twon24* z));
   179 	    iq[i] =  (int32_t)(z-two24*fw);
   180 	    z     =  q[j-1]+fw;
   181 	}
   182 
   183     /* compute n */
   184 	z  = scalbn(z,q0);		/* actual value of z */
   185 	z -= 8.0*floor(z*0.125);		/* trim off integer >= 8 */
   186 	n  = (int32_t) z;
   187 	z -= (double)n;
   188 	ih = 0;
   189 	if(q0>0) {	/* need iq[jz-1] to determine n */
   190 	    i  = (iq[jz-1]>>(24-q0)); n += i;
   191 	    iq[jz-1] -= i<<(24-q0);
   192 	    ih = iq[jz-1]>>(23-q0);
   193 	}
   194 	else if(q0==0) ih = iq[jz-1]>>23;
   195 	else if(z>=0.5) ih=2;
   196 
   197 	if(ih>0) {	/* q > 0.5 */
   198 	    n += 1; carry = 0;
   199 	    for(i=0;i<jz ;i++) {	/* compute 1-q */
   200 		j = iq[i];
   201 		if(carry==0) {
   202 		    if(j!=0) {
   203 			carry = 1; iq[i] = 0x1000000- j;
   204 		    }
   205 		} else  iq[i] = 0xffffff - j;
   206 	    }
   207 	    if(q0>0) {		/* rare case: chance is 1 in 12 */
   208 	        switch(q0) {
   209 	        case 1:
   210 	    	   iq[jz-1] &= 0x7fffff; break;
   211 	    	case 2:
   212 	    	   iq[jz-1] &= 0x3fffff; break;
   213 	        }
   214 	    }
   215 	    if(ih==2) {
   216 		z = one - z;
   217 		if(carry!=0) z -= scalbn(one,q0);
   218 	    }
   219 	}
   220 
   221     /* check if recomputation is needed */
   222 	if(z==zero) {
   223 	    j = 0;
   224 	    for (i=jz-1;i>=jk;i--) j |= iq[i];
   225 	    if(j==0) { /* need recomputation */
   226 		for(k=1;iq[jk-k]==0;k++);   /* k = no. of terms needed */
   227 
   228 		for(i=jz+1;i<=jz+k;i++) {   /* add q[jz+1] to q[jz+k] */
   229 		    f[jx+i] = (double) ipio2[jv+i];
   230 		    for(j=0,fw=0.0;j<=jx;j++) fw += x[j]*f[jx+i-j];
   231 		    q[i] = fw;
   232 		}
   233 		jz += k;
   234 		goto recompute;
   235 	    }
   236 	}
   237 
   238     /* chop off zero terms */
   239 	if(z==0.0) {
   240 	    jz -= 1; q0 -= 24;
   241 	    while(iq[jz]==0) { jz--; q0-=24;}
   242 	} else { /* break z into 24-bit if necessary */
   243 	    z = scalbn(z,-q0);
   244 	    if(z>=two24) {
   245 		fw = (double)((int32_t)(twon24*z));
   246 		iq[jz] = (int32_t)(z-two24*fw);
   247 		jz += 1; q0 += 24;
   248 		iq[jz] = (int32_t) fw;
   249 	    } else iq[jz] = (int32_t) z ;
   250 	}
   251 
   252     /* convert integer "bit" chunk to floating-point value */
   253 	fw = scalbn(one,q0);
   254 	for(i=jz;i>=0;i--) {
   255 	    q[i] = fw*(double)iq[i]; fw*=twon24;
   256 	}
   257 
   258     /* compute PIo2[0,...,jp]*q[jz,...,0] */
   259 	for(i=jz;i>=0;i--) {
   260 	    for(fw=0.0,k=0;k<=jp&&k<=jz-i;k++) fw += PIo2[k]*q[i+k];
   261 	    fq[jz-i] = fw;
   262 	}
   263 
   264     /* compress fq[] into y[] */
   265 	switch(prec) {
   266 	    case 0:
   267 		fw = 0.0;
   268 		for (i=jz;i>=0;i--) fw += fq[i];
   269 		y[0] = (ih==0)? fw: -fw;
   270 		break;
   271 	    case 1:
   272 	    case 2:
   273 		fw = 0.0;
   274 		for (i=jz;i>=0;i--) fw += fq[i];
   275 		y[0] = (ih==0)? fw: -fw;
   276 		fw = fq[0]-fw;
   277 		for (i=1;i<=jz;i++) fw += fq[i];
   278 		y[1] = (ih==0)? fw: -fw;
   279 		break;
   280 	    case 3:	/* painful */
   281 		for (i=jz;i>0;i--) {
   282 		    fw      = fq[i-1]+fq[i];
   283 		    fq[i]  += fq[i-1]-fw;
   284 		    fq[i-1] = fw;
   285 		}
   286 		for (i=jz;i>1;i--) {
   287 		    fw      = fq[i-1]+fq[i];
   288 		    fq[i]  += fq[i-1]-fw;
   289 		    fq[i-1] = fw;
   290 		}
   291 		for (fw=0.0,i=jz;i>=2;i--) fw += fq[i];
   292 		if(ih==0) {
   293 		    y[0] =  fq[0]; y[1] =  fq[1]; y[2] =  fw;
   294 		} else {
   295 		    y[0] = -fq[0]; y[1] = -fq[1]; y[2] = -fw;
   296 		}
   297 	}
   298 	return n&7;
   299 }