workbench/libs/mathieeedoubtrans/ieeedpatan.c

   1 /*
   2     Copyright © 1995-2003, The AROS Development Team. All rights reserved.
   3     $Id$
   4 */
   5
   6 #include "mathieeedoubtrans_intern.h"
   7
   8 /*
   9     FUNCTION
  10       Calculates the angle of a given number representing the tangent
  11       of that angle. The angle will be in radians.
  12
  13     RESULT
  14       IEEE double precision floating point number
  15
  16     NOTES
  17
  18     EXAMPLE
  19
  20     BUGS
  21
  22     SEE ALSO
  23
  24     INTERNALS
  25
  26     HISTORY
  27 */
  28
  29 AROS_LHQUAD1(double, IEEEDPAtan,
  30     AROS_LHAQUAD(double, y, D0, D1),
  31     struct MathIeeeDoubTransBase *, MathIeeeDoubTransBase, 5, MathIeeeDoubTrans
  32 )
  33 {
  34     AROS_LIBFUNC_INIT
  35
  36     QUAD yabs;
  37     QUAD ysquared, ycubed, ypow5, one;
  38
  39     AND64C
  40     (
  41         yabs, y,
  42         (IEEEDPMantisse_Mask_Hi | IEEEDPExponent_Mask_Hi),
  43         (IEEEDPMantisse_Mask_Lo | IEEEDPExponent_Mask_Lo)
  44     );
  45
  46     /* check for +- infinity -> output: +-pi/2 */
  47     if (is_eqC(yabs, IEEEDPPInfty_Hi, IEEEDPPInfty_Lo))
  48     {
  49         AND64QC(y, IEEEDPSign_Mask_Hi, IEEEDPSign_Mask_Lo);
  50         OR64QC(y, pio2_hi_Hi, pio2_hi_Lo);
  51         return y;
  52     }
  53     /* atan(0) = 0 */
  54     if (is_eqC(yabs, 0, 0) /* 0 == yabs */ )
  55     {
  56         SetSR(Zero_Bit, Zero_Bit | Negative_Bit | Overflow_Bit);
  57         return y;
  58     }
  59
  60     ysquared = IEEEDPMul(yabs, yabs);
  61     ycubed = IEEEDPMul(yabs, ysquared);
  62
  63     Set_Value64C(one, one_Hi, one_Lo);
  64
  65     /* atan(x >= 860) = pi/2 - 1/x + 1/(3*x^3) */
  66     if (Get_High32of64(yabs)  >= 0x408ae000)
  67     {
  68         QUAD tmp1, onethird;
  69         Set_Value64C(tmp1, pio2_hi_Hi, pio2_hi_Lo);
  70         Set_Value64C(onethird, onethird_Hi, onethird_Lo);
  71
  72         tmp1 = IEEEDPAdd(tmp1,IEEEDPDiv(IEEEDPSub(onethird, ysquared),ycubed));
  73
  74         if (is_eq(yabs,y)) /* arg has positive sign */
  75         {
  76             SetSR(0, Zero_Bit | Negative_Bit | Overflow_Bit);
  77             return tmp1;
  78         }
  79         else
  80         {
  81             SetSR(Negative_Bit, Zero_Bit | Negative_Bit | Overflow_Bit);
  82             OR64QC(tmp1, IEEEDPSign_Mask_Hi, IEEEDPSign_Mask_Lo);
  83             return tmp1;
  84         }
  85     }
  86
  87     /*
  88           atan(x >= 128) = pi/2 - 1/x + 1/(3*x^3) -1/(5*x^5)
  89         = pi/2 + (1/3*x^2 - x^4 - 1/5) / x^5
  90     */
  91     if (Get_High32of64(yabs) >= 0x40600000)
  92     {
  93         QUAD tmp1, onethird, onefifth;
  94         Set_Value64C(tmp1, pio2_hi_Hi, pio2_hi_Lo);
  95         Set_Value64C(onethird, onethird_Hi, onethird_Lo);
  96         Set_Value64C(onefifth, onefifth_Hi, onefifth_Lo);
  97         ypow5 = IEEEDPMul(ycubed, ysquared);
  98
  99         tmp1 = IEEEDPAdd
 100         (
 101             tmp1,
 102             IEEEDPDiv(
 103             IEEEDPSub(
 104             IEEEDPSub(
 105             IEEEDPMul(onethird,ysquared),
 106             IEEEDPMul(ysquared,ysquared)),onefifth),ypow5 )
 107         );
 108
 109         if (is_eq(yabs,y)) /* arg has positive sign */
 110         {
 111             SetSR(0, Zero_Bit | Negative_Bit | Overflow_Bit);
 112             return tmp1;
 113         }
 114         else
 115         {
 116             SetSR(Negative_Bit, Zero_Bit | Negative_Bit | Overflow_Bit);
 117             OR64QC(tmp1, IEEEDPSign_Mask_Hi, IEEEDPSign_Mask_Lo);
 118             return tmp1;
 119         }
 120     }
 121
 122     /* atan(x) = asin(x/sqrt(1+x^2)) */
 123
 124     return IEEEDPAsin(IEEEDPDiv(y,IEEEDPSqrt(IEEEDPAdd(one,ysquared))));
 125
 126     AROS_LIBFUNC_EXIT
 127 }