fpu/fpatan.cc

   1 /*============================================================================
   2 This source file is an extension to the SoftFloat IEC/IEEE Floating-point
   3 Arithmetic Package, Release 2b, written for Bochs (x86 achitecture simulator)
   4 floating point emulation.
   5
   6 THIS SOFTWARE IS DISTRIBUTED AS IS, FOR FREE.  Although reasonable effort has
   7 been made to avoid it, THIS SOFTWARE MAY CONTAIN FAULTS THAT WILL AT TIMES
   8 RESULT IN INCORRECT BEHAVIOR.  USE OF THIS SOFTWARE IS RESTRICTED TO PERSONS
   9 AND ORGANIZATIONS WHO CAN AND WILL TAKE FULL RESPONSIBILITY FOR ALL LOSSES,
  10 COSTS, OR OTHER PROBLEMS THEY INCUR DUE TO THE SOFTWARE, AND WHO FURTHERMORE
  11 EFFECTIVELY INDEMNIFY JOHN HAUSER AND THE INTERNATIONAL COMPUTER SCIENCE
  12 INSTITUTE (possibly via similar legal warning) AGAINST ALL LOSSES, COSTS, OR
  13 OTHER PROBLEMS INCURRED BY THEIR CUSTOMERS AND CLIENTS DUE TO THE SOFTWARE.
  14
  15 Derivative works are acceptable, even for commercial purposes, so long as
  16 (1) the source code for the derivative work includes prominent notice that
  17 the work is derivative, and (2) the source code includes prominent notice with
  18 these four paragraphs for those parts of this code that are retained.
  19 =============================================================================*/
  20
  21 /*============================================================================
  22  * Written for Bochs (x86 achitecture simulator) by
  23  *            Stanislav Shwartsman [sshwarts at sourceforge net]
  24  * ==========================================================================*/
  25
  26 #define FLOAT128
  27
  28 #include "softfloatx80.h"
  29 #include "softfloat-round-pack.h"
  30 #include "fpu_constant.h"
  31
  32 #define FPATAN_ARR_SIZE 11
  33
  34 static const float128 float128_one =
  35         packFloat128(BX_CONST64(0x3fff000000000000), BX_CONST64(0x0000000000000000));
  36 static const float128 float128_sqrt3 =
  37         packFloat128(BX_CONST64(0x3fffbb67ae8584ca), BX_CONST64(0xa73b25742d7078b8));
  38 static const floatx80 floatx80_pi  =
  39         packFloatx80(0, 0x4000, BX_CONST64(0xc90fdaa22168c235));
  40
  41 static const float128 float128_pi2 =
  42         packFloat128(BX_CONST64(0x3fff921fb54442d1), BX_CONST64(0x8469898CC5170416));
  43 static const float128 float128_pi4 =
  44         packFloat128(BX_CONST64(0x3ffe921fb54442d1), BX_CONST64(0x8469898CC5170416));
  45 static const float128 float128_pi6 =
  46         packFloat128(BX_CONST64(0x3ffe0c152382d736), BX_CONST64(0x58465BB32E0F580F));
  47
  48 static float128 atan_arr[FPATAN_ARR_SIZE] =
  49 {
  50     PACK_FLOAT_128(0x3fff000000000000, 0x0000000000000000), /*  1 */
  51     PACK_FLOAT_128(0xbffd555555555555, 0x5555555555555555), /*  3 */
  52     PACK_FLOAT_128(0x3ffc999999999999, 0x999999999999999a), /*  5 */
  53     PACK_FLOAT_128(0xbffc249249249249, 0x2492492492492492), /*  7 */
  54     PACK_FLOAT_128(0x3ffbc71c71c71c71, 0xc71c71c71c71c71c), /*  9 */
  55     PACK_FLOAT_128(0xbffb745d1745d174, 0x5d1745d1745d1746), /* 11 */
  56     PACK_FLOAT_128(0x3ffb3b13b13b13b1, 0x3b13b13b13b13b14), /* 13 */
  57     PACK_FLOAT_128(0xbffb111111111111, 0x1111111111111111), /* 15 */
  58     PACK_FLOAT_128(0x3ffae1e1e1e1e1e1, 0xe1e1e1e1e1e1e1e2), /* 17 */
  59     PACK_FLOAT_128(0xbffaaf286bca1af2, 0x86bca1af286bca1b), /* 19 */
  60     PACK_FLOAT_128(0x3ffa861861861861, 0x8618618618618618)  /* 21 */
  61 };
  62
  63 extern float128 OddPoly(float128 x, float128 *arr, unsigned n, float_status_t &status);
  64
  65 /* |x| < 1/4 */
  66 static float128 poly_atan(float128 x1, float_status_t &status)
  67 {
  68 /*
  69     //                 3     5     7     9     11     13     15     17
  70     //                x     x     x     x     x      x      x      x
  71     // atan(x) ~ x - --- + --- - --- + --- - ---- + ---- - ---- + ----
  72     //                3     5     7     9     11     13     15     17
  73     //
  74     //                 2     4     6     8     10     12     14     16
  75     //                x     x     x     x     x      x      x      x
  76     //   = x * [ 1 - --- + --- - --- + --- - ---- + ---- - ---- + ---- ]
  77     //                3     5     7     9     11     13     15     17
  78     //
  79     //           5                          5
  80     //          --       4k                --        4k+2
  81     //   p(x) = >  C  * x           q(x) = >  C   * x
  82     //          --  2k                     --  2k+1
  83     //          k=0                        k=0
  84     //
  85     //                            2
  86     //    atan(x) ~ x * [ p(x) + x * q(x) ]
  87     //
  88 */
  89     return OddPoly(x1, atan_arr, FPATAN_ARR_SIZE, status);
  90 }
  91
  92 // =================================================
  93 // FPATAN                  Compute y * log (x)
  94 //                                        2
  95 // =================================================
  96
  97 //
  98 // Uses the following identities:
  99 //
 100 // 1. ----------------------------------------------------------
 101 //
 102 //   atan(-x) = -atan(x)
 103 //
 104 // 2. ----------------------------------------------------------
 105 //
 106 //                             x + y
 107 //   atan(x) + atan(y) = atan -------, xy < 1
 108 //                             1-xy
 109 //
 110 //                             x + y
 111 //   atan(x) + atan(y) = atan ------- + PI, x > 0, xy > 1
 112 //                             1-xy
 113 //
 114 //                             x + y
 115 //   atan(x) + atan(y) = atan ------- - PI, x < 0, xy > 1
 116 //                             1-xy
 117 //
 118 // 3. ----------------------------------------------------------
 119 //
 120 //   atan(x) = atan(INF) + atan(- 1/x)
 121 //
 122 //                           x-1
 123 //   atan(x) = PI/4 + atan( ----- )
 124 //                           x+1
 125 //
 126 //                           x * sqrt(3) - 1
 127 //   atan(x) = PI/6 + atan( ----------------- )
 128 //                             x + sqrt(3)
 129 //
 130 // 4. ----------------------------------------------------------
 131 //                   3     5     7     9                 2n+1
 132 //                  x     x     x     x              n  x
 133 //   atan(x) = x - --- + --- - --- + --- - ... + (-1)  ------ + ...
 134 //                  3     5     7     9                 2n+1
 135 //
 136
 137 floatx80 fpatan(floatx80 a, floatx80 b, float_status_t &status)
 138 {
 139     // handle unsupported extended double-precision floating encodings
 140     if (floatx80_is_unsupported(a)) {
 141         float_raise(status, float_flag_invalid);
 142         return floatx80_default_nan;
 143     }
 144
 145     Bit64u aSig = extractFloatx80Frac(a);
 146     Bit32s aExp = extractFloatx80Exp(a);
 147     int aSign = extractFloatx80Sign(a);
 148     Bit64u bSig = extractFloatx80Frac(b);
 149     Bit32s bExp = extractFloatx80Exp(b);
 150     int bSign = extractFloatx80Sign(b);
 151
 152     int zSign = aSign ^ bSign;
 153
 154     if (bExp == 0x7FFF)
 155     {
 156         if ((Bit64u) (bSig<<1))
 157             return propagateFloatx80NaN(a, b, status);
 158
 159         if (aExp == 0x7FFF) {
 160             if ((Bit64u) (aSig<<1))
 161                 return propagateFloatx80NaN(a, b, status);
 162
 163             if (aSign) {   /* return 3PI/4 */
 164                 return roundAndPackFloatx80(80, bSign,
 165                         FLOATX80_3PI4_EXP, FLOAT_3PI4_HI, FLOAT_3PI4_LO, status);
 166             }
 167             else {         /* return  PI/4 */
 168                 return roundAndPackFloatx80(80, bSign,
 169                         FLOATX80_PI4_EXP, FLOAT_PI_HI, FLOAT_PI_LO, status);
 170             }
 171         }
 172
 173         if (aSig && (aExp == 0))
 174             float_raise(status, float_flag_denormal);
 175
 176         /* return PI/2 */
 177         return roundAndPackFloatx80(80, bSign, FLOATX80_PI2_EXP, FLOAT_PI_HI, FLOAT_PI_LO, status);
 178     }
 179     if (aExp == 0x7FFF)
 180     {
 181         if ((Bit64u) (aSig<<1))
 182             return propagateFloatx80NaN(a, b, status);
 183
 184         if (bSig && (bExp == 0))
 185             float_raise(status, float_flag_denormal);
 186
 187 return_PI_or_ZERO:
 188
 189         if (aSign) {   /* return PI */
 190             return roundAndPackFloatx80(80, bSign, FLOATX80_PI_EXP, FLOAT_PI_HI, FLOAT_PI_LO, status);
 191         } else {       /* return  0 */
 192             return packFloatx80(bSign, 0, 0);
 193         }
 194     }
 195     if (bExp == 0)
 196     {
 197         if (bSig == 0) {
 198              if (aSig && (aExp == 0)) float_raise(status, float_flag_denormal);
 199              goto return_PI_or_ZERO;
 200         }
 201
 202         float_raise(status, float_flag_denormal);
 203         normalizeFloatx80Subnormal(bSig, &bExp, &bSig);
 204     }
 205     if (aExp == 0)
 206     {
 207         if (aSig == 0)   /* return PI/2 */
 208             return roundAndPackFloatx80(80, bSign, FLOATX80_PI2_EXP, FLOAT_PI_HI, FLOAT_PI_LO, status);
 209
 210         float_raise(status, float_flag_denormal);
 211         normalizeFloatx80Subnormal(aSig, &aExp, &aSig);
 212     }
 213
 214     float_raise(status, float_flag_inexact);
 215
 216     /* |a| = |b| ==> return PI/4 */
 217     if (aSig == bSig && aExp == bExp)
 218         return roundAndPackFloatx80(80, bSign, FLOATX80_PI4_EXP, FLOAT_PI_HI, FLOAT_PI_LO, status);
 219
 220     /* ******************************** */
 221     /* using float128 for approximation */
 222     /* ******************************** */
 223
 224     float128 a128 = normalizeRoundAndPackFloat128(0, aExp-0x10, aSig, 0, status);
 225     float128 b128 = normalizeRoundAndPackFloat128(0, bExp-0x10, bSig, 0, status);
 226     float128 x;
 227     int swap = 0, add_pi6 = 0, add_pi4 = 0;
 228
 229     if (aExp > bExp || (aExp == bExp && aSig > bSig))
 230     {
 231         x = float128_div(b128, a128, status);
 232     }
 233     else {
 234         x = float128_div(a128, b128, status);
 235         swap = 1;
 236     }
 237
 238     Bit32s xExp = extractFloat128Exp(x);
 239
 240     if (xExp <= EXP_BIAS-40)
 241         goto approximation_completed;
 242
 243     if (x.hi >= BX_CONST64(0x3ffe800000000000)) // 3/4 < x < 1
 244     {
 245         /*
 246         arctan(x) = arctan((x-1)/(x+1)) + pi/4
 247         */
 248         float128 t1 = float128_sub(x, float128_one, status);
 249         float128 t2 = float128_add(x, float128_one, status);
 250         x = float128_div(t1, t2, status);
 251         add_pi4 = 1;
 252     }
 253     else
 254     {
 255         /* argument correction */
 256         if (xExp >= 0x3FFD)                     // 1/4 < x < 3/4
 257         {
 258             /*
 259             arctan(x) = arctan((x*sqrt(3)-1)/(x+sqrt(3))) + pi/6
 260             */
 261             float128 t1 = float128_mul(x, float128_sqrt3, status);
 262             float128 t2 = float128_add(x, float128_sqrt3, status);
 263             x = float128_sub(t1, float128_one, status);
 264             x = float128_div(x, t2, status);
 265             add_pi6 = 1;
 266         }
 267     }
 268
 269     x = poly_atan(x, status);
 270     if (add_pi6) x = float128_add(x, float128_pi6, status);
 271     if (add_pi4) x = float128_add(x, float128_pi4, status);
 272
 273 approximation_completed:
 274     if (swap) x = float128_sub(float128_pi2, x, status);
 275     floatx80 result = float128_to_floatx80(x, status);
 276     if (zSign) floatx80_chs(result);
 277     int rSign = extractFloatx80Sign(result);
 278     if (!bSign && rSign)
 279         return floatx80_add(result, floatx80_pi, status);
 280     if (bSign && !rSign)
 281         return floatx80_sub(result, floatx80_pi, status);
 282     return result;
 283 }