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