Back to index

glibc  2.9
Classes | Defines | Functions
mpa.h File Reference
This graph shows which files directly or indirectly include this file:

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)

Class Documentation

struct mp_no

Definition at line 44 of file mpa.h.

Class Members
double d
int e
union number

Definition at line 57 of file mpa.h.

Class Members
double d
int i

Define Documentation

#define ABS (   x)    ((x) < 0 ? -(x) : (x))

Definition at line 66 of file mpa.h.

#define EX   x->e

Definition at line 62 of file mpa.h.

#define EY   y->e

Definition at line 63 of file mpa.h.

#define EZ   z->e

Definition at line 64 of file mpa.h.

#define X   x->d

Definition at line 59 of file mpa.h.

#define Y   y->d

Definition at line 60 of file mpa.h.

#define Z   z->d

Definition at line 61 of file mpa.h.


Function Documentation

int __acr ( const mp_no ,
const 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;
}

Here is the call graph for this function:

Here is the caller graph for this function:

void __add ( const mp_no ,
const mp_no ,
mp_no ,
int   
)

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

Here is the call graph for this function:

Here is the caller graph for this function:

void __c32 ( mp_no ,
mp_no ,
mp_no ,
int   
)

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.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);
  }
  __sub(&one,&c,y,p);
  __cpy(&s,z,p);
}

Here is the caller graph for this function:

void __cpy ( const mp_no ,
mp_no ,
int   
)

Definition at line 99 of file mpa.c.

                                            {
  int i;

  EY = EX;
  for (i=0; i <= p; i++)    Y[i] = X[i];

  return;
}

Here is the caller graph for this function:

void __cpymn ( const mp_no ,
int  ,
mp_no ,
int   
)

Definition at line 115 of file mpa.c.

                                                     {

  int i,k;

  EY = EX;     k=MIN(m,n);
  for (i=0; i <= k; i++)    Y[i] = X[i];
  for (   ; i <= n; i++)    Y[i] = ZERO;

  return;
}
int __cr ( const mp_no ,
const mp_no ,
int   
)

Definition at line 86 of file mpa.c.

                                                 {
  int i;

  if      (X[0] > Y[0])  i= 1;
  else if (X[0] < Y[0])  i=-1;
  else if (X[0] < ZERO ) i= __acr(y,x,p);
  else                   i= __acr(x,y,p);

  return i;
}

Here is the call graph for this function:

void __dbl_mp ( double  ,
mp_no ,
int   
)

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

Here is the caller graph for this function:

void __dvd ( const mp_no ,
const mp_no ,
mp_no ,
int   
)

Definition at line 501 of file mpa.c.

                                                            {

  mp_no w;

  if (X[0] == ZERO)    Z[0] = ZERO;
  else                {__inv(y,&w,p);   __mul(x,&w,z,p);}
  return;
}

Here is the call graph for this function:

Here is the caller graph for this function:

void __inv ( const mp_no ,
mp_no ,
int   
)

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

Here is the call graph for this function:

Here is the caller graph for this function:

void __mp_dbl ( const mp_no ,
double *  ,
int   
)

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

Here is the call graph for this function:

Here is the caller graph for this function:

void __mpatan ( mp_no ,
mp_no ,
int   
)

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

Here is the caller graph for this function:

void __mpatan2 ( mp_no ,
mp_no ,
mp_no ,
int   
)

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

Here is the call graph for this function:

void __mpexp ( mp_no ,
mp_no __y,
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;
}

Here is the caller graph for this function:

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.e-5;
    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<=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);
    }
    else __mul(&c,&hp,y,p);
    n = (int) t;
    if (x < 0) { y->d[0] = - y->d[0]; n = -n; }
    return (n&3);
  }
}

Here is the caller graph for this function:

void __mpsqrt ( mp_no ,
mp_no ,
int   
)

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

Here is the caller graph for this function:

void __mul ( const mp_no ,
const mp_no ,
mp_no ,
int   
)

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=k-p; i2=p+1; }
    else        {i1=1;   i2=k;   }
    for (i=i1,j=i2-1; 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;
}

Here is the caller graph for this function:

void __sub ( const mp_no ,
const mp_no ,
mp_no ,
int   
)

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

Here is the call graph for this function:

Here is the caller graph for this function: