Back to index

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

Go to the source code of this file.

Functions

double fastiroot (double)
void __mpsqrt (mp_no *x, mp_no *y, int p)

Function Documentation

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

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 call graph for this function:

Here is the caller graph for this function:

double fastiroot ( double  x)

Definition at line 85 of file mpsqrt.c.

                           {
  union {int i[2]; double d;} p,q;
  double y,z, t;
  int n;
  static const double c0 = 0.99674, c1 = -0.53380, c2 = 0.45472, c3 = -0.21553;

  p.d = x;
  p.i[HIGH_HALF] = (p.i[HIGH_HALF] & 0x3FFFFFFF ) | 0x3FE00000 ;
  q.d = x;
  y = p.d;
  z = y -1.0;
  n = (q.i[HIGH_HALF] - p.i[HIGH_HALF])>>1;
  z = ((c3*z + c2)*z + c1)*z + c0;            /* 2**-7         */
  z = z*(1.5 - 0.5*y*z*z);                    /* 2**-14        */
  p.d = z*(1.5 - 0.5*y*z*z);                  /* 2**-28        */
  p.i[HIGH_HALF] -= n;
  t = x*p.d;
  return p.d*(1.5 - 0.5*p.d*t);
}

Here is the caller graph for this function: