/* 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 __branred(double x, double *a, double *aa)
{
int i,k;
#if 0
*
* 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
- * the Free Software Foundation; either version 2 of the License, or
+ * the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
- *
+ *
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
*
* You should have received a copy of the GNU Lesser General Public License
* along with this program; if not, write to the Free Software
- * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+ * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
*/
/************************************************************************/
/* */
-/* MODULE_NAME:halfulp.c */
-/* */
+/* MODULE_NAME:halfulp.c */
+/* */
/* FUNCTIONS:halfulp */
/* FILES NEEDED: mydefs.h dla.h endian.h */
/* uroot.c */
/*3. if x can be represented by x=2**n for some integer n. */
/************************************************************************/
-#include "endian.h"
+#include "endian.h"
#include "mydefs.h"
#include "dla.h"
3, 3, 3, 3, 3, 3, 3, 3 };
-double halfulp(double x, double y)
+double __halfulp(double x, double y)
{
mynumber v;
double z,u,uu,j1,j2,j3,j4,j5;
int4 k,l,m,n;
if (y <= 0) { /*if power is negative or zero */
v.x = y;
- if (v.i[LOW_HALF] != 0) return -10.0;
+ if (v.i[LOW_HALF] != 0) return -10.0;
v.x = x;
- if (v.i[LOW_HALF] != 0) return -10.0;
- if ((v.i[HIGH_HALF]&0x000fffff) != 0) return -10; /* if x =2 ^ n */
+ if (v.i[LOW_HALF] != 0) return -10.0;
+ if ((v.i[HIGH_HALF]&0x000fffff) != 0) return -10; /* if x =2 ^ n */
k = ((v.i[HIGH_HALF]&0x7fffffff)>>20)-1023; /* find this n */
z = (double) k;
return (z*y == -1075.0)?0: -10.0;
/* if y > 0 */
v.x = y;
if (v.i[LOW_HALF] != 0) return -10.0;
-
+
v.x=x;
- /* case where x = 2**n for some integer n */
+ /* case where x = 2**n for some integer n */
if (((v.i[HIGH_HALF]&0x000fffff)|v.i[LOW_HALF]) == 0) {
k=(v.i[HIGH_HALF]>>20)-1023;
return (((double) k)*y == -1075.0)?0:-10.0;
- }
-
+ }
+
v.x = y;
k = v.i[HIGH_HALF];
m = k<<12;
l = 0;
- while (m)
+ while (m)
{m = m<<1; l++; }
n = (k&0x000fffff)|0x00100000;
n = n>>(20-l); /* n is the odd integer of y */
if (n > 34) return -10.0;
k = -k;
if (k>5) return -10.0;
-
+
/* now treat x */
while (k>0) {
- z = usqrt(x);
+ z = __ieee754_sqrt(x);
EMULV(z,z,u,uu,j1,j2,j3,j4,j5);
if (((u-x)+uu) != 0) break;
x = z;
k--;
}
- if (k) return -10.0;
-
+ if (k) return -10.0;
+
/* it is impossible that n == 2, so the mantissa of x must be short */
-
+
v.x = x;
if (v.i[LOW_HALF]) return -10.0;
k = v.i[HIGH_HALF];
while (m) {m = m<<1; l++; }
m = (k&0x000fffff)|0x00100000;
m = m>>(20-l); /* m is the odd integer of x */
-
+
/* now check whether the length of m**n is at most 54 bits */
-
+
if (m > tab54[n-3]) return -10.0;
-
+
/* yes, it is - now compute x**n by simple multiplications */
-
+
u = x;
for (k=1;k<n;k++) u = u*x;
return u;
}
-
-
#include "endian.h"
#include "mpa.h"
-int mpranred(double, mp_no *, int);
+int __mpranred(double, mp_no *, int);
void __c32(mp_no *, mp_no *, mp_no *, int);
void __mptan(double x, mp_no *mpy, int p) {
int n;
mp_no mpw, mpc, mps;
- n = mpranred(x, &mpw, p) & 0x00000001; /* negative or positive result */
+ n = __mpranred(x, &mpw, p) & 0x00000001; /* negative or positive result */
__c32(&mpw, &mpc, &mps, p); /* computing sin(x) and cos(x) */
if (n) /* second or fourth quarter of unit circle */
- { dvd(&mpc,&mps,mpy,p);
+ { __dvd(&mpc,&mps,mpy,p);
mpy->d[0] *= MONE;
} /* tan is negative in this area */
- else dvd(&mps,&mpc,mpy,p);
+ else __dvd(&mps,&mpc,mpy,p);
return;
}
static double bsloww(double x, double dx, double orig, int n);
static double bsloww1(double x, double dx, double orig, int n);
static double bsloww2(double x, double dx, double orig, int n);
-int branred(double x, double *a, double *aa);
+int __branred(double x, double *a, double *aa);
static double cslow2(double x);
static double csloww(double x, double dx, double orig);
static double csloww1(double x, double dx, double orig);
/* -----------------281474976710656 <|x| <2^1024----------------------------*/
else if (k < 0x7ff00000) {
- n = branred(x,&a,&da);
+ n = __branred(x,&a,&da);
switch (n) {
case 0:
if (a*a < 0.01588) return bsloww(a,da,x,n);
else if (k < 0x7ff00000) {/* 281474976710656 <|x| <2^1024 */
- n = branred(x,&a,&da);
+ n = __branred(x,&a,&da);
switch (n) {
case 1:
if (a*a < 0.01588) return bsloww(-a,-da,x,n);
mp_no mpy;
#endif
- int branred(double, double *, double *);
- int mpranred(double, mp_no *, int);
+ int __branred(double, double *, double *);
+ int __mpranred(double, mp_no *, int);
/* x=+-INF, x=NaN */
num.d = x; ux = num.i[HIGH_HALF];
/* (---) The case 1e8 < abs(x) < 2**1024 */
/* Range reduction by algorithm iii */
- n = (branred(x,&a,&da)) & 0x00000001;
+ n = (__branred(x,&a,&da)) & 0x00000001;
EADD(a,da,t1,t2) a=t1; da=t2;
if (a<ZERO) {ya=-a; yya=-da; sy=MONE;}
else {ya= a; yya= da; sy= ONE;}
/* Second stage */
/* Reduction by algorithm iv */
- p=10; n = (mpranred(x,&mpa,p)) & 0x00000001;
+ p=10; n = (__mpranred(x,&mpa,p)) & 0x00000001;
__mp_dbl(&mpa,&a,p); __dbl_mp(a,&mpt1,p);
__sub(&mpa,&mpt1,&mpt2,p); __mp_dbl(&mpt2,&da,p);
#endif
for (i=1;i<=p;i++) mpk.d[i]=0;
- mul(x,x,&x2,p);
- cpy(&oofac27,&gor,p);
- cpy(&gor,&sum,p);
+ __mul(x,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(&gor,&mpk,&mpt1,p);
+ __cpy(&mpt1,&gor,p);
+ __mul(&x2,&sum,&mpt1,p);
+ __sub(&gor,&mpt1,&sum,p);
}
- mul(x,&sum,y,p);
+ __mul(x,&sum,y,p);
}
/**********************************************************************/
#endif
for (i=1;i<=p;i++) mpk.d[i]=0;
- mul(x,x,&x2,p);
+ __mul(x,x,&x2,p);
mpk.d[1]=27.0;
- mul(&oofac27,&mpk,&gor,p);
- cpy(&gor,&sum,p);
+ __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(&gor,&mpk,&mpt1,p);
+ __cpy(&mpt1,&gor,p);
+ __mul(&x2,&sum,&mpt1,p);
+ __sub(&gor,&mpt1,&sum,p);
}
- mul(&x2,&sum,y,p);
+ __mul(&x2,&sum,y,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;
- cpy(x,&u,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(&mpt,&c,&t1,p);
- mul(&t1,&c,&t2,p);
- add(&t2,&t2,&c,p);
+ __mul(&c,&s,&t,p);
+ __sub(&s,&t,&t1,p);
+ __add(&t1,&t1,&s,p);
+ __sub(&mpt,&c,&t1,p);
+ __mul(&t1,&c,&t2,p);
+ __add(&t2,&t2,&c,p);
}
- sub(&one,&c,y,p);
- cpy(&s,z,p);
+ __sub(&one,&c,y,p);
+ __cpy(&s,z,p);
}
/************************************************************************/
p=32;
__dbl_mp(res,&a,p);
__dbl_mp(0.5*(res1-res),&b,p);
- add(&a,&b,&c,p);
+ __add(&a,&b,&c,p);
if (x>0.8)
- { sub(&hp,&c,&a,p);
+ { __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);
+ __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;
p=32;
__dbl_mp(res,&a,p);
__dbl_mp(0.5*(res1-res),&b,p);
- add(&a,&b,&c,p);
+ __add(&a,&b,&c,p);
if (x>2.4)
- { sub(&pi,&c,&a,p);
+ { __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);
+ { __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);
+ __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;
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); }
+ __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;
p=32;
__dbl_mp(x,&a,p);
__dbl_mp(dx,&b,p);
- add(&a,&b,&c,p);
+ __add(&a,&b,&c,p);
if (x>0.8)
- { sub(&hp,&c,&b,p);
+ { __sub(&hp,&c,&b,p);
__c32(&b,&a,&c,p);
}
else __c32(&c,&a,&b,p); /* a = cos(x+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 __mpranred(double x, mp_no *y, int p)
{
number v;
double t,xn;
v.d = t;
n =v.i[LOW_HALF]&3;
__dbl_mp(xn,&a,p);
- mul(&a,&hp,&b,p);
+ __mul(&a,&hp,&b,p);
__dbl_mp(x,&c,p);
- sub(&c,&b,y,p);
+ __sub(&c,&b,y,p);
return n;
}
else { /* if x is very big more precision required */
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);
+ __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] >= 8388608.0)
{ t +=1.0;
- sub(&c,&one,&b,p);
- mul(&b,&hp,y,p);
+ __sub(&c,&one,&b,p);
+ __mul(&b,&hp,y,p);
}
- else mul(&c,&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);
mp_no u,s,c;
double y;
p=32;
- n=mpranred(x,&u,p); /* n is 0, 1, 2 or 3 */
+ 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:
double y;
p=32;
- n=mpranred(x,&u,p); /* n is 0, 1, 2 or 3 */
+ 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*/
/* 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);
- add(&mpy,&mpcor,&mpw,p);
- sub(&mpy,&mpcor,&mpz,p);
+ __mul(&mpeps,&mpy,&mpcor,p);
+ __add(&mpy,&mpcor,&mpw,p);
+ __sub(&mpy,&mpcor,&mpz,p);
__mp_dbl(&mpw, &w, p);
__mp_dbl(&mpz, &z, p);
if (w == z) return w;
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 __halfulp(double x,double y);
double slowpow(double x, double y, double z) {
double res,res1;
static const mp_no eps = {-3,{1.0,4.0}};
int p;
- res = halfulp(x,y); /* halfulp() returns -10 or x^y */
+ 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 */
p = 10; /* p=precision */
__dbl_mp(y,&mpy,p);
__dbl_mp(z,&mpz,p);
__mplog(&mpx, &mpz, p); /* log(x) = z */
- mul(&mpy,&mpz,&mpw,p); /* y * z =w */
+ __mul(&mpy,&mpz,&mpw,p); /* y * z =w */
__mpexp(&mpw, &mpp, p); /* e^w =pp */
- add(&mpp,&eps,&mpr,p); /* pp+eps =r */
+ __add(&mpp,&eps,&mpr,p); /* pp+eps =r */
__mp_dbl(&mpr, &res, p);
- sub(&mpp,&eps,&mpr1,p); /* pp -eps =r1 */
+ __sub(&mpp,&eps,&mpr1,p); /* pp -eps =r1 */
__mp_dbl(&mpr1, &res1, p); /* converting into double precision */
if (res == res1) return res;
__dbl_mp(y,&mpy,p);
__dbl_mp(z,&mpz,p);
__mplog(&mpx, &mpz, p); /* log(c)=z */
- mul(&mpy,&mpz,&mpw,p); /* y*z =w */
+ __mul(&mpy,&mpz,&mpw,p); /* y*z =w */
__mpexp(&mpw, &mpp, p); /* e^w=pp */
__mp_dbl(&mpp, &res, p); /* converting into double precision */
return res;
projects / kopensolaris-gnu / glibc.git / commitdiff