trig.c

操作系统SunOS 4.1.3版本的源码

#ifndef lint
static  char sccsid[] = "@(#)trig.c 1.1 92/07/30 SMI"; 
#endif

/*
 * Copyright (c) 1988 by Sun Microsystems, Inc.
 */

/* sin(x), cos(x), tan(x), sincos(x,s,c)
 * Coefficients obtained from Peter Tang's "Some Software 
 * Implementations of the Functions Sin and Cos." 
 * Code by K.C. Ng for SUN 4.0, Feb 3, 1987. Revised on April, 1988.
 *
 * Static kernel function:
 *	trig()		... sin, cos, tan, and sincos
 *      argred()        ... huge argument reduction
 *
 * Method.
 *      Let S and C denote the polynomial approximations to sin and cos 
 *      respectively on [-PI/4, +PI/4].
 *      1. Assume the argument x is reduced to y1+y2 = x-k*pi/2 in
 *	   [-pi/2 , +pi/2], and let n = k mod 4.
 *	2. Let S=S(y1+y2), C=C(y1+y2). Depending on n, we have
 *
 *          n        sin(x)      cos(x)        tan(x)
 *     ----------------------------------------------------------
 *	    0	       S	   C		 S/C
 *	    1	       C	  -S		-C/S
 *	    2	      -S	  -C		 S/C
 *	    3	      -C	   S		-C/S
 *     ----------------------------------------------------------
 *
 *   Notes:
 *      1. S = S(y1+y2) is computed by:
 *		y1 + ((y2-0.5*z*y2) - (y1*z)*ss)
 *         where z=y1*y1 and ss is an approximation of (y1-sin(y1))/y1^3.
 *
 *	2. C = C(y1+y2) is computed by:
 *		z  = y1*y1
 *		hz = z*0.5
 *		t  = y1*y2
 *		if(z >= thresh1) 	c = 5/8   - ((hz-3/8 )-(cc-t));
 *		else if (z >= thresh2)  c = 13/16 - ((hz-3/16)-(cc-t));
 *		else 			c = 1     - ( hz      -(cc-t));
 *         where
 *	        cc is an approximation of cos(y1)-1+y1*y1/2 
 *              thresh1 = (acos(3/4)**2)  = 2^-1 * Hex 1.0B70C6D604DD5 
 *              thresh2 = (acos(7/8)**2)  = 2^-2 * Hex 1.0584C22231902 
 *
 *	3. Since S(y1+y2)/C(y1+y2) = tan(y1+y2) ~ tan(y1) + y2 + y1*y1*y2, we 
 *	   have
 *	      S(y1+y2)/C(y1+y2) = S(y1)/C(y1) + y2 + y1*y1*y2 
 *				= [S(y1)+C(y1)*(y2+y1*y1*y2)]/C(y1)
 *	   and hence 
 *	      C(y1+y2)/S(y1+y2) = C(y1)/[S(y1)+C(y1)*(y2+y1*y1*y2)]
 *
 *         For better accuracy, S/C and -C/S are computed as follow. 
 *	   Using the same notation as 1 and 2, 
 *		z  = y1*y1
 *		hz = ysq*0.5
 *		if(z >= thresh1) 	c = 5/8   - ((hz-3/8 )-cc);
 *		else if (z >= thresh2)  c = 13/16 - ((hz-3/16)-cc);
 *		else 			c = 1     - ( hz      -cc);
 *		
 *				     hz - (cc+ss)
 *		 S/C  =  y1 + ( y1*(--------------- +y1*y2) + y2)
 *					 c 
 *
 *		-C/S  = -c/(y1+(y1*ss+c*(y2+z*y2)))
 *
 * Special cases:
 *      Let trig be any of sin, cos, or tan.
 *      trig(+-INF)  is NaN, with signals;
 *      trig(NaN)    is that NaN;
 *
 * Accuracy:
 *	According to fp_pi, computer TRIG(x) returns 
 *		trig(x)    	nearly rounded if fp_pi == 0 =fp_pi_infinte;
 *		trig(x*pi/PI66)	nearly rounded if fp_pi == 1 == fp_pi_66;
 *		trig(x*pi/PI53)	nearly rounded if fp_pi == 2 == fp_pi_53;
 *
 *      In decimal:
 *			pi   = 3.141592653589793 23846264338327 ..... 
 *    		66 bits PI66 = 3.141592653589793 23845859885078 ..... ,
 *    		53 bits PI53 = 3.141592653589793 115997963 ..... ,
 *
 *      In hexadecimal:
 *              	pi   = 3.243F6A8885A308D313198A2E....
 *    		66 bits PI66 = 3.243F6A8885A308D3 =  2 * 1.921FB54442D184698
 *    		53 bits PI53 = 3.243F6A8885A30    =  2 * 1.921FB54442D18    
 *
 * 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.
 */

#include <math.h>

enum fp_pi_type fp_pi = fp_pi_66; 	/* initialize to 66 bit PI */
static int argred(),dummy();
static double trig(),sincos_c;

double sin(x)
double x;
{
	return trig(x,0);
}

double cos(x)
double x;
{
	return trig(x,1);
}

double tan(x)
double x;
{
	return trig(x,2);
}

void sincos(x,s,c)
double x,*s,*c;
{
	*s = trig(x,3);
	*c = sincos_c;
}


/*
 * thresh1:  (acos(3/4))**2
 * thresh1:  (acos(7/8))**2
 * invpio2:  53 bits of 2/pi
 * pio4:     53 bits of pi/4
 * pio2:     53 bits of pi/2
 * pio2_1:   first  33 bit of pi/2
 * pio2_1t:  pi/2 - pio2_1
 * pio2_1t5: (53 bit pi)/2 - pio2_1
 * pio2_2:   second 33 bit of pi/2
 * pio2_2t:  pi/2 - pio2_2
 * pio2_3:   third  33 bit of pi/2
 * pio2_3t:  pi/2 - pio2_3
 */
static double   /* see above for the definition of the values */
thresh1 = 0.52234479296242375401249,    /* 2^ -1  * Hex  1.0B70C6D604DD5 */
thresh2 = 0.25538924535466389631466,    /* 2^ -2  * Hex  1.0584C22231902 */
invpio2 = 0.636619772367581343075535,   /* 2^ -1  * Hex  1.45F306DC9C883 */
pio4    = 0.78539816339744830961566,    /* 2^ -1  * Hex  1.921FB54442D18 */
pio2    = 1.57079632679489661923,       /* 2^  0  * Hex  1.921FB54442D18 */
pio2_1  = 1.570796326734125614166,      /* 2^  0  * Hex  1.921FB54400000 */
pio2_1t = 6.077100506506192601475e-11,  /* 2^-34  * Hex  1.0B4611A626331 */
pio2_1t5= 6.077094383272196864709e-11,  /* 2^-34  * Hex  1.0B46000000000 */
pio2_2  = 6.077100506303965976596e-11,  /* 2^-34  * Hex  1.0B4611A600000 */
pio2_2t = 2.022266248795950732400e-21,  /* 2^-69  * Hex  1.3198A2E037073 */
pio2_3  = 2.022266248711166455796e-21,  /* 2^-69  * Hex  1.3198A2E000000 */
pio2_3t = 8.478427660368899643959e-32;  /* 2^-104 * Hex  1.B839A252049C1 */

static double C[] = {
	/* C[n] = coefficients for cos(x)-1+x*x/2 with pi, n=0,...,5 */
 4.16666666666666019E-2    , /*C0 2^ -5 * Hex  1.555555555554C */
-1.38888888888744744E-3    , /*C1 2^-10 * Hex -1.6C16C16C1521F */
 2.48015872896717822E-5    , /*C2 2^-16 * Hex  1.A01A019CBF671 */
-2.75573144009911252E-7    , /*C3 2^-22 * Hex -1.27E4F812B495B */
 2.08757292566166689E-9    , /*C4 2^-29 * Hex  1.1EE9F14CDC590 */
-1.13599319556004135E-11   , /*C5 2^-37 * Hex -1.8FB12BAF59D4B */
 0.0,
 0.0,
	/* C[8+n] = coefficients for cos(x)-1+x*x/2 with PI66, n=0,...,5 */
 4.16666666666666019E-2    , /*C0 2^ -5 * Hex  1.555555555554C */
-1.38888888888744744E-3    , /*C1 2^-10 * Hex -1.6C16C16C1521F */
 2.48015872896717822E-5    , /*C2 2^-16 * Hex  1.A01A019CBF671 */
-2.75573144009911252E-7    , /*C3 2^-22 * Hex -1.27E4F812B495B */
 2.08757292566166689E-9    , /*C4 2^-29 * Hex  1.1EE9F14CDC590 */
-1.13599319556004135E-11   , /*C5 2^-37 * Hex -1.8FB12BAF59D4B */
 0.0,
 0.0,
	/* C[16+n] = coefficients for cos(x)-1+x*x/2 with PI53, n=0,...,5 */
 4.1666666666666504759E-2    , /*C0 2^ -5 * Hex  1.555555555553E */
-1.3888888888865301516E-3    , /*C1 2^-10 * Hex -1.6C16C16C14199 */
 2.4801587269650015769E-5    , /*C2 2^-16 * Hex  1.A01A01971CAEB */
-2.7557304623183959811E-7    , /*C3 2^-22 * Hex -1.27E4F1314AD1A */
 2.0873958177697780076E-9    , /*C4 2^-29 * Hex  1.1EE3B60DDDC8C */
-1.1250289076471311557E-11   , /*C5 2^-37 * Hex -1.8BD5986B2A52E */
 0.0,
 0.0,
};

static double S[] = {
	/* S[n] = coefficients for (x-sin(x))/x with pi, n=0,...,5 */
 1.66666666666666796E-1    , /* S0 2^ -3 * Hex  1.555555555555A */
-8.33333333333178931E-3    , /* S1 2^ -7 * Hex -1.1111111110D97 */
 1.98412698361250482E-4    , /* S2 2^-13 * Hex  1.A01A019E4A9AC */
-2.75573156035895542E-6    , /* S3 2^-19 * Hex -1.71DE37262ECF8 */
 2.50510254394993115E-8    , /* S4 2^-26 * Hex  1.AE5F933569B98 */
-1.59108690260756780E-10   , /* S5 2^-33 * Hex -1.5DE23AD2495F2 */
 0.0,
 0.0,
	/* S[8+n] = coefficients for (x-sin(x))/x with PI66, n=0,...,5 */
 1.66666666666666796E-1    , /* S0 2^ -3 * Hex  1.555555555555A */
-8.33333333333178931E-3    , /* S1 2^ -7 * Hex -1.1111111110D97 */
 1.98412698361250482E-4    , /* S2 2^-13 * Hex  1.A01A019E4A9AC */
-2.75573156035895542E-6    , /* S3 2^-19 * Hex -1.71DE37262ECF8 */
 2.50510254394993115E-8    , /* S4 2^-26 * Hex  1.AE5F933569B98 */
-1.59108690260756780E-10   , /* S5 2^-33 * Hex -1.5DE23AD2495F2 */
 0.0,
 0.0,
	/* S[16+n] = coefficients for (x-sin(x))/x with PI53, n=0,...,5 */
 1.6666666666666463126E-1    , /* S0 2^ -3 * Hex  1.555555555550C */
-8.3333333332992771264E-3    , /* S1 2^ -7 * Hex -1.111111110C461 */
 1.9841269816180999116E-4    , /* S2 2^-13 * Hex  1.A01A019746345 */
-2.7557309793219876880E-6    , /* S3 2^-19 * Hex -1.71DE3209CDCD9 */
 2.5050225177523807003E-8    , /* S4 2^-26 * Hex  1.AE5C0E319A4EF */
-1.5868926979889205164E-10   , /* S5 2^-33 * Hex -1.5CF61DF672B13 */
 0.0,
 0.0,
};

static double			
medium  = 1647099.0,    /* 2^19*(pi/2): threshold for medium arg red */
tiny   	= 1.0e-10,	/* fl(1.0 + tiny*tiny) == 1.0 */
big	= 1.0e10,	/* 1/tiny */
zero	= 0.0,
one 	= 1.0,
half	= 0.5,
c5_8	= 0.625,
c3_8	= 0.375,
c13_16	= 0.8125,
c3_16	= 0.1875;

static double trig(x,k) 
double x; int k;
{
        double y1,y2,t,z,hz,ss,w,cc,fn;
	register j,m;
        int n,signx;
	if(!finite(x))  return sincos_c=x-x;	/* trig(NaN or INF) is NaN */
	signx = signbit(x); t = fabs(x);
	if(t<0.002246) {
	    if (t<tiny) {
		dummy(t+big);	/* create inexact flag if x!=0 */
		if(k==3) {sincos_c = one; return x;}
		return (k==1)? one:x; 
	    }
	    z = x*x;
	    if (t<0.000085) {
		switch(k) {
		case 0: return x - z*x*0.1666666666666666666;
		case 1: return 1 - 0.5*z;
		case 2: return x + z*x*0.3333333333333333333;
		case 3: sincos_c = 1 - 0.5*z;
			return x - z*x*0.1666666666666666666;
		}
	    }
	    switch(k) {
	    case 0:	return 
		x-z*x*(0.1666666666666666666-0.008333333333333333333*z);
	    case 1:	return 
		1-z*(0.5-0.04166666666666666666*z);
	    case 2:	return 
		x+z*x*(0.3333333333333333333+0.133333333333333333333*z);
	    case 3:	sincos_c = 1-z*(0.5-0.04166666666666666666*z); return
		x-z*x*(0.1666666666666666666-0.008333333333333333333*z);

	    }
	}
	n  = j = 0; 	/* j will be used to indicate whether y2=0 */
	y1 = x;
	if(t > pio4) {
	    switch((int)fp_pi) {
		case 0: j = 1;
			if(t<=pio2) {
			    n  = 1;
			    z  = t - pio2_1;
			    y1 = z - pio2_1t;
			    if(fabs(y1)>0.00390625) y2=(z-y1)-pio2_1t;
			    else {
				z -= pio2_2;
				y1 = z - pio2_2t;
				if(fabs(y1)>9.094947e-13) y2=(z-y1)-pio2_2t;
				else {
				    z -= pio2_3;
				    y1 = z - pio2_3t;
				    y2 = (z-y1)-pio2_3t;
				}
			    }
			    if(signx==1) {n= -n; y1= -y1; y2= -y2;}
			} else if(t<=medium) {
			    n  = (int) (t*invpio2+half);
			    fn = n;
			    z  = t - fn*pio2_1;
			    w  = fn*pio2_1t;
			    y1 = z - w;
			    if(fabs(y1)>0.00390625) y2=(z-y1)-w;
			    else {
				z -= fn*pio2_2;
			        w  = fn*pio2_2t;
				y1 = z - w;
				if(fabs(y1)>9.094947e-13) y2=(z-y1)-w;
				else {
				    z -= fn*pio2_3;
				    w  = fn*pio2_3t;
				    y1 = z - w;
				    y2 = (z-y1)-w;
				}
			    }
			    if(signx==1) {n= -n; y1= -y1; y2= -y2;}
			} else n = argred(x,&y1,&y2);
			break;
		case 1: j = 1;
			if(t<=pio2) {
			    n  = 1;
			    z  = t - pio2_1;
			    y1 = z - pio2_2;
			    y2 = (z-y1) - pio2_2;
			    if(signx==1) {n= -n; y1= -y1; y2= -y2;}
			} else if(t<=medium) {
			    n  = (int) (t*invpio2+half);
			    z  = (t - n*pio2_1);
			    w  = n*pio2_2;
			    y1 = z-w;
			    y2 = (z-y1)-w;
			    if(signx==1) {n= -n; y1= -y1; y2= -y2;}
			} else n = argred(x,&y1,&y2);
			break;
		case 2: if(t<=pio2) {
			    n  = 1;
			    y1 = t - pio2;
			    if(signx==1) {n= -1; y1= -y1;}
			} else if(t<=medium) {
			    n  = (int) (t*invpio2+half);
			    y1 = (t - n*pio2_1)-n*pio2_1t5;
			    if(signx==1) {n= -n; y1= -y1;}
			} else {j=1; n = argred(x,&y1,&y2);}
			break;