Format sincos32.c

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;
 }
-/******************************************************************/