home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / LIBSRC.ZOO / libsrc / math / log1p.c < prev    next >
C/C++ Source or Header  |  1992-02-03  |  6KB  |  171 lines
  1. /*
  2.  * Copyright (c) 1985 Regents of the University of California.
  3.  * All rights reserved.
  4.  *
  5.  * Redistribution and use in source and binary forms, with or without
  6.  * modification, are permitted provided that the following conditions
  7.  * are met:
  8.  * 1. Redistributions of source code must retain the above copyright
  9.  *    notice, this list of conditions and the following disclaimer.
  10.  * 2. Redistributions in binary form must reproduce the above copyright
  11.  *    notice, this list of conditions and the following disclaimer in the
  12.  *    documentation and/or other materials provided with the distribution.
  13.  * 3. All advertising materials mentioning features or use of this software
  14.  *    must display the following acknowledgement:
  15.  *    This product includes software developed by the University of
  16.  *    California, Berkeley and its contributors.
  17.  * 4. Neither the name of the University nor the names of its contributors
  18.  *    may be used to endorse or promote products derived from this software
  19.  *    without specific prior written permission.
  20.  *
  21.  * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
  22.  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  23.  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  24.  * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
  25.  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  26.  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
  27.  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  28.  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  29.  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  30.  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  31.  * SUCH DAMAGE.
  32.  */
  33.  
  34. #ifndef lint
  35. static char sccsid[] = "@(#)log1p.c    5.6 (Berkeley) 10/9/90";
  36. #endif /* not lint */
  37.  
  38. /* LOG1P(x) 
  39.  * RETURN THE LOGARITHM OF 1+x
  40.  * DOUBLE PRECISION (VAX D FORMAT 56 bits, IEEE DOUBLE 53 BITS)
  41.  * CODED IN C BY K.C. NG, 1/19/85; 
  42.  * REVISED BY K.C. NG on 2/6/85, 3/7/85, 3/24/85, 4/16/85.
  43.  * 
  44.  * Required system supported functions:
  45.  *    scalb(x,n) 
  46.  *    copysign(x,y)
  47.  *    logb(x)    
  48.  *    finite(x)
  49.  *
  50.  * Required kernel function:
  51.  *    log__L(z)
  52.  *
  53.  * Method :
  54.  *    1. Argument Reduction: find k and f such that 
  55.  *            1+x  = 2^k * (1+f), 
  56.  *       where  sqrt(2)/2 < 1+f < sqrt(2) .
  57.  *
  58.  *    2. Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
  59.  *         = 2s + 2/3 s**3 + 2/5 s**5 + .....,
  60.  *       log(1+f) is computed by
  61.  *
  62.  *                 log(1+f) = 2s + s*log__L(s*s)
  63.  *       where
  64.  *        log__L(z) = z*(L1 + z*(L2 + z*(... (L6 + z*L7)...)))
  65.  *
  66.  *       See log__L() for the values of the coefficients.
  67.  *
  68.  *    3. Finally,  log(1+x) = k*ln2 + log(1+f).  
  69.  *
  70.  *    Remarks 1. In step 3 n*ln2 will be stored in two floating point numbers
  71.  *           n*ln2hi + n*ln2lo, where ln2hi is chosen such that the last 
  72.  *           20 bits (for VAX D format), or the last 21 bits ( for IEEE 
  73.  *           double) is 0. This ensures n*ln2hi is exactly representable.
  74.  *        2. In step 1, f may not be representable. A correction term c
  75.  *            for f is computed. It follows that the correction term for
  76.  *           f - t (the leading term of log(1+f) in step 2) is c-c*x. We
  77.  *           add this correction term to n*ln2lo to attenuate the error.
  78.  *
  79.  *
  80.  * Special cases:
  81.  *    log1p(x) is NaN with signal if x < -1; log1p(NaN) is NaN with no signal;
  82.  *    log1p(INF) is +INF; log1p(-1) is -INF with signal;
  83.  *    only log1p(0)=0 is exact for finite argument.
  84.  *
  85.  * Accuracy:
  86.  *    log1p(x) returns the exact log(1+x) nearly rounded. In a test run 
  87.  *    with 1,536,000 random arguments on a VAX, the maximum observed
  88.  *    error was .846 ulps (units in the last place).
  89.  *
  90.  * Constants:
  91.  * The hexadecimal values are the intended ones for the following constants.
  92.  * The decimal values may be used, provided that the compiler will convert
  93.  * from decimal to binary accurately enough to produce the hexadecimal values
  94.  * shown.
  95.  */
  96.  
  97. #include <errno.h>
  98. #include "mathimpl.h"
  99.  
  100. vc(ln2hi, 6.9314718055829871446E-1  ,7217,4031,0000,f7d0,   0, .B17217F7D00000)
  101. vc(ln2lo, 1.6465949582897081279E-12 ,bcd5,2ce7,d9cc,e4f1, -39, .E7BCD5E4F1D9CC)
  102. vc(sqrt2, 1.4142135623730950622E0   ,04f3,40b5,de65,33f9,   1, .B504F333F9DE65)
  103.  
  104. ic(ln2hi, 6.9314718036912381649E-1,   -1, 1.62E42FEE00000)
  105. ic(ln2lo, 1.9082149292705877000E-10, -33, 1.A39EF35793C76)
  106. ic(sqrt2, 1.4142135623730951455E0,     0, 1.6A09E667F3BCD)
  107.  
  108. #ifdef vccast
  109. #define    ln2hi    vccast(ln2hi)
  110. #define    ln2lo    vccast(ln2lo)
  111. #define    sqrt2    vccast(sqrt2)
  112. #endif
  113.  
  114. double log1p(x)
  115. double x;
  116. {
  117.     const static double zero=0.0, negone= -1.0, one=1.0, 
  118.               half=1.0/2.0, small=1.0E-20;   /* 1+small == 1 */
  119.     double z,s,t,c;
  120.     int k;
  121.  
  122. #if !defined(vax)&&!defined(tahoe)
  123.     if(x!=x) return(x);    /* x is NaN */
  124. #endif    /* !defined(vax)&&!defined(tahoe) */
  125.  
  126.     if(finite(x)) {
  127.        if( x > negone ) {
  128.  
  129.        /* argument reduction */
  130.           if(copysign(x,one)<small) return(x);
  131.           k=logb(one+x); z=scalb(x,-k); t=scalb(one,-k);
  132.           if(z+t >= sqrt2 ) 
  133.           { k += 1 ; z *= half; t *= half; }
  134.           t += negone; x = z + t;
  135.           c = (t-x)+z ;        /* correction term for x */
  136.  
  137.         /* compute log(1+x)  */
  138.               s = x/(2+x); t = x*x*half;
  139.           c += (k*ln2lo-c*x);
  140.           z = c+s*(t+log__L(s*s));
  141.           x += (z - t) ;
  142.  
  143.           return(k*ln2hi+x);
  144.        }
  145.     /* end of if (x > negone) */
  146.  
  147.         else {
  148. #if defined(vax)||defined(tahoe)
  149.         if ( x == negone )
  150.             return (infnan(-ERANGE));    /* -INF */
  151.         else
  152.             return (infnan(EDOM));    /* NaN */
  153. #else    /* defined(vax)||defined(tahoe) */
  154.         /* x = -1, return -INF with signal */
  155.         if ( x == negone ) return( negone/zero );
  156.  
  157.         /* negative argument for log, return NaN with signal */
  158.             else return ( zero / zero );
  159. #endif    /* defined(vax)||defined(tahoe) */
  160.         }
  161.     }
  162.     /* end of if (finite(x)) */
  163.  
  164.     /* log(-INF) is NaN */
  165.     else if(x<0) 
  166.          return(zero/zero);
  167.  
  168.     /* log(+INF) is INF */
  169.     else return(x);      
  170. }
  171.