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