home *** CD-ROM | disk | FTP | other *** search
/ Fresh Fish 4 / FreshFish_May-June1994.bin / bsd / src / libm / libm-amiga / common_source / log.c < prev    next >
C/C++ Source or Header  |  1993-09-23  |  5KB  |  164 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[] = "@(#)log.c    5.6 (Berkeley) 10/9/90";
  36. #endif /* not lint */
  37.  
  38. /* LOG(X)
  39.  * RETURN THE LOGARITHM OF x 
  40.  * DOUBLE PRECISION (VAX D FORMAT 56 bits or IEEE DOUBLE 53 BITS)
  41.  * CODED IN C BY K.C. NG, 1/19/85;
  42.  * REVISED BY K.C. NG on 2/7/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.  *            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(x) = k*ln2 + log(1+f).  (Here n*ln2 will be stored
  69.  *       in two floating point number: n*ln2hi + n*ln2lo, n*ln2hi is exact
  70.  *       since the last 20 bits of ln2hi is 0.)
  71.  *
  72.  * Special cases:
  73.  *    log(x) is NaN with signal if x < 0 (including -INF) ; 
  74.  *    log(+INF) is +INF; log(0) is -INF with signal;
  75.  *    log(NaN) is that NaN with no signal.
  76.  *
  77.  * Accuracy:
  78.  *    log(x) returns the exact log(x) nearly rounded. In a test run with
  79.  *    1,536,000 random arguments on a VAX, the maximum observed error was
  80.  *    .826 ulps (units in the last place).
  81.  *
  82.  * Constants:
  83.  * The hexadecimal values are the intended ones for the following constants.
  84.  * The decimal values may be used, provided that the compiler will convert
  85.  * from decimal to binary accurately enough to produce the hexadecimal values
  86.  * shown.
  87.  */
  88.  
  89. #include <errno.h>
  90. #include "mathimpl.h"
  91.  
  92. vc(ln2hi, 6.9314718055829871446E-1  ,7217,4031,0000,f7d0,   0, .B17217F7D00000)
  93. vc(ln2lo, 1.6465949582897081279E-12 ,bcd5,2ce7,d9cc,e4f1, -39, .E7BCD5E4F1D9CC)
  94. vc(sqrt2, 1.4142135623730950622E0   ,04f3,40b5,de65,33f9,   1, .B504F333F9DE65)
  95.  
  96. ic(ln2hi, 6.9314718036912381649E-1,   -1, 1.62E42FEE00000)
  97. ic(ln2lo, 1.9082149292705877000E-10, -33, 1.A39EF35793C76)
  98. ic(sqrt2, 1.4142135623730951455E0,     0, 1.6A09E667F3BCD)
  99.  
  100. #ifdef vccast
  101. #define    ln2hi    vccast(ln2hi)
  102. #define    ln2lo    vccast(ln2lo)
  103. #define    sqrt2    vccast(sqrt2)
  104. #endif
  105.  
  106.  
  107. double log(x)
  108. double x;
  109. {
  110.     const static double zero=0.0, negone= -1.0, half=1.0/2.0;
  111.     double s,z,t;
  112.     int k,n;
  113.  
  114. #if !defined(vax)&&!defined(tahoe)
  115.     if(x!=x) return(x);    /* x is NaN */
  116. #endif    /* !defined(vax)&&!defined(tahoe) */
  117.     if(finite(x)) {
  118.        if( x > zero ) {
  119.  
  120.        /* argument reduction */
  121.           k=logb(x);   x=scalb(x,-k);
  122.           if(k == -1022) /* subnormal no. */
  123.            {n=logb(x); x=scalb(x,-n); k+=n;} 
  124.           if(x >= sqrt2 ) {k += 1; x *= half;}
  125.           x += negone ;
  126.  
  127.        /* compute log(1+x)  */
  128.               s=x/(2+x); t=x*x*half;
  129.           z=k*ln2lo+s*(t+log__L(s*s));
  130.           x += (z - t) ;
  131.  
  132.           return(k*ln2hi+x);
  133.        }
  134.     /* end of if (x > zero) */
  135.  
  136.        else {
  137. #if defined(vax)||defined(tahoe)
  138.         if ( x == zero )
  139.             return (infnan(-ERANGE));    /* -INF */
  140.         else
  141.             return (infnan(EDOM));    /* NaN */
  142. #else    /* defined(vax)||defined(tahoe) */
  143.         /* zero argument, return -INF with signal */
  144.         if ( x == zero )
  145.             return( negone/zero );
  146.  
  147.         /* negative argument, return NaN with signal */
  148.         else 
  149.             return ( zero / zero );
  150. #endif    /* defined(vax)||defined(tahoe) */
  151.         }
  152.     }
  153.     /* end of if (finite(x)) */
  154.     /* NOTREACHED if defined(vax)||defined(tahoe) */
  155.  
  156.     /* log(-INF) is NaN with signal */
  157.     else if (x<0) 
  158.         return(zero/zero);      
  159.  
  160.     /* log(+INF) is +INF */
  161.     else return(x);      
  162.  
  163. }
  164.