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
31 static const floatx80 floatx80_negone = packFloatx80(1, 0x3fff, BX_CONST64(0x8000000000000000));
32 static const floatx80 floatx80_neghalf = packFloatx80(1, 0x3ffe, BX_CONST64(0x8000000000000000));
33 static const float128 float128_ln2 =
34 packFloat128(BX_CONST64(0x3ffe62e42fefa39e), BX_CONST64(0xf35793c7673007e6));
35
36 #ifdef BETTER_THAN_PENTIUM
37
38 #define LN2_SIG_HI BX_CONST64(0xb17217f7d1cf79ab)
39 #define LN2_SIG_LO BX_CONST64(0xc9e3b39800000000) /* 96 bit precision */
40
41 #else
42
43 #define LN2_SIG_HI BX_CONST64(0xb17217f7d1cf79ab)
44 #define LN2_SIG_LO BX_CONST64(0xc000000000000000) /* 67-bit precision */
45
46 #endif
47
48 #define EXP_ARR_SIZE 15
49
50 static float128 exp_arr[EXP_ARR_SIZE] =
51 {
52 PACK_FLOAT_128(0x3fff000000000000, 0x0000000000000000), /* 1 */
53 PACK_FLOAT_128(0x3ffe000000000000, 0x0000000000000000), /* 2 */
54 PACK_FLOAT_128(0x3ffc555555555555, 0x5555555555555555), /* 3 */
55 PACK_FLOAT_128(0x3ffa555555555555, 0x5555555555555555), /* 4 */
56 PACK_FLOAT_128(0x3ff8111111111111, 0x1111111111111111), /* 5 */
57 PACK_FLOAT_128(0x3ff56c16c16c16c1, 0x6c16c16c16c16c17), /* 6 */
58 PACK_FLOAT_128(0x3ff2a01a01a01a01, 0xa01a01a01a01a01a), /* 7 */
59 PACK_FLOAT_128(0x3fefa01a01a01a01, 0xa01a01a01a01a01a), /* 8 */
60 PACK_FLOAT_128(0x3fec71de3a556c73, 0x38faac1c88e50017), /* 9 */
61 PACK_FLOAT_128(0x3fe927e4fb7789f5, 0xc72ef016d3ea6679), /* 10 */
62 PACK_FLOAT_128(0x3fe5ae64567f544e, 0x38fe747e4b837dc7), /* 11 */
63 PACK_FLOAT_128(0x3fe21eed8eff8d89, 0x7b544da987acfe85), /* 12 */
64 PACK_FLOAT_128(0x3fde6124613a86d0, 0x97ca38331d23af68), /* 13 */
65 PACK_FLOAT_128(0x3fda93974a8c07c9, 0xd20badf145dfa3e5), /* 14 */
66 PACK_FLOAT_128(0x3fd6ae7f3e733b81, 0xf11d8656b0ee8cb0) /* 15 */
67 };
68
69 extern float128 EvalPoly(float128 x, float128 *arr, int n, float_status_t &status);
70
71 /* required -1 < x < 1 */
poly_exp(float128 x,float_status_t & status)72 static float128 poly_exp(float128 x, float_status_t &status)
73 {
74 /*
75 // 2 3 4 5 6 7 8 9
76 // x x x x x x x x x
77 // e - 1 ~ x + --- + --- + --- + --- + --- + --- + --- + --- + ...
78 // 2! 3! 4! 5! 6! 7! 8! 9!
79 //
80 // 2 3 4 5 6 7 8
81 // x x x x x x x x
82 // = x [ 1 + --- + --- + --- + --- + --- + --- + --- + --- + ... ]
83 // 2! 3! 4! 5! 6! 7! 8! 9!
84 //
85 // 8 8
86 // -- 2k -- 2k+1
87 // p(x) = > C * x q(x) = > C * x
88 // -- 2k -- 2k+1
89 // k=0 k=0
90 //
91 // x
92 // e - 1 ~ x * [ p(x) + x * q(x) ]
93 //
94 */
95 float128 t = EvalPoly(x, exp_arr, EXP_ARR_SIZE, status);
96 return float128_mul(t, x, status);
97 }
98
99 // =================================================
100 // x
101 // FX2M1 Compute 2 - 1
102 // =================================================
103
104 //
105 // Uses the following identities:
106 //
107 // 1. ----------------------------------------------------------
108 // x x*ln(2)
109 // 2 = e
110 //
111 // 2. ----------------------------------------------------------
112 // 2 3 4 5 n
113 // x x x x x x x
114 // e = 1 + --- + --- + --- + --- + --- + ... + --- + ...
115 // 1! 2! 3! 4! 5! n!
116 //
117
f2xm1(floatx80 a,float_status_t & status)118 floatx80 f2xm1(floatx80 a, float_status_t &status)
119 {
120 Bit64u zSig0, zSig1, zSig2;
121
122 // handle unsupported extended double-precision floating encodings
123 if (floatx80_is_unsupported(a))
124 {
125 float_raise(status, float_flag_invalid);
126 return floatx80_default_nan;
127 }
128
129 Bit64u aSig = extractFloatx80Frac(a);
130 Bit32s aExp = extractFloatx80Exp(a);
131 int aSign = extractFloatx80Sign(a);
132
133 if (aExp == 0x7FFF) {
134 if ((Bit64u) (aSig<<1))
135 return propagateFloatx80NaN(a, status);
136
137 return (aSign) ? floatx80_negone : a;
138 }
139
140 if (aExp == 0) {
141 if (aSig == 0) return a;
142 float_raise(status, float_flag_denormal | float_flag_inexact);
143 normalizeFloatx80Subnormal(aSig, &aExp, &aSig);
144
145 tiny_argument:
146 mul128By64To192(LN2_SIG_HI, LN2_SIG_LO, aSig, &zSig0, &zSig1, &zSig2);
147 if (0 < (Bit64s) zSig0) {
148 shortShift128Left(zSig0, zSig1, 1, &zSig0, &zSig1);
149 --aExp;
150 }
151 return
152 roundAndPackFloatx80(80, aSign, aExp, zSig0, zSig1, status);
153 }
154
155 float_raise(status, float_flag_inexact);
156
157 if (aExp < 0x3FFF)
158 {
159 if (aExp < FLOATX80_EXP_BIAS-68)
160 goto tiny_argument;
161
162 /* ******************************** */
163 /* using float128 for approximation */
164 /* ******************************** */
165
166 float128 x = floatx80_to_float128(a, status);
167 x = float128_mul(x, float128_ln2, status);
168 x = poly_exp(x, status);
169 return float128_to_floatx80(x, status);
170 }
171 else
172 {
173 if (a.exp == 0xBFFF && ! (aSig<<1))
174 return floatx80_neghalf;
175
176 return a;
177 }
178 }
179