backend: add double version of cos

cp from fdlibm and pass the cft after refined

Signed-off-by: rander <rander.wang@intel.com>
Tested-by: Yang Rong <rong.r.yang@intel.com>

2 files changed, 409 insertions, 0 deletions

diff --git a/backend/src/libocl/tmpl/ocl_math_common.tmpl.cl b/backend/src/libocl/tmpl/ocl_math_common.tmpl.cl
index 819dc269..4f0fb90d 100644
--- a/backend/src/libocl/tmpl/ocl_math_common.tmpl.cl
+++ b/backend/src/libocl/tmpl/ocl_math_common.tmpl.cl
@@ -463,6 +463,414 @@ OVERLOADABLE double copysign(double x, double y)
 	return as_double(ux);
 }
 
+double __scalbn (double x, int n)
+{
+	double two54   =  1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
+	twom54  =  5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */
+	huge   = 1.0e+300,
+	tiny   = 1.0e-300;
+
+	int  k,hx,lx;
+	hx = __HI(x);
+	lx = __LO(x);
+		k = (hx&0x7ff00000)>>20;		/* extract exponent */
+		if (k==0) {				/* 0 or subnormal x */
+			if ((lx|(hx&0x7fffffff))==0) return x; /* +-0 */
+		x *= two54;
+		hx = __HI(x);
+		k = ((hx&0x7ff00000)>>20) - 54;
+			if (n< -50000) return tiny*x; 	/*underflow*/
+		}
+		if (k==0x7ff) return x+x;		/* NaN or Inf */
+		k = k+n;
+		if (k >  0x7fe) return huge*copysign(huge,x); /* overflow  */
+		if (k > 0) 				/* normal result */
+		{__setHigh(&x, (hx&0x800fffff)|(k<<20)); return x;}
+		if (k <= -54)
+			{
+			if (n > 50000) 	/* in case integer overflow in n+k */
+		return huge*copysign(huge,x);	/*overflow*/
+		else return tiny*copysign(tiny,x); 	/*underflow*/
+		}
+		k += 54;				/* subnormal result */
+		__setHigh(&x, (hx&0x800fffff)|(k<<20));
+		return x*twom54;
+}
+
+int __kernel_rem_pio2(double *x, double *y, int e0, int nx, int prec, const int *ipio2)
+{
+	double zero   = 0.0,
+	one	= 1.0,
+	two24   =  1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
+	twon24  =  5.96046447753906250000e-08; /* 0x3E700000, 0x00000000 */
+
+	double PIo2[] = {
+	  1.57079625129699707031e+00, /* 0x3FF921FB, 0x40000000 */
+	  7.54978941586159635335e-08, /* 0x3E74442D, 0x00000000 */
+	  5.39030252995776476554e-15, /* 0x3CF84698, 0x80000000 */
+	  3.28200341580791294123e-22, /* 0x3B78CC51, 0x60000000 */
+	  1.27065575308067607349e-29, /* 0x39F01B83, 0x80000000 */
+	  1.22933308981111328932e-36, /* 0x387A2520, 0x40000000 */
+	  2.73370053816464559624e-44, /* 0x36E38222, 0x80000000 */
+	  2.16741683877804819444e-51, /* 0x3569F31D, 0x00000000 */
+	};
+
+	const int init_jk[] = {2,3,4,6}; /* initial value for jk */
+
+	int jz,jx,jv,jp,jk,carry,n,iq[20],i,j,k,m,q0,ih;
+	double z,fw,f[20],fq[20],q[20];
+
+	/* initialize jk*/
+	jk = init_jk[prec];
+	jp = jk;
+
+	/* determine jx,jv,q0, note that 3>q0 */
+	jx =  nx-1;
+	jv = (e0-3)/24; if(jv<0) jv=0;
+	q0 =  e0-24*(jv+1);
+
+	/* set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk] */
+	j = jv-jx; m = jx+jk;
+	for(i=0;i<=m;i++,j++) f[i] = (j<0)? zero : (double) ipio2[j];
+
+	/* compute q[0],q[1],...q[jk] */
+	for (i=0;i<=jk;i++) {
+		for(j=0,fw=0.0;j<=jx;j++) fw += x[j]*f[jx+i-j]; q[i] = fw;
+	}
+
+	jz = jk;
+recompute:
+	/* distill q[] into iq[] reversingly */
+	for(i=0,j=jz,z=q[jz];j>0;i++,j--) {
+		fw	=  (double)((int)(twon24* z));
+		iq[i] =  (int)(z-two24*fw);
+		z	 =  q[j-1]+fw;
+	}
+
+	/* compute n */
+	z  = __scalbn(z,q0);		/* actual value of z */
+	z -= 8.0*floor(z*0.125);		/* trim off integer >= 8 */
+	n  = (int) z;
+	z -= (double)n;
+	ih = 0;
+	if(q0>0) {	/* need iq[jz-1] to determine n */
+		i  = (iq[jz-1]>>(24-q0)); n += i;
+		iq[jz-1] -= i<<(24-q0);
+		ih = iq[jz-1]>>(23-q0);
+	}
+	else if(q0==0) ih = iq[jz-1]>>23;
+	else if(z>=0.5) ih=2;
+
+	if(ih>0) {	/* q > 0.5 */
+		n += 1; carry = 0;
+		for(i=0;i<jz ;i++) {	/* compute 1-q */
+		j = iq[i];
+		if(carry==0) {
+			if(j!=0) {
+			carry = 1; iq[i] = 0x1000000- j;
+			}
+		} else  iq[i] = 0xffffff - j;
+		}
+		if(q0>0) {		/* rare case: chance is 1 in 12 */
+			switch(q0) {
+			case 1:
+			   iq[jz-1] &= 0x7fffff; break;
+			case 2:
+			   iq[jz-1] &= 0x3fffff; break;
+			}
+		}
+		if(ih==2) {
+		z = one - z;
+		if(carry!=0) z -= __scalbn(one,q0);
+		}
+	}
+
+	/* check if recomputation is needed */
+	if(z==zero) {
+		j = 0;
+		for (i=jz-1;i>=jk;i--) j |= iq[i];
+		if(j==0) { /* need recomputation */
+		for(k=1;iq[jk-k]==0;k++);   /* k = no. of terms needed */
+
+		for(i=jz+1;i<=jz+k;i++) {   /* add q[jz+1] to q[jz+k] */
+			f[jx+i] = (double) ipio2[jv+i];
+			for(j=0,fw=0.0;j<=jx;j++) fw += x[j]*f[jx+i-j];
+			q[i] = fw;
+		}
+		jz += k;
+		goto recompute;
+		}
+	}
+
+	/* chop off zero terms */
+	if(z==0.0) {
+		jz -= 1; q0 -= 24;
+		while(iq[jz]==0) { jz--; q0-=24;}
+	} else { /* break z into 24-bit if necessary */
+		z = __scalbn(z,-q0);
+		if(z>=two24) {
+		fw = (double)((int)(twon24*z));
+		iq[jz] = (int)(z-two24*fw);
+		jz += 1; q0 += 24;
+		iq[jz] = (int) fw;
+		} else iq[jz] = (int) z ;
+	}
+
+	/* convert integer "bit" chunk to floating-point value */
+	fw = __scalbn(one,q0);
+	for(i=jz;i>=0;i--) {
+		q[i] = fw*(double)iq[i]; fw*=twon24;
+	}
+
+	/* compute PIo2[0,...,jp]*q[jz,...,0] */
+	for(i=jz;i>=0;i--) {
+		for(fw=0.0,k=0;k<=jp&&k<=jz-i;k++) fw += PIo2[k]*q[i+k];
+		fq[jz-i] = fw;
+	}
+
+	/* compress fq[] into y[] */
+	switch(prec) {
+		case 0:
+		fw = 0.0;
+		for (i=jz;i>=0;i--) fw += fq[i];
+		y[0] = (ih==0)? fw: -fw;
+		break;
+		case 1:
+		case 2:
+		fw = 0.0;
+		for (i=jz;i>=0;i--) fw += fq[i];
+		y[0] = (ih==0)? fw: -fw;
+		fw = fq[0]-fw;
+		for (i=1;i<=jz;i++) fw += fq[i];
+		y[1] = (ih==0)? fw: -fw;
+		break;
+		case 3:	/* painful */
+		for (i=jz;i>0;i--) {
+			fw	  = fq[i-1]+fq[i];
+			fq[i]  += fq[i-1]-fw;
+			fq[i-1] = fw;
+		}
+		for (i=jz;i>1;i--) {
+			fw	  = fq[i-1]+fq[i];
+			fq[i]  += fq[i-1]-fw;
+			fq[i-1] = fw;
+		}
+		for (fw=0.0,i=jz;i>=2;i--) fw += fq[i];
+		if(ih==0) {
+			y[0] =  fq[0]; y[1] =  fq[1]; y[2] =  fw;
+		} else {
+			y[0] = -fq[0]; y[1] = -fq[1]; y[2] = -fw;
+		}
+	}
+	return n&7;
+}
+
+int __ieee754_rem_pio2(double x, double *y)
+{
+	int two_over_pi[] = {
+	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,
+	};
+
+	int npio2_hw[] = {
+	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,
+	};
+
+	double zero =  0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
+	halfD =  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 */
+
+	double z,w,t,r,fn;
+	double tx[3];
+	int e0,i,j,nx,n,ix,hx;
+
+	hx = __HI(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  = (int) (t*invpio2+halfD);
+		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 {
+			j  = ix>>20;
+			y[0] = r-w;
+			i = j-(((__HI(y[0]))>>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;
+			i = j-(((__HI(y[0]))>>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) */
+	__setLow(&z,  __LO(x));
+	e0 	= (ix>>20)-1046;	/* e0 = ilogb(z)-23; */
+	__setHigh(&z, ix - (e0<<20));
+	for(i=0;i<2;i++) {
+		tx[i] = (double)((int)(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;
+}
+
+double __kernel_cos(double x, double y)
+{
+	double one =  1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
+	C1  =  4.16666666666666019037e-02, /* 0x3FA55555, 0x5555554C */
+	C2  = -1.38888888888741095749e-03, /* 0xBF56C16C, 0x16C15177 */
+	C3  =  2.48015872894767294178e-05, /* 0x3EFA01A0, 0x19CB1590 */
+	C4  = -2.75573143513906633035e-07, /* 0xBE927E4F, 0x809C52AD */
+	C5  =  2.08757232129817482790e-09, /* 0x3E21EE9E, 0xBDB4B1C4 */
+	C6  = -1.13596475577881948265e-11; /* 0xBDA8FAE9, 0xBE8838D4 */
+
+	double a,hz,z,r,qx;
+	int ix;
+	ix = __HI(x)&0x7fffffff;	/* ix = |x|'s high word*/
+	if(ix<0x3e400000) {			/* if x < 2**27 */
+		if(((int)x)==0) return one;		/* generate inexact */
+	}
+	z  = x*x;
+	r  = z*(C1+z*(C2+z*(C3+z*(C4+z*(C5+z*C6)))));
+	if(ix < 0x3FD33333) 			/* if |x| < 0.3 */
+		return one - (0.5*z - (z*r - x*y));
+	else {
+		if(ix > 0x3fe90000) {		/* x > 0.78125 */
+		qx = 0.28125;
+		} else {
+			__setHigh(&qx,  ix-0x00200000);	/* x/4 */
+			__setLow(&qx, 0);
+		}
+		hz = 0.5*z-qx;
+		a  = one-qx;
+		return a - (hz - (z*r-x*y));
+	}
+}
+
+double __kernel_sin(double x, double y, int iy)
+{
+	double halfD =  5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
+	S1  = -1.66666666666666324348e-01, /* 0xBFC55555, 0x55555549 */
+	S2  =  8.33333333332248946124e-03, /* 0x3F811111, 0x1110F8A6 */
+	S3  = -1.98412698298579493134e-04, /* 0xBF2A01A0, 0x19C161D5 */
+	S4  =  2.75573137070700676789e-06, /* 0x3EC71DE3, 0x57B1FE7D */
+	S5  = -2.50507602534068634195e-08, /* 0xBE5AE5E6, 0x8A2B9CEB */
+	S6  =  1.58969099521155010221e-10; /* 0x3DE5D93A, 0x5ACFD57C */
+
+	double z,r,v;
+	int ix;
+	ix = __HI(x)&0x7fffffff;	/* high word of x */
+	if(ix<0x3e400000)			/* |x| < 2**-27 */
+	   {if((int)x==0) return x;}		/* generate inexact */
+	z	=  x*x;
+	v	=  z*x;
+	r	=  S2+z*(S3+z*(S4+z*(S5+z*S6)));
+	if(iy==0) return x+v*(S1+z*r);
+	else	  return x-((z*(halfD*y-v*r)-y)-v*S1);
+}
+
+OVERLOADABLE double cos(double x)
+{
+	double y[2],z=0.0;
+	int n, ix;
+
+	/* High word of x. */
+	ix = __HI(x);
+
+	/* |x| ~< pi/4 */
+	ix &= 0x7fffffff;
+	if(ix <= 0x3fe921fb) return __kernel_cos(x,z);
+
+	/* cos(Inf or NaN) is NaN */
+	else if (ix>=0x7ff00000) return x-x;
+
+	/* argument reduction needed */
+	else {
+		n = __ieee754_rem_pio2(x,y);
+		switch(n&3) {
+		case 0: return  __kernel_cos(y[0],y[1]);
+		case 1: return -__kernel_sin(y[0],y[1],1);
+		case 2: return -__kernel_cos(y[0],y[1]);
+		default:
+				return  __kernel_sin(y[0],y[1],1);
+		}
+	}
+}
+
 OVERLOADABLE double fabs(double x)
 {
     long  qw = as_ulong(x);
diff --git a/backend/src/libocl/tmpl/ocl_math_common.tmpl.h b/backend/src/libocl/tmpl/ocl_math_common.tmpl.h
index 76ce952e..1802533f 100644
--- a/backend/src/libocl/tmpl/ocl_math_common.tmpl.h
+++ b/backend/src/libocl/tmpl/ocl_math_common.tmpl.h
@@ -33,6 +33,7 @@ OVERLOADABLE double atan2pi(double x, double y);
 OVERLOADABLE double atanh(double x);
 OVERLOADABLE double ceil(double x);
 OVERLOADABLE double copysign(double x, double y);
+OVERLOADABLE double cos(double x);
 OVERLOADABLE double fabs(double x);
 OVERLOADABLE double fdim(double x, double y);
 OVERLOADABLE double floor(double x);