workbench/libs/mathieeedoubtrans/ieeedplog.c

   1 /*
   2     Copyright © 1995-2008, The AROS Development Team. All rights reserved.
   3     $Id$
   4 */
   5
   6 #include "mathieeedoubtrans_intern.h"
   7
   8 /*
   9     FUNCTION
  10       Calculate logarithm (base e) of the given IEEE double precision number
  11
  12     RESULT
  13       IEEE double precision number
  14
  15       flags:
  16         zero     : result is zero
  17         negative : result is negative
  18         overflow : argument was negative
  19
  20     NOTES
  21
  22     EXAMPLE
  23
  24     BUGS
  25
  26     SEE ALSO
  27
  28     INTERNALS
  29       ALGORITHM:
  30
  31       If the Argument is negative set overflow-flag and return NAN.
  32       If the Argument is 0 return 0xFFF0000000000000.
  33       If the Argument is pos. Infinity return pos. Infinity.
  34
  35       All other cases:
  36
  37       (ld is the logarithm with base 2)
  38       (ln is the logarithm with base e)
  39       y = M * 2^E
  40
  41       ln y = ln ( M * 2^E ) =
  42
  43            = ln M + ln 2^E =
  44
  45            = ln M + E * ln (2) =
  46
  47              ld M        ld 2
  48            = ----- + E * ----- =      [ld 2 = 1]
  49              ld e        ld e
  50
  51              ld M + E
  52            = --------
  53              ld e
  54
  55       ld e can be precalculated, of course.
  56       For calculating ld M see file intern_ieeespld.c
  57
  58     HISTORY
  59 */
  60
  61 AROS_LHQUAD1(double, IEEEDPLog,
  62     AROS_LHAQUAD(double, y, D0, D1),
  63     struct MathIeeeDoubTransBase *, MathIeeeDoubTransBase, 14, MathIeeeDoubTrans
  64 )
  65 {
  66     AROS_LIBFUNC_INIT
  67
  68     QUAD Res, tmp, Exponent64, ld_M;
  69     LONG Exponent;
  70     /* check for negative sign */
  71     if ( is_lessSC(y, 0x0, 0x0) /* y<0 */)
  72     {
  73         SetSR(Overflow_Bit, Zero_Bit | Negative_Bit | Overflow_Bit);
  74         Set_Value64C(Res, IEEEDPNAN_Hi, IEEEDPNAN_Lo);
  75         return Res;
  76     }
  77
  78     if ( is_eqC(y, 0x0, 0x0) )
  79     {
  80         Set_Value64C
  81         (
  82             Res,
  83             (IEEEDPSign_Mask_Hi + IEEEDPExponent_Mask_Hi),
  84             (IEEEDPSign_Mask_Lo + IEEEDPExponent_Mask_Lo)
  85         );
  86         return Res;
  87     }
  88     /* check for argument == 0 or argument == +infinity */
  89     if
  90     (
  91            is_eqC(y, IEEEDPPInfty_Hi, IEEEDPPInfty_Lo)
  92         || is_eqC(y, IEEEDPExponent_Mask_Hi, IEEEDPExponent_Mask_Lo)
  93     )
  94     {
  95         return y;
  96     }
  97
  98     /* convert the Exponent of the argument (y) to the ieeedp-format */
  99     Exponent = ((Get_High32of64(y) & IEEEDPExponent_Mask_Hi) >> 20) - 0x3fe;
 100     Exponent64 = IEEEDPFlt(Exponent);
 101
 102     Set_Value64(tmp, y);
 103     AND64QC(tmp, IEEEDPMantisse_Mask_Hi, IEEEDPMantisse_Mask_Lo);
 104     OR64QC(tmp, 0x3fe00000, 0x0);
 105     ld_M = intern_IEEEDPLd
 106     (
 107         (struct MathIeeeDoubTransBase *) MathIeeeDoubTransBase,
 108         tmp
 109     );
 110     /*
 111                       ld M + E
 112         ln(fnum1) =  --------
 113                         ld e
 114     */
 115     Set_Value64C(tmp, 0x3fe62e42, 0xfefa39ef ); /* ln 2*/
 116     return IEEEDPMul( IEEEDPAdd(ld_M, Exponent64), tmp);
 117
 118     AROS_LIBFUNC_EXIT
 119 }