workbench/libs/mathieeesingtrans/ieeespasin.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 arcussin of the given number
  21
  22     RESULT
  23       IEEE single precision floating point number
  24
  25
  26       flags:
  27         zero     : result is zero
  28         negative : result is negative
  29         overflow : fnum &lt; -1  or  fnum &gt; 1
  30
  31     NOTES
  32
  33     EXAMPLE
  34
  35     BUGS
  36
  37     SEE ALSO
  38
  39     INTERNALS
  40
  41     HISTORY
  42 */
  43
  44 AROS_LH1(float, IEEESPAsin,
  45     AROS_LHA(float, y, D0),
  46     struct Library *, MathIeeeSingTransBase, 19, MathIeeeSingTrans
  47 )
  48 {
  49     AROS_LIBFUNC_INIT
  50
  51     /* 1> |x| >= 0.5 */
  52     LONG t,w,p,q,c,r,s,ix;
  53     ix = y & (IEEESPMantisse_Mask | IEEESPExponent_Mask); /* ix = |fnum| */
  54
  55     if (one == ix) /* |y| = 1 -> result = +-(pi/2) */
  56     {
  57         return (pio2 | (y & IEEESPSign_Mask  ));
  58     }
  59
  60     if (1 == IEEESPCmp(ix,one)) /* |y| > 1 */
  61     {
  62         SetSR(Overflow_Bit, Zero_Bit | Overflow_Bit | Negative_Bit);
  63         return 0x7fffffff;
  64     }
  65
  66     /* error: 1 ulp (unit in the last place)*/
  67     if ( -1 == IEEESPCmp(ix, onehalf)) /* |y| < 0.5 */
  68     {
  69         t = IEEESPMul(y, y);
  70         p = IEEESPMul(t, IEEESPAdd(pS0,
  71             IEEESPMul(t, IEEESPAdd(pS1,
  72             IEEESPMul(t, IEEESPAdd(pS2,
  73             IEEESPMul(t, IEEESPAdd(pS3,
  74             IEEESPMul(t, IEEESPAdd(pS4,
  75             IEEESPMul(t, pS5) ))))))))));
  76         q = IEEESPAdd(one,
  77             IEEESPMul(t, IEEESPAdd(qS1,
  78             IEEESPMul(t, IEEESPAdd(qS2,
  79             IEEESPMul(t, IEEESPAdd(qS3,
  80             IEEESPMul(t, qS4))))))));
  81         w = IEEESPDiv(p, q);
  82
  83         return IEEESPAdd(y, IEEESPMul(y, w));
  84     }
  85
  86
  87     w = IEEESPSub(one, ix) ; /* w = 1 - fnum ; y = 1-x */
  88     t = IEEESPMul(w, onehalf);  /* t = w / 2    ; z = y/2 */
  89     p = IEEESPMul(t,IEEESPAdd(pS0,
  90         IEEESPMul(t,IEEESPAdd(pS1,
  91         IEEESPMul(t,IEEESPAdd(pS2,
  92         IEEESPMul(t,IEEESPAdd(pS3,
  93         IEEESPMul(t,IEEESPAdd(pS4,
  94         IEEESPMul(t,pS5)))))))))));
  95     q = IEEESPAdd(one,
  96         IEEESPMul(t,IEEESPAdd(qS1,
  97         IEEESPMul(t,IEEESPAdd(qS2,
  98         IEEESPMul(t,IEEESPAdd(qS3,
  99         IEEESPMul(t,qS4))))))));
 100     s = IEEESPSqrt(t);      /* s = sqrt(t)  ; s = sqrt(z) */
 101
 102     if(1 == IEEESPCmp(ix, 0x3f79999a /*0.975*/ )) /* |fnum| > 0.975 */
 103     {
 104         /*error: 2 ulp (4 ulp when |fnum| close to 1) */
 105         w = IEEESPDiv(p, q); /* w = p / q; */
 106         /* res = pi/2-(2*(s+s*w)) */
 107         t = IEEESPSub(pio2, IEEESPMul(two,IEEESPAdd(s,IEEESPMul(s,w))));
 108     }
 109     else
 110     {
 111         /* error: 2 ulp */
 112         w = s;
 113         c = IEEESPDiv(IEEESPSub(t,IEEESPMul(w,w)),IEEESPAdd(s,w)); /* c=(t-w*w)/(s+w) */
 114         r = IEEESPDiv(p,q);
 115         p = IEEESPAdd(IEEESPAdd(c,c),IEEESPMul(IEEESPAdd(s,s),r));
 116         q = IEEESPSub(pio4, IEEESPAdd(w,w));
 117         t = IEEESPSub(pio4, IEEESPSub(p,q));
 118     }
 119
 120     return (t | (y & IEEESPSign_Mask )) ;
 121
 122     AROS_LIBFUNC_EXIT
 123 }