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