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