workbench/libs/mathieeesingtrans/ieeespacos.c

   1 /*
   2     Copyright © 1995-2003, The AROS Development Team. All rights reserved.
   3     $Id$
   4 */
   5 /*
   6  * ====================================================
   7  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   8  *
   9  * Developed at SunSoft, a Sun Microsystems, Inc. business.
  10  * Permission to use, copy, modify, and distribute this
  11  * software is freely granted, provided that this notice
  12  * is preserved.
  13  * ====================================================
  14  */
  15
  16 #include "mathieeesingtrans_intern.h"
  17
  18 /*
  19     FUNCTION
  20       Calculate arcuscos of the given number
  21
  22     RESULT
  23
  24       IEEE single precision floating point number
  25
  26
  27       flags:
  28       zero     : result is zero
  29       negative : 0 (not possible)
  30       overflow : y &lt; -1  or  y &gt; 1
  31
  32     NOTES
  33
  34     EXAMPLE
  35
  36     BUGS
  37
  38     SEE ALSO
  39
  40     INTERNALS
  41
  42     HISTORY
  43 */
  44
  45 AROS_LH1(float, IEEESPAcos,
  46     AROS_LHA(float, y, D0),
  47     struct Library *, MathIeeeSingTransBase, 20, MathIeeeSingTrans
  48 )
  49 {
  50     AROS_LIBFUNC_INIT
  51
  52     /* 1> |x| >= 0.5 */
  53     LONG z,p,q,r,w,s,c,ix,df;
  54     ix = y & (IEEESPMantisse_Mask | IEEESPExponent_Mask); /* ix = |y|  */
  55
  56     z = IEEESPCmp(ix,one);
  57
  58     if (1==z) /* |y| > 1 */
  59     {
  60         SetSR(Overflow_Bit, Zero_Bit | Overflow_Bit | Negative_Bit);
  61         return 0x7fffffff;
  62     }
  63
  64     if (0==z) /* |y| = 1 */
  65     {
  66         if (y & IEEESPSign_Mask) /* y = -1 */ return pi;
  67         SetSR(Zero_Bit, Zero_Bit | Overflow_Bit | Negative_Bit);
  68         return 0;
  69     }
  70
  71     /* error: 1 ulp (unit in the last place) */
  72     if (-1 == IEEESPCmp(ix,onehalf)) /* |fnum1| < 0.5 */
  73     {
  74         z = IEEESPMul(y, y);
  75         p = IEEESPMul(z, IEEESPAdd(pS0,
  76             IEEESPMul(z, IEEESPAdd(pS1,
  77             IEEESPMul(z, IEEESPAdd(pS2,
  78             IEEESPMul(z, IEEESPAdd(pS3,
  79             IEEESPMul(z, IEEESPAdd(pS4,
  80             IEEESPMul(z, pS5)))))))))));
  81         q = IEEESPAdd(one,
  82             IEEESPMul(z, IEEESPAdd(qS1,
  83             IEEESPMul(z, IEEESPAdd(qS2,
  84             IEEESPMul(z, IEEESPAdd(qS3,
  85             IEEESPMul(z, qS4))))))));
  86         r = IEEESPDiv(p, q);
  87
  88         return (IEEESPSub(pio2,IEEESPAdd(y,IEEESPMul(y,r))));
  89     }
  90
  91     /* error: 1 ulp */
  92     if (y & IEEESPSign_Mask) /* y < -0.5 */
  93     {
  94         z = IEEESPMul(onehalf, IEEESPAdd(one, y));
  95         p = IEEESPMul(z, IEEESPAdd(pS0,
  96             IEEESPMul(z, IEEESPAdd(pS1,
  97             IEEESPMul(z, IEEESPAdd(pS2,
  98             IEEESPMul(z, IEEESPAdd(pS3,
  99             IEEESPMul(z, IEEESPAdd(pS4,
 100             IEEESPMul(z, pS5)))))))))));
 101         q = IEEESPAdd(one,
 102             IEEESPMul(z, IEEESPAdd(qS1,
 103             IEEESPMul(z, IEEESPAdd(qS2,
 104             IEEESPMul(z, IEEESPAdd(qS3,
 105             IEEESPMul(z, qS4))))))));
 106         s = IEEESPSqrt(z);
 107         r = IEEESPDiv(p,q); /* r = p/q; */
 108         w = IEEESPMul(r,s);
 109
 110         return IEEESPSub(pi, IEEESPMul(two, IEEESPAdd(s, w)));
 111     }
 112
 113     /* fnum1 > 0.5 */
 114     /* error : 1 ulp */
 115     z = IEEESPMul(onehalf, IEEESPSub(one, y));
 116     s = IEEESPSqrt(z);
 117     df = s;
 118     //df = df & 0xfff000ff;
 119     c = IEEESPDiv(IEEESPSub(z, IEEESPMul(df,df)), IEEESPAdd(df,s) );
 120     p = IEEESPMul(z, IEEESPAdd(pS0,
 121         IEEESPMul(z, IEEESPAdd(pS1,
 122         IEEESPMul(z, IEEESPAdd(pS2,
 123         IEEESPMul(z, IEEESPAdd(pS3,
 124         IEEESPMul(z, IEEESPAdd(pS4,
 125         IEEESPMul(z, pS5)))))))))));
 126     q = IEEESPAdd(one,
 127         IEEESPMul(z, IEEESPAdd(qS1,
 128         IEEESPMul(z, IEEESPAdd(qS2,
 129         IEEESPMul(z, IEEESPAdd(qS3,
 130         IEEESPMul(z, qS4))))))));
 131     r = IEEESPDiv(p, q);
 132     w = IEEESPAdd(c, IEEESPMul(r,s));
 133
 134     return IEEESPMul(two, IEEESPAdd(df,w));
 135
 136     AROS_LIBFUNC_EXIT
 137 }