workbench/libs/mathieeedoubtrans/ieeedpasin.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 "mathieeedoubtrans_intern.h"
  17
  18 /*
  19     FUNCTION
  20       Calculate the arcus sine of the IEEE double precision number
  21
  22     RESULT
  23       IEEE double precision floating point number
  24
  25       flags:
  26         zero     : result is zero
  27         negative : result is negative
  28         overflow : argument is out of range
  29
  30     NOTES
  31
  32     EXAMPLE
  33
  34     BUGS
  35
  36     SEE ALSO
  37
  38     INTERNALS
  39 */
  40
  41 AROS_LHQUAD1(double, IEEEDPAsin,
  42     AROS_LHAQUAD(double, x, D0, D1),
  43     struct MathIeeeDoubTransBase *, MathIeeeDoubTransBase, 19, MathIeeeDoubTrans
  44 )
  45 {
  46     AROS_LIBFUNC_INIT
  47
  48     QUAD t,w,p,q,c,r,s;
  49     int hx,ix;
  50
  51     /*
  52         if this is a 32bit -compiler we need to define some variables, otherwise
  53         these are available as 64bit constants
  54     */
  55 #   if defined AROS_64BIT_TYPE || defined __GNUC__
  56 #   else
  57 #       undef two_64
  58 #       undef pio2_lo_64
  59 #       undef one_64
  60         const QUAD two_64     = two,
  61                    pio2_lo_64 = pio2_lo,
  62                    one_64     = oneC;
  63 #   endif
  64
  65     hx = Get_High32of64(x);
  66     ix = hx & 0x7fffffff;
  67     if(ix >= 0x3ff00000)
  68     {    /* |x| >= 1 */
  69         if(( (ix-0x3ff00000) | Get_Low32of64(x) ) == 0)
  70         {  /* |x|==1 -> asin(1) = pi/2 */
  71             /*
  72                 if this is a 32bit -compiler we need to define some variables, otherwise
  73                 these are available as 64bit constants
  74             */
  75             if (hx > 0)
  76             {
  77 #               if defined AROS_64BIT_TYPE || defined __GNUC__
  78 #               else
  79 #                   undef pio2_hi_64
  80                     QUAD pio2_hi_64 = pio2_hi;
  81 #               endif
  82
  83                 return pio2_hi_64;
  84             }
  85             else
  86             {
  87 #               if defined AROS_64BIT_TYPE || defined __GNUC__
  88 #               else
  89 #                   undef neg_pio2_hi_64
  90                     QUAD neg_pio2_hi_64 = neg_pio2_hi;
  91 #               endif
  92
  93                 return neg_pio2_hi_64;
  94             }
  95         }
  96         else  /* asin(>1) = NAN  */
  97         {
  98             /*
  99                 if this is a 32bit -compiler we need to define some variables,
 100                 otherwise these are available as 64bit constants
 101             */
 102 #           if defined AROS_64BIT_TYPE || defined __GNUC__
 103 #           else
 104 #               undef IEEEDPNAN_64
 105                 QUAD IEEEDPNAN_64;
 106                 Set_Value64C(IEEEDPNAN_64, IEEEDPNAN_Hi, IEEEDPNAN_Lo);
 107 #           endif
 108
 109             return IEEEDPNAN_64;         /* acos(>1)= NAN */
 110         } /* else */
 111     } /* if */
 112
 113     {
 114         /*
 115             if this is a 32bit -compiler we need to define some variables,
 116             otherwise these are available as 64bit constants
 117         */
 118
 119 #       if defined AROS_64BIT_TYPE || defined __GNUC__
 120 #       else
 121 #           undef pS0_64
 122 #           undef pS1_64
 123 #           undef pS2_64
 124 #           undef pS3_64
 125 #           undef pS4_64
 126 #           undef pS5_64
 127 #           undef qS1_64
 128 #           undef qS2_64
 129 #           undef qS3_64
 130 #           undef qS4_64
 131
 132 #           undef pio2_hi_64
 133
 134             /* Ok, let's define some constants */
 135             const QUAD pS0_64 = pS0,
 136                        pS1_64 = pS1,
 137                        pS2_64 = pS2,
 138                        pS3_64 = pS3,
 139                        pS4_64 = pS4,
 140                        pS5_64 = pS5;
 141             const QUAD qS1_64 = qS1,
 142                        qS2_64 = qS2,
 143                        qS3_64 = qS3,
 144                        qS4_64 = qS4;
 145             const QUAD pio2_hi_64 = pio2_hi;
 146 #       endif
 147
 148         if(ix < 0x3fe00000)
 149         {         /* |x| < 0.5 */
 150             if(ix <= 0x3e400000) return x;/* if|x|<2**-27 -> asin(x)=x */
 151
 152             Set_Value64(t, IEEEDPMul(x, x));
 153             Set_Value64
 154             (
 155                 p,
 156                 IEEEDPMul(t, IEEEDPAdd(pS0_64,
 157                 IEEEDPMul(t, IEEEDPAdd(pS1_64,
 158                 IEEEDPMul(t, IEEEDPAdd(pS2_64,
 159                 IEEEDPMul(t, IEEEDPAdd(pS3_64,
 160                 IEEEDPMul(t, IEEEDPAdd(pS4_64,
 161                 IEEEDPMul(t, pS5_64)))))))))))
 162             );
 163
 164             Set_Value64
 165             (
 166                 q,
 167                 IEEEDPAdd(one_64,
 168                 IEEEDPMul(t, IEEEDPAdd(qS1_64,
 169                 IEEEDPMul(t, IEEEDPAdd(qS2_64,
 170                 IEEEDPMul(t, IEEEDPAdd(qS3_64,
 171                 IEEEDPMul(t, qS4_64))))))))
 172             );
 173             Set_Value64(w, IEEEDPDiv(p, q));
 174             Set_Value64(t, IEEEDPAdd(x, IEEEDPMul(x, w)));
 175
 176             return t;
 177         } /* if */
 178         {
 179             /*
 180                 if this is a 32bit -compiler we need to define some variables,
 181                 otherwise these are available as 64bit constants
 182             */
 183 #           if defined AROS_64BIT_TYPE || defined __GNUC__
 184 #           else
 185 #               undef one_64
 186 #               undef onehalf_64
 187 #               undef two_64
 188 #               undef pio2_hi_64
 189 #               undef pio2_lo_64
 190 #               undef pio4_hi_64
 191                 const QUAD one_64     = oneC,
 192                            onehalf_64 = onehalf,
 193                            two_64     = two,
 194                            pio2_hi_64 = pio2_hi,
 195                            pio2_lo_64 = pio2_lo,
 196                            pio4_hi_64 = pio4_hi;
 197 #           endif
 198             /* 1 > |x| >= 0.5 */
 199             AND64QC
 200             (
 201                 x,
 202                 (IEEEDPMantisse_Mask_Hi | IEEEDPExponent_Mask_Hi),
 203                 (IEEEDPMantisse_Mask_Lo | IEEEDPExponent_Mask_Lo)
 204             );
 205             Set_Value64(w, IEEEDPSub(one_64, x));
 206             Set_Value64(t, IEEEDPMul(w, onehalf_64));
 207             Set_Value64
 208             (
 209                 p,
 210                 IEEEDPMul(t, IEEEDPAdd(pS0_64,
 211                 IEEEDPMul(t, IEEEDPAdd(pS1_64,
 212                 IEEEDPMul(t, IEEEDPAdd(pS2_64,
 213                 IEEEDPMul(t, IEEEDPAdd(pS3_64,
 214                 IEEEDPMul(t, IEEEDPAdd(pS4_64,
 215                 IEEEDPMul(t, pS5_64)))))))))))
 216             );
 217
 218             Set_Value64
 219             (
 220                 q,
 221                 IEEEDPAdd(one_64,
 222                 IEEEDPMul(t, IEEEDPAdd(qS1_64,
 223                 IEEEDPMul(t, IEEEDPAdd(qS2_64,
 224                 IEEEDPMul(t, IEEEDPAdd(qS3_64,
 225                 IEEEDPMul(t, qS4_64))))))))
 226             );
 227             Set_Value64(s, IEEEDPSqrt(t));
 228
 229             if(ix >= 0x3fef3333) /* if |x| > 0.975 */
 230             {
 231                 Set_Value64(w, IEEEDPDiv(p,q));
 232                 Set_Value64
 233                 (
 234                     t,
 235                     IEEEDPSub
 236                     (
 237                         pio2_hi_64, IEEEDPSub
 238                         (
 239                             IEEEDPMul
 240                             (
 241                                 two_64, IEEEDPAdd
 242                                 (
 243                                     s, IEEEDPMul(s,w)
 244                                 )
 245                             ),
 246                             pio2_lo_64
 247                         )
 248                     )
 249                 );
 250             }
 251             else
 252             {
 253                 Set_Value64(w,s);
 254                 AND64QC(w,0xffffffff,0x0);
 255                 Set_Value64
 256                 (
 257                     c, IEEEDPDiv(IEEEDPSub(t, IEEEDPMul(w, w)), IEEEDPAdd(s,w))
 258                 );
 259                 Set_Value64(r, IEEEDPDiv(p,q));
 260                 Set_Value64
 261                 (
 262                     p,
 263                     IEEEDPSub
 264                     (
 265                         IEEEDPMul(IEEEDPMul(two_64, s), r),
 266                         IEEEDPSub(pio2_lo_64, IEEEDPMul(two_64, c))
 267                     )
 268                 );
 269                 Set_Value64(q, IEEEDPSub(pio4_hi_64,IEEEDPMul(two_64,w)));
 270                 Set_Value64(t, IEEEDPSub(pio4_hi_64,IEEEDPSub(p,q)));
 271             } /* else */
 272             if (hx > 0)
 273             {
 274                 return t;
 275             }
 276             else
 277             {
 278                 OR64QC(t, IEEEDPSign_Mask_Hi, IEEEDPSign_Mask_Lo);
 279                 return t;
 280             } /* else */
 281         }
 282     }
 283
 284     AROS_LIBFUNC_EXIT
 285 }