workbench/libs/mathtrans/spasin.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 "mathtrans_intern.h"
  17
  18 /*
  19     FUNCTION
  20       Calculate arcussin of the given number
  21
  22     RESULT
  23       Motorola fast floating point number
  24
  25       flags:
  26         zero     : result is zero
  27         negative : result is negative
  28         overflow : fnum &lt; -1  or  fnum &gt; 1
  29
  30     NOTES
  31
  32     EXAMPLE
  33
  34     BUGS
  35
  36     SEE ALSO
  37
  38     INTERNALS
  39
  40     HISTORY
  41 */
  42
  43 AROS_LH1(float, SPAsin,
  44     AROS_LHA(float, fnum1, D0),
  45     struct Library *, MathTransBase, 19, MathTrans
  46 )
  47 {
  48     AROS_LIBFUNC_INIT
  49
  50     /* 1> |x| >= 0.5 */
  51     LONG t,w,p,q,c,r,s,ix;
  52     ix = fnum1 & (FFPMantisse_Mask | FFPExponent_Mask); /* ix = |fnum| */
  53
  54     if ((LONG)one == ix) /* |fnum1| = 1 -> result = +-(pi/2) */
  55     {
  56         return (pio2 | (fnum1 & FFPSign_Mask  ));
  57     }
  58
  59     if (1 == SPCmp(ix,one)) /* |fnum1| > 1 */
  60     {
  61         SetSR(Overflow_Bit, Zero_Bit | Overflow_Bit | Negative_Bit);
  62         return -1;
  63     }
  64
  65     /* error: 1 ulp (unit in the last place)*/
  66     if (-1 == SPCmp(ix,onehalf)) /* |fnum1| < 0.5 */
  67     {
  68         if (-1 == SPCmp(ix,0xb89b6736)) /* |fnum1| < 70422/10000000 */
  69         {
  70             return fnum1;
  71         }
  72
  73         t = SPMul(fnum1, fnum1);
  74         p = SPMul(t, SPAdd(pS0,
  75             SPMul(t, SPAdd(pS1,
  76             SPMul(t, SPAdd(pS2,
  77             SPMul(t, SPAdd(pS3,
  78             SPMul(t, SPAdd(pS4,
  79             SPMul(t, pS5)))))))))));
  80         q = SPAdd(one,
  81             SPMul(t, SPAdd(qS1,
  82             SPMul(t, SPAdd(qS2,
  83             SPMul(t, SPAdd(qS3,
  84             SPMul(t, qS4))))))));
  85         w = SPDiv(q, p);
  86
  87         return SPAdd(fnum1, SPMul(fnum1, w));
  88     }
  89
  90
  91     w = SPSub(ix, one) ; /* w = 1 - fnum ; y = 1-x */
  92     t = SPMul(w, onehalf);  /* t = w / 2    ; z = y/2 */
  93     p = SPMul(t,SPAdd(pS0,
  94         SPMul(t,SPAdd(pS1,
  95         SPMul(t,SPAdd(pS2,
  96         SPMul(t,SPAdd(pS3,
  97         SPMul(t,SPAdd(pS4,
  98         SPMul(t,pS5)))))))))));
  99     q = SPAdd(one,
 100         SPMul(t,SPAdd(qS1,
 101         SPMul(t,SPAdd(qS2,
 102         SPMul(t,SPAdd(qS3,
 103         SPMul(t,qS4))))))));
 104     s = SPSqrt(t);      /* s = sqrt(t)  ; s = sqrt(z) */
 105
 106     if(1 == SPCmp(ix, 0xf9999a40 /*0.975*/ )) /* |fnum| > 0.975 */
 107     {
 108         /*error: 2 ulp (4 ulp when |fnum| close to 1) */
 109         w = SPDiv(q,p); /* w = p / q; */
 110         /* res = pi/2-(2*(s+s*w)) */
 111         t = SPSub(SPMul(two,SPAdd(s,SPMul(s,w))),pio2);
 112         t = SPAdd(t, 0x8000002b); /* for better accuracy */
 113     }
 114     else
 115     {
 116         /* error: 2 ulp */
 117         w = s;
 118         c = SPDiv(SPAdd(s,w),SPSub(SPMul(w,w),t)); /* c=(t-w*w)/(s+w) */
 119         r = SPDiv(q,p);
 120         p = SPAdd(SPAdd(c,c),SPMul(SPAdd(s,s),r));
 121         q = SPSub(SPAdd(w,w) ,pio4);
 122         t = SPSub(SPSub(q,p) ,pio4);
 123     }
 124
 125     return (t | (fnum1 & FFPSign_Mask )) ;
 126
 127     AROS_LIBFUNC_EXIT
 128 }