glibc
2.9

Go to the source code of this file.
Classes  
struct  mp_no 
union  number 
Defines  
#define  X x>d 
#define  Y y>d 
#define  Z z>d 
#define  EX x>e 
#define  EY y>e 
#define  EZ z>e 
#define  ABS(x) ((x) < 0 ? (x) : (x)) 
Functions  
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) 
void  __mp_dbl (const mp_no *, double *, int) 
void  __dbl_mp (double, mp_no *, int) 
void  __add (const mp_no *, const mp_no *, mp_no *, int) 
void  __sub (const mp_no *, const mp_no *, mp_no *, int) 
void  __mul (const mp_no *, const mp_no *, mp_no *, int) 
void  __inv (const mp_no *, mp_no *, int) 
void  __dvd (const mp_no *, const mp_no *, mp_no *, int) 
void  __mpatan (mp_no *, mp_no *, int) 
void  __mpatan2 (mp_no *, mp_no *, mp_no *, int) 
void  __mpsqrt (mp_no *, mp_no *, int) 
void  __mpexp (mp_no *, mp_no *__y, int) 
void  __c32 (mp_no *, mp_no *, mp_no *, int) 
int  __mpranred (double, mp_no *, int) 
Definition at line 67 of file mpa.c.
{ int i; if (X[0] == ZERO) { if (Y[0] == ZERO) i= 0; else i=1; } else if (Y[0] == ZERO) i= 1; else { if (EX > EY) i= 1; else if (EX < EY) i=1; else i= mcr(x,y,p); } return i; }
Definition at line 380 of file mpa.c.
{ int n; if (X[0] == ZERO) {__cpy(y,z,p); return; } else if (Y[0] == ZERO) {__cpy(x,z,p); return; } if (X[0] == Y[0]) { if (__acr(x,y,p) > 0) {add_magnitudes(x,y,z,p); Z[0] = X[0]; } else {add_magnitudes(y,x,z,p); Z[0] = Y[0]; } } else { if ((n=__acr(x,y,p)) == 1) {sub_magnitudes(x,y,z,p); Z[0] = X[0]; } else if (n == 1) {sub_magnitudes(y,x,z,p); Z[0] = Y[0]; } else Z[0] = ZERO; } return; }
Definition at line 112 of file sincos32.c.
{ 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); u.e=u.e1; 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); } __sub(&one,&c,y,p); __cpy(&s,z,p); }
Definition at line 251 of file mpa.c.
{ int i,n; double u; /* Sign */ if (x == ZERO) {Y[0] = ZERO; return; } else if (x > ZERO) Y[0] = ONE; else {Y[0] = MONE; x=x; } /* Exponent */ for (EY=ONE; x >= RADIX; EY += ONE) x *= RADIXI; for ( ; x < ONE; EY = ONE) x *= RADIX; /* Digits */ n=MIN(p,4); for (i=1; i<=n; i++) { u = (x + TWO52)  TWO52; if (u>x) u = ONE; Y[i] = u; x = u; x *= RADIX; } for ( ; i<=p; i++) Y[i] = ZERO; return; }
Definition at line 469 of file mpa.c.
{ int i; #if 0 int l; #endif double t; mp_no z,w; static const int np1[] = {0,0,0,0,1,2,2,2,2,3,3,3,3,3,3,3,3,3, 4,4,4,4,4,4,4,4,4,4,4,4,4,4,4}; const mp_no mptwo = {1,{1.0,2.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}}; __cpy(x,&z,p); z.e=0; __mp_dbl(&z,&t,p); t=ONE/t; __dbl_mp(t,y,p); EY = EX; for (i=0; i<np1[p]; i++) { __cpy(y,&w,p); __mul(x,&w,y,p); __sub(&mptwo,y,&z,p); __mul(&w,&z,y,p); } return; }
Definition at line 233 of file mpa.c.
{ #if 0 int i,k; double a,c,u,v,z[5]; #endif if (X[0] == ZERO) {*y = ZERO; return; } if (EX> 42) norm(x,y,p); else if (EX==42 && X[1]>=TWO10) norm(x,y,p); else denorm(x,y,p); }
Definition at line 39 of file mpatan.c.
{ #include "mpatan.h" int i,m,n; double dx; 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}}, mptwo = {0,{0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}}, mptwoim1 = {0,{0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, 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,mpsm,mpt,mpt1,mpt2,mpt3; /* 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;} } 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 */ __mul(x,x,&mpsm,p); if (m==0) __cpy(x,&mps,p); else { for (i=0; i<m; i++) { __add(&mpone,&mpsm,&mpt1,p); __mpsqrt(&mpt1,&mpt2,p); __add(&mpt2,&mpt2,&mpt1,p); __add(&mptwo,&mpsm,&mpt2,p); __add(&mpt1,&mpt2,&mpt3,p); __dvd(&mpsm,&mpt3,&mpt1,p); __cpy(&mpt1,&mpsm,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; __dvd(&mpsm,&mptwoim1,&mpt,p); for (i=n1; i>1; i) { mptwoim1.d[1] = TWO; __dvd(&mpsm,&mptwoim1,&mpt1,p); __mul(&mpsm,&mpt,&mpt2,p); __sub(&mpt1,&mpt2,&mpt,p); } __mul(&mps,&mpt,&mpt1,p); __sub(&mps,&mpt1,&mpt,p); /* Compute Atan(x) */ mptwoim1.d[1] = twom[m].d; __mul(&mptwoim1,&mpt,y,p); return; }
Definition at line 46 of file mpatan2.c.
{ 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}}; mp_no mpt1,mpt2,mpt3; if (X[0] <= ZERO) { mpone.e = 1; mpone.d[0] = mpone.d[1] = ONE; __dvd(x,y,&mpt1,p); __mul(&mpt1,&mpt1,&mpt2,p); if (mpt1.d[0] != ZERO) mpt1.d[0] = ONE; __add(&mpt2,&mpone,&mpt3,p); __mpsqrt(&mpt3,&mpt2,p); __add(&mpt1,&mpt2,&mpt3,p); mpt3.d[0]=Y[0]; __mpatan(&mpt3,&mpt1,p); __add(&mpt1,&mpt1,z,p); } else { __dvd(y,x,&mpt1,p); __mpatan(&mpt1,z,p); } return; }
Definition at line 39 of file mpexp.c.
{ 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}; 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}; 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}}; 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}}; 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}}; mp_no mps,mpak,mpt1,mpt2; /* Choose m,n and compute a=2**(m) */ n = np[p]; m1 = m1p[p]; a = twomm1[p].d; for (i=0; i<EX; i++) a *= RADIXI; for ( ; i>EX; i) a *= RADIX; b = X[1]*RADIXI; m2 = 24*EX; for (; b<HALF; m2) { a *= TWO; b *= TWO; } if (b == HALF) { for (i=2; i<=p; i++) { if (X[i]!=ZERO) break; } if (i==p+1) { m2; a *= TWO; } } if ((m=m1+m2) <= 0) { m=0; a=ONE; for (i=n1; i>0; i,n) { if (m1np[i][p]+m2>0) break; } } /* Compute s=x*2**(m). Put result in mps */ __dbl_mp(a,&mpt1,p); __mul(x,&mpt1,&mps,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; __dvd(&mps,&mpk,&mpt1,p); __add(&mpone,&mpt1,&mpak,p); for (k=n1; k>1; k) { __mul(&mps,&mpak,&mpt1,p); mpk.d[1]=nn[k].d; __dvd(&mpt1,&mpk,&mpt2,p); __add(&mpone,&mpt2,&mpak,p); } __mul(&mps,&mpak,&mpt1,p); __add(&mpone,&mpt1,&mpt2,p); /* Raise polynomial value to the power of 2**m. Put result in y */ for (k=0,j=0; k<m; ) { __mul(&mpt2,&mpt2,&mpt1,p); k++; if (k==m) { j=1; break; } __mul(&mpt1,&mpt1,&mpt2,p); k++; } if (j) __cpy(&mpt1,y,p); else __cpy(&mpt2,y,p); return; }
int __mpranred  (  double  , 
mp_no *  ,  
int  
) 
Definition at line 230 of file sincos32.c.
{ number v; double t,xn; int i,k,n; static const mp_no one = {1,{1.0,1.0}}; 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.e5; 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<=pc.e;i++) c.d[i]=c.d[i+c.e]; for (i=p+1c.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); } else __mul(&c,&hp,y,p); n = (int) t; if (x < 0) { y>d[0] =  y>d[0]; n = n; } return (n&3); } }
Definition at line 46 of file mpsqrt.c.
{ #include "mpsqrt.h" int i,m,ex,ey; double dx,dy; mp_no mphalf = {0,{0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}}, mp3halfs = {0,{0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, 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 mpxn,mpz,mpu,mpt1,mpt2; /* Prepare multiprecision 1/2 and 3/2 */ mphalf.e =0; mphalf.d[0] =ONE; mphalf.d[1] =HALFRAD; mp3halfs.e=1; mp3halfs.d[0]=ONE; mp3halfs.d[1]=ONE; mp3halfs.d[2]=HALFRAD; ex=EX; ey=EX/2; __cpy(x,&mpxn,p); mpxn.e = (ey+ey); __mp_dbl(&mpxn,&dx,p); dy=fastiroot(dx); __dbl_mp(dy,&mpu,p); __mul(&mpxn,&mphalf,&mpz,p); m=mp[p]; for (i=0; i<m; i++) { __mul(&mpu,&mpu,&mpt1,p); __mul(&mpt1,&mpz,&mpt2,p); __sub(&mp3halfs,&mpt2,&mpt1,p); __mul(&mpu,&mpt1,&mpt2,p); __cpy(&mpt2,&mpu,p); } __mul(&mpxn,&mpu,y,p); EY += ey; return; }
Definition at line 429 of file mpa.c.
{ int i, i1, i2, j, k, k2; double u; /* Is z=0? */ if (X[0]*Y[0]==ZERO) { Z[0]=ZERO; return; } /* Multiply, add and carry */ k2 = (p<3) ? p+p : p+3; Z[k2]=ZERO; for (k=k2; k>1; ) { if (k > p) {i1=kp; i2=p+1; } else {i1=1; i2=k; } for (i=i1,j=i21; i<i2; i++,j) Z[k] += X[i]*Y[j]; u = (Z[k] + CUTTER)CUTTER; if (u > Z[k]) u = RADIX; Z[k] = u; Z[k] = u*RADIXI; } /* Is there a carry beyond the most significant digit? */ if (Z[1] == ZERO) { for (i=1; i<=p; i++) Z[i]=Z[i+1]; EZ = EX + EY  1; } else EZ = EX + EY; Z[0] = X[0] * Y[0]; return; }
Definition at line 404 of file mpa.c.
{ int n; if (X[0] == ZERO) {__cpy(y,z,p); Z[0] = Z[0]; return; } else if (Y[0] == ZERO) {__cpy(x,z,p); return; } if (X[0] != Y[0]) { if (__acr(x,y,p) > 0) {add_magnitudes(x,y,z,p); Z[0] = X[0]; } else {add_magnitudes(y,x,z,p); Z[0] = Y[0]; } } else { if ((n=__acr(x,y,p)) == 1) {sub_magnitudes(x,y,z,p); Z[0] = X[0]; } else if (n == 1) {sub_magnitudes(y,x,z,p); Z[0] = Y[0]; } else Z[0] = ZERO; } return; }