diff options
Diffstat (limited to 'sysdeps/ieee754')
| -rw-r--r-- | sysdeps/ieee754/dbl-64/sincos32.c | 490 |
1 files changed, 254 insertions, 236 deletions
diff --git a/sysdeps/ieee754/dbl-64/sincos32.c b/sysdeps/ieee754/dbl-64/sincos32.c index 954db66d6b..ac27fcda83 100644 --- a/sysdeps/ieee754/dbl-64/sincos32.c +++ b/sysdeps/ieee754/dbl-64/sincos32.c @@ -48,312 +48,330 @@ # define SECTION #endif -/****************************************************************/ -/* Compute Multi-Precision sin() function for given p. Receive */ -/* Multi Precision number x and result stored at y */ -/****************************************************************/ +/* Compute Multi-Precision sin() function for given p. Receive Multi Precision + number x and result stored at y. */ static void SECTION -ss32(mp_no *x, mp_no *y, int p) { +ss32 (mp_no *x, mp_no *y, int p) +{ int i; double a; - mp_no mpt1,x2,gor,sum ,mpk={1,{1.0}}; - for (i=1;i<=p;i++) mpk.d[i]=0; + mp_no mpt1, x2, gor, sum, mpk = {1, {1.0}}; + for (i = 1; i <= p; i++) + mpk.d[i] = 0; - __sqr(x,&x2,p); - __cpy(&oofac27,&gor,p); - __cpy(&gor,&sum,p); - for (a=27.0;a>1.0;a-=2.0) { - mpk.d[1]=a*(a-1.0); - __mul(&gor,&mpk,&mpt1,p); - __cpy(&mpt1,&gor,p); - __mul(&x2,&sum,&mpt1,p); - __sub(&gor,&mpt1,&sum,p); - } - __mul(x,&sum,y,p); + __sqr (x, &x2, p); + __cpy (&oofac27, &gor, p); + __cpy (&gor, &sum, p); + for (a = 27.0; a > 1.0; a -= 2.0) + { + mpk.d[1] = a * (a - 1.0); + __mul (&gor, &mpk, &mpt1, p); + __cpy (&mpt1, &gor, p); + __mul (&x2, &sum, &mpt1, p); + __sub (&gor, &mpt1, &sum, p); + } + __mul (x, &sum, y, p); } -/**********************************************************************/ -/* Compute Multi-Precision cos() function for given p. Receive Multi */ -/* Precision number x and result stored at y */ -/**********************************************************************/ +/* Compute Multi-Precision cos() function for given p. Receive Multi Precision + number x and result stored at y. */ static void SECTION -cc32(mp_no *x, mp_no *y, int p) { +cc32 (mp_no *x, mp_no *y, int p) +{ int i; double a; - mp_no mpt1,x2,gor,sum ,mpk={1,{1.0}}; - for (i=1;i<=p;i++) mpk.d[i]=0; + mp_no mpt1, x2, gor, sum, mpk = {1, {1.0}}; + for (i = 1; i <= p; i++) + mpk.d[i] = 0; - __sqr(x,&x2,p); - mpk.d[1]=27.0; - __mul(&oofac27,&mpk,&gor,p); - __cpy(&gor,&sum,p); - for (a=26.0;a>2.0;a-=2.0) { - mpk.d[1]=a*(a-1.0); - __mul(&gor,&mpk,&mpt1,p); - __cpy(&mpt1,&gor,p); - __mul(&x2,&sum,&mpt1,p); - __sub(&gor,&mpt1,&sum,p); - } - __mul(&x2,&sum,y,p); + __sqr (x, &x2, p); + mpk.d[1] = 27.0; + __mul (&oofac27, &mpk, &gor, p); + __cpy (&gor, &sum, p); + for (a = 26.0; a > 2.0; a -= 2.0) + { + mpk.d[1] = a * (a - 1.0); + __mul (&gor, &mpk, &mpt1, p); + __cpy (&mpt1, &gor, p); + __mul (&x2, &sum, &mpt1, p); + __sub (&gor, &mpt1, &sum, p); + } + __mul (&x2, &sum, y, p); } -/***************************************************************************/ -/* c32() computes both sin(x), cos(x) as Multi precision numbers */ -/***************************************************************************/ +/* Compute both sin(x), cos(x) as Multi precision numbers. */ void SECTION -__c32(mp_no *x, mp_no *y, mp_no *z, int p) { - mp_no u,t,t1,t2,c,s; +__c32 (mp_no *x, mp_no *y, mp_no *z, int p) +{ + mp_no u, t, t1, t2, c, s; int i; - __cpy(x,&u,p); - u.e=u.e-1; - cc32(&u,&c,p); - ss32(&u,&s,p); - for (i=0;i<24;i++) { - __mul(&c,&s,&t,p); - __sub(&s,&t,&t1,p); - __add(&t1,&t1,&s,p); - __sub(&mptwo,&c,&t1,p); - __mul(&t1,&c,&t2,p); - __add(&t2,&t2,&c,p); - } - __sub(&mpone,&c,y,p); - __cpy(&s,z,p); + __cpy (x, &u, p); + u.e = u.e - 1; + cc32 (&u, &c, p); + ss32 (&u, &s, p); + for (i = 0; i < 24; i++) + { + __mul (&c, &s, &t, p); + __sub (&s, &t, &t1, p); + __add (&t1, &t1, &s, p); + __sub (&mptwo, &c, &t1, p); + __mul (&t1, &c, &t2, p); + __add (&t2, &t2, &c, p); + } + __sub (&mpone, &c, y, p); + __cpy (&s, z, p); } -/************************************************************************/ -/*Routine receive double x and two double results of sin(x) and return */ -/*result which is more accurate */ -/*Computing sin(x) with multi precision routine c32 */ -/************************************************************************/ +/* Receive double x and two double results of sin(x) and return result which is + more accurate, computing sin(x) with multi precision routine c32. */ double SECTION -__sin32(double x, double res, double res1) { +__sin32 (double x, double res, double res1) +{ int p; - mp_no a,b,c; - p=32; - __dbl_mp(res,&a,p); - __dbl_mp(0.5*(res1-res),&b,p); - __add(&a,&b,&c,p); - if (x>0.8) - { __sub(&hp,&c,&a,p); - __c32(&a,&b,&c,p); - } - else __c32(&c,&a,&b,p); /* b=sin(0.5*(res+res1)) */ - __dbl_mp(x,&c,p); /* c = x */ - __sub(&b,&c,&a,p); - /* if a>0 return min(res,res1), otherwise return max(res,res1) */ - if (a.d[0]>0) return (res<res1)?res:res1; - else return (res>res1)?res:res1; + mp_no a, b, c; + p = 32; + __dbl_mp (res, &a, p); + __dbl_mp (0.5 * (res1 - res), &b, p); + __add (&a, &b, &c, p); + if (x > 0.8) + { + __sub (&hp, &c, &a, p); + __c32 (&a, &b, &c, p); + } + else + __c32 (&c, &a, &b, p); /* b=sin(0.5*(res+res1)) */ + __dbl_mp (x, &c, p); /* c = x */ + __sub (&b, &c, &a, p); + /* if a > 0 return min (res, res1), otherwise return max (res, res1). */ + if (a.d[0] > 0) + return (res < res1) ? res : res1; + else + return (res > res1) ? res : res1; } -/************************************************************************/ -/*Routine receive double x and two double results of cos(x) and return */ -/*result which is more accurate */ -/*Computing cos(x) with multi precision routine c32 */ -/************************************************************************/ +/* Receive double x and two double results of cos(x) and return result which is + more accurate, computing cos(x) with multi precision routine c32. */ double SECTION -__cos32(double x, double res, double res1) { +__cos32 (double x, double res, double res1) +{ int p; - mp_no a,b,c; - p=32; - __dbl_mp(res,&a,p); - __dbl_mp(0.5*(res1-res),&b,p); - __add(&a,&b,&c,p); - if (x>2.4) - { __sub(&pi,&c,&a,p); - __c32(&a,&b,&c,p); - b.d[0]=-b.d[0]; - } - else if (x>0.8) - { __sub(&hp,&c,&a,p); - __c32(&a,&c,&b,p); - } - else __c32(&c,&b,&a,p); /* b=cos(0.5*(res+res1)) */ - __dbl_mp(x,&c,p); /* c = x */ - __sub(&b,&c,&a,p); - /* if a>0 return max(res,res1), otherwise return min(res,res1) */ - if (a.d[0]>0) return (res>res1)?res:res1; - else return (res<res1)?res:res1; + mp_no a, b, c; + p = 32; + __dbl_mp (res, &a, p); + __dbl_mp (0.5 * (res1 - res), &b, p); + __add (&a, &b, &c, p); + if (x > 2.4) + { + __sub (&pi, &c, &a, p); + __c32 (&a, &b, &c, p); + b.d[0] = -b.d[0]; + } + else if (x > 0.8) + { + __sub (&hp, &c, &a, p); + __c32 (&a, &c, &b, p); + } + else + __c32 (&c, &b, &a, p); /* b=cos(0.5*(res+res1)) */ + __dbl_mp (x, &c, p); /* c = x */ + __sub (&b, &c, &a, p); + /* if a > 0 return max (res, res1), otherwise return min (res, res1). */ + if (a.d[0] > 0) + return (res > res1) ? res : res1; + else + return (res < res1) ? res : res1; } -/*******************************************************************/ -/*Compute sin(x+dx) as Multi Precision number and return result as */ -/* double */ -/*******************************************************************/ +/* Compute sin(x+dx) as Multi Precision number and return result as double. */ double SECTION -__mpsin(double x, double dx) { +__mpsin (double x, double dx) +{ int p; double y; - mp_no a,b,c; - p=32; - __dbl_mp(x,&a,p); - __dbl_mp(dx,&b,p); - __add(&a,&b,&c,p); - if (x>0.8) { __sub(&hp,&c,&a,p); __c32(&a,&b,&c,p); } - else __c32(&c,&a,&b,p); /* b = sin(x+dx) */ - __mp_dbl(&b,&y,p); + mp_no a, b, c; + p = 32; + __dbl_mp (x, &a, p); + __dbl_mp (dx, &b, p); + __add (&a, &b, &c, p); + if (x > 0.8) + { + __sub (&hp, &c, &a, p); + __c32 (&a, &b, &c, p); + } + else + __c32 (&c, &a, &b, p); /* b = sin(x+dx) */ + __mp_dbl (&b, &y, p); return y; } -/*******************************************************************/ -/* Compute cos()of double-length number (x+dx) as Multi Precision */ -/* number and return result as double */ -/*******************************************************************/ +/* Compute cos() of double-length number (x+dx) as Multi Precision number and + return result as double. */ double SECTION -__mpcos(double x, double dx) { +__mpcos (double x, double dx) +{ int p; double y; - mp_no a,b,c; - p=32; - __dbl_mp(x,&a,p); - __dbl_mp(dx,&b,p); - __add(&a,&b,&c,p); - if (x>0.8) - { __sub(&hp,&c,&b,p); - __c32(&b,&c,&a,p); - } - else __c32(&c,&a,&b,p); /* a = cos(x+dx) */ - __mp_dbl(&a,&y,p); + mp_no a, b, c; + p = 32; + __dbl_mp (x, &a, p); + __dbl_mp (dx, &b, p); + __add (&a, &b, &c, p); + if (x > 0.8) + { + __sub (&hp, &c, &b, p); + __c32 (&b, &c, &a, p); + } + else + __c32 (&c, &a, &b, p); /* a = cos(x+dx) */ + __mp_dbl (&a, &y, p); return y; } -/******************************************************************/ -/* mpranred() performs range reduction of a double number x into */ -/* multi precision number y, such that y=x-n*pi/2, abs(y)<pi/4, */ -/* n=0,+-1,+-2,.... */ -/* Return int which indicates in which quarter of circle x is */ -/******************************************************************/ +/* Perform range reduction of a double number x into multi precision number y, + such that y = x - n * pi / 2, abs (y) < pi / 4, n = 0, +-1, +-2, ... + Return int which indicates in which quarter of circle x is. */ int SECTION -__mpranred(double x, mp_no *y, int p) +__mpranred (double x, mp_no *y, int p) { number v; - double t,xn; - int i,k,n; - mp_no a,b,c; + double t, xn; + int i, k, n; + mp_no a, b, c; - if (ABS(x) < 2.8e14) { - t = (x*hpinv.d + toint.d); - xn = t - toint.d; - v.d = t; - n =v.i[LOW_HALF]&3; - __dbl_mp(xn,&a,p); - __mul(&a,&hp,&b,p); - __dbl_mp(x,&c,p); - __sub(&c,&b,y,p); - return n; - } - else { /* if x is very big more precision required */ - __dbl_mp(x,&a,p); - a.d[0]=1.0; - k = a.e-5; - if (k < 0) k=0; - b.e = -k; - b.d[0] = 1.0; - for (i=0;i<p;i++) b.d[i+1] = toverp[i+k]; - __mul(&a,&b,&c,p); - t = c.d[c.e]; - for (i=1;i<=p-c.e;i++) c.d[i]=c.d[i+c.e]; - for (i=p+1-c.e;i<=p;i++) c.d[i]=0; - c.e=0; - if (c.d[1] >= HALFRAD) - { t +=1.0; - __sub(&c,&mpone,&b,p); - __mul(&b,&hp,y,p); + if (ABS (x) < 2.8e14) + { + t = (x * hpinv.d + toint.d); + xn = t - toint.d; + v.d = t; + n = v.i[LOW_HALF] & 3; + __dbl_mp (xn, &a, p); + __mul (&a, &hp, &b, p); + __dbl_mp (x, &c, p); + __sub (&c, &b, y, p); + return n; + } + else + { + /* If x is very big more precision required. */ + __dbl_mp (x, &a, p); + a.d[0] = 1.0; + k = a.e - 5; + if (k < 0) + k = 0; + b.e = -k; + b.d[0] = 1.0; + for (i = 0; i < p; i++) + b.d[i + 1] = toverp[i + k]; + __mul (&a, &b, &c, p); + t = c.d[c.e]; + for (i = 1; i <= p - c.e; i++) + c.d[i] = c.d[i + c.e]; + for (i = p + 1 - c.e; i <= p; i++) + c.d[i] = 0; + c.e = 0; + if (c.d[1] >= HALFRAD) + { + t += 1.0; + __sub (&c, &mpone, &b, p); + __mul (&b, &hp, y, p); + } + else + __mul (&c, &hp, y, p); + n = (int) t; + if (x < 0) + { + y->d[0] = -y->d[0]; + n = -n; + } + return (n & 3); } - else __mul(&c,&hp,y,p); - n = (int) t; - if (x < 0) { y->d[0] = - y->d[0]; n = -n; } - return (n&3); - } } -/*******************************************************************/ -/* Multi-Precision sin() function subroutine, for p=32. It is */ -/* based on the routines mpranred() and c32(). */ -/*******************************************************************/ +/* Multi-Precision sin() function subroutine, for p = 32. It is based on the + routines mpranred() and c32(). */ double SECTION -__mpsin1(double x) +__mpsin1 (double x) { int p; int n; - mp_no u,s,c; + mp_no u, s, c; double y; - p=32; - n=__mpranred(x,&u,p); /* n is 0, 1, 2 or 3 */ - __c32(&u,&c,&s,p); - switch (n) { /* in which quarter of unit circle y is*/ - case 0: - __mp_dbl(&s,&y,p); - return y; - break; + p = 32; + n = __mpranred (x, &u, p); /* n is 0, 1, 2 or 3. */ + __c32 (&u, &c, &s, p); + /* Convert result based on which quarter of unit circle y is in. */ + switch (n) + { + case 0: + __mp_dbl (&s, &y, p); + return y; + break; - case 2: - __mp_dbl(&s,&y,p); - return -y; - break; + case 2: + __mp_dbl (&s, &y, p); + return -y; + break; - case 1: - __mp_dbl(&c,&y,p); - return y; - break; + case 1: + __mp_dbl (&c, &y, p); + return y; + break; - case 3: - __mp_dbl(&c,&y,p); - return -y; - break; - - } - return 0; /* unreachable, to make the compiler happy */ + case 3: + __mp_dbl (&c, &y, p); + return -y; + break; + } + /* Unreachable, to make the compiler happy. */ + return 0; } -/*****************************************************************/ -/* Multi-Precision cos() function subroutine, for p=32. It is */ -/* based on the routines mpranred() and c32(). */ -/*****************************************************************/ - +/* Multi-Precision cos() function subroutine, for p = 32. It is based on the + routines mpranred() and c32(). */ double SECTION -__mpcos1(double x) +__mpcos1 (double x) { int p; int n; - mp_no u,s,c; + mp_no u, s, c; double y; - p=32; - n=__mpranred(x,&u,p); /* n is 0, 1, 2 or 3 */ - __c32(&u,&c,&s,p); - switch (n) { /* in what quarter of unit circle y is*/ + p = 32; + n = __mpranred (x, &u, p); /* n is 0, 1, 2 or 3. */ + __c32 (&u, &c, &s, p); + /* Convert result based on which quarter of unit circle y is in. */ + switch (n) + { + case 0: + __mp_dbl (&c, &y, p); + return y; + break; - case 0: - __mp_dbl(&c,&y,p); - return y; - break; + case 2: + __mp_dbl (&c, &y, p); + return -y; + break; - case 2: - __mp_dbl(&c,&y,p); - return -y; - break; + case 1: + __mp_dbl (&s, &y, p); + return -y; + break; - case 1: - __mp_dbl(&s,&y,p); - return -y; - break; - - case 3: - __mp_dbl(&s,&y,p); - return y; - break; - - } - return 0; /* unreachable, to make the compiler happy */ + case 3: + __mp_dbl (&s, &y, p); + return y; + break; + } + /* Unreachable, to make the compiler happy. */ + return 0; } -/******************************************************************/ |