diff options
Diffstat (limited to 'sysdeps/ieee754/dbl-64')
24 files changed, 418 insertions, 197 deletions
diff --git a/sysdeps/ieee754/dbl-64/branred.c b/sysdeps/ieee754/dbl-64/branred.c index 76015f0c5c..c8483034af 100644 --- a/sysdeps/ieee754/dbl-64/branred.c +++ b/sysdeps/ieee754/dbl-64/branred.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001 Free Software Foundation, Inc. + * Copyright (C) 2001, 2011 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -38,6 +38,10 @@ #include "branred.h" #include "math_private.h" +#ifndef SECTION +# define SECTION +#endif + /*******************************************************************/ /* Routine branred() performs range reduction of a double number */ @@ -45,7 +49,9 @@ /* x=n*pi/2+(a+aa), abs(a+aa)<pi/4, n=0,+-1,+-2,.... */ /* Routine return integer (n mod 4) */ /*******************************************************************/ -int __branred(double x, double *a, double *aa) +int +SECTION +__branred(double x, double *a, double *aa) { int i,k; #if 0 diff --git a/sysdeps/ieee754/dbl-64/doasin.c b/sysdeps/ieee754/dbl-64/doasin.c index 64abc3cbb1..14958b5ca2 100644 --- a/sysdeps/ieee754/dbl-64/doasin.c +++ b/sysdeps/ieee754/dbl-64/doasin.c @@ -34,11 +34,17 @@ #include <dla.h> #include "math_private.h" +#ifndef SECTION +# define SECTION +#endif + /********************************************************************/ /* Compute arcsin(x,dx,v) of double-length number (x+dx) the result */ /* stored in v where v= v[0]+v[1] =arcsin(x+dx) */ /********************************************************************/ -void __doasin(double x, double dx, double v[]) { +void +SECTION +__doasin(double x, double dx, double v[]) { #include "doasin.h" diff --git a/sysdeps/ieee754/dbl-64/dosincos.c b/sysdeps/ieee754/dbl-64/dosincos.c index 712d585b9e..e8890ff8de 100644 --- a/sysdeps/ieee754/dbl-64/dosincos.c +++ b/sysdeps/ieee754/dbl-64/dosincos.c @@ -39,6 +39,10 @@ #include "dosincos.h" #include "math_private.h" +#ifndef SECTION +# define SECTION +#endif + extern const union { int4 i[880]; @@ -52,7 +56,9 @@ extern const union /*(x+dx) between 0 and PI/4 */ /***********************************************************************/ -void __dubsin(double x, double dx, double v[]) { +void +SECTION +__dubsin(double x, double dx, double v[]) { double r,s,c,cc,d,dd,d2,dd2,e,ee, sn,ssn,cs,ccs,ds,dss,dc,dcc; #ifndef DLA_FMS @@ -106,7 +112,9 @@ void __dubsin(double x, double dx, double v[]) { /*(x+dx) between 0 and PI/4 */ /**********************************************************************/ -void __dubcos(double x, double dx, double v[]) { +void +SECTION +__dubcos(double x, double dx, double v[]) { double r,s,c,cc,d,dd,d2,dd2,e,ee, sn,ssn,cs,ccs,ds,dss,dc,dcc; #ifndef DLA_FMS @@ -172,7 +180,9 @@ void __dubcos(double x, double dx, double v[]) { /* Routine receive Double-Length number (x+dx) and computes cos(x+dx) */ /* as Double-Length number and store it in array v */ /**********************************************************************/ -void __docos(double x, double dx, double v[]) { +void +SECTION +__docos(double x, double dx, double v[]) { double y,yy,p,w[2]; if (x>0) {y=x; yy=dx;} else {y=-x; yy=-dx;} diff --git a/sysdeps/ieee754/dbl-64/e_asin.c b/sysdeps/ieee754/dbl-64/e_asin.c index cd4cc2e2c2..65319c0b58 100644 --- a/sysdeps/ieee754/dbl-64/e_asin.c +++ b/sysdeps/ieee754/dbl-64/e_asin.c @@ -42,6 +42,10 @@ #include "uasncs.h" #include "math_private.h" +#ifndef SECTION +# define SECTION +#endif + void __doasin(double x, double dx, double w[]); void __dubsin(double x, double dx, double v[]); void __dubcos(double x, double dx, double v[]); @@ -53,7 +57,9 @@ double __cos32(double x, double res, double res1); /* An ultimate asin routine. Given an IEEE double machine number x */ /* it computes the correctly rounded (to nearest) value of arcsin(x) */ /***************************************************************************/ -double __ieee754_asin(double x){ +double +SECTION +__ieee754_asin(double x){ double x1,x2,xx,s1,s2,res1,p,t,res,r,cor,cc,y,c,z,w[2]; mynumber u,v; int4 k,m,n; @@ -334,7 +340,9 @@ strong_alias (__ieee754_asin, __asin_finite) /* */ /*******************************************************************/ -double __ieee754_acos(double x) +double +SECTION +__ieee754_acos(double x) { double x1,x2,xx,s1,s2,res1,p,t,res,r,cor,cc,y,c,z,w[2],eps; #if 0 diff --git a/sysdeps/ieee754/dbl-64/e_atan2.c b/sysdeps/ieee754/dbl-64/e_atan2.c index 9caacccf4c..64dae3e8d5 100644 --- a/sysdeps/ieee754/dbl-64/e_atan2.c +++ b/sysdeps/ieee754/dbl-64/e_atan2.c @@ -44,6 +44,10 @@ #include "atnat2.h" #include "math_private.h" +#ifndef SECTION +# define SECTION +#endif + /************************************************************************/ /* An ultimate atan2 routine. Given two IEEE double machine numbers y,x */ /* it computes the correctly rounded (to nearest) value of atan2(y,x). */ @@ -59,7 +63,9 @@ static double signArctan2(double y,double z) static double normalized(double ,double,double ,double); void __mpatan2(mp_no *,mp_no *,mp_no *,int); -double __ieee754_atan2(double y,double x) { +double +SECTION +__ieee754_atan2(double y,double x) { int i,de,ux,dx,uy,dy; #if 0 @@ -384,7 +390,9 @@ strong_alias (__ieee754_atan2, __atan2_finite) #endif /* Treat the Denormalized case */ -static double normalized(double ax,double ay,double y, double z) +static double +SECTION +normalized(double ax,double ay,double y, double z) { int p; mp_no mpx,mpy,mpz,mperr,mpz2,mpt1; p=6; @@ -394,7 +402,9 @@ static double normalized(double ax,double ay,double y, double z) return signArctan2(y,z); } /* Stage 3: Perform a multi-Precision computation */ -static double atan2Mp(double x,double y,const int pr[]) +static double +SECTION +atan2Mp(double x,double y,const int pr[]) { double z1,z2; int i,p; diff --git a/sysdeps/ieee754/dbl-64/e_exp.c b/sysdeps/ieee754/dbl-64/e_exp.c index 48bbb05ed8..e7a839d42e 100644 --- a/sysdeps/ieee754/dbl-64/e_exp.c +++ b/sysdeps/ieee754/dbl-64/e_exp.c @@ -40,13 +40,19 @@ #include "uexp.tbl" #include "math_private.h" +#ifndef SECTION +# define SECTION +#endif + double __slowexp(double); /***************************************************************************/ /* An ultimate exp routine. Given an IEEE double machine number x */ /* it computes the correctly rounded (to nearest) value of e^x */ /***************************************************************************/ -double __ieee754_exp(double x) { +double +SECTION +__ieee754_exp(double x) { double bexp, t, eps, del, base, y, al, bet, res, rem, cor; mynumber junk1, junk2, binexp = {{0,0}}; #if 0 @@ -156,7 +162,9 @@ strong_alias (__ieee754_exp, __exp_finite) /*else return e^(x + xx) (always positive ) */ /************************************************************************/ -double __exp1(double x, double xx, double error) { +double +SECTION +__exp1(double x, double xx, double error) { double bexp, t, eps, del, base, y, al, bet, res, rem, cor; mynumber junk1, junk2, binexp = {{0,0}}; #if 0 diff --git a/sysdeps/ieee754/dbl-64/e_log.c b/sysdeps/ieee754/dbl-64/e_log.c index 7a0a26f251..e45520eba8 100644 --- a/sysdeps/ieee754/dbl-64/e_log.c +++ b/sysdeps/ieee754/dbl-64/e_log.c @@ -41,13 +41,19 @@ #include "MathLib.h" #include "math_private.h" +#ifndef SECTION +# define SECTION +#endif + void __mplog(mp_no *, mp_no *, int); /*********************************************************************/ /* An ultimate log routine. Given an IEEE double machine number x */ /* it computes the correctly rounded (to nearest) value of log(x). */ /*********************************************************************/ -double __ieee754_log(double x) { +double +SECTION +__ieee754_log(double x) { #define M 4 static const int pr[M]={8,10,18,32}; int i,j,n,ux,dx,p; diff --git a/sysdeps/ieee754/dbl-64/e_pow.c b/sysdeps/ieee754/dbl-64/e_pow.c index 94b1ab8961..350e93986d 100644 --- a/sysdeps/ieee754/dbl-64/e_pow.c +++ b/sysdeps/ieee754/dbl-64/e_pow.c @@ -43,6 +43,10 @@ #include "upow.tbl" #include "math_private.h" +#ifndef SECTION +# define SECTION +#endif + double __exp1(double x, double xx, double error); static double log1(double x, double *delta, double *error); @@ -55,7 +59,9 @@ static int checkint(double x); /* An ultimate power routine. Given two IEEE double machine numbers y,x */ /* it computes the correctly rounded (to nearest) value of X^y. */ /***************************************************************************/ -double __ieee754_pow(double x, double y) { +double +SECTION +__ieee754_pow(double x, double y) { double z,a,aa,error, t,a1,a2,y1,y2; #if 0 double gor=1.0; @@ -160,7 +166,9 @@ strong_alias (__ieee754_pow, __pow_finite) /**************************************************************************/ /* Computing x^y using more accurate but more slow log routine */ /**************************************************************************/ -static double power1(double x, double y) { +static double +SECTION +power1(double x, double y) { double z,a,aa,error, t,a1,a2,y1,y2; z = my_log2(x,&aa,&error); t = y*134217729.0; @@ -183,7 +191,9 @@ static double power1(double x, double y) { /* + the parameter delta. */ /* The result is bounded by error (rightmost argument) */ /****************************************************************************/ -static double log1(double x, double *delta, double *error) { +static double +SECTION +log1(double x, double *delta, double *error) { int i,j,m; #if 0 int n; @@ -275,7 +285,9 @@ static double log1(double x, double *delta, double *error) { /* Computing log(x)(x is left argument).The result is return double + delta.*/ /* The result is bounded by error (right argument) */ /****************************************************************************/ -static double my_log2(double x, double *delta, double *error) { +static double +SECTION +my_log2(double x, double *delta, double *error) { int i,j,m; #if 0 int n; @@ -369,7 +381,9 @@ static double my_log2(double x, double *delta, double *error) { /* Routine receives a double x and checks if it is an integer. If not */ /* it returns 0, else it returns 1 if even or -1 if odd. */ /**********************************************************************/ -static int checkint(double x) { +static int +SECTION +checkint(double x) { union {int4 i[2]; double x;} u; int k,m,n; #if 0 diff --git a/sysdeps/ieee754/dbl-64/halfulp.c b/sysdeps/ieee754/dbl-64/halfulp.c index 31bd2daf4d..6018309427 100644 --- a/sysdeps/ieee754/dbl-64/halfulp.c +++ b/sysdeps/ieee754/dbl-64/halfulp.c @@ -40,6 +40,10 @@ #include <dla.h> #include "math_private.h" +#ifndef SECTION +# define SECTION +#endif + static const int4 tab54[32] = { 262143, 11585, 1782, 511, 210, 107, 63, 42, 30, 22, 17, 14, 12, 10, 9, 7, @@ -47,7 +51,9 @@ static const int4 tab54[32] = { 3, 3, 3, 3, 3, 3, 3, 3 }; -double __halfulp(double x, double y) +double +SECTION +__halfulp(double x, double y) { mynumber v; double z,u,uu; diff --git a/sysdeps/ieee754/dbl-64/mpa.c b/sysdeps/ieee754/dbl-64/mpa.c index ad5a639c4b..39c640882b 100644 --- a/sysdeps/ieee754/dbl-64/mpa.c +++ b/sysdeps/ieee754/dbl-64/mpa.c @@ -47,11 +47,18 @@ #include "mpa.h" #include "mpa2.h" #include <sys/param.h> /* For MIN() */ + +#ifndef SECTION +# define SECTION +#endif + +#ifndef NO___ACR /* mcr() compares the sizes of the mantissas of two multiple precision */ /* numbers. Mantissas are compared regardless of the signs of the */ /* numbers, even if x->d[0] or y->d[0] are zero. Exponents are also */ /* disregarded. */ -static int mcr(const mp_no *x, const mp_no *y, int p) { +static int +mcr(const mp_no *x, const mp_no *y, int p) { int i; for (i=1; i<=p; i++) { if (X[i] == Y[i]) continue; @@ -61,9 +68,9 @@ static int mcr(const mp_no *x, const mp_no *y, int p) { } - /* acr() compares the absolute values of two multiple precision numbers */ -static int __acr(const mp_no *x, const mp_no *y, int p) { +int +__acr(const mp_no *x, const mp_no *y, int p) { int i; if (X[0] == ZERO) { @@ -79,10 +86,11 @@ static int __acr(const mp_no *x, const mp_no *y, int p) { return i; } +#endif #if 0 -/* cr90 compares the values of two multiple precision numbers */ +/* cr() compares the values of two multiple precision numbers */ static int __cr(const mp_no *x, const mp_no *y, int p) { int i; @@ -119,8 +127,6 @@ static void __cpymn(const mp_no *x, int m, mp_no *y, int n) { EY = EX; k=MIN(m,n); for (i=0; i <= k; i++) Y[i] = X[i]; for ( ; i <= n; i++) Y[i] = ZERO; - - return; } #endif @@ -177,7 +183,6 @@ static void norm(const mp_no *x, double *y, int p) for (i=1; i>EX; i--) c *= RADIXI; *y = c; - return; #undef R } @@ -225,8 +230,6 @@ static void denorm(const mp_no *x, double *y, int p) c = X[0]*((z[1] + R*(z[2] + R*z[3])) - TWO10); *y = c*TWOM1032; - return; - #undef R } @@ -252,7 +255,9 @@ void __mp_dbl(const mp_no *x, double *y, int p) { /* number *y. If the precision p is too small the result is truncated. x is */ /* left unchanged. */ -void __dbl_mp(double x, mp_no *y, int p) { +void +SECTION +__dbl_mp(double x, mp_no *y, int p) { int i,n; double u; @@ -273,7 +278,6 @@ void __dbl_mp(double x, mp_no *y, int p) { if (u>x) u -= ONE; Y[i] = u; x -= u; x *= RADIX; } for ( ; i<=p; i++) Y[i] = ZERO; - return; } @@ -283,7 +287,9 @@ void __dbl_mp(double x, mp_no *y, int p) { /* No guard digit is used. The result equals the exact sum, truncated. */ /* *x & *y are left unchanged. */ -static void add_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) { +static void +SECTION +add_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) { int i,j,k; @@ -325,7 +331,9 @@ static void add_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) { /* or y&z. One guard digit is used. The error is less than one ulp. */ /* *x & *y are left unchanged. */ -static void sub_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) { +static void +SECTION +sub_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) { int i,j,k; @@ -372,8 +380,6 @@ static void sub_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) { Z[k++] = Z[i++]; for (; k <= p; ) Z[k++] = ZERO; - - return; } @@ -381,7 +387,9 @@ static void sub_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) { /* but not x&z or y&z. One guard digit is used. The error is less than */ /* one ulp. *x & *y are left unchanged. */ -void __add(const mp_no *x, const mp_no *y, mp_no *z, int p) { +void +SECTION +__add(const mp_no *x, const mp_no *y, mp_no *z, int p) { int n; @@ -397,7 +405,6 @@ void __add(const mp_no *x, const mp_no *y, mp_no *z, int p) { else if (n == -1) {sub_magnitudes(y,x,z,p); Z[0] = Y[0]; } else Z[0] = ZERO; } - return; } @@ -405,7 +412,9 @@ void __add(const mp_no *x, const mp_no *y, mp_no *z, int p) { /* overlap but not x&z or y&z. One guard digit is used. The error is */ /* less than one ulp. *x & *y are left unchanged. */ -void __sub(const mp_no *x, const mp_no *y, mp_no *z, int p) { +void +SECTION +__sub(const mp_no *x, const mp_no *y, mp_no *z, int p) { int n; @@ -421,7 +430,6 @@ void __sub(const mp_no *x, const mp_no *y, mp_no *z, int p) { else if (n == -1) {sub_magnitudes(y,x,z,p); Z[0] = -Y[0]; } else Z[0] = ZERO; } - return; } @@ -430,7 +438,9 @@ void __sub(const mp_no *x, const mp_no *y, mp_no *z, int p) { /* truncated to p digits. In case p>3 the error is bounded by 1.001 ulp. */ /* *x & *y are left unchanged. */ -void __mul(const mp_no *x, const mp_no *y, mp_no *z, int p) { +void +SECTION +__mul(const mp_no *x, const mp_no *y, mp_no *z, int p) { int i, i1, i2, j, k, k2; double u; @@ -461,7 +471,6 @@ void __mul(const mp_no *x, const mp_no *y, mp_no *z, int p) { EZ = EX + EY; Z[0] = X[0] * Y[0]; - return; } @@ -470,7 +479,9 @@ void __mul(const mp_no *x, const mp_no *y, mp_no *z, int p) { /* 2.001*r**(1-p) for p>3. */ /* *x=0 is not permissible. *x is left unchanged. */ -static void __inv(const mp_no *x, mp_no *y, int p) { +static +SECTION +void __inv(const mp_no *x, mp_no *y, int p) { int i; #if 0 int l; @@ -493,7 +504,6 @@ static void __inv(const mp_no *x, mp_no *y, int p) { __sub(&mptwo,y,&z,p); __mul(&w,&z,y,p); } - return; } @@ -502,11 +512,12 @@ static void __inv(const mp_no *x, mp_no *y, int p) { /* Relative error bound = 2.001*r**(1-p) for p=2, 2.063*r**(1-p) for p=3 */ /* and 3.001*r**(1-p) for p>3. *y=0 is not permissible. */ -void __dvd(const mp_no *x, const mp_no *y, mp_no *z, int p) { +void +SECTION +__dvd(const mp_no *x, const mp_no *y, mp_no *z, int p) { mp_no w; if (X[0] == ZERO) Z[0] = ZERO; else {__inv(y,&w,p); __mul(x,&w,z,p);} - return; } diff --git a/sysdeps/ieee754/dbl-64/mpa.h b/sysdeps/ieee754/dbl-64/mpa.h index 3ca0ca5f00..5647ab7b4f 100644 --- a/sysdeps/ieee754/dbl-64/mpa.h +++ b/sysdeps/ieee754/dbl-64/mpa.h @@ -64,7 +64,7 @@ typedef union { int i[2]; double d; } number; #define ABS(x) ((x) < 0 ? -(x) : (x)) -// int __acr(const mp_no *, const mp_no *, int); +int __acr(const mp_no *, const mp_no *, int); // int __cr(const mp_no *, const mp_no *, int); void __cpy(const mp_no *, mp_no *, int); // void __cpymn(const mp_no *, int, mp_no *, int); diff --git a/sysdeps/ieee754/dbl-64/mpatan.c b/sysdeps/ieee754/dbl-64/mpatan.c index ee21c25138..f40873ea59 100644 --- a/sysdeps/ieee754/dbl-64/mpatan.c +++ b/sysdeps/ieee754/dbl-64/mpatan.c @@ -1,8 +1,7 @@ - /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001 Free Software Foundation + * Copyright (C) 2001, 2011 Free Software Foundation * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -34,11 +33,19 @@ #include "endian.h" #include "mpa.h" -void __mpsqrt(mp_no *, mp_no *, int); -void __mpatan(mp_no *x, mp_no *y, int p) { +#ifndef SECTION +# define SECTION +#endif + #include "mpatan.h" +void __mpsqrt(mp_no *, mp_no *, int); + +void +SECTION +__mpatan(mp_no *x, mp_no *y, int p) { + int i,m,n; double dx; mp_no @@ -54,19 +61,19 @@ void __mpatan(mp_no *x, mp_no *y, int p) { mp_no mps,mpsm,mpt,mpt1,mpt2,mpt3; - /* Choose m and initiate mpone, mptwo & mptwoim1 */ + /* Choose m and initiate mpone, mptwo & mptwoim1 */ if (EX>0) m=7; else if (EX<0) m=0; else { __mp_dbl(x,&dx,p); dx=ABS(dx); for (m=6; m>0; m--) - {if (dx>xm[m].d) break;} + {if (dx>__atan_xm[m].d) break;} } mpone.e = mptwo.e = mptwoim1.e = 1; mpone.d[0] = mpone.d[1] = mptwo.d[0] = mptwoim1.d[0] = ONE; mptwo.d[1] = TWO; - /* Reduce x m times */ + /* Reduce x m times */ __mul(x,x,&mpsm,p); if (m==0) __cpy(x,&mps,p); else { @@ -82,8 +89,8 @@ void __mpatan(mp_no *x, mp_no *y, int p) { __mpsqrt(&mpsm,&mps,p); mps.d[0] = X[0]; } - /* Evaluate a truncated power series for Atan(s) */ - n=np[p]; mptwoim1.d[1] = twonm1[p].d; + /* Evaluate a truncated power series for Atan(s) */ + n=__atan_np[p]; mptwoim1.d[1] = __atan_twonm1[p].d; __dvd(&mpsm,&mptwoim1,&mpt,p); for (i=n-1; i>1; i--) { mptwoim1.d[1] -= TWO; @@ -94,8 +101,8 @@ void __mpatan(mp_no *x, mp_no *y, int p) { __mul(&mps,&mpt,&mpt1,p); __sub(&mps,&mpt1,&mpt,p); - /* Compute Atan(x) */ - mptwoim1.d[1] = twom[m].d; + /* Compute Atan(x) */ + mptwoim1.d[1] = __atan_twom[m].d; __mul(&mptwoim1,&mpt,y,p); return; diff --git a/sysdeps/ieee754/dbl-64/mpatan.h b/sysdeps/ieee754/dbl-64/mpatan.h index d420ff3408..003b06c695 100644 --- a/sysdeps/ieee754/dbl-64/mpatan.h +++ b/sysdeps/ieee754/dbl-64/mpatan.h @@ -1,8 +1,7 @@ - /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001 Free Software Foundation, Inc. + * Copyright (C) 2001, 2011 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -29,9 +28,18 @@ #ifndef MPATAN_H #define MPATAN_H +extern const number __atan_xm[8] attribute_hidden; +extern const number __atan_twonm1[33] attribute_hidden; +extern const number __atan_twom[8] attribute_hidden; +extern const number __atan_one attribute_hidden; +extern const number __atan_two attribute_hidden; +extern const int __atan_np[33] attribute_hidden; + + +#ifndef AVOID_MPATAN_H #ifdef BIG_ENDI - static const number - xm[8] = { /* x[m] */ + const number + __atan_xm[8] = { /* x[m] */ /**/ {{0x00000000, 0x00000000} }, /* 0.0 */ /**/ {{0x3f8930be, 0x00000000} }, /* 0.0123 */ /**/ {{0x3f991687, 0x00000000} }, /* 0.0245 */ @@ -40,9 +48,9 @@ /**/ {{0x3fc95810, 0x00000000} }, /* 0.198 */ /**/ {{0x3fda7ef9, 0x00000000} }, /* 0.414 */ /**/ {{0x3ff00000, 0x00000000} }, /* 1.0 */ - }; - static const number - twonm1[33] = { /* 2n-1 */ + }; + const number + __atan_twonm1[33] = { /* 2n-1 */ /**/ {{0x00000000, 0x00000000} }, /* 0 */ /**/ {{0x00000000, 0x00000000} }, /* 0 */ /**/ {{0x00000000, 0x00000000} }, /* 0 */ @@ -76,10 +84,10 @@ /**/ {{0x405b4000, 0x00000000} }, /* 109 */ /**/ {{0x405c4000, 0x00000000} }, /* 113 */ /**/ {{0x405d4000, 0x00000000} }, /* 117 */ - }; + }; - static const number - twom[8] = { /* 2**m */ + const number + __atan_twom[8] = { /* 2**m */ /**/ {{0x3ff00000, 0x00000000} }, /* 1.0 */ /**/ {{0x40000000, 0x00000000} }, /* 2.0 */ /**/ {{0x40100000, 0x00000000} }, /* 4.0 */ @@ -88,17 +96,17 @@ /**/ {{0x40400000, 0x00000000} }, /* 32.0 */ /**/ {{0x40500000, 0x00000000} }, /* 64.0 */ /**/ {{0x40600000, 0x00000000} }, /* 128.0 */ - }; + }; - static const number -/**/ one = {{0x3ff00000, 0x00000000} }, /* 1 */ -/**/ two = {{0x40000000, 0x00000000} }; /* 2 */ + const number +/**/ __atan_one = {{0x3ff00000, 0x00000000} }, /* 1 */ +/**/ __atan_two = {{0x40000000, 0x00000000} }; /* 2 */ #else #ifdef LITTLE_ENDI - static const number - xm[8] = { /* x[m] */ + const number + __atan_xm[8] = { /* x[m] */ /**/ {{0x00000000, 0x00000000} }, /* 0.0 */ /**/ {{0x00000000, 0x3f8930be} }, /* 0.0123 */ /**/ {{0x00000000, 0x3f991687} }, /* 0.0245 */ @@ -107,9 +115,9 @@ /**/ {{0x00000000, 0x3fc95810} }, /* 0.198 */ /**/ {{0x00000000, 0x3fda7ef9} }, /* 0.414 */ /**/ {{0x00000000, 0x3ff00000} }, /* 1.0 */ - }; - static const number - twonm1[33] = { /* 2n-1 */ + }; + const number +__atan_twonm1[33] = { /* 2n-1 */ /**/ {{0x00000000, 0x00000000} }, /* 0 */ /**/ {{0x00000000, 0x00000000} }, /* 0 */ /**/ {{0x00000000, 0x00000000} }, /* 0 */ @@ -143,10 +151,10 @@ /**/ {{0x00000000, 0x405b4000} }, /* 109 */ /**/ {{0x00000000, 0x405c4000} }, /* 113 */ /**/ {{0x00000000, 0x405d4000} }, /* 117 */ - }; + }; - static const number - twom[8] = { /* 2**m */ + const number + __atan_twom[8] = { /* 2**m */ /**/ {{0x00000000, 0x3ff00000} }, /* 1.0 */ /**/ {{0x00000000, 0x40000000} }, /* 2.0 */ /**/ {{0x00000000, 0x40100000} }, /* 4.0 */ @@ -155,20 +163,21 @@ /**/ {{0x00000000, 0x40400000} }, /* 32.0 */ /**/ {{0x00000000, 0x40500000} }, /* 64.0 */ /**/ {{0x00000000, 0x40600000} }, /* 128.0 */ - }; + }; - static const number -/**/ one = {{0x00000000, 0x3ff00000} }, /* 1 */ -/**/ two = {{0x00000000, 0x40000000} }; /* 2 */ + const number +/**/ __atan_one = {{0x00000000, 0x3ff00000} }, /* 1 */ +/**/ __atan_two = {{0x00000000, 0x40000000} }; /* 2 */ #endif #endif -#define ONE one.d -#define TWO two.d - - static const int - np[33] = { 0, 0, 0, 0, 6, 8,10,11,13,15,17,19,21,23,25,27,28, - 30,32,34,36,38,40,42,43,45,47,49,51,53,55,57,59}; + const int + __atan_np[33] = { 0, 0, 0, 0, 6, 8,10,11,13,15,17,19,21,23,25,27,28, + 30,32,34,36,38,40,42,43,45,47,49,51,53,55,57,59}; #endif +#endif + +#define ONE __atan_one.d +#define TWO __atan_two.d diff --git a/sysdeps/ieee754/dbl-64/mpatan2.c b/sysdeps/ieee754/dbl-64/mpatan2.c index 8977ec9042..1deb056417 100644 --- a/sysdeps/ieee754/dbl-64/mpatan2.c +++ b/sysdeps/ieee754/dbl-64/mpatan2.c @@ -1,8 +1,7 @@ - /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001 Free Software Foundation + * Copyright (C) 2001, 2011 Free Software Foundation * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -38,18 +37,24 @@ #include "mpa.h" +#ifndef SECTION +# define SECTION +#endif + void __mpsqrt(mp_no *, mp_no *, int); void __mpatan(mp_no *, mp_no *, int); /* Multi-Precision Atan2(y,x) function subroutine, for p >= 4. */ /* y=0 is not permitted if x<=0. No error messages are given. */ -void __mpatan2(mp_no *y, mp_no *x, mp_no *z, int p) { +void +SECTION +__mpatan2(mp_no *y, mp_no *x, mp_no *z, int p) { static const double ZERO = 0.0, ONE = 1.0; mp_no mpone = {0,{0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, - 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, - 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}}; + 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, + 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}}; mp_no mpt1,mpt2,mpt3; diff --git a/sysdeps/ieee754/dbl-64/mpexp.c b/sysdeps/ieee754/dbl-64/mpexp.c index e2ab71b2cc..b0cffe2fe5 100644 --- a/sysdeps/ieee754/dbl-64/mpexp.c +++ b/sysdeps/ieee754/dbl-64/mpexp.c @@ -1,8 +1,7 @@ - /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001 Free Software Foundation + * Copyright (C) 2001, 2011 Free Software Foundation * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -34,34 +33,40 @@ #include "mpa.h" #include "mpexp.h" +#ifndef SECTION +# define SECTION +#endif + /* Multi-Precision exponential function subroutine (for p >= 4, */ /* 2**(-55) <= abs(x) <= 1024). */ -void __mpexp(mp_no *x, mp_no *y, int p) { +void +SECTION +__mpexp(mp_no *x, mp_no *y, int p) { int i,j,k,m,m1,m2,n; double a,b; static const int np[33] = {0,0,0,0,3,3,4,4,5,4,4,5,5,5,6,6,6,6,6,6, - 6,6,6,6,7,7,7,7,8,8,8,8,8}; + 6,6,6,6,7,7,7,7,8,8,8,8,8}; static const int m1p[33]= {0,0,0,0,17,23,23,28,27,38,42,39,43,47,43,47,50,54, - 57,60,64,67,71,74,68,71,74,77,70,73,76,78,81}; + 57,60,64,67,71,74,68,71,74,77,70,73,76,78,81}; static const int m1np[7][18] = { - { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, - { 0, 0, 0, 0,36,48,60,72, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, - { 0, 0, 0, 0,24,32,40,48,56,64,72, 0, 0, 0, 0, 0, 0, 0}, - { 0, 0, 0, 0,17,23,29,35,41,47,53,59,65, 0, 0, 0, 0, 0}, - { 0, 0, 0, 0, 0, 0,23,28,33,38,42,47,52,57,62,66, 0, 0}, - { 0, 0, 0, 0, 0, 0, 0, 0,27, 0, 0,39,43,47,51,55,59,63}, - { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,43,47,50,54}}; + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + { 0, 0, 0, 0,36,48,60,72, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + { 0, 0, 0, 0,24,32,40,48,56,64,72, 0, 0, 0, 0, 0, 0, 0}, + { 0, 0, 0, 0,17,23,29,35,41,47,53,59,65, 0, 0, 0, 0, 0}, + { 0, 0, 0, 0, 0, 0,23,28,33,38,42,47,52,57,62,66, 0, 0}, + { 0, 0, 0, 0, 0, 0, 0, 0,27, 0, 0,39,43,47,51,55,59,63}, + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,43,47,50,54}}; mp_no mpone = {0,{0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, - 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, - 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}}; + 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, + 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}}; mp_no mpk = {0,{0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, - 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, - 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}}; + 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, + 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}}; mp_no mps,mpak,mpt1,mpt2; /* Choose m,n and compute a=2**(-m) */ - n = np[p]; m1 = m1p[p]; a = twomm1[p].d; + n = np[p]; m1 = m1p[p]; a = __mpexp_twomm1[p].d; for (i=0; i<EX; i++) a *= RADIXI; for ( ; i>EX; i--) a *= RADIX; b = X[1]*RADIXI; m2 = 24*EX; @@ -81,12 +86,12 @@ void __mpexp(mp_no *x, mp_no *y, int p) { /* Evaluate the polynomial. Put result in mpt2 */ mpone.e=1; mpone.d[0]=ONE; mpone.d[1]=ONE; - mpk.e = 1; mpk.d[0] = ONE; mpk.d[1]=nn[n].d; + mpk.e = 1; mpk.d[0] = ONE; mpk.d[1]=__mpexp_nn[n].d; __dvd(&mps,&mpk,&mpt1,p); __add(&mpone,&mpt1,&mpak,p); for (k=n-1; k>1; k--) { __mul(&mps,&mpak,&mpt1,p); - mpk.d[1]=nn[k].d; + mpk.d[1]=__mpexp_nn[k].d; __dvd(&mpt1,&mpk,&mpt2,p); __add(&mpone,&mpt2,&mpak,p); } diff --git a/sysdeps/ieee754/dbl-64/mpexp.h b/sysdeps/ieee754/dbl-64/mpexp.h index a0b08ccb41..7985060a8c 100644 --- a/sysdeps/ieee754/dbl-64/mpexp.h +++ b/sysdeps/ieee754/dbl-64/mpexp.h @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001 Free Software Foundation, Inc. + * Copyright (C) 2001, 2011 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -28,9 +28,20 @@ #ifndef MPEXP_H #define MPEXP_H +extern const number __mpexp_twomm1[33] attribute_hidden; +extern const number __mpexp_nn[9] attribute_hidden; +extern const number __mpexp_radix attribute_hidden; +extern const number __mpexp_radixi attribute_hidden; +extern const number __mpexp_zero attribute_hidden; +extern const number __mpexp_one attribute_hidden; +extern const number __mpexp_two attribute_hidden; +extern const number __mpexp_half attribute_hidden; + + +#ifndef AVOID_MPEXP_H #ifdef BIG_ENDI - static const number - twomm1[33] = { /* 2**-m1 */ + const number + __mpexp_twomm1[33] = { /* 2**-m1 */ /**/ {{0x00000000, 0x00000000} }, /* 0 */ /**/ {{0x00000000, 0x00000000} }, /* 0 */ /**/ {{0x00000000, 0x00000000} }, /* 0 */ @@ -65,8 +76,8 @@ /**/ {{0x3b100000, 0x00000000} }, /* 2**-78 */ /**/ {{0x3ae00000, 0x00000000} }, /* 2**-81 */ }; - static const number - nn[9]={ /* n */ + const number + __mpexp_nn[9]={ /* n */ /**/ {{0x00000000, 0x00000000} }, /* 0 */ /**/ {{0x3ff00000, 0x00000000} }, /* 1 */ /**/ {{0x40000000, 0x00000000} }, /* 2 */ @@ -78,18 +89,18 @@ /**/ {{0x40200000, 0x00000000} }, /* 8 */ }; - static const number -/**/ radix = {{0x41700000, 0x00000000} }, /* 2**24 */ -/**/ radixi = {{0x3e700000, 0x00000000} }, /* 2**-24 */ -/**/ zero = {{0x00000000, 0x00000000} }, /* 0 */ -/**/ one = {{0x3ff00000, 0x00000000} }, /* 1 */ -/**/ two = {{0x40000000, 0x00000000} }, /* 2 */ -/**/ half = {{0x3fe00000, 0x00000000} }; /* 1/2 */ + const number +/**/ __mpexp_radix = {{0x41700000, 0x00000000} }, /* 2**24 */ +/**/ __mpexp_radixi = {{0x3e700000, 0x00000000} }, /* 2**-24 */ +/**/ __mpexp_zero = {{0x00000000, 0x00000000} }, /* 0 */ +/**/ __mpexp_one = {{0x3ff00000, 0x00000000} }, /* 1 */ +/**/ __mpexp_two = {{0x40000000, 0x00000000} }, /* 2 */ +/**/ __mpexp_half = {{0x3fe00000, 0x00000000} }; /* 1/2 */ #else #ifdef LITTLE_ENDI - static const number - twomm1[33] = { /* 2**-m1 */ + const number + __mpexp_twomm1[33] = { /* 2**-m1 */ /**/ {{0x00000000, 0x00000000} }, /* 0 */ /**/ {{0x00000000, 0x00000000} }, /* 0 */ /**/ {{0x00000000, 0x00000000} }, /* 0 */ @@ -124,8 +135,8 @@ /**/ {{0x00000000, 0x3b100000} }, /* 2**-78 */ /**/ {{0x00000000, 0x3ae00000} }, /* 2**-81 */ }; - static const number - nn[9]={ /* n */ + const number + __mpexp_nn[9]={ /* n */ /**/ {{0x00000000, 0x00000000} }, /* 0 */ /**/ {{0x00000000, 0x3ff00000} }, /* 1 */ /**/ {{0x00000000, 0x40000000} }, /* 2 */ @@ -137,22 +148,23 @@ /**/ {{0x00000000, 0x40200000} }, /* 8 */ }; - static const number -/**/ radix = {{0x00000000, 0x41700000} }, /* 2**24 */ -/**/ radixi = {{0x00000000, 0x3e700000} }, /* 2**-24 */ -/**/ zero = {{0x00000000, 0x00000000} }, /* 0 */ -/**/ one = {{0x00000000, 0x3ff00000} }, /* 1 */ -/**/ two = {{0x00000000, 0x40000000} }, /* 2 */ -/**/ half = {{0x00000000, 0x3fe00000} }; /* 1/2 */ + const number +/**/ __mpexp_radix = {{0x00000000, 0x41700000} }, /* 2**24 */ +/**/ __mpexp_radixi = {{0x00000000, 0x3e700000} }, /* 2**-24 */ +/**/ __mpexp_zero = {{0x00000000, 0x00000000} }, /* 0 */ +/**/ __mpexp_one = {{0x00000000, 0x3ff00000} }, /* 1 */ +/**/ __mpexp_two = {{0x00000000, 0x40000000} }, /* 2 */ +/**/ __mpexp_half = {{0x00000000, 0x3fe00000} }; /* 1/2 */ #endif #endif +#endif -#define RADIX radix.d -#define RADIXI radixi.d -#define ZERO zero.d -#define ONE one.d -#define TWO two.d -#define HALF half.d +#define RADIX __mpexp_radix.d +#define RADIXI __mpexp_radixi.d +#define ZERO __mpexp_zero.d +#define ONE __mpexp_one.d +#define TWO __mpexp_two.d +#define HALF __mpexp_half.d #endif diff --git a/sysdeps/ieee754/dbl-64/mpsqrt.c b/sysdeps/ieee754/dbl-64/mpsqrt.c index bea623296b..d1a80f9091 100644 --- a/sysdeps/ieee754/dbl-64/mpsqrt.c +++ b/sysdeps/ieee754/dbl-64/mpsqrt.c @@ -32,6 +32,12 @@ #include "endian.h" #include "mpa.h" +#ifndef SECTION +# define SECTION +#endif + +#include "mpsqrt.h" + /****************************************************************************/ /* Multi-Precision square root function subroutine for precision p >= 4. */ /* The relative error is bounded by 3.501*r**(1-p), where r=2**24. */ @@ -42,9 +48,9 @@ static double fastiroot(double); -void __mpsqrt(mp_no *x, mp_no *y, int p) { -#include "mpsqrt.h" - +void +SECTION +__mpsqrt(mp_no *x, mp_no *y, int p) { int i,m,ex,ey; double dx,dy; mp_no @@ -64,7 +70,7 @@ void __mpsqrt(mp_no *x, mp_no *y, int p) { __mp_dbl(&mpxn,&dx,p); dy=fastiroot(dx); __dbl_mp(dy,&mpu,p); __mul(&mpxn,&mphalf,&mpz,p); - m=mp[p]; + m=__mpsqrt_mp[p]; for (i=0; i<m; i++) { __mul(&mpu,&mpu,&mpt1,p); __mul(&mpt1,&mpz,&mpt2,p); @@ -81,7 +87,9 @@ void __mpsqrt(mp_no *x, mp_no *y, int p) { /* Compute a double precision approximation for 1/sqrt(x) */ /* with the relative error bounded by 2**-51. */ /***********************************************************/ -static double fastiroot(double x) { +static double +SECTION +fastiroot(double x) { union {int i[2]; double d;} p,q; double y,z, t; int n; diff --git a/sysdeps/ieee754/dbl-64/mpsqrt.h b/sysdeps/ieee754/dbl-64/mpsqrt.h index 729d57af2c..86fa397b22 100644 --- a/sysdeps/ieee754/dbl-64/mpsqrt.h +++ b/sysdeps/ieee754/dbl-64/mpsqrt.h @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001 Free Software Foundation, Inc. + * Copyright (C) 2001, 2011 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -28,24 +28,31 @@ #ifndef MPSQRT_H #define MPSQRT_H +extern const number __mpsqrt_one attribute_hidden; +extern const number __mpsqrt_halfrad attribute_hidden; +extern const int __mpsqrt_mp[33] attribute_hidden; + + +#ifndef AVOID_MPSQRT_H #ifdef BIG_ENDI - static const number -/**/ one = {{0x3ff00000, 0x00000000} }, /* 1 */ -/**/ halfrad = {{0x41600000, 0x00000000} }; /* 2**23 */ + const number +/**/ __mpsqrt_one = {{0x3ff00000, 0x00000000} }, /* 1 */ +/**/ __mpsqrt_halfrad = {{0x41600000, 0x00000000} }; /* 2**23 */ #else #ifdef LITTLE_ENDI - static const number -/**/ one = {{0x00000000, 0x3ff00000} }, /* 1 */ -/**/ halfrad = {{0x00000000, 0x41600000} }; /* 2**23 */ + const number +/**/ __mpsqrt_one = {{0x00000000, 0x3ff00000} }, /* 1 */ +/**/ __mpsqrt_halfrad = {{0x00000000, 0x41600000} }; /* 2**23 */ #endif #endif -#define ONE one.d -#define HALFRAD halfrad.d + const int __mpsqrt_mp[33] = {0,0,0,0,1,2,2,2,2,3,3,3,3,3,3,3,3,4,4,4,4,4,4,4, + 4,4,4,4,4,4,4,4,4}; +#endif - static const int mp[33] = {0,0,0,0,1,2,2,2,2,3,3,3,3,3,3,3,3,4,4,4,4,4,4,4, - 4,4,4,4,4,4,4,4,4}; +#define ONE __mpsqrt_one.d +#define HALFRAD __mpsqrt_halfrad.d #endif diff --git a/sysdeps/ieee754/dbl-64/mptan.c b/sysdeps/ieee754/dbl-64/mptan.c index 267445a19e..e1e5d9b925 100644 --- a/sysdeps/ieee754/dbl-64/mptan.c +++ b/sysdeps/ieee754/dbl-64/mptan.c @@ -1,8 +1,7 @@ - /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001 Free Software Foundation + * Copyright (C) 2001, 2011 Free Software Foundation * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -38,10 +37,16 @@ #include "endian.h" #include "mpa.h" +#ifndef SECTION +# define SECTION +#endif + int __mpranred(double, mp_no *, int); void __c32(mp_no *, mp_no *, mp_no *, int); -void __mptan(double x, mp_no *mpy, int p) { +void +SECTION +__mptan(double x, mp_no *mpy, int p) { static const double MONE = -1.0; diff --git a/sysdeps/ieee754/dbl-64/s_sin.c b/sysdeps/ieee754/dbl-64/s_sin.c index 02d428ca03..6f19f158f1 100644 --- a/sysdeps/ieee754/dbl-64/s_sin.c +++ b/sysdeps/ieee754/dbl-64/s_sin.c @@ -55,6 +55,10 @@ #include "MathLib.h" #include "math_private.h" +#ifndef SECTION +# define SECTION +#endif + extern const union { int4 i[880]; @@ -92,7 +96,9 @@ static double csloww2(double x, double dx, double orig, int n); /* An ultimate sin routine. Given an IEEE double machine number x */ /* it computes the correctly rounded (to nearest) value of sin(x) */ /*******************************************************************/ -double __sin(double x){ +double +SECTION +__sin(double x){ double xx,res,t,cor,y,s,c,sn,ssn,cs,ccs,xn,a,da,db,eps,xn1,xn2; #if 0 double w[2]; @@ -349,7 +355,9 @@ double __sin(double x){ /* it computes the correctly rounded (to nearest) value of cos(x) */ /*******************************************************************/ -double __cos(double x) +double +SECTION +__cos(double x) { double y,xx,res,t,cor,s,c,sn,ssn,cs,ccs,xn,a,da,db,eps,xn1,xn2; mynumber u,v; @@ -596,7 +604,9 @@ double __cos(double x) /* precision and if still doesn't accurate enough by mpsin or dubsin */ /************************************************************************/ -static double slow(double x) { +static double +SECTION +slow(double x) { static const double th2_36 = 206158430208.0; /* 1.5*2**37 */ double y,x1,x2,xx,r,t,res,cor,w[2]; x1=(x+th2_36)-th2_36; @@ -620,7 +630,9 @@ static const double th2_36 = 206158430208.0; /* 1.5*2**37 */ /* and if result still doesn't accurate enough by mpsin or dubsin */ /*******************************************************************************/ -static double slow1(double x) { +static double +SECTION +slow1(double x) { mynumber u; double sn,ssn,cs,ccs,s,c,w[2],y,y1,y2,c1,c2,xx,cor,res; static const double t22 = 6291456.0; @@ -656,7 +668,9 @@ static double slow1(double x) { /* Routine compute sin(x) for 0.855469 <|x|<2.426265 by __sincostab.tbl */ /* and if result still doesn't accurate enough by mpsin or dubsin */ /**************************************************************************/ -static double slow2(double x) { +static double +SECTION +slow2(double x) { mynumber u; double sn,ssn,cs,ccs,s,c,w[2],y,y1,y2,e1,e2,xx,cor,res,del; static const double t22 = 6291456.0; @@ -708,7 +722,9 @@ static double slow2(double x) { /* result.And if result not accurate enough routine calls mpsin1 or dubsin */ /***************************************************************************/ -static double sloww(double x,double dx, double orig) { +static double +SECTION +sloww(double x,double dx, double orig) { static const double th2_36 = 206158430208.0; /* 1.5*2**37 */ double y,x1,x2,xx,r,t,res,cor,w[2],a,da,xn; union {int4 i[2]; double x;} v; @@ -755,7 +771,9 @@ static double sloww(double x,double dx, double orig) { /* accurate enough routine calls mpsin1 or dubsin */ /***************************************************************************/ -static double sloww1(double x, double dx, double orig) { +static double +SECTION +sloww1(double x, double dx, double orig) { mynumber u; double sn,ssn,cs,ccs,s,c,w[2],y,y1,y2,c1,c2,xx,cor,res; static const double t22 = 6291456.0; @@ -797,7 +815,9 @@ static double sloww1(double x, double dx, double orig) { /* accurate enough routine calls mpsin1 or dubsin */ /***************************************************************************/ -static double sloww2(double x, double dx, double orig, int n) { +static double +SECTION +sloww2(double x, double dx, double orig, int n) { mynumber u; double sn,ssn,cs,ccs,s,c,w[2],y,y1,y2,e1,e2,xx,cor,res; static const double t22 = 6291456.0; @@ -841,7 +861,9 @@ static double sloww2(double x, double dx, double orig, int n) { /* result.And if result not accurate enough routine calls other routines */ /***************************************************************************/ -static double bsloww(double x,double dx, double orig,int n) { +static double +SECTION +bsloww(double x,double dx, double orig,int n) { static const double th2_36 = 206158430208.0; /* 1.5*2**37 */ double y,x1,x2,xx,r,t,res,cor,w[2]; #if 0 @@ -874,7 +896,9 @@ static double bsloww(double x,double dx, double orig,int n) { /* And if result not accurate enough routine calls other routines */ /***************************************************************************/ -static double bsloww1(double x, double dx, double orig,int n) { +static double +SECTION +bsloww1(double x, double dx, double orig,int n) { mynumber u; double sn,ssn,cs,ccs,s,c,w[2],y,y1,y2,c1,c2,xx,cor,res; static const double t22 = 6291456.0; @@ -917,7 +941,9 @@ mynumber u; /* And if result not accurate enough routine calls other routines */ /***************************************************************************/ -static double bsloww2(double x, double dx, double orig, int n) { +static double +SECTION +bsloww2(double x, double dx, double orig, int n) { mynumber u; double sn,ssn,cs,ccs,s,c,w[2],y,y1,y2,e1,e2,xx,cor,res; static const double t22 = 6291456.0; @@ -959,7 +985,9 @@ mynumber u; /* precision and if still doesn't accurate enough by mpcos or docos */ /************************************************************************/ -static double cslow2(double x) { +static double +SECTION +cslow2(double x) { mynumber u; double sn,ssn,cs,ccs,s,c,w[2],y,y1,y2,e1,e2,xx,cor,res; static const double t22 = 6291456.0; @@ -1002,7 +1030,9 @@ static double cslow2(double x) { /***************************************************************************/ -static double csloww(double x,double dx, double orig) { +static double +SECTION +csloww(double x,double dx, double orig) { static const double th2_36 = 206158430208.0; /* 1.5*2**37 */ double y,x1,x2,xx,r,t,res,cor,w[2],a,da,xn; union {int4 i[2]; double x;} v; @@ -1051,7 +1081,9 @@ static double csloww(double x,double dx, double orig) { /* accurate enough routine calls other routines */ /***************************************************************************/ -static double csloww1(double x, double dx, double orig) { +static double +SECTION +csloww1(double x, double dx, double orig) { mynumber u; double sn,ssn,cs,ccs,s,c,w[2],y,y1,y2,c1,c2,xx,cor,res; static const double t22 = 6291456.0; @@ -1095,7 +1127,9 @@ static double csloww1(double x, double dx, double orig) { /* accurate enough routine calls other routines */ /***************************************************************************/ -static double csloww2(double x, double dx, double orig, int n) { +static double +SECTION +csloww2(double x, double dx, double orig, int n) { mynumber u; double sn,ssn,cs,ccs,s,c,w[2],y,y1,y2,e1,e2,xx,cor,res; static const double t22 = 6291456.0; diff --git a/sysdeps/ieee754/dbl-64/s_tan.c b/sysdeps/ieee754/dbl-64/s_tan.c index 4089e0d65f..f0fcd677ae 100644 --- a/sysdeps/ieee754/dbl-64/s_tan.c +++ b/sysdeps/ieee754/dbl-64/s_tan.c @@ -41,10 +41,16 @@ #include "MathLib.h" #include "math.h" +#ifndef SECTION +# define SECTION +#endif + static double tanMp(double); void __mptan(double, mp_no *, int); -double tan(double x) { +double +SECTION +tan(double x) { #include "utan.h" #include "utan.tbl" @@ -486,7 +492,9 @@ double tan(double x) { /* multiple precision stage */ /* Convert x to multi precision number,compute tan(x) by mptan() routine */ /* and converts result back to double */ -static double tanMp(double x) +static double +SECTION +tanMp(double x) { int p; double y; diff --git a/sysdeps/ieee754/dbl-64/sincos32.c b/sysdeps/ieee754/dbl-64/sincos32.c index a4f896a465..e39aaeea06 100644 --- a/sysdeps/ieee754/dbl-64/sincos32.c +++ b/sysdeps/ieee754/dbl-64/sincos32.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001 Free Software Foundation + * Copyright (C) 2001, 2011 Free Software Foundation * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -45,11 +45,17 @@ #include "sincos32.h" #include "math_private.h" +#ifndef SECTION +# define SECTION +#endif + /****************************************************************/ /* Compute Multi-Precision sin() function for given p. Receive */ /* Multi Precision number x and result stored at y */ /****************************************************************/ -static void ss32(mp_no *x, mp_no *y, int p) { +static void +SECTION +ss32(mp_no *x, mp_no *y, int p) { int i; double a; #if 0 @@ -79,7 +85,9 @@ static void ss32(mp_no *x, mp_no *y, int p) { /* Compute Multi-Precision cos() function for given p. Receive Multi */ /* Precision number x and result stored at y */ /**********************************************************************/ -static void cc32(mp_no *x, mp_no *y, int p) { +static void +SECTION +cc32(mp_no *x, mp_no *y, int p) { int i; double a; #if 0 @@ -109,7 +117,9 @@ static void cc32(mp_no *x, mp_no *y, int p) { /***************************************************************************/ /* c32() computes both sin(x), cos(x) as Multi precision numbers */ /***************************************************************************/ -void __c32(mp_no *x, mp_no *y, mp_no *z, int p) { +void +SECTION +__c32(mp_no *x, mp_no *y, mp_no *z, int p) { static const mp_no mpt={1,{1.0,2.0}}, one={1,{1.0,1.0}}; mp_no u,t,t1,t2,c,s; int i; @@ -134,7 +144,9 @@ void __c32(mp_no *x, mp_no *y, mp_no *z, int p) { /*result which is more accurate */ /*Computing sin(x) with multi precision routine c32 */ /************************************************************************/ -double __sin32(double x, double res, double res1) { +double +SECTION +__sin32(double x, double res, double res1) { int p; mp_no a,b,c; p=32; @@ -158,7 +170,9 @@ double __sin32(double x, double res, double res1) { /*result which is more accurate */ /*Computing cos(x) with multi precision routine c32 */ /************************************************************************/ -double __cos32(double x, double res, double res1) { +double +SECTION +__cos32(double x, double res, double res1) { int p; mp_no a,b,c; p=32; @@ -172,12 +186,12 @@ double __cos32(double x, double res, double res1) { } else if (x>0.8) { __sub(&hp,&c,&a,p); - __c32(&a,&c,&b,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>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; } @@ -186,7 +200,9 @@ double __cos32(double x, double res, double res1) { /*Compute sin(x+dx) as Multi Precision number and return result as */ /* double */ /*******************************************************************/ -double __mpsin(double x, double dx) { +double +SECTION +__mpsin(double x, double dx) { int p; double y; mp_no a,b,c; @@ -204,7 +220,9 @@ double __mpsin(double x, double dx) { /* Compute cos()of double-length number (x+dx) as Multi Precision */ /* number and return result as double */ /*******************************************************************/ -double __mpcos(double x, double dx) { +double +SECTION +__mpcos(double x, double dx) { int p; double y; mp_no a,b,c; @@ -227,7 +245,9 @@ double __mpcos(double x, double dx) { /* n=0,+-1,+-2,.... */ /* Return int which indicates in which quarter of circle x is */ /******************************************************************/ -int __mpranred(double x, mp_no *y, int p) +int +SECTION +__mpranred(double x, mp_no *y, int p) { number v; double t,xn; @@ -275,7 +295,9 @@ int __mpranred(double x, mp_no *y, int p) /* Multi-Precision sin() function subroutine, for p=32. It is */ /* based on the routines mpranred() and c32(). */ /*******************************************************************/ -double __mpsin1(double x) +double +SECTION +__mpsin1(double x) { int p; int n; @@ -314,7 +336,9 @@ double __mpsin1(double x) /* based on the routines mpranred() and c32(). */ /*****************************************************************/ -double __mpcos1(double x) +double +SECTION +__mpcos1(double x) { int p; int n; diff --git a/sysdeps/ieee754/dbl-64/slowexp.c b/sysdeps/ieee754/dbl-64/slowexp.c index 78c107f709..6a6bce31ba 100644 --- a/sysdeps/ieee754/dbl-64/slowexp.c +++ b/sysdeps/ieee754/dbl-64/slowexp.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001 Free Software Foundation + * Copyright (C) 2001, 2011 Free Software Foundation * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -31,10 +31,16 @@ #include "mpa.h" #include "math_private.h" +#ifndef SECTION +# define SECTION +#endif + void __mpexp(mp_no *x, mp_no *y, int p); /*Converting from double precision to Multi-precision and calculating e^x */ -double __slowexp(double x) { +double +SECTION +__slowexp(double x) { double w,z,res,eps=3.0e-26; #if 0 double y; @@ -47,7 +53,7 @@ double __slowexp(double x) { p=6; __dbl_mp(x,&mpx,p); /* Convert a double precision number x */ - /* into a multiple precision number mpx with prec. p. */ + /* into a multiple precision number mpx with prec. p. */ __mpexp(&mpx, &mpy, p); /* Multi-Precision exponential function */ __dbl_mp(eps,&mpeps,p); __mul(&mpeps,&mpy,&mpcor,p); diff --git a/sysdeps/ieee754/dbl-64/slowpow.c b/sysdeps/ieee754/dbl-64/slowpow.c index e11a532bf8..0c57e6d4f2 100644 --- a/sysdeps/ieee754/dbl-64/slowpow.c +++ b/sysdeps/ieee754/dbl-64/slowpow.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001 Free Software Foundation + * Copyright (C) 2001, 2011 Free Software Foundation * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -35,12 +35,18 @@ #include "mpa.h" #include "math_private.h" +#ifndef SECTION +# define SECTION +#endif + void __mpexp(mp_no *x, mp_no *y, int p); void __mplog(mp_no *x, mp_no *y, int p); double ulog(double); double __halfulp(double x,double y); -double __slowpow(double x, double y, double z) { +double +SECTION +__slowpow(double x, double y, double z) { double res,res1; mp_no mpx, mpy, mpz,mpw,mpp,mpr,mpr1; static const mp_no eps = {-3,{1.0,4.0}}; @@ -48,7 +54,7 @@ double __slowpow(double x, double y, double z) { res = __halfulp(x,y); /* halfulp() returns -10 or x^y */ if (res >= 0) return res; /* if result was really computed by halfulp */ - /* else, if result was not really computed by halfulp */ + /* else, if result was not really computed by halfulp */ p = 10; /* p=precision */ __dbl_mp(x,&mpx,p); __dbl_mp(y,&mpy,p); |