summaryrefslogtreecommitdiff
path: root/sysdeps/ieee754/dbl-64/sincos32.c
diff options
context:
space:
mode:
Diffstat (limited to 'sysdeps/ieee754/dbl-64/sincos32.c')
-rw-r--r--sysdeps/ieee754/dbl-64/sincos32.c523
1 files changed, 266 insertions, 257 deletions
diff --git a/sysdeps/ieee754/dbl-64/sincos32.c b/sysdeps/ieee754/dbl-64/sincos32.c
index 954db66d6b..e42fd27250 100644
--- a/sysdeps/ieee754/dbl-64/sincos32.c
+++ b/sysdeps/ieee754/dbl-64/sincos32.c
@@ -43,317 +43,326 @@
#include "mpa.h"
#include "sincos32.h"
#include <math_private.h>
+#include <stap-probe.h>
#ifndef SECTION
# 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 && res >= res1) || (a.d[0] <= 0 && res <= res1))
+ res = res1;
+ LIBC_PROBE (slowasin, 2, &res, &x);
+ return res;
}
-/************************************************************************/
-/*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 && res <= res1) || (a.d[0] <= 0 && res >= res1))
+ res = res1;
+ LIBC_PROBE (slowacos, 2, &res, &x);
+ return res;
}
-/*******************************************************************/
-/*Compute sin(x+dx) as Multi Precision number and return result as */
-/* double */
-/*******************************************************************/
+/* Compute sin() of double-length number (X + DX) as Multi Precision number and
+ return result as double. If REDUCE_RANGE is true, X is assumed to be the
+ original input and DX is ignored. */
double
SECTION
-__mpsin(double x, double dx) {
- int p;
+__mpsin (double x, double dx, bool reduce_range)
+{
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);
- return y;
-}
+ mp_no a, b, c, s;
+ int n;
+ int p = 32;
-/*******************************************************************/
-/* Compute cos()of double-length number (x+dx) as Multi Precision */
-/* number and return result as double */
-/*******************************************************************/
-double
-SECTION
-__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);
- return y;
-}
+ if (reduce_range)
+ {
+ n = __mpranred (x, &a, p); /* n is 0, 1, 2 or 3. */
+ __c32 (&a, &c, &s, p);
+ }
+ else
+ {
+ n = -1;
+ __dbl_mp (x, &b, p);
+ __dbl_mp (dx, &c, p);
+ __add (&b, &c, &a, p);
+ if (x > 0.8)
+ {
+ __sub (&hp, &a, &b, p);
+ __c32 (&b, &s, &c, p);
+ }
+ else
+ __c32 (&a, &c, &s, p); /* b = sin(x+dx) */
+ }
-/******************************************************************/
-/* 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 */
-/******************************************************************/
-int
-SECTION
-__mpranred(double x, mp_no *y, int p)
-{
- number v;
- double t,xn;
- int i,k,n;
- mp_no a,b,c;
+ /* Convert result based on which quarter of unit circle y is in. */
+ switch (n)
+ {
+ case 1:
+ __mp_dbl (&c, &y, p);
+ break;
+
+ case 3:
+ __mp_dbl (&c, &y, p);
+ y = -y;
+ break;
+
+ case 2:
+ __mp_dbl (&s, &y, p);
+ y = -y;
+ break;
- 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);
+ /* Quadrant not set, so the result must be sin (X + DX), which is also in
+ S. */
+ case 0:
+ default:
+ __mp_dbl (&s, &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);
- }
+ LIBC_PROBE (slowsin, 3, &x, &dx, &y);
+ return y;
}
-/*******************************************************************/
-/* Multi-Precision sin() function subroutine, for p=32. It is */
-/* based on the routines mpranred() and c32(). */
-/*******************************************************************/
+/* Compute cos() of double-length number (X + DX) as Multi Precision number and
+ return result as double. If REDUCE_RANGE is true, X is assumed to be the
+ original input and DX is ignored. */
double
SECTION
-__mpsin1(double x)
+__mpcos (double x, double dx, bool reduce_range)
{
- int p;
- int n;
- 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;
+ mp_no a, b, c, s;
+ int n;
+ int p = 32;
- case 2:
- __mp_dbl(&s,&y,p);
- return -y;
- break;
+ if (reduce_range)
+ {
+ n = __mpranred (x, &a, p); /* n is 0, 1, 2 or 3. */
+ __c32 (&a, &c, &s, p);
+ }
+ else
+ {
+ n = -1;
+ __dbl_mp (x, &b, p);
+ __dbl_mp (dx, &c, p);
+ __add (&b, &c, &a, p);
+ if (x > 0.8)
+ {
+ __sub (&hp, &a, &b, p);
+ __c32 (&b, &s, &c, p);
+ }
+ else
+ __c32 (&a, &c, &s, p); /* a = cos(x+dx) */
+ }
- case 1:
- __mp_dbl(&c,&y,p);
- return y;
- break;
+ /* Convert result based on which quarter of unit circle y is in. */
+ switch (n)
+ {
+ case 1:
+ __mp_dbl (&s, &y, p);
+ y = -y;
+ break;
- case 3:
- __mp_dbl(&c,&y,p);
- return -y;
- break;
+ case 3:
+ __mp_dbl (&s, &y, p);
+ break;
- }
- return 0; /* unreachable, to make the compiler happy */
-}
+ case 2:
+ __mp_dbl (&c, &y, p);
+ y = -y;
+ break;
-/*****************************************************************/
-/* Multi-Precision cos() function subroutine, for p=32. It is */
-/* based on the routines mpranred() and c32(). */
-/*****************************************************************/
+ /* Quadrant not set, so the result must be cos (X + DX), which is also
+ stored in C. */
+ case 0:
+ default:
+ __mp_dbl (&c, &y, p);
+ }
+ LIBC_PROBE (slowcos, 3, &x, &dx, &y);
+ return y;
+}
-double
+/* 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
-__mpcos1(double x)
+__mpranred (double x, mp_no *y, int p)
{
- int p;
- int n;
- 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*/
-
- case 0:
- __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 3:
- __mp_dbl(&s,&y,p);
- return y;
- break;
+ number v;
+ double t, xn;
+ int i, k, n;
+ mp_no a, b, c;
- }
- return 0; /* unreachable, to make the compiler happy */
+ 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);
+ }
}
-/******************************************************************/
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-05 15:32:35 +0000