home *** CD-ROM | disk | FTP | other *** search
/ GEMini Atari / GEMini_Atari_CD-ROM_Walnut_Creek_December_1993.iso / files / language / icon8_fx / log.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-10-23  |  7.4 KB  |  286 lines
  1. /************************************************************************
  2.  *                                                                      *
  3.  * A replacement log() routine for PML library.  It really computes     *
  4.  * base 2 logarithm which can be multiplied by a proper constant        *
  5.  * to provide final answer.  A rational approximation of an original    *
  6.  * is replaced by a polynomial one.  In practice, with a software       *
  7.  * floating point, this gives a better precision.                       *
  8.  *                                                                      *
  9.  * Michal Jaegermann, December 1990                                     *
  10.  *                                                                      *
  11.  * If GCC_HACK is defined then we are folding log and log10 routines    *
  12.  * by making in assembler an extra entry point.  Do not define that     *
  13.  * for portable routine!!                                               *
  14.  *                                                                      *
  15.  ************************************************************************
  16.  */
  17. /* #define GCC_HACK */
  18. /************************************************************************
  19.  *                                    *
  20.  *                N O T I C E                *
  21.  *                                    *
  22.  *            Copyright Abandoned, 1987, Fred Fish        *
  23.  *                                    *
  24.  *    This previously copyrighted work has been placed into the    *
  25.  *    public domain by the author (Fred Fish) and may be freely used    *
  26.  *    for any purpose, private or commercial.  I would appreciate    *
  27.  *    it, as a courtesy, if this notice is left in all copies and    *
  28.  *    derivative works.  Thank you, and enjoy...            *
  29.  *                                    *
  30.  *    The author makes no warranty of any kind with respect to this    *
  31.  *    product and explicitly disclaims any implied warranties of    *
  32.  *    merchantability or fitness for any particular purpose.        *
  33.  *                                    *
  34.  ************************************************************************
  35.  */
  36.  
  37.  
  38. /*
  39.  *  FUNCTION
  40.  *
  41.  *    log   double precision natural log
  42.  *
  43.  *  KEY WORDS
  44.  *
  45.  *    log
  46.  *    machine independent routines
  47.  *    math libraries
  48.  *
  49.  *  DESCRIPTION
  50.  *
  51.  *    Returns double precision natural log of double precision
  52.  *    floating point argument.
  53.  *
  54.  *  USAGE
  55.  *
  56.  *    double log (x)
  57.  *    double x;
  58.  *
  59.  *  REFERENCES
  60.  *
  61.  *    Computer Approximations, J.F. Hart et al, John Wiley & Sons,
  62.  *    1968, pp. 105-111.
  63.  *
  64.  *  RESTRICTIONS
  65.  *
  66.  *    The absolute error in the approximating rational function is
  67.  *    10**(-19.38).  Note that contrary to DEC's assertion
  68.  *    in the F4P user's guide, the error is absolute and not
  69.  *    relative.
  70.  *      ** modified ** theoretical precision is only 10**(-16.5);
  71.  *      it works better in practice. 
  72.  *
  73.  *    This error bound assumes exact arithmetic
  74.  *    in the polynomial evaluation.  Additional rounding and
  75.  *    truncation errors may occur as the argument is reduced
  76.  *    to the range over which the polynomial approximation
  77.  *    is valid, and as the polynomial is evaluated using
  78.  *    finite-precision arithmetic.
  79.  *    
  80.  *  PROGRAMMER
  81.  *
  82.  *    Fred Fish
  83.  *      Modifications - Michal Jaegermann
  84.  *
  85.  *  INTERNALS
  86.  *
  87.  *    Computes log(X) from:
  88.  *
  89.  *      (1)    If argument is zero then flag an error
  90.  *        and return minus infinity (or rather its
  91.  *        machine representation).
  92.  *
  93.  *      (2)    If argument is negative then flag an
  94.  *        error and return minus infinity.
  95.  *
  96.  *      (3)    Given that x = m * 2**k then extract
  97.  *        mantissa m and exponent k.
  98.  *        (3a)  If m = 0.5 then we know what is the final
  99.  *              result and calculations can be shortened. 
  100.  *
  101.  *      (4)    Transform m which is in range [0.5,1.0]
  102.  *        to range [1/sqrt(2), sqrt(2)] by:
  103.  *
  104.  *            s = m * sqrt(2)
  105.  *
  106.  *      (5)    Compute z = (s - 1) / (s + 1)
  107.  *        (5a)  Modified - combine steps (4) and (5) computing
  108.  *              z = (m - 1/sqrt(2)) / (m + 1/sqrt(2))
  109.  *
  110.  *      (6)    Now use the approximation from HART
  111.  *        page 111 to find log(s):
  112.  *
  113.  *        log(s) = z * ( P(z**2) / Q(z**2) )
  114.  *
  115.  *        Where:
  116.  *
  117.  *        P(z**2) =  SUM [ Pj * (z**(2*j)) ]
  118.  *        over j = {0,1,2,3}
  119.  *
  120.  *        Q(z**2) =  SUM [ Qj * (z**(2*j)) ]
  121.  *        over j = {0,1,2,3}
  122.  *
  123.  *        P0 =  -0.240139179559210509868484e2
  124.  *        P1 =  0.30957292821537650062264e2
  125.  *        P2 =  -0.96376909336868659324e1
  126.  *        P3 =  0.4210873712179797145
  127.  *        Q0 =  -0.120069589779605254717525e2
  128.  *        Q1 =  0.19480966070088973051623e2
  129.  *        Q2 =  -0.89111090279378312337e1
  130.  *        Q3 =  1.0000
  131.  *
  132.  *        (coefficients from HART table #2705 pg 229)
  133.   *      (6a)    Modified - compute, as a polynomial of z, an
  134.  *              approximation of log2(m * sqrt(2)). Coefficients
  135.  *              for this polynomial are not given explicitely in HART
  136.  *              but can be computed from table #2665, for example.
  137.  *
  138.  *      (7)    Finally, compute log(x) from:
  139.  *
  140.  *        log(x) = (k * log(2)) - log(sqrt(2)) + log(s)
  141.  *
  142.  *        (7a)  log2(x) = k - 1/2 + log2(m * sqrt(2)).  Multiply
  143.  *              by a constant to adjust a logarithm base.
  144.  *
  145.  */
  146.  
  147. #include <stdio.h>
  148. #include <pmluser.h>
  149. #include "pml.h"
  150.  
  151. #define HALFSQRT2        0.70710678118654752440
  152.  
  153. static double log_coeffs[] = {
  154.     2.88539008177793058026,
  155.     0.9617966939212138784,
  156.     0.577078018095541030,
  157.     0.4121983030496934,
  158.     0.32062199711811,
  159.     0.2612917312170,
  160.     0.24451102108
  161. };
  162.  
  163. #ifdef GCC_HACK
  164. static double log_of2[] = {
  165.     LN2, 0.30102999566398119521 /* log10(2) */
  166. };
  167.  
  168. static char funcname[] = "log";
  169. #else
  170. static char funcname[] = "log2";
  171. #endif
  172.  
  173.  
  174.  
  175. #ifdef GCC_HACK
  176. /*
  177.  * This fake function header may change even from version to a version
  178.  * of gcc compiler.  Ensure that first four assembler instructions in
  179.  * log and log10 are the same!
  180.  */
  181. __asm__("
  182.     .text
  183.     .even
  184.     .globl _log10
  185. _log10:");
  186. #ifdef __MSHORT__
  187. __asm__("
  188.     link a6,#-24
  189.     moveml #0x3e30,sp@-");
  190. #else
  191. __asm__("
  192.     link a6,#-28
  193.     moveml #0x3f30,sp@-");
  194. #endif  /* __MSHORT__ */
  195. __asm__("
  196.     movel a6@(8),d2
  197.     movel a6@(12),d3
  198.     moveq #1,d6
  199.     jra   lgentry");
  200.  
  201. double log (x)
  202. #else
  203. double log2 (x)
  204. #endif  /* GCC_HACK */
  205.  
  206. double x;
  207. {
  208.     auto int k;
  209.     auto double s;
  210.     auto struct exception xcpt;
  211.     auto char * xcptstr;
  212. #ifdef GCC_HACK
  213.     auto int index;
  214. #endif
  215.     extern double frexp ();
  216.     extern double poly ();
  217.  
  218.  
  219.     DBUG_ENTER (funcname);
  220.     DBUG_3 ("login", "arg %le", x);
  221. #ifdef GCC_HACK
  222.     index = 0;
  223.  
  224. __asm__("
  225. lgentry:");
  226. #endif
  227.  
  228.     if (x <= 0.0) {
  229.     xcpt.retval = -MAXDOUBLE;
  230.     xcpt.name = funcname;
  231.     xcpt.arg1 = x;
  232.         if (x == 0.0) {
  233.         xcpt.type = SING;
  234.         xcptstr = "SINGULARITY";
  235.     }
  236.     else {
  237.         xcpt.type = DOMAIN;
  238.         xcptstr = "DOMAIN";
  239.     }
  240.     if (!matherr (&xcpt)) {
  241.         fprintf (stderr, "%s: %s error\n", xcptstr, funcname);
  242.         errno = EDOM;
  243.     }
  244.     }
  245.     else {
  246.     s = frexp (x, &k);
  247.     if (0.5 == s ) {
  248.         s = -1.0;
  249.     }
  250.     else {
  251.         /* range reduction with s multiplied by SQRT2 */
  252.         s = (s - HALFSQRT2) / (s + HALFSQRT2);
  253.         /* polynomial approximation */
  254.         s *= poly ((sizeof(log_coeffs)/sizeof(double)) - 1, 
  255.                    log_coeffs, s * s);
  256.         /* and subtract log2 of SQRT2 */
  257.         s -= 0.5;
  258.     }
  259.     /* add the original binary exponent */
  260.     s += k;
  261.     /* multiply to get a requested logarithm */
  262. #ifdef GCC_HACK
  263.     xcpt.retval = s * log_of2[index];
  264. #else
  265.     xcpt.retval = s;
  266. #endif
  267.     }
  268.     DBUG_3 ("logout", "result %le", xcpt.retval);
  269.     DBUG_RETURN (xcpt.retval);
  270. }
  271.  
  272. #ifndef GCC_HACK
  273. double log (x)
  274. double x;
  275. {
  276.     return (LN2 * log2(x));
  277. }
  278.  
  279. double log10 (x)
  280. double x;
  281. {
  282.     return (0.30102999566398119521 * log2(x));
  283. }
  284.  
  285. #endif /* GCC_HACK */
  286.