mirror of
https://sourceware.org/git/glibc.git
synced 2024-11-22 13:00:06 +00:00
Format sincos32.c
This commit is contained in:
Siddhesh Poyarekar
parent
11ca09e932
commit
97a0650b8a
@ -1,3 +1,7 @@
|
||||
2013-09-18 Siddhesh Poyarekar <siddhesh@redhat.com>
|
||||
|
||||
* sysdeps/ieee754/dbl-64/sincos32.c: Fix code formatting.
|
||||
|
||||
2013-09-17 Joseph Myers <joseph@codesourcery.com>
|
||||
|
||||
[BZ #15966]
|
||||
|
||||
@ -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;
|
||||
}
|
||||
/******************************************************************/
|
||||
|
||||
Loading…
Reference in New Issue
Block a user