Back to index

glibc  2.9
Defines | Functions
e_log.c File Reference
#include "endian.h"
#include "dla.h"
#include "mpa.h"
#include "MathLib.h"
#include "math_private.h"
#include "ulog.tbl"
#include "ulog.h"

Go to the source code of this file.

Defines

#define M   4

Functions

void __mplog (mp_no *, mp_no *, int)
double __ieee754_log (double x)

Define Documentation

#define M   4

Function Documentation

double __ieee754_log ( double  x)

Definition at line 50 of file e_log.c.

                               {
#define M 4
  static const int pr[M]={8,10,18,32};
  int i,j,n,ux,dx,p;
#if 0
  int k;
#endif
  double dbl_n,u,p0,q,r0,w,nln2a,luai,lubi,lvaj,lvbj,
         sij,ssij,ttij,A,B,B0,y,y1,y2,polI,polII,sa,sb,
         t1,t2,t3,t4,t5,t6,t7,t8,t,ra,rb,ww,
         a0,aa0,s1,s2,ss2,s3,ss3,a1,aa1,a,aa,b,bb,c;
  number num;
  mp_no mpx,mpy,mpy1,mpy2,mperr;

#include "ulog.tbl"
#include "ulog.h"

  /* Treating special values of x ( x<=0, x=INF, x=NaN etc.). */

  num.d = x;  ux = num.i[HIGH_HALF];  dx = num.i[LOW_HALF];
  n=0;
  if (ux < 0x00100000) {
    if (((ux & 0x7fffffff) | dx) == 0)  return MHALF/ZERO; /* return -INF */
    if (ux < 0) return (x-x)/ZERO;                         /* return NaN  */
    n -= 54;    x *= two54.d;                              /* scale x     */
    num.d = x;
  }
  if (ux >= 0x7ff00000) return x+x;                        /* INF or NaN  */

  /* Regular values of x */

  w = x-ONE;
  if (ABS(w) > U03) { goto case_03; }


  /*--- Stage I, the case abs(x-1) < 0.03 */

  t8 = MHALF*w;
  EMULV(t8,w,a,aa,t1,t2,t3,t4,t5)
  EADD(w,a,b,bb)

  /* Evaluate polynomial II */
  polII = (b0.d+w*(b1.d+w*(b2.d+w*(b3.d+w*(b4.d+
          w*(b5.d+w*(b6.d+w*(b7.d+w*b8.d))))))))*w*w*w;
  c = (aa+bb)+polII;

  /* End stage I, case abs(x-1) < 0.03 */
  if ((y=b+(c+b*E2)) == b+(c-b*E2))  return y;

  /*--- Stage II, the case abs(x-1) < 0.03 */

  a = d11.d+w*(d12.d+w*(d13.d+w*(d14.d+w*(d15.d+w*(d16.d+
            w*(d17.d+w*(d18.d+w*(d19.d+w*d20.d))))))));
  EMULV(w,a,s2,ss2,t1,t2,t3,t4,t5)
  ADD2(d10.d,dd10.d,s2,ss2,s3,ss3,t1,t2)
  MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
  ADD2(d9.d,dd9.d,s2,ss2,s3,ss3,t1,t2)
  MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
  ADD2(d8.d,dd8.d,s2,ss2,s3,ss3,t1,t2)
  MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
  ADD2(d7.d,dd7.d,s2,ss2,s3,ss3,t1,t2)
  MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
  ADD2(d6.d,dd6.d,s2,ss2,s3,ss3,t1,t2)
  MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
  ADD2(d5.d,dd5.d,s2,ss2,s3,ss3,t1,t2)
  MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
  ADD2(d4.d,dd4.d,s2,ss2,s3,ss3,t1,t2)
  MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
  ADD2(d3.d,dd3.d,s2,ss2,s3,ss3,t1,t2)
  MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
  ADD2(d2.d,dd2.d,s2,ss2,s3,ss3,t1,t2)
  MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
  MUL2(w,ZERO,s2,ss2,s3,ss3,t1,t2,t3,t4,t5,t6,t7,t8)
  ADD2(w,ZERO,    s3,ss3, b, bb,t1,t2)

  /* End stage II, case abs(x-1) < 0.03 */
  if ((y=b+(bb+b*E4)) == b+(bb-b*E4))  return y;
  goto stage_n;

  /*--- Stage I, the case abs(x-1) > 0.03 */
  case_03:

  /* Find n,u such that x = u*2**n,   1/sqrt(2) < u < sqrt(2)  */
  n += (num.i[HIGH_HALF] >> 20) - 1023;
  num.i[HIGH_HALF] = (num.i[HIGH_HALF] & 0x000fffff) | 0x3ff00000;
  if (num.d > SQRT_2) { num.d *= HALF;  n++; }
  u = num.d;  dbl_n = (double) n;

  /* Find i such that ui=1+(i-75)/2**8 is closest to u (i= 0,1,2,...,181) */
  num.d += h1.d;
  i = (num.i[HIGH_HALF] & 0x000fffff) >> 12;

  /* Find j such that vj=1+(j-180)/2**16 is closest to v=u/ui (j= 0,...,361) */
  num.d = u*Iu[i].d + h2.d;
  j = (num.i[HIGH_HALF] & 0x000fffff) >> 4;

  /* Compute w=(u-ui*vj)/(ui*vj) */
  p0=(ONE+(i-75)*DEL_U)*(ONE+(j-180)*DEL_V);
  q=u-p0;   r0=Iu[i].d*Iv[j].d;   w=q*r0;

  /* Evaluate polynomial I */
  polI = w+(a2.d+a3.d*w)*w*w;

  /* Add up everything */
  nln2a = dbl_n*LN2A;
  luai  = Lu[i][0].d;   lubi  = Lu[i][1].d;
  lvaj  = Lv[j][0].d;   lvbj  = Lv[j][1].d;
  EADD(luai,lvaj,sij,ssij)
  EADD(nln2a,sij,A  ,ttij)
  B0 = (((lubi+lvbj)+ssij)+ttij)+dbl_n*LN2B;
  B  = polI+B0;

  /* End stage I, case abs(x-1) >= 0.03 */
  if ((y=A+(B+E1)) == A+(B-E1))  return y;


  /*--- Stage II, the case abs(x-1) > 0.03 */

  /* Improve the accuracy of r0 */
  EMULV(p0,r0,sa,sb,t1,t2,t3,t4,t5)
  t=r0*((ONE-sa)-sb);
  EADD(r0,t,ra,rb)

  /* Compute w */
  MUL2(q,ZERO,ra,rb,w,ww,t1,t2,t3,t4,t5,t6,t7,t8)

  EADD(A,B0,a0,aa0)

  /* Evaluate polynomial III */
  s1 = (c3.d+(c4.d+c5.d*w)*w)*w;
  EADD(c2.d,s1,s2,ss2)
  MUL2(s2,ss2,w,ww,s3,ss3,t1,t2,t3,t4,t5,t6,t7,t8)
  MUL2(s3,ss3,w,ww,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
  ADD2(s2,ss2,w,ww,s3,ss3,t1,t2)
  ADD2(s3,ss3,a0,aa0,a1,aa1,t1,t2)

  /* End stage II, case abs(x-1) >= 0.03 */
  if ((y=a1+(aa1+E3)) == a1+(aa1-E3)) return y;


  /* Final stages. Use multi-precision arithmetic. */
  stage_n:

  for (i=0; i<M; i++) {
    p = pr[i];
    __dbl_mp(x,&mpx,p);  __dbl_mp(y,&mpy,p);
    __mplog(&mpx,&mpy,p);
    __dbl_mp(e[i].d,&mperr,p);
    __add(&mpy,&mperr,&mpy1,p);  __sub(&mpy,&mperr,&mpy2,p);
    __mp_dbl(&mpy1,&y1,p);       __mp_dbl(&mpy2,&y2,p);
    if (y1==y2)   return y1;
  }
  return y1;
}

Here is the call graph for this function:

Here is the caller graph for this function:

void __mplog ( mp_no ,
mp_no ,
int   
)

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