glibc
2.9

Go to the source code of this file.
Functions  
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 __halfulp  (  double  x, 
double  y  
) 
Definition at line 52 of file halfulp.c.
{ 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; 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 */ 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 */ 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) {m = m<<1; l++; } n = (k&0x000fffff)0x00100000; n = n>>(20l); /* n is the odd integer of y */ k = ((k>>20) 1023)l; /* y = n*2**k */ if (k>5) return 10.0; if (k>0) for (;k>0;k) n *= 2; if (n > 34) return 10.0; k = k; if (k>5) return 10.0; /* now treat x */ while (k>0) { z = __ieee754_sqrt(x); EMULV(z,z,u,uu,j1,j2,j3,j4,j5); if (((ux)+uu) != 0) break; x = z; k; } 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]; m = k<<12; l = 0; while (m) {m = m<<1; l++; } m = (k&0x000fffff)0x00100000; m = m>>(20l); /* m is the odd integer of x */ /* now check whether the length of m**n is at most 54 bits */ if (m > tab54[n3]) 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; }
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; }
Definition at line 43 of file mplog.c.
{ #include "mplog.h" int i,m; #if 0 int j,k,m1,m2,n; double a,b; #endif static const int mp[33] = {0,0,0,0,0,1,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}; 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; /* Choose m and initiate mpone */ m = mp[p]; mpone.e = 1; mpone.d[0]=mpone.d[1]=ONE; /* Perform m newton iterations to solve for y: exp(y)x=0. */ /* The iterations formula is: y(n+1)=y(n)+(x*exp(y(n))1). */ __cpy(y,&mpt1,p); for (i=0; i<m; i++) { mpt1.d[0]=mpt1.d[0]; __mpexp(&mpt1,&mpt2,p); __mul(x,&mpt2,&mpt1,p); __sub(&mpt1,&mpone,&mpt2,p); __add(y,&mpt2,&mpt1,p); __cpy(&mpt1,y,p); } return; }
double __slowpow  (  double  x, 
double  y,  
double  z  
) 
Definition at line 44 of file slowpow.c.
{ double res, res1; long double ldw, ldz, ldpp; static const long double ldeps = 0x4.0p96; 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 */ /* Compute pow as long double, 106 bits */ ldz = __ieee754_logl ((long double) x); ldw = (long double) y *ldz; ldpp = __ieee754_expl (ldw); res = (double) (ldpp + ldeps); res1 = (double) (ldpp  ldeps); if (res != res1) /* if result still not accurate enough */ { /* use mpa for higher persision. */ mp_no mpx, mpy, mpz, mpw, mpp, mpr, mpr1; static const mp_no eps = { 3, {1.0, 4.0} }; int p; p = 10; /* p=precision 240 bits */ __dbl_mp (x, &mpx, p); __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 */ __mpexp (&mpw, &mpp, p); /* e^w =pp */ __add (&mpp, &eps, &mpr, p); /* pp+eps =r */ __mp_dbl (&mpr, &res, p); __sub (&mpp, &eps, &mpr1, p); /* pp eps =r1 */ __mp_dbl (&mpr1, &res1, p); /* converting into double precision */ if (res == res1) return res; /* if we get here result wasn't calculated exactly, continue for more exact calculation using 768 bits. */ p = 32; __dbl_mp (x, &mpx, p); __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 */ __mpexp (&mpw, &mpp, p); /* e^w=pp */ __mp_dbl (&mpp, &res, p); /* converting into double precision */ } return res; }
double ulog  (  double  ) 