[picodrive.git] / platform / gp2x / code940 / uClibc / e_rem_pio2.c

/* @(#)e_rem_pio2.c 5.1 93/09/24 */
/*
 * ====================================================
 * 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_pio2.c,v 1.8 1995/05/10 20:46:02 jtc Exp $";
#endif

/* __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
 */
#ifdef __STDC__
static const int32_t two_over_pi[] = {
#else
static int32_t two_over_pi[] = {
#endif
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,
};

#ifdef __STDC__
static const int32_t npio2_hw[] = {
#else
static int32_t npio2_hw[] = {
#endif
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)
 */

#ifdef __STDC__
static const double
#else
static double
#endif
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 */

#ifdef __STDC__
	int32_t __ieee754_rem_pio2(double x, double *y)
#else
	int32_t __ieee754_rem_pio2(x,y)
	double x,y[];
#endif
{
	double z=0.0,w,t,r,fn;
	double tx[3];
	int32_t e0,i,j,nx,n,ix,hx;
	u_int32_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 {
	        u_int32_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;
}
Commit	Line	Data
	1	/* @(#)e_rem_pio2.c 5.1 93/09/24 */
	2	/*
	3	* ====================================================
	4	* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
	5	*
	6	* Developed at SunPro, a Sun Microsystems, Inc. business.
	7	* Permission to use, copy, modify, and distribute this
	8	* software is freely granted, provided that this notice
	9	* is preserved.
	10	* ====================================================
	11	*/
	12
	13	#if defined(LIBM_SCCS) && !defined(lint)
	14	static char rcsid[] = "$NetBSD: e_rem_pio2.c,v 1.8 1995/05/10 20:46:02 jtc Exp $";
	15	#endif
	16
	17	/* __ieee754_rem_pio2(x,y)
	18	*
	19	* return the remainder of x rem pi/2 in y[0]+y[1]
	20	* use __kernel_rem_pio2()
	21	*/
	22
	23	#include "math.h"
	24	#include "math_private.h"
	25
	26	/*
	27	* Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi
	28	*/
	29	#ifdef __STDC__
	30	static const int32_t two_over_pi[] = {
	31	#else
	32	static int32_t two_over_pi[] = {
	33	#endif
	34	0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62,
	35	0x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A,
	36	0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129,
	37	0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41,
	38	0x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8,
	39	0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF,
	40	0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5,
	41	0xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08,
	42	0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3,
	43	0x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880,
	44	0x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B,
	45	};
	46
	47	#ifdef __STDC__
	48	static const int32_t npio2_hw[] = {
	49	#else
	50	static int32_t npio2_hw[] = {
	51	#endif
	52	0x3FF921FB, 0x400921FB, 0x4012D97C, 0x401921FB, 0x401F6A7A, 0x4022D97C,
	53	0x4025FDBB, 0x402921FB, 0x402C463A, 0x402F6A7A, 0x4031475C, 0x4032D97C,
	54	0x40346B9C, 0x4035FDBB, 0x40378FDB, 0x403921FB, 0x403AB41B, 0x403C463A,
	55	0x403DD85A, 0x403F6A7A, 0x40407E4C, 0x4041475C, 0x4042106C, 0x4042D97C,
	56	0x4043A28C, 0x40446B9C, 0x404534AC, 0x4045FDBB, 0x4046C6CB, 0x40478FDB,
	57	0x404858EB, 0x404921FB,
	58	};
	59
	60	/*
	61	* invpio2: 53 bits of 2/pi
	62	* pio2_1: first 33 bit of pi/2
	63	* pio2_1t: pi/2 - pio2_1
	64	* pio2_2: second 33 bit of pi/2
	65	* pio2_2t: pi/2 - (pio2_1+pio2_2)
	66	* pio2_3: third 33 bit of pi/2
	67	* pio2_3t: pi/2 - (pio2_1+pio2_2+pio2_3)
	68	*/
	69
	70	#ifdef __STDC__
	71	static const double
	72	#else
	73	static double
	74	#endif
	75	zero = 0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
	76	half = 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
	77	two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
	78	invpio2 = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
	79	pio2_1 = 1.57079632673412561417e+00, /* 0x3FF921FB, 0x54400000 */
	80	pio2_1t = 6.07710050650619224932e-11, /* 0x3DD0B461, 0x1A626331 */
	81	pio2_2 = 6.07710050630396597660e-11, /* 0x3DD0B461, 0x1A600000 */
	82	pio2_2t = 2.02226624879595063154e-21, /* 0x3BA3198A, 0x2E037073 */
	83	pio2_3 = 2.02226624871116645580e-21, /* 0x3BA3198A, 0x2E000000 */
	84	pio2_3t = 8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */
	85
	86	#ifdef __STDC__
	87	int32_t __ieee754_rem_pio2(double x, double *y)
	88	#else
	89	int32_t __ieee754_rem_pio2(x,y)
	90	double x,y[];
	91	#endif
	92	{
	93	double z=0.0,w,t,r,fn;
	94	double tx[3];
	95	int32_t e0,i,j,nx,n,ix,hx;
	96	u_int32_t low;
	97
	98	GET_HIGH_WORD(hx,x); /* high word of x */
	99	ix = hx&0x7fffffff;
	100	if(ix<=0x3fe921fb) /* \|x\| ~<= pi/4 , no need for reduction */
	101	{y[0] = x; y[1] = 0; return 0;}
	102	if(ix<0x4002d97c) { /* \|x\| < 3pi/4, special case with n=+-1 */
	103	if(hx>0) {
	104	z = x - pio2_1;
	105	if(ix!=0x3ff921fb) { /* 33+53 bit pi is good enough */
	106	y[0] = z - pio2_1t;
	107	y[1] = (z-y[0])-pio2_1t;
	108	} else { /* near pi/2, use 33+33+53 bit pi */
	109	z -= pio2_2;
	110	y[0] = z - pio2_2t;
	111	y[1] = (z-y[0])-pio2_2t;
	112	}
	113	return 1;
	114	} else { /* negative x */
	115	z = x + pio2_1;
	116	if(ix!=0x3ff921fb) { /* 33+53 bit pi is good enough */
	117	y[0] = z + pio2_1t;
	118	y[1] = (z-y[0])+pio2_1t;
	119	} else { /* near pi/2, use 33+33+53 bit pi */
	120	z += pio2_2;
	121	y[0] = z + pio2_2t;
	122	y[1] = (z-y[0])+pio2_2t;
	123	}
	124	return -1;
	125	}
	126	}
	127	if(ix<=0x413921fb) { /* \|x\| ~<= 2^19(pi/2), medium size /
	128	t = fabs(x);
	129	n = (int32_t) (t*invpio2+half);
	130	fn = (double)n;
	131	r = t-fn*pio2_1;
	132	w = fnpio2_1t; / 1st round good to 85 bit */
	133	if(n<32&&ix!=npio2_hw[n-1]) {
	134	y[0] = r-w; /* quick check no cancellation */
	135	} else {
	136	u_int32_t high;
	137	j = ix>>20;
	138	y[0] = r-w;
	139	GET_HIGH_WORD(high,y[0]);
	140	i = j-((high>>20)&0x7ff);
	141	if(i>16) { /* 2nd iteration needed, good to 118 */
	142	t = r;
	143	w = fn*pio2_2;
	144	r = t-w;
	145	w = fn*pio2_2t-((t-r)-w);
	146	y[0] = r-w;
	147	GET_HIGH_WORD(high,y[0]);
	148	i = j-((high>>20)&0x7ff);
	149	if(i>49) { /* 3rd iteration need, 151 bits acc */
	150	t = r; /* will cover all possible cases */
	151	w = fn*pio2_3;
	152	r = t-w;
	153	w = fn*pio2_3t-((t-r)-w);
	154	y[0] = r-w;
	155	}
	156	}
	157	}
	158	y[1] = (r-y[0])-w;
	159	if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;}
	160	else return n;
	161	}
	162	/*
	163	* all other (large) arguments
	164	*/
	165	if(ix>=0x7ff00000) { /* x is inf or NaN */
	166	y[0]=y[1]=x-x; return 0;
	167	}
	168	/* set z = scalbn(\|x\|,ilogb(x)-23) */
	169	GET_LOW_WORD(low,x);
	170	SET_LOW_WORD(z,low);
	171	e0 = (ix>>20)-1046; /* e0 = ilogb(z)-23; */
	172	SET_HIGH_WORD(z, ix - ((int32_t)(e0<<20)));
	173	for(i=0;i<2;i++) {
	174	tx[i] = (double)((int32_t)(z));
	175	z = (z-tx[i])*two24;
	176	}
	177	tx[2] = z;
	178	nx = 3;
	179	while(tx[nx-1]==zero) nx--; /* skip zero term */
	180	n = __kernel_rem_pio2(tx,y,e0,nx,2,two_over_pi);
	181	if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;}
	182	return n;
	183	}