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