home *** CD-ROM | disk | FTP | other *** search
/ Die Ultimative Software-P…i Collection 1996 & 1997 / Die Ultimative Software-Pakete CD-ROM fur Atari Collection 1996 & 1997.iso / g / gnu_c / pmlsrc23.zoo / pmlsrc / log.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-03-19  |  9.8 KB  |  423 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.  * Do not define that with gcc -O2 !!                    *
  15.  *                                                                      *
  16.  * 68881 support added by Michael Ritzert, November 1991        *
  17.  ************************************************************************
  18.  */
  19.  
  20. /************************************************************************
  21.  *                                    *
  22.  *                N O T I C E                *
  23.  *                                    *
  24.  *            Copyright Abandoned, 1987, Fred Fish        *
  25.  *                                    *
  26.  *    This previously copyrighted work has been placed into the    *
  27.  *    public domain by the author (Fred Fish) and may be freely used    *
  28.  *    for any purpose, private or commercial.  I would appreciate    *
  29.  *    it, as a courtesy, if this notice is left in all copies and    *
  30.  *    derivative works.  Thank you, and enjoy...            *
  31.  *                                    *
  32.  *    The author makes no warranty of any kind with respect to this    *
  33.  *    product and explicitly disclaims any implied warranties of    *
  34.  *    merchantability or fitness for any particular purpose.        *
  35.  *                                    *
  36.  ************************************************************************
  37.  */
  38.  
  39.  
  40. /*
  41.  *  FUNCTION
  42.  *
  43.  *    log   double precision natural log
  44.  *
  45.  *  KEY WORDS
  46.  *
  47.  *    log
  48.  *    machine independent routines
  49.  *    math libraries
  50.  *
  51.  *  DESCRIPTION
  52.  *
  53.  *    Returns double precision natural log of double precision
  54.  *    floating point argument.
  55.  *
  56.  *  USAGE
  57.  *
  58.  *    double log (x)
  59.  *    double x;
  60.  *
  61.  *  REFERENCES
  62.  *
  63.  *    Computer Approximations, J.F. Hart et al, John Wiley & Sons,
  64.  *    1968, pp. 105-111.
  65.  *
  66.  *  RESTRICTIONS
  67.  *
  68.  *    The absolute error in the approximating rational function is
  69.  *    10**(-19.38).  Note that contrary to DEC's assertion
  70.  *    in the F4P user's guide, the error is absolute and not
  71.  *    relative.
  72.  *      ** modified ** theoretical precision is only 10**(-16.5);
  73.  *      it works better in practice.
  74.  *
  75.  *    This error bound assumes exact arithmetic
  76.  *    in the polynomial evaluation.  Additional rounding and
  77.  *    truncation errors may occur as the argument is reduced
  78.  *    to the range over which the polynomial approximation
  79.  *    is valid, and as the polynomial is evaluated using
  80.  *    finite-precision arithmetic.
  81.  *    
  82.  *  PROGRAMMER
  83.  *
  84.  *    Fred Fish
  85.  *      Modifications - Michal Jaegermann
  86.  *
  87.  *  INTERNALS
  88.  *
  89.  *    Computes log(X) from:
  90.  *
  91.  *      (1)    If argument is zero then flag an error
  92.  *        and return minus infinity (or rather its
  93.  *        machine representation).
  94.  *
  95.  *      (2)    If argument is negative then flag an
  96.  *        error and return minus infinity.
  97.  *
  98.  *      (3)    Given that x = m * 2**k then extract
  99.  *        mantissa m and exponent k.
  100.  *        (3a)  If m = 0.5 then we know what is the final
  101.  *              result and calculations can be shortened.
  102.  *
  103.  *      (4)    Transform m which is in range [0.5,1.0]
  104.  *        to range [1/sqrt(2), sqrt(2)] by:
  105.  *
  106.  *            s = m * sqrt(2)
  107.  *
  108.  *      (5)    Compute z = (s - 1) / (s + 1)
  109.  *        (5a)  Modified - combine steps (4) and (5) computing
  110.  *              z = (m - 1/sqrt(2)) / (m + 1/sqrt(2))
  111.  *
  112.  *      (6)    Now use the approximation from HART
  113.  *        page 111 to find log(s):
  114.  *
  115.  *        log(s) = z * ( P(z**2) / Q(z**2) )
  116.  *
  117.  *        Where:
  118.  *
  119.  *        P(z**2) =  SUM [ Pj * (z**(2*j)) ]
  120.  *        over j = {0,1,2,3}
  121.  *
  122.  *        Q(z**2) =  SUM [ Qj * (z**(2*j)) ]
  123.  *        over j = {0,1,2,3}
  124.  *
  125.  *        P0 =  -0.240139179559210509868484e2
  126.  *        P1 =  0.30957292821537650062264e2
  127.  *        P2 =  -0.96376909336868659324e1
  128.  *        P3 =  0.4210873712179797145
  129.  *        Q0 =  -0.120069589779605254717525e2
  130.  *        Q1 =  0.19480966070088973051623e2
  131.  *        Q2 =  -0.89111090279378312337e1
  132.  *        Q3 =  1.0000
  133.  *
  134.  *        (coefficients from HART table #2705 pg 229)
  135.  *      (6a)    Modified - compute, as a polynomial of z, an
  136.  *              approximation of log2(m * sqrt(2)). Coefficients
  137.  *              for this polynomial are not given explicitely in HART
  138.  *              but can be computed from table #2665, for example.
  139.  *
  140.  *      (7)    Finally, compute log(x) from:
  141.  *
  142.  *        log(x) = (k * log(2)) - log(sqrt(2)) + log(s)
  143.  *
  144.  *        (7a)  log2(x) = k - 1/2 + log(m * sqrt(2)).  Multiply
  145.  *              by a constant to adjust a logarithm base.
  146.  *
  147.  */
  148.  
  149. #if !defined (__M68881__) && !defined (sfp004)
  150.  
  151. #include <stdio.h>
  152. #include <math.h>
  153. #include "pml.h"
  154.  
  155.             /* mjr++                */
  156.             /* use a different routine instead    */
  157.             /* see end of listing            */
  158.  
  159. # define HALFSQRT2        0.70710678118654752440
  160.  
  161. static double log_coeffs[] = {
  162.     2.88539008177793058026,
  163.     0.9617966939212138784,
  164.     0.577078018095541030,
  165.     0.4121983030496934,
  166.     0.32062199711811,
  167.     0.2612917312170,
  168.     0.24451102108
  169. };
  170.  
  171. #ifdef __GCC_HACK__
  172. static double log_of2[] = {
  173.     LN2, 0.30102999566398119521 /* log10(2) */
  174. };
  175. #endif
  176.  
  177. static char funcname[] = "log";
  178.  
  179. #ifdef __GCC_HACK__
  180. /*
  181.  * This fake function header may change even from version to a version
  182.  * of gcc compiler.  Ensure that first four assembler instructions in
  183.  * log and log10 are the same!
  184.  */
  185. __asm__("
  186.     .text
  187.     .even
  188.     .globl _log10
  189. _log10:");
  190. #ifdef __MSHORT__
  191. __asm__("
  192.     link a6,#-32
  193.     moveml #0x3e30,sp@-");
  194. #else
  195. __asm__("
  196.     link a6,#-36
  197.     moveml #0x3f30,sp@-");
  198. #endif  /* __MSHORT__ */
  199. __asm__("
  200.     movel a6@(8),d4
  201.     movel a6@(12),d5
  202.     moveq #1,d6
  203.     jra   lgentry");
  204. #endif  /* __GCC_HACK__ */
  205.  
  206. double log (x)
  207. double x;
  208. {
  209.     int k;
  210.     double s;
  211.     struct exception xcpt;
  212.     char * xcptstr;
  213. #ifdef __GCC_HACK__
  214.     int index;
  215.  
  216.     index = 0;
  217. __asm__("
  218. lgentry:");
  219. #endif
  220.  
  221.     if (x <= 0.0) {
  222.     xcpt.name = funcname;
  223.     xcpt.arg1 = x;
  224.         if (x == 0.0) {
  225.         xcpt.retval = -MAXDOUBLE;
  226.         xcpt.type = SING;
  227.         xcptstr = "SINGULARITY";
  228.     }
  229.     else {
  230.         xcpt.retval = MAXDOUBLE;
  231.         xcpt.type = DOMAIN;
  232.         xcptstr = "DOMAIN";
  233.     }
  234.     if (!matherr (&xcpt)) {
  235.         fprintf (stderr, "%s: %s error\n", xcptstr, funcname);
  236.         errno = EDOM;
  237.     }
  238.     }
  239.     else {
  240.     s = frexp (x, &k);
  241.     if (0.5 == s ) {
  242.         s = -1.0;
  243.     }
  244.     else {
  245.         /* range reduction with s multiplied by SQRT2 */
  246.         s = (s - HALFSQRT2) / (s + HALFSQRT2);
  247.         /* polynomial approximation */
  248.         s *= poly ((sizeof(log_coeffs)/sizeof(double)) - 1,
  249.                    log_coeffs, s * s);
  250.         /* and subtract log2 of SQRT2 */
  251.         s -= 0.5;
  252.     }
  253.     /* add the original binary exponent */
  254.     s += k;
  255.     /* multiply to get a requested logarithm */
  256. #ifdef __GCC_HACK__
  257.     xcpt.retval = s * log_of2[index];
  258. #else
  259.     xcpt.retval = s * LN2;
  260. #endif
  261.     }
  262.     return (xcpt.retval);
  263. }
  264. #endif /* !__M68881__ && ! sfp004    */
  265.  
  266. /*****************************************************************/
  267.  
  268. #ifdef    sfp004
  269.  
  270. __asm("
  271.  
  272. comm =     -6
  273. resp =    -16
  274. zahl =      0
  275.  
  276. ");    /* end asm    */
  277.  
  278. #endif    sfp004
  279. #if defined (__M68881__) || defined (sfp004)
  280.     __asm(".text; .even");
  281.  
  282. # ifdef    ERROR_CHECK
  283.  
  284.     __asm("
  285.  
  286. _Overflow:
  287.     .ascii \"OVERFLOW\\0\"
  288. _Domain:
  289.     .ascii \"DOMAIN\\0\"
  290. _Error_String:
  291.     .ascii \"log: %s error\\n\\0\"
  292. .even
  293.  
  294. | pml compatible log
  295. | m.ritzert 7.12.1991
  296. | ritzert@dfg.dbp.de
  297. |
  298. |    /* NAN  = {7fffffff,ffffffff}        */
  299. |    /* +Inf = {7ff00000,00000000}        */
  300. |    /* -Inf = {fff00000,00000000}        */
  301. |    /* MAX_D= {7fee42d1,30773b76}        */
  302. |    /* MIN_D= {ffee42d1,30773b76}        */
  303.  
  304. .even
  305. double_max:
  306.     .long    0x7fee42d1
  307.     .long    0x30273b76
  308. double_min:
  309.     .long    0xffee42d1
  310.     .long    0x30273b76
  311. NaN:
  312.     .long    0x7fffffff
  313.     .long    0xffffffff
  314. p_Inf:
  315.     .long    0x7ff00000
  316.     .long    0x00000000
  317. m_Inf:
  318.     .long    0xfff00000
  319.     .long    0x00000000
  320. ");
  321. # endif    ERROR_CHECK
  322.  
  323.     __asm("
  324. .even
  325.     .globl _log
  326. _log:
  327.     ");    /* end asm    */
  328.  
  329. #endif    /* __M68881__ || sfp004    */
  330. #ifdef    __M68881__
  331.  
  332.     __asm("
  333.     flognd    a7@(4), fp0    | log
  334.     fmoved    fp0,a7@-    | push result
  335.     moveml    a7@+,d0-d1    | return_value
  336.     ");    /* end asm    */
  337.  
  338. #endif    __M68881__
  339. #ifdef    sfp004
  340.     __asm("
  341.     lea    0xfffa50,a0
  342.     movew    #0x5414,a0@(comm)    | specify function
  343.     cmpiw    #0x8900,a0@(resp)    | check
  344.     movel    a7@(4),a0@        | load arg_hi
  345.     movel    a7@(8),a0@        | load arg_low
  346.     movew    #0x7400,a0@(comm)    | result to d0
  347.     .long    0x0c688900, 0xfff067f8    | wait
  348.     movel    a0@,d0
  349.     movel    a0@,d1
  350.     ");    /* end asm    */
  351.  
  352. #endif    sfp004
  353. #if defined (__M68881__) || defined (sfp004)
  354. # ifdef    ERROR_CHECK
  355.     __asm("
  356.     lea    double_max,a0    |
  357.     swap    d0        | exponent into lower word
  358.     cmpw    a0@(16),d0    | == NaN ?
  359.     beq    error_nan    |
  360.     cmpw    a0@(24),d0    | == + Infinity ?
  361.     beq    error_plus    |
  362.     cmpw    a0@(32),d0    | == - Infinity ?
  363.     beq    error_minus    |
  364.     swap    d0        | result ok,
  365.     rts            | restore d0
  366. ");
  367. #ifndef    __MSHORT__
  368. __asm("
  369. error_minus:
  370.     swap    d0
  371.     moveml    d0-d1,a7@-
  372.     movel    #63,_errno    | errno = ERANGE
  373.     pea    _Overflow    | for printf
  374.     bra    error_exit    |
  375. error_plus:
  376.     swap    d0
  377.     moveml    d0-d1,a7@-
  378.     movel    #63,_errno    | NAN => errno = EDOM
  379.     pea    _Overflow    | for printf
  380.     bra    error_exit    |
  381. error_nan:
  382.     moveml    a0@(24),d0-d1    | result = +inf
  383.     moveml    d0-d1,a7@-
  384.     movel    #62,_errno    | NAN => errno = EDOM
  385.     pea    _Domain        | for printf
  386. ");
  387. #else    __MSHORT__
  388. __asm("
  389. error_minus:
  390.     swap    d0
  391.     moveml    d0-d1,a7@-
  392.     movew    #63,_errno    | errno = ERANGE
  393.     pea    _Overflow    | for printf
  394.     bra    error_exit    |
  395. error_plus:
  396.     swap    d0
  397.     moveml    d0-d1,a7@-
  398.     movew    #63,_errno    | NAN => errno = EDOM
  399.     pea    _Overflow    | for printf
  400.     bra    error_exit    |
  401. error_nan:
  402.     moveml    a0@(24),d0-d1    | result = +inf
  403.     moveml    d0-d1,a7@-
  404.     movew    #62,_errno    | NAN => errno = EDOM
  405.     pea    _Domain        | for printf
  406. ");
  407. #endif    __MSHORT__
  408. __asm("
  409. error_exit:
  410.     pea    _Error_String    |
  411.     pea    __iob+52    |
  412.     jbsr    _fprintf    |
  413.     addl    #12,a7        |
  414.     moveml    a7@+,d0-d1
  415.     rts
  416.     ");
  417. # else    ERROR_CHECK
  418.  
  419. __asm("rts");
  420.  
  421. # endif    ERROR_CHECK
  422. #endif /* __M68881__ || sfp004    */
  423.