4.3BSD - /usr/src/usr.lib/libm/log1p.c

   1: /*
   2:  * Copyright (c) 1985 Regents of the University of California.
   3:  *
   4:  * Use and reproduction of this software are granted  in  accordance  with
   5:  * the terms and conditions specified in  the  Berkeley  Software  License
   6:  * Agreement (in particular, this entails acknowledgement of the programs'
   7:  * source, and inclusion of this notice) with the additional understanding
   8:  * that  all  recipients  should regard themselves as participants  in  an
   9:  * ongoing  research  project and hence should  feel  obligated  to report
  10:  * their  experiences (good or bad) with these elementary function  codes,
  11:  * using "sendbug 4bsd-bugs@BERKELEY", to the authors.
  12:  */
  13: 
  14: #ifndef lint
  15: static char sccsid[] = "@(#)log1p.c	1.3 (Berkeley) 8/21/85";
  16: #endif not lint
  17: 
  18: /* LOG1P(x)
  19:  * RETURN THE LOGARITHM OF 1+x
  20:  * DOUBLE PRECISION (VAX D FORMAT 56 bits, IEEE DOUBLE 53 BITS)
  21:  * CODED IN C BY K.C. NG, 1/19/85;
  22:  * REVISED BY K.C. NG on 2/6/85, 3/7/85, 3/24/85, 4/16/85.
  23:  *
  24:  * Required system supported functions:
  25:  *	scalb(x,n)
  26:  *	copysign(x,y)
  27:  *	logb(x)
  28:  *	finite(x)
  29:  *
  30:  * Required kernel function:
  31:  *	log__L(z)
  32:  *
  33:  * Method :
  34:  *	1. Argument Reduction: find k and f such that
  35:  *			1+x  = 2^k * (1+f),
  36:  *	   where  sqrt(2)/2 < 1+f < sqrt(2) .
  37:  *
  38:  *	2. Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
  39:  *		 = 2s + 2/3 s**3 + 2/5 s**5 + .....,
  40:  *	   log(1+f) is computed by
  41:  *
  42:  *	     		log(1+f) = 2s + s*log__L(s*s)
  43:  *	   where
  44:  *		log__L(z) = z*(L1 + z*(L2 + z*(... (L6 + z*L7)...)))
  45:  *
  46:  *	   See log__L() for the values of the coefficients.
  47:  *
  48:  *	3. Finally,  log(1+x) = k*ln2 + log(1+f).
  49:  *
  50:  *	Remarks 1. In step 3 n*ln2 will be stored in two floating point numbers
  51:  *		   n*ln2hi + n*ln2lo, where ln2hi is chosen such that the last
  52:  *		   20 bits (for VAX D format), or the last 21 bits ( for IEEE
  53:  *		   double) is 0. This ensures n*ln2hi is exactly representable.
  54:  *		2. In step 1, f may not be representable. A correction term c
  55:  *	 	   for f is computed. It follows that the correction term for
  56:  *		   f - t (the leading term of log(1+f) in step 2) is c-c*x. We
  57:  *		   add this correction term to n*ln2lo to attenuate the error.
  58:  *
  59:  *
  60:  * Special cases:
  61:  *	log1p(x) is NaN with signal if x < -1; log1p(NaN) is NaN with no signal;
  62:  *	log1p(INF) is +INF; log1p(-1) is -INF with signal;
  63:  *	only log1p(0)=0 is exact for finite argument.
  64:  *
  65:  * Accuracy:
  66:  *	log1p(x) returns the exact log(1+x) nearly rounded. In a test run
  67:  *	with 1,536,000 random arguments on a VAX, the maximum observed
  68:  *	error was .846 ulps (units in the last place).
  69:  *
  70:  * Constants:
  71:  * The hexadecimal values are the intended ones for the following constants.
  72:  * The decimal values may be used, provided that the compiler will convert
  73:  * from decimal to binary accurately enough to produce the hexadecimal values
  74:  * shown.
  75:  */
  76: 
  77: #ifdef VAX  /* VAX D format */
  78: #include <errno.h>
  79: 
  80: /* double static */
  81: /* ln2hi  =  6.9314718055829871446E-1    , Hex  2^  0   *  .B17217F7D00000 */
  82: /* ln2lo  =  1.6465949582897081279E-12   , Hex  2^-39   *  .E7BCD5E4F1D9CC */
  83: /* sqrt2  =  1.4142135623730950622E0     ; Hex  2^  1   *  .B504F333F9DE65 */
  84: static long     ln2hix[] = { 0x72174031, 0x0000f7d0};
  85: static long     ln2lox[] = { 0xbcd52ce7, 0xd9cce4f1};
  86: static long     sqrt2x[] = { 0x04f340b5, 0xde6533f9};
  87: #define    ln2hi    (*(double*)ln2hix)
  88: #define    ln2lo    (*(double*)ln2lox)
  89: #define    sqrt2    (*(double*)sqrt2x)
  90: #else   /* IEEE double */
  91: double static
  92: ln2hi  =  6.9314718036912381649E-1    , /*Hex  2^ -1   *  1.62E42FEE00000 */
  93: ln2lo  =  1.9082149292705877000E-10   , /*Hex  2^-33   *  1.A39EF35793C76 */
  94: sqrt2  =  1.4142135623730951455E0     ; /*Hex  2^  0   *  1.6A09E667F3BCD */
  95: #endif
  96: 
  97: double log1p(x)
  98: double x;
  99: {
 100:     static double zero=0.0, negone= -1.0, one=1.0,
 101:               half=1.0/2.0, small=1.0E-20;   /* 1+small == 1 */
 102:     double logb(),copysign(),scalb(),log__L(),z,s,t,c;
 103:     int k,finite();
 104: 
 105: #ifndef VAX
 106:     if(x!=x) return(x); /* x is NaN */
 107: #endif
 108: 
 109:     if(finite(x)) {
 110:        if( x > negone ) {
 111: 
 112:        /* argument reduction */
 113:           if(copysign(x,one)<small) return(x);
 114:           k=logb(one+x); z=scalb(x,-k); t=scalb(one,-k);
 115:           if(z+t >= sqrt2 )
 116:           { k += 1 ; z *= half; t *= half; }
 117:           t += negone; x = z + t;
 118:           c = (t-x)+z ;     /* correction term for x */
 119: 
 120:        /* compute log(1+x)  */
 121:               s = x/(2+x); t = x*x*half;
 122:           c += (k*ln2lo-c*x);
 123:           z = c+s*(t+log__L(s*s));
 124:           x += (z - t) ;
 125: 
 126:           return(k*ln2hi+x);
 127:        }
 128:     /* end of if (x > negone) */
 129: 
 130:         else {
 131: #ifdef VAX
 132:         extern double infnan();
 133:         if ( x == negone )
 134:             return (infnan(-ERANGE));   /* -INF */
 135:         else
 136:             return (infnan(EDOM));  /* NaN */
 137: #else   /* IEEE double */
 138:         /* x = -1, return -INF with signal */
 139:         if ( x == negone ) return( negone/zero );
 140: 
 141:         /* negative argument for log, return NaN with signal */
 142:             else return ( zero / zero );
 143: #endif
 144:         }
 145:     }
 146:     /* end of if (finite(x)) */
 147: 
 148:     /* log(-INF) is NaN */
 149:     else if(x<0)
 150:          return(zero/zero);
 151: 
 152:     /* log(+INF) is INF */
 153:     else return(x);
 154: }

Defined functions

Defined variables

Defined macros