Back to index

glibc  2.9
Functions
mplog.c File Reference
#include "endian.h"
#include "mpa.h"
#include "mplog.h"

Go to the source code of this file.

Functions

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

Function Documentation

void __mpexp ( mp_no ,
mp_no ,
int   
)

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

Here is the caller graph for this function: