![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to fyl2x.c CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [Previous NeXT emulator] / previous / src / softfloat|
Return to fyl2x.c CVS log | Up to [Previous NeXT emulator] / previous / src / softfloat |
1.1 root 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: float_raise(float_flag_invalid)
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: * Adapted for lib/softfloat in MESS by Hans Ostermeyer (03/2012)
25: * ==========================================================================*/
26:
27: #define FLOAT128
28:
29: #define USE_estimateDiv128To64
30: #include "mamesf.h"
31: #include "softfloat.h"
32: #include "fpu_constant.h"
33:
34: static const floatx80 floatx80_log10_2 = { 0x3ffd, 0x9a209a84fbcff798U };
35: static const floatx80 floatx80_ln_2 = { 0x3ffe, 0xb17217f7d1cf79acU };
36: static const floatx80 floatx80_one = { 0x3fff, 0x8000000000000000U };
37: static const floatx80 floatx80_default_nan = { floatx80_default_nan_high, floatx80_default_nan_low };
38:
39: #define packFloat_128(zHi, zLo) {(zHi), (zLo)}
40: #define PACK_FLOAT_128(hi,lo) packFloat_128(LIT64(hi),LIT64(lo))
41:
42: #define EXP_BIAS 0x3FFF
43:
44:
45: static const float128 float128_one =
46: packFloat_128(0x3fff000000000000U, 0x0000000000000000U);
47: static const float128 float128_two =
48: packFloat_128(0x4000000000000000U, 0x0000000000000000U);
49:
50: static const float128 float128_ln2inv2 =
51: packFloat_128(0x400071547652b82fU, 0xe1777d0ffda0d23aU);
52:
53: #define SQRT2_HALF_SIG 0xb504f333f9de6484U
54:
55: extern float128 OddPoly(float128 x, float128 *arr, unsigned n);
56:
57: #define L2_ARR_SIZE 9
58:
59: static float128 ln_arr[L2_ARR_SIZE] =
60: {
61: PACK_FLOAT_128(0x3fff000000000000, 0x0000000000000000), /* 1 */
62: PACK_FLOAT_128(0x3ffd555555555555, 0x5555555555555555), /* 3 */
63: PACK_FLOAT_128(0x3ffc999999999999, 0x999999999999999a), /* 5 */
64: PACK_FLOAT_128(0x3ffc249249249249, 0x2492492492492492), /* 7 */
65: PACK_FLOAT_128(0x3ffbc71c71c71c71, 0xc71c71c71c71c71c), /* 9 */
66: PACK_FLOAT_128(0x3ffb745d1745d174, 0x5d1745d1745d1746), /* 11 */
67: PACK_FLOAT_128(0x3ffb3b13b13b13b1, 0x3b13b13b13b13b14), /* 13 */
68: PACK_FLOAT_128(0x3ffb111111111111, 0x1111111111111111), /* 15 */
69: PACK_FLOAT_128(0x3ffae1e1e1e1e1e1, 0xe1e1e1e1e1e1e1e2) /* 17 */
70: };
71:
72: static float128 poly_ln(float128 x1)
73: {
74: /*
75: //
76: // 3 5 7 9 11 13 15
77: // 1+u u u u u u u u
78: // 1/2 ln --- ~ u + --- + --- + --- + --- + ---- + ---- + ---- =
79: // 1-u 3 5 7 9 11 13 15
80: //
81: // 2 4 6 8 10 12 14
82: // u u u u u u u
83: // = u * [ 1 + --- + --- + --- + --- + ---- + ---- + ---- ] =
84: // 3 5 7 9 11 13 15
85: //
86: // 3 3
87: // -- 4k -- 4k+2
88: // p(u) = > C * u q(u) = > C * u
89: // -- 2k -- 2k+1
90: // k=0 k=0
91: //
92: // 1+u 2
93: // 1/2 ln --- ~ u * [ p(u) + u * q(u) ]
94: // 1-u
95: //
96: */
97: return OddPoly(x1, ln_arr, L2_ARR_SIZE);
98: }
99:
100: /* required sqrt(2)/2 < x < sqrt(2) */
101: static float128 poly_l2(float128 x)
102: {
103: /* using float128 for approximation */
104: float128 x_p1 = float128_add(x, float128_one);
105: float128 x_m1 = float128_sub(x, float128_one);
106: x = float128_div(x_m1, x_p1);
107: x = poly_ln(x);
108: x = float128_mul(x, float128_ln2inv2);
109: return x;
110: }
111:
112: static float128 poly_l2p1(float128 x)
113: {
114: /* using float128 for approximation */
115: float128 x_p2 = float128_add(x, float128_two);
116: x = float128_div(x, x_p2);
117: x = poly_ln(x);
118: x = float128_mul(x, float128_ln2inv2);
119: return x;
120: }
121:
122: // =================================================
123: // FYL2X Compute y * log (x)
124: // 2
125: // =================================================
126:
127: //
128: // Uses the following identities:
129: //
130: // 1. ----------------------------------------------------------
131: // ln(x)
132: // log (x) = -------, ln (x*y) = ln(x) + ln(y)
133: // 2 ln(2)
134: //
135: // 2. ----------------------------------------------------------
136: // 1+u x-1
137: // ln (x) = ln -----, when u = -----
138: // 1-u x+1
139: //
140: // 3. ----------------------------------------------------------
141: // 3 5 7 2n+1
142: // 1+u u u u u
143: // ln ----- = 2 [ u + --- + --- + --- + ... + ------ + ... ]
144: // 1-u 3 5 7 2n+1
145: //
146:
147: static floatx80 fyl2x(floatx80 a, floatx80 b)
148: {
149: uint64_t aSig = extractFloatx80Frac(a);
150: int32_t aExp = extractFloatx80Exp(a);
151: int aSign = extractFloatx80Sign(a);
152: uint64_t bSig = extractFloatx80Frac(b);
153: int32_t bExp = extractFloatx80Exp(b);
154: int bSign = extractFloatx80Sign(b);
155:
156: int zSign = bSign ^ 1;
157:
158: if (aExp == 0x7FFF) {
159: if ((uint64_t) (aSig<<1)
160: || ((bExp == 0x7FFF) && (uint64_t) (bSig<<1)))
161: {
162: return propagateFloatx80NaN(a, b);
163: }
164: if (aSign)
165: {
166: invalid:
167: float_raise(float_flag_invalid);
168: return floatx80_default_nan;
169: }
170: else {
171: if (bExp == 0) {
172: if (bSig == 0) goto invalid;
173: float_raise(float_flag_denormal);
174: }
175: return packFloatx80(bSign, 0x7FFF, floatx80_default_infinity_low);
176: }
177: }
178: if (bExp == 0x7FFF)
179: {
180: if ((uint64_t) (bSig<<1)) return propagateFloatx80NaN(a, b);
181: if (aSign && (uint64_t)(aExp | aSig)) goto invalid;
182: if (aSig && (aExp == 0))
183: float_raise(float_flag_denormal);
184: if (aExp < 0x3FFF) {
185: return packFloatx80(zSign, 0x7FFF, floatx80_default_infinity_low);
186: }
187: if (aExp == 0x3FFF && ((uint64_t) (aSig<<1) == 0)) goto invalid;
188: return packFloatx80(bSign, 0x7FFF, floatx80_default_infinity_low);
189: }
190: if (aExp == 0) {
191: if (aSig == 0) {
192: if ((bExp | bSig) == 0) goto invalid;
193: float_raise(float_flag_divbyzero);
194: return packFloatx80(zSign, 0x7FFF, floatx80_default_infinity_low);
195: }
196: if (aSign) goto invalid;
197: float_raise(float_flag_denormal);
198: normalizeFloatx80Subnormal(aSig, &aExp, &aSig);
199: }
200: if (aSign) goto invalid;
201: if (bExp == 0) {
202: if (bSig == 0) {
203: if (aExp < 0x3FFF) return packFloatx80(zSign, 0, 0);
204: return packFloatx80(bSign, 0, 0);
205: }
206: float_raise(float_flag_denormal);
207: normalizeFloatx80Subnormal(bSig, &bExp, &bSig);
208: }
209: if (aExp == 0x3FFF && ((uint64_t) (aSig<<1) == 0))
210: return packFloatx80(bSign, 0, 0);
211:
212: float_raise(float_flag_inexact);
213:
214: int ExpDiff = aExp - 0x3FFF;
215: aExp = 0;
216: if (aSig >= SQRT2_HALF_SIG) {
217: ExpDiff++;
218: aExp--;
219: }
220:
221: /* ******************************** */
222: /* using float128 for approximation */
223: /* ******************************** */
224:
225: uint64_t zSig0, zSig1;
226: shift128Right(aSig<<1, 0, 16, &zSig0, &zSig1);
227: float128 x = packFloat128(0, aExp+0x3FFF, zSig0, zSig1);
228: x = poly_l2(x);
229: x = float128_add(x, int64_to_float128((int64_t) ExpDiff));
230: return floatx80_mul(b, float128_to_floatx80(x));
231: }
232:
233: // =================================================
234: // FYL2XP1 Compute y * log (x + 1)
235: // 2
236: // =================================================
237:
238: //
239: // Uses the following identities:
240: //
241: // 1. ----------------------------------------------------------
242: // ln(x)
243: // log (x) = -------
244: // 2 ln(2)
245: //
246: // 2. ----------------------------------------------------------
247: // 1+u x
248: // ln (x+1) = ln -----, when u = -----
249: // 1-u x+2
250: //
251: // 3. ----------------------------------------------------------
252: // 3 5 7 2n+1
253: // 1+u u u u u
254: // ln ----- = 2 [ u + --- + --- + --- + ... + ------ + ... ]
255: // 1-u 3 5 7 2n+1
256: //
257:
258: floatx80 fyl2xp1(floatx80 a, floatx80 b)
259: {
260: int32_t aExp, bExp;
261: uint64_t aSig, bSig, zSig0, zSig1, zSig2;
262: int aSign, bSign;
263:
264: aSig = extractFloatx80Frac(a);
265: aExp = extractFloatx80Exp(a);
266: aSign = extractFloatx80Sign(a);
267: bSig = extractFloatx80Frac(b);
268: bExp = extractFloatx80Exp(b);
269: bSign = extractFloatx80Sign(b);
270: int zSign = aSign ^ bSign;
271:
272: if (aExp == 0x7FFF) {
273: if ((uint64_t) (aSig<<1)
274: || ((bExp == 0x7FFF) && (uint64_t) (bSig<<1)))
275: {
276: return propagateFloatx80NaN(a, b);
277: }
278: if (aSign)
279: {
280: invalid:
281: float_raise(float_flag_invalid);
282: return floatx80_default_nan;
283: }
284: else {
285: if (bExp == 0) {
286: if (bSig == 0) goto invalid;
287: float_raise(float_flag_denormal);
288: }
289: return packFloatx80(bSign, 0x7FFF, floatx80_default_infinity_low);
290: }
291: }
292: if (bExp == 0x7FFF)
293: {
294: if ((uint64_t) (bSig<<1))
295: return propagateFloatx80NaN(a, b);
296:
297: if (aExp == 0) {
298: if (aSig == 0) goto invalid;
299: float_raise(float_flag_denormal);
300: }
301:
302: return packFloatx80(zSign, 0x7FFF, floatx80_default_infinity_low);
303: }
304: if (aExp == 0) {
305: if (aSig == 0) {
306: if (bSig && (bExp == 0)) float_raise(float_flag_denormal);
307: return packFloatx80(zSign, 0, 0);
308: }
309: float_raise(float_flag_denormal);
310: normalizeFloatx80Subnormal(aSig, &aExp, &aSig);
311: }
312: if (bExp == 0) {
313: if (bSig == 0) return packFloatx80(zSign, 0, 0);
314: float_raise(float_flag_denormal);
315: normalizeFloatx80Subnormal(bSig, &bExp, &bSig);
316: }
317:
318: float_raise(float_flag_inexact);
319:
320: if (aSign && aExp >= 0x3FFF)
321: return a;
322:
323: if (aExp >= 0x3FFC) // big argument
324: {
325: return fyl2x(floatx80_add(a, floatx80_one), b);
326: }
327:
328: // handle tiny argument
329: if (aExp < EXP_BIAS-70)
330: {
331: // first order approximation, return (a*b)/ln(2)
332: int32_t zExp = aExp + FLOAT_LN2INV_EXP - 0x3FFE;
333:
334: mul128By64To192(FLOAT_LN2INV_HI, FLOAT_LN2INV_LO, aSig, &zSig0, &zSig1, &zSig2);
335: if (0 < (int64_t) zSig0) {
336: shortShift128Left(zSig0, zSig1, 1, &zSig0, &zSig1);
337: --zExp;
338: }
339:
340: zExp = zExp + bExp - 0x3FFE;
341: mul128By64To192(zSig0, zSig1, bSig, &zSig0, &zSig1, &zSig2);
342: if (0 < (int64_t) zSig0) {
343: shortShift128Left(zSig0, zSig1, 1, &zSig0, &zSig1);
344: --zExp;
345: }
346:
347: return
348: roundAndPackFloatx80(80, aSign ^ bSign, zExp, zSig0, zSig1);
349: }
350:
351: /* ******************************** */
352: /* using float128 for approximation */
353: /* ******************************** */
354:
355: shift128Right(aSig<<1, 0, 16, &zSig0, &zSig1);
356: float128 x = packFloat128(aSign, aExp, zSig0, zSig1);
357: x = poly_l2p1(x);
358: return floatx80_mul(b, float128_to_floatx80(x));
359: }
360:
361: floatx80 floatx80_flognp1(floatx80 a)
362: {
363: return fyl2xp1(a, floatx80_ln_2);
364: }
365:
366: floatx80 floatx80_flogn(floatx80 a)
367: {
368: return fyl2x(a, floatx80_ln_2);
369: }
370:
371: floatx80 floatx80_flog2(floatx80 a)
372: {
373: return fyl2x(a, floatx80_one);
374: }
375:
376: floatx80 floatx80_flog10(floatx80 a)
377: {
378: return fyl2x(a, floatx80_log10_2);
379: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.