trig.c

操作系统SunOS 4.1.3版本的源码

	    }
	}

    /* compute sin, cos, tan, sincos */ 
	z = y1*y1;
	m  = ((int) fp_pi)<<3;
	if(k<2) { 		/* sin, cos */
	    n += k;
	    if((n&1)==0) {	/* sin on primary range */
		t = y1*z;
		ss=S[m]+z*(S[m+1]+z*(S[m+2]+z*(S[m+3]+z*(S[m+4]+z*S[m+5]))));
		if(j==0) w = y1 - t*ss; else w = y1 - (t*ss -(y2-0.5*(z*y2)));
		return ((n&2)==0)? w: -w;		/* n=0,2 */
	    } else {
		cc = (z*z)*(C[m]+z*(C[m+1]+z*(C[m+2]+
			z*(C[m+3]+z*(C[m+4]+z*C[m+5])))));
		hz = z*half;
		if(j!=0) cc -= y1*y2;
		if     (z>=thresh1) 	w = c5_8  - ((hz-c3_8 )-cc);
		else if(z>=thresh2) 	w = c13_16- ((hz-c3_16)-cc);
		else 			w = one   - (hz        -cc); 
		return ((n&2)==0)? w: -w;		/* n=1,3 */
	    }
	} else {		/* tan,sincos */
	    t =z*z;
	    ss=(S[m]+z*(S[m+1]+z*(S[m+2]+z*(S[m+3]+z*(S[m+4]+z*S[m+5])))));
	    cc=t*(C[m]+z*(C[m+1]+z*(C[m+2]+z*(C[m+3]+z*(C[m+4]+z*C[m+5])))));
	    hz = z*half;

	    if (k==2) {		/* tan */
	        ss *= z;
	        if     (z>=thresh1) 	w = c5_8  - ((hz-c3_8 )-cc);
	        else if(z>=thresh2) 	w = c13_16- ((hz-c3_16)-cc);
	        else 			w = one   - (hz        -cc); 
	        if ((n&1)==0) { 	/*  tan =  sin/cos */
		    if(j==0) return y1+y1*((hz-(ss+cc))/w);	else

    /* return y1+(y1*(y1*y2+(hz-(ss+cc))/w)+y2): rearrange to force ordering */
		    return y1-(y1*(((ss+cc)-hz)/w-y1*y2)-y2);
	        } else {		/*  tan = -cos/sin */
		    if(j==0) return -w/(y1-y1*ss); else
		    return -w/(y1-(y1*ss-w*(y2+y2*z)));
	        }
	    } else {	/* sincos: make sure to get same answers as sin,cos */
	        if(j!=0) cc -= y1*y2;
	        if     (z>=thresh1) 	w = c5_8  - ((hz-c3_8 )-cc);
	        else if(z>=thresh2) 	w = c13_16- ((hz-c3_16)-cc);
	        else 			w = one   - (hz        -cc); 
	        t = y1*z;
	        t = (j==0)? y1 - t*ss: y1 - (t*ss -(y2-hz*y2));
	        switch (n&3) {
		    case 0:  sincos_c = w;  return t;
		    case 1:  sincos_c = -t; return w;
		    case 2:  sincos_c = -w; return -t;
		    default: sincos_c = t;  return -w;
	        }
	    }
	}
}

static int dummy(x)
double x;
{
	return 1;
}

/*
 * huge argument reduction routine
 */
static long p_53[] = {
0xA2F983, 0x6E4E44, 0x16F3C4, 0xEA69B5, 0xD3E131, 0x60E1D2,
0xD7982A, 0xC031F5, 0xD67BCC, 0xFD1375, 0x60919B, 0x3FA0BB,
0x612ABB, 0x714F9B, 0x03DA8A, 0xC05948, 0xD023F4, 0x5AFA37,
0x51631D, 0xCD7A90, 0xC0474A, 0xF6A6F3, 0x1A52E1, 0x5C3927,
0x3ADA45, 0x4E2DB5, 0x64E8C4, 0x274A5B, 0xB74ADC, 0x1E6591,
0x2822BE, 0x4771F5, 0x12A63F, 0x83BD35, 0x2488CA, 0x1FE1BE,
0x42C21A, 0x682569, 0x2AFB91, 0x68ADE1, 0x4A42E5, 0x9BE357,
0xB79675, 0xCE998A, 0x83AF8B, 0xE645E6, 0xDF0789, 0x9E9747,
0xAA15FF, 0x358C3F, 0xAF3141, 0x72A3F7, 0x2BF1D4, 0xF3AD96,
0x7D759F, 0x257FCE, 0x29FB69, 0xB1B42C, 0xC32DE1, 0x8C0BBD,
0x31EC2F, 0x942026, 0x85DCE7, 0x653FF3, 0x136FA7, 0x0D7A5F,
};	/* 396 Hex digits (476 decimal) of 2/PI, PI = 53 bits of pi */

static long p_66[] = {
0xA2F983, 0x6E4E44, 0x152A00, 0x062BC4, 0x0DA276, 0xBED4C1,
0xFDF905, 0x5CD5BA, 0x767CEC, 0x1F80D6, 0xC26053, 0x3A0070,
0x107C2A, 0xF68EE9, 0x687B7A, 0xB990AA, 0x38DE4B, 0x96CFF3,
0x92735E, 0x8B34F6, 0x195BFC, 0x27F88E, 0xA93EC5, 0x3958A5,
0x3E5D13, 0x1C55A8, 0x5B4A8B, 0xA42E04, 0x12D105, 0x35580D,
0xF62347, 0x450900, 0xB98BCA, 0xF7E8A4, 0xA2E5D5, 0x69BC52,
0xF0381D, 0x1A0A88, 0xFE8714, 0x7F6735, 0xBB7D4D, 0xC6F642,
0xB27E80, 0x6191BF, 0xB6B750, 0x52776E, 0xD60FD0, 0x607DCC,
0x68BFAF, 0xED69FC, 0x6EB305, 0xD2557D, 0x25BDFB, 0x3E4AA1,
0x84472D, 0x8B0376, 0xF77740, 0xD290DF, 0x15EC8C, 0x45A5C3,
0x6181EF, 0xC5E7E8, 0xD8909C, 0xF62144, 0x298428, 0x6E5D9D,
};	/* 396 Hex digits (476 decimal) of 2/PI, PI = 66 bits of pi */

static long p_inf[] = {
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, 
};	/* 396 Hex digits (476 decimal) of 2/pi */

	/* p_##[i] contains the (24*i)-th to (24*i+23)-th bit
	 * of 2/pi or 2/PI after binary point. The corresponding 
	 * floating value fv[i] = fv(pr[i]) is
	 *	fv[i] = scalbn((double)pr[i],-24*(i+1)).
	 * Note that unless fv[i] = 0, we have
	 *	2**(-24*(i+1)) <= fv[i] < 2**(-24*i).
	 *	2**(-24*(i+1)) <= ulp(fv[i]) .
	 */

/*
 * Constants:
 * The hexadecimal values are the intended ones for the following constants.
 * The decimal values may be used, provided that the compiler will convert
 * from decimal to binary accurately enough to produce the hexadecimal values
 * shown.
 */
static double			
two24  = 16777216.0,
twon24 = 5.9604644775390625E-8,
p1     = 1.5707962512969970703E0,       /* first  24 bit of pi/2 */
p2     = 7.5497894158615963534E-8,      /* second 24 bit of pi/2 */
p3_inf = 5.3903025299577647655E-15,     /* third  24 bit of pi/2 */
p3_66  = 5.3903008358918702569E-15,     /* third  18 bit of pi/2 */
p3_53  = 5.3290705182007513940E-15;     /* third   5 bit of pi/2 */

static int argred(x,y1,y2) 	/* return the last three bits of N */
double x, *y1, *y2;
{
	long *pr;
	double x1,x2,x3,t,t1,t2,t3,p3,y,z,fvm1,fvm2;
	double q[10],fv[10];		
	double q1,q2,q3;
	int signx,k,n,mm=0,n1=1;
	unsigned long *py = (unsigned long *) &y;

    /* determine double ordering */
	if((*(1+(long *)&one))==0x3ff00000) n1=0;

    /* save sign bit and replace x by its absolute value */
	signx = signbit(x); 
	x = fabs(x);
	n = ilogb(x);
    /* Let  x <=| y stands for either y is zero or x is less than the least
     * non-zero significant bit of y. 
     * We will break x into three 24-bit pieces x1+x2+x3.
     * Here we choose k so that 
     * 		2**3 <=| x1*fv[k-1] 
     * where fv[k-1] is the floating value of pr[k-1].
     * Analysis: since 2**(-24*k) <=| fv[k-1] < 2**(-24*(k-1)) and
     *		       2**(n-23)  <=| x1      < 2**(n+1),
     *                 2**(n-23)  <=| x2*2**24< 2**(n+1),
     *                 2**(n-23)  <=| x3*2**48< 2**(n+1),
     * we have
     *		    2**(n-24k-23) <=| x1*fv[k-1] < 2**(n-24k+25).
     *		    2**(n-24k-23) <=| x2*fv[k-2] < 2**(n-24k+25).
     *		    2**(n-24k-23) <=| x3*fv[k-3] < 2**(n-24k+25).
     * Thus 
     *		    2**(n-24k-28) <=| x1*fv[k-1]+x2*fv[k-2]+x3*fv[k-3]
     * Thus n-24k-23 >= 3 or (n-26)/24 >= k.
     */
	k = (n-26)/24;
	if(k<=0) k = 0;			/* k must be >= 0 */
	else x = scalbn(x,-24*k);
    /* break x in three 24-bit pieces */
	y = x;
	py[n1] &= 0xffffffe0; 	/* 48 bit of x */
	x3 = y;
	py[n1] &= 0xe0000000; 	/* 24 bit of x */
	x1 = y;
	x2 = x3-y;
	x3 = x-x3;

    /* convert pr to fv */
	pr = p_inf; p3 = p3_inf;
	if(fp_pi==fp_pi_66) {pr = p_66;p3 = p3_66;}
	if(fp_pi==fp_pi_53) {pr = p_53;p3 = p3_53;}
	if(k>=1) {
	    fvm1 = (double)pr[k-1];
	    if(k> 1) fvm2 = (double)pr[k-2]*two24;
	}
	z     = twon24;
	fv[0] = ((double)pr[k])*z;
	z    *= twon24;
	fv[1] = ((double)pr[k+1])*z;
	z    *= twon24;
	fv[2] = ((double)pr[k+2])*z;
	z    *= twon24;
	fv[3] = ((double)pr[k+3])*z;
	z    *= twon24;
	fv[4] = ((double)pr[k+4])*z;

    /* compute q[0],q[1],...q[4] */
	if (x3==zero) {
	    if (x2==zero) {
		q[0] = x1*fv[0];
		q[1] = x1*fv[1];
		q[2] = x1*fv[2];
		q[3] = x1*fv[3];
		q[4] = x1*fv[4];
	    } else {
		if(k>0) q[0] = x1*fv[0] + x2*fvm1; else q[0] = x1*fv[0];
		q[1] = x1*fv[1]+x2*fv[0];
		q[2] = x1*fv[2]+x2*fv[1];
		q[3] = x1*fv[3]+x2*fv[2];
		q[4] = x1*fv[4]+x2*fv[3];
	    }
	} else {
	    if(k>=1) {
		if(k>1) q[0] = x1*fv[0]+x2*fvm1 +x3*fvm2;
		else    q[0] = x1*fv[0]+x2*fvm1 ;
	    	q[1] = x1*fv[1]+x2*fv[0]+x3*fvm1;
	    } else {
		q[0] = x1*fv[0];
	    	q[1] = x1*fv[1]+x2*fv[0];
	    }
	    q[2] = x1*fv[2]+x2*fv[1]+x3*fv[0];
	    q[3] = x1*fv[3]+x2*fv[2]+x3*fv[1];
	    q[4] = x1*fv[4]+x2*fv[3]+x3*fv[2];
	}

    /* trim off integer part that is bigger than 8 */
	q[0] -= 8.0*aint(q[0]*0.125);
	q[1] -= 8.0*aint(q[1]*0.125);
	n = q[2]*0.125; 
	if(n!=0) q[2] -= (double) (n<<3);

    /* default situation sum q[0] to q[4] */
    /* sum q[i] into t1 and t2 with integer part in mm and t1 < 0.5 */
	t   = q[4]; 
	t  += q[3]; 
	t  += q[2]; 
	t  += q[1]; 
	t  += q[0];
	t2  = q[0]-t; 
	t2 += q[1];
	t2 += q[2];
	t2 += q[3];
	t2 += q[4];
	mm  = (int) t; t -= (double) mm;
	t1  = t+t2;
	t2 += t - t1;
	if(t1>0.5||(t1==0.5&&t2>zero)) {
	    mm += 1;
	    t   = t1-1.0;
	    t1  = t+t2;
	    t2 += t-t1;
	}

    /* may need to re-calculate t1 and t2 */
	if( ilogb(t1)-ilogb(q[3]) < 20 )  {
	    t2  = q[1];
	    t   = q[0] - (double) mm;
	    t1  = t2+t;
	    t2  = q[2]-((t1-t)-t2);
	    t   = t1;
	    t1  = t2+t;
	    t2  = q[3]-((t1-t)-t2);
	    t   = t1;
	    t1  = t2+t;
	    t2  = q[4]-((t1-t)-t2);
	    z  *= twon24;
	    fv[5] = ((double)pr[k+5])*z;
	    if (x3==zero) {
	    	if (x2==zero) 	q[5] = x1*fv[5];
	        else 		q[5] = x1*fv[5]+x2*fv[4];
	    } else   		q[5] = x1*fv[5]+x2*fv[4]+x3*fv[3];
	    if( ilogb(t1)-ilogb(q[4]) >= 20 ) {
		t2 += q[5]; goto final_adjust;
	    } else {
	    	t   = t1;
	    	t1  = t2+t;
	    	t2  = q[5]-((t1-t)-t2);
	    }
	    z  *= twon24;
	    fv[6] = ((double)pr[k+6])*z;
	    if (x3==zero) {
	    	if (x2==zero) 	q[6] = x1*fv[6];
	        else 		q[6] = x1*fv[6]+x2*fv[5];
	    } else   		q[6] = x1*fv[6]+x2*fv[5]+x3*fv[4];
	    if( ilogb(t1)-ilogb(q[5]) >= 20 ) { 
		t2 += q[6]; goto final_adjust;
	    } else {
	    	t   = t1;
	    	t1  = t2+t;
	    	t2  = q[6]-((t1-t)-t2);
	    }
	    z  *= twon24;
	    fv[7] = ((double)pr[k+7])*z;
	    if (x3==zero) {
	    	if (x2==zero) 	q[7] = x1*fv[7];
	        else 		q[7] = x1*fv[7]+x2*fv[6];
	    } else   		q[7] = x1*fv[7]+x2*fv[6]+x3*fv[5];
	    t2 += q[7]; 
	final_adjust:
	    t   = t1;
	    t1  = t2+t;
	    t2  = t2 - (t1-t);
	}

    /* break t1+t2 into 3 pieces t1+t2+t3 with t1 and t2 24-bit */
	z  	 = y  = t1;
	t  	 = t2;
	py[n1]  &= 0xffffffe0;
	t3 	 = y;
	py[n1]  &= 0xe0000000;
	t1 	 = y;
	t2	 = t3 - y ;
	t3 	 = t - (t3 - z);


    /* compute (p1+p2+p3)*(t1+t2+t3) to y1 and y2 */
			/*  p1*t1 p1*t2 p1*t3 */
			/*        p2*t1 p2*t2 */
			/*		p3*t1 */
	q3  = p2*t2+p3*t1;
	q3 += p1*t3;
	q2  = p1*t2+p2*t1;
	q1  = p1*t1;
	t1  = q1+q2;
	t2  = q3 - ((t1-q1)-q2);
	*y1 = t1+t2;
	*y2 = t2 - (*y1 - t1);
	if (signx==1)  {*y1= -(*y1); *y2= -(*y2); mm= -mm;}
	return mm&7;
}