workbench/libs/mathieeedoubtrans/ieeedpsin.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       Calculate the sine of a given IEEE double precision number in radians
  11
  12     RESULT
  13       IEEE double precision floating point number
  14
  15       flags:
  16         zero     : result is zero
  17         negative : result is negative
  18         overflow : 0
  19
  20     NOTES
  21
  22     EXAMPLE
  23
  24     BUGS
  25
  26     SEE ALSO
  27
  28     INTERNALS
  29       Algorithm for Calculation of sin(y):
  30       <code>
  31          z    = floor ( |y| / pi );
  32          y_1  = |y| - z * pi;        => 0 <= y_1 < pi
  33
  34          if (y_1 > pi/2 ) then y_1 = pi - y_1;
  35
  36          => 0 <= y_1 < pi/2
  37
  38          Res = y - y^3/3! + y^5/5! - y^7/7! + y^9/9! - y^11/11! =
  39              = y(1+y^2(-1/3!+y^2(1/5!+y^2(-1/7!+y^2(1/9!-1/11!y^2)))));
  40
  41          if (y < 0)
  42            Res = -Res;
  43
  44          if (z was an odd number)
  45            Res = -Res;
  46       </code>
  47
  48     HISTORY
  49 */
  50
  51 AROS_LHQUAD1(double, IEEEDPSin,
  52     AROS_LHAQUAD(double, y, D0, D1),
  53     struct MathIeeeDoubTransBase *, MathIeeeDoubTransBase, 6, MathIeeeDoubTrans
  54 )
  55 {
  56     AROS_LIBFUNC_INIT
  57
  58     QUAD z,Res,ysquared,yabs, Qtmp;
  59     AND64C
  60     (
  61         yabs, y,
  62         (IEEEDPMantisse_Mask_Hi | IEEEDPExponent_Mask_Hi),
  63         (IEEEDPMantisse_Mask_Lo | IEEEDPExponent_Mask_Lo)
  64     );
  65
  66     if (is_eqC(yabs, IEEEDPPInfty_Hi, IEEEDPPInfty_Lo))
  67     {
  68         SetSR(Overflow_Bit, Zero_Bit | Negative_Bit | Overflow_Bit);
  69         Set_Value64C(Res, IEEEDPNAN_Hi, IEEEDPNAN_Lo);
  70         return Res;
  71     }
  72
  73     z = IEEEDPFloor(IEEEDPDiv(yabs, pi));
  74     Qtmp  = IEEEDPMul(z,pi);
  75     OR64QC(Qtmp, IEEEDPSign_Mask_Hi, IEEEDPSign_Mask_Lo) /* Qtmp=-Qtmp */
  76     yabs = IEEEDPAdd(yabs, Qtmp);
  77     if (is_greaterC(yabs, pio2_Hi, pio2_Lo))
  78     {
  79         OR64QC(Qtmp, IEEEDPSign_Mask_Hi, IEEEDPSign_Mask_Lo) /* Qtmp=-Qtmp */
  80         yabs  =IEEEDPAdd(pi, yabs);
  81     }
  82     ysquared = IEEEDPMul(yabs,yabs);
  83     Res = IEEEDPMul(yabs,
  84           IEEEDPAdd(sinf1,
  85           IEEEDPMul(ysquared,
  86           IEEEDPAdd(sinf2,
  87           IEEEDPMul(ysquared,
  88           IEEEDPAdd(sinf3,
  89           IEEEDPMul(ysquared,
  90           IEEEDPAdd(sinf4,
  91           IEEEDPMul(ysquared,
  92           IEEEDPAdd(sinf5,
  93           IEEEDPMul(ysquared,
  94           IEEEDPAdd(sinf6,
  95           IEEEDPMul(ysquared,
  96           IEEEDPAdd(sinf7,
  97           IEEEDPMul(ysquared, sinf8)))))))))))))));
  98
  99     if (is_lessSC(y, 0x0, 0x0) )
 100     {
 101         XOR64QC(Res, IEEEDPSign_Mask_Hi, IEEEDPSign_Mask_Lo);
 102     }
 103
 104     if (TRUE == intern_IEEEDPisodd(z))
 105     {
 106         XOR64QC(Res, IEEEDPSign_Mask_Hi, IEEEDPSign_Mask_Lo);
 107     }
 108
 109     if (is_eqC(Res, 0x0, 0x0))
 110     {
 111         SetSR(Zero_Bit, Zero_Bit | Negative_Bit | Overflow_Bit);
 112         return Res;
 113     }
 114
 115     if (is_lessSC(Res, 0x0 ,0x0))
 116     {
 117         SetSR(Negative_Bit, Zero_Bit | Negative_Bit | Overflow_Bit);
 118     }
 119
 120     return Res;
 121
 122     AROS_LIBFUNC_EXIT
 123 }