Back to index

glibc  2.9
Functions
slowpow.c File Reference
#include "mpa.h"
#include "math_private.h"

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)

Function Documentation

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>>(20-l);                       /*   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 (((u-x)+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>>(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;
}

Here is the call graph for this function:

void __mpexp ( mp_no x,
mp_no y,
int  p 
)

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=n-1; 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=n-1; 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;
}

Here is the call graph for this function:

void __mplog ( mp_no x,
mp_no y,
int  p 
)

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;
}

Here is the call graph for this function:

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.0p-96;

  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;
}

Here is the call graph for this function:

Here is the caller graph for this function:

double ulog ( double  )