![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to softfloat_decimal.c CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [HATARI the Atari ST Emulator] / hatari / src / cpu / softfloat|
Return to softfloat_decimal.c CVS log | Up to [HATARI the Atari ST Emulator] / hatari / src / cpu / softfloat |
1.1 root 1: /*============================================================================
2:
3: This C source file is an extension to the SoftFloat IEC/IEEE Floating-point
4: Arithmetic Package, Release 2a.
5:
6: =============================================================================*/
7:
8: #include <stdint.h>
9:
10: #include "sysconfig.h"
11: #include "sysdeps.h"
12:
13: #define DECIMAL_LOG 0
14:
15: #if DECIMAL_LOG
16: #define decimal_log write_log
17: #else
18: #define decimal_log(fmt, ...)
19: #endif
20:
21: #include "softfloat.h"
22: #include "softfloat-macros.h"
23: #include "softfloat/softfloat-specialize.h"
24:
25: /*----------------------------------------------------------------------------
26: | Methods for converting decimal floats to binary extended precision floats.
27: *----------------------------------------------------------------------------*/
28:
29: static void round128to64(flag aSign, int32_t *aExp, uint64_t *aSig0, uint64_t *aSig1, float_status *status)
30: {
31: flag increment;
32: int32_t zExp;
33: uint64_t zSig0, zSig1;
34:
35: zExp = *aExp;
36: zSig0 = *aSig0;
37: zSig1 = *aSig1;
38:
39: increment = ( (int64_t) zSig1 < 0 );
40: if (status->float_rounding_mode != float_round_nearest_even) {
41: if (status->float_rounding_mode == float_round_to_zero) {
42: increment = 0;
43: } else {
44: if (aSign) {
45: increment = (status->float_rounding_mode == float_round_down) && zSig1;
46: } else {
47: increment = (status->float_rounding_mode == float_round_up) && zSig1;
48: }
49: }
50: }
51:
52: if (increment) {
53: ++zSig0;
54: if (zSig0 == 0) {
55: ++zExp;
56: zSig0 = LIT64(0x8000000000000000);
57: } else {
58: zSig0 &= ~ (((uint64_t) (zSig1<<1) == 0) & (status->float_rounding_mode == float_round_nearest_even));
59: }
60: } else {
61: if ( zSig0 == 0 ) zExp = 0;
62: }
63:
64: *aExp = zExp;
65: *aSig0 = zSig0;
66: *aSig1 = 0;
67: }
68:
69: static void mul128by128round(int32_t *aExp, uint64_t *aSig0, uint64_t *aSig1, int32_t bExp, uint64_t bSig0, uint64_t bSig1, float_status *status)
70: {
71: int32_t zExp;
72: uint64_t zSig0, zSig1, zSig2, zSig3;
73:
74: zExp = *aExp;
75: zSig0 = *aSig0;
76: zSig1 = *aSig1;
77:
78: round128to64(0, &bExp, &bSig0, &bSig1, status);
79:
80: zExp += bExp - 0x3FFE;
81: mul128To256(zSig0, zSig1, bSig0, bSig1, &zSig0, &zSig1, &zSig2, &zSig3);
82: zSig1 |= (zSig2 | zSig3) != 0;
83: if ( 0 < (int64_t) zSig0 ) {
84: shortShift128Left( zSig0, zSig1, 1, &zSig0, &zSig1 );
85: --zExp;
86: }
87: *aExp = zExp;
88: *aSig0 = zSig0;
89: *aSig1 = zSig1;
90:
91: round128to64(0, aExp, aSig0, aSig1, status);
92: }
93:
94: static void mul128by128(int32_t *aExp, uint64_t *aSig0, uint64_t *aSig1, int32_t bExp, uint64_t bSig0, uint64_t bSig1)
95: {
96: int32_t zExp;
97: uint64_t zSig0, zSig1, zSig2, zSig3;
98:
99: zExp = *aExp;
100: zSig0 = *aSig0;
101: zSig1 = *aSig1;
102:
103: zExp += bExp - 0x3FFE;
104: mul128To256(zSig0, zSig1, bSig0, bSig1, &zSig0, &zSig1, &zSig2, &zSig3);
105: zSig1 |= (zSig2 | zSig3) != 0;
106: if ( 0 < (int64_t) zSig0 ) {
107: shortShift128Left( zSig0, zSig1, 1, &zSig0, &zSig1 );
108: --zExp;
109: }
110: *aExp = zExp;
111: *aSig0 = zSig0;
112: *aSig1 = zSig1;
113: }
114:
115: static void div128by128(int32_t *paExp, uint64_t *paSig0, uint64_t *paSig1, int32_t bExp, uint64_t bSig0, uint64_t bSig1)
116: {
117: int32_t zExp, aExp;
118: uint64_t zSig0, zSig1, aSig0, aSig1;
119: uint64_t rem0, rem1, rem2, rem3, term0, term1, term2, term3;
120:
121: aExp = *paExp;
122: aSig0 = *paSig0;
123: aSig1 = *paSig1;
124:
125: zExp = aExp - bExp + 0x3FFE;
126: if ( le128( bSig0, bSig1, aSig0, aSig1 ) ) {
127: shift128Right( aSig0, aSig1, 1, &aSig0, &aSig1 );
128: ++zExp;
129: }
130: zSig0 = estimateDiv128To64( aSig0, aSig1, bSig0 );
131: mul128By64To192( bSig0, bSig1, zSig0, &term0, &term1, &term2 );
132: sub192( aSig0, aSig1, 0, term0, term1, term2, &rem0, &rem1, &rem2 );
133: while ( (int64_t) rem0 < 0 ) {
134: --zSig0;
135: add192( rem0, rem1, rem2, 0, bSig0, bSig1, &rem0, &rem1, &rem2 );
136: }
137: zSig1 = estimateDiv128To64( rem1, rem2, bSig0 );
138: if ( ( zSig1 & 0x3FFF ) <= 4 ) {
139: mul128By64To192( bSig0, bSig1, zSig1, &term1, &term2, &term3 );
140: sub192( rem1, rem2, 0, term1, term2, term3, &rem1, &rem2, &rem3 );
141: while ( (int64_t) rem1 < 0 ) {
142: --zSig1;
143: add192( rem1, rem2, rem3, 0, bSig0, bSig1, &rem1, &rem2, &rem3 );
144: }
145: zSig1 |= ( ( rem1 | rem2 | rem3 ) != 0 );
146: }
147:
148: *paExp = zExp;
149: *paSig0 = zSig0;
150: *paSig1 = zSig1;
151: }
152:
153: #if 0
154:
155: void tentoint128(flag mSign, flag eSign, int32_t *aExp, uint64_t *aSig0, uint64_t *aSig1, int32_t scale, float_status *status)
156: {
157: int32_t mExp;
158: uint64_t mSig0, mSig1;
159:
160:
161: *aExp = 0x3FFF;
162: *aSig0 = LIT64(0x8000000000000000);
163: *aSig1 = 0;
164:
165: mExp = 0x4002;
166: mSig0 = LIT64(0xA000000000000000);
167: mSig1 = 0;
168:
169:
170: while (scale) {
171: if (scale & 1) {
172: mul128by128round(aExp, aSig0, aSig1, mExp, mSig0, mSig1, status);
173: }
174: mul128by128(&mExp, &mSig0, &mSig1, mExp, mSig0, mSig1);
175: scale >>= 1;
176: }
177: }
178:
179: #else
180:
181: static void tentoint128(flag mSign, flag eSign, int32_t *aExp, uint64_t *aSig0, uint64_t *aSig1, int32_t scale, float_status *status)
182: {
183: int8_t save_rounding_mode;
184: int32_t mExp;
185: uint64_t mSig0, mSig1;
186:
187: save_rounding_mode = status->float_rounding_mode;
188: switch (status->float_rounding_mode) {
189: case float_round_nearest_even:
190: break;
191: case float_round_down:
192: if (mSign != eSign) {
193: set_float_rounding_mode(float_round_up, status);
194: }
195: break;
196: case float_round_up:
197: if (mSign != eSign) {
198: set_float_rounding_mode(float_round_down, status);
199: }
200: break;
201: case float_round_to_zero:
202: if (eSign == 0) {
203: set_float_rounding_mode(float_round_down, status);
204: } else {
205: set_float_rounding_mode(float_round_up, status);
206: }
207: break;
208: default:
209: break;
210: }
211:
212: *aExp = 0x3FFF;
213: *aSig0 = LIT64(0x8000000000000000);
214: *aSig1 = 0;
215:
216: mExp = 0x4002;
217: mSig0 = LIT64(0xA000000000000000);
218: mSig1 = 0;
219:
220: while (scale) {
221: if (scale & 1) {
222: mul128by128round(aExp, aSig0, aSig1, mExp, mSig0, mSig1, status);
223: }
224: mul128by128(&mExp, &mSig0, &mSig1, mExp, mSig0, mSig1);
225: scale >>= 1;
226: }
227:
228: set_float_rounding_mode(save_rounding_mode, status);
229: }
230:
231: #endif
232:
233: static int64_t tentointdec(int32_t scale)
234: {
235: uint64_t decM, decX;
236:
237: decX = 1;
238: decM = 10;
239:
240: while (scale) {
241: if (scale & 1) {
242: decX *= decM;
243: }
244: decM *= decM;
245: scale >>= 1;
246: }
247:
248: return decX;
249: }
250:
251:
252: static int64_t float128toint64(flag zSign, int32_t zExp, uint64_t zSig0, uint64_t zSig1, float_status *status)
253: {
254: int8_t roundingMode;
255: flag roundNearestEven, increment;
256: int64_t z;
257:
258: shift128RightJamming(zSig0, zSig1, 0x403E - zExp, &zSig0, &zSig1);
259:
260: roundingMode = status->float_rounding_mode;
261: roundNearestEven = (roundingMode == float_round_nearest_even);
262: increment = ((int64_t)zSig1 < 0);
263: if (!roundNearestEven) {
264: if (roundingMode == float_round_to_zero) {
265: increment = 0;
266: } else {
267: if (zSign) {
268: increment = (roundingMode == float_round_down ) && zSig1;
269: } else {
270: increment = (roundingMode == float_round_up ) && zSig1;
271: }
272: }
273: }
274: if (increment) {
275: ++zSig0;
276: zSig0 &= ~ (((uint64_t)(zSig1<<1) == 0) & roundNearestEven);
277: }
278: z = zSig0;
279: if (zSig1) float_raise(float_flag_inexact, status);
280: return z;
281: }
282:
283: static int32_t getDecimalExponent(int32_t aExp, uint64_t aSig)
284: {
285: flag zSign;
286: int32_t zExp, shiftCount;
287: uint64_t zSig0, zSig1;
288:
289: if (aSig == 0 || aExp == 0x3FFF) {
290: return 0;
291: }
292: if (aExp < 0) {
293: return -4932;
294: }
295:
296: aSig ^= LIT64(0x8000000000000000);
297: aExp -= 0x3FFF;
298: zSign = (aExp < 0);
299: aExp = zSign ? -aExp : aExp;
300: shiftCount = 31 - countLeadingZeros32(aExp);
301: zExp = 0x3FFF + shiftCount;
302:
303: if (shiftCount < 0) {
304: shortShift128Left(aSig, 0, -shiftCount, &zSig0, &zSig1);
305: } else {
306: shift128Right(aSig, 0, shiftCount, &zSig0, &zSig1);
307: aSig = (uint64_t)aExp << (63 - shiftCount);
308: if (zSign) {
309: sub128(aSig, 0, zSig0, zSig1, &zSig0, &zSig1);
310: } else {
311: add128(aSig, 0, zSig0, zSig1, &zSig0, &zSig1);
312: }
313: }
314:
315: shiftCount = countLeadingZeros64(zSig0);
316: shortShift128Left(zSig0, zSig1, shiftCount, &zSig0, &zSig1);
317: zExp -= shiftCount;
318: mul128by128(&zExp, &zSig0, &zSig1, 0x3FFD, LIT64(0x9A209A84FBCFF798), LIT64(0x8F8959AC0B7C9178));
319:
320: shiftCount = 0x403E - zExp;
321: shift128RightJamming(zSig0, zSig1, shiftCount, &zSig0, &zSig1);
322:
323: if ((int64_t)zSig1 < 0) {
324: ++zSig0;
325: zSig0 &= ~(((int64_t)(zSig1<<1) == 0) & 1);
326: }
327:
328: zExp = zSign ? -zSig0 : zSig0;
329:
330: return zExp;
331: }
332:
333: /*----------------------------------------------------------------------------
334: | Decimal to binary
335: *----------------------------------------------------------------------------*/
336:
337: floatx80 floatdecimal_to_floatx80(floatx80 a, float_status *status)
338: {
339: flag decSign, zSign, decExpSign;
340: int32_t decExp, zExp, xExp, shiftCount;
341: uint64_t decSig, zSig0, zSig1, xSig0, xSig1;
342:
343: decSign = extractFloatx80Sign(a);
344: decExp = extractFloatx80Exp(a);
345: decSig = extractFloatx80Frac(a);
346:
347: if (decExp == 0x7FFF) return a;
348:
349: if (decExp == 0 && decSig == 0) return a;
350:
351: decExpSign = (decExp >> 14) & 1;
352: decExp &= 0x3FFF;
353:
354: shiftCount = countLeadingZeros64( decSig );
355: zExp = 0x403E - shiftCount;
356: zSig0 = decSig << shiftCount;
357: zSig1 = 0;
358: zSign = decSign;
359:
360: tentoint128(decSign, decExpSign, &xExp, &xSig0, &xSig1, decExp, status);
361:
362: if (decExpSign) {
363: div128by128(&zExp, &zSig0, &zSig1, xExp, xSig0, xSig1);
364: } else {
365: mul128by128(&zExp, &zSig0, &zSig1, xExp, xSig0, xSig1);
366: }
367:
368: if (zSig1) float_raise(float_flag_decimal, status);
369: round128to64(zSign, &zExp, &zSig0, &zSig1, status);
370:
371: return packFloatx80( zSign, zExp, zSig0 );
372:
373: }
374:
375: /*----------------------------------------------------------------------------
376: | Binary to decimal
377: *----------------------------------------------------------------------------*/
378:
379: floatx80 floatx80_to_floatdecimal(floatx80 a, int32_t *k, float_status *status)
380: {
381: flag aSign, decSign;
382: int32_t aExp, decExp, zExp, xExp;
383: uint64_t aSig, decSig, decX, zSig0, zSig1, xSig0, xSig1;
384: flag ictr, lambda;
385: int32_t kfactor, ilog, iscale, len;
386:
387: aSign = extractFloatx80Sign(a);
388: aExp = extractFloatx80Exp(a);
389: aSig = extractFloatx80Frac(a);
390:
391: if (aExp == 0x7FFF) {
392: if ((uint64_t) (aSig<<1)) return propagateFloatx80NaNOneArg(a, status);
393: return a;
394: }
395:
396: if (aExp == 0) {
397: if (aSig == 0) return packFloatx80(aSign, 0, 0);
398: normalizeFloatx80Subnormal(aSig, &aExp, &aSig);
399: }
400:
401: kfactor = *k;
402:
403: ilog = getDecimalExponent(aExp, aSig);
404:
405: ictr = 0;
406:
407: try_again:
408: decimal_log(_T("ILOG = %i\n"), ilog);
409:
410: if (kfactor > 0) {
411: if (kfactor > 17) {
412: kfactor = 17;
413: float_raise(float_flag_invalid, status);
414: }
415: len = kfactor;
416: } else {
417: len = ilog + 1 - kfactor;
418: if (len > 17) {
419: len = 17;
420: }
421: if (len < 1) {
422: len = 1;
423: }
424: if (kfactor > ilog) {
425: ilog = kfactor;
426: decimal_log(_T("ILOG is kfactor = %i\n"), ilog);
427: }
428: }
429:
430: decimal_log(_T("LEN = %i\n"),len);
431:
432: lambda = 0;
433: iscale = ilog + 1 - len;
434:
435: if (iscale < 0) {
436: lambda = 1;
437: iscale = -iscale;
438: }
439:
440: decimal_log(_T("ISCALE = %i, LAMBDA = %i\n"),iscale, lambda);
441:
442: tentoint128(lambda, 0, &xExp, &xSig0, &xSig1, iscale, status);
443:
444: decimal_log(_T("AFTER tentoint128: zExp = %04x, zSig0 = %16llx, zSig1 = %16llx\n"), xExp, xSig0, xSig1);
445:
446: zExp = aExp;
447: zSig0 = aSig;
448: zSig1 = 0;
449:
450: if (lambda) {
451: mul128by128(&zExp, &zSig0, &zSig1, xExp, xSig0, xSig1);
452: } else {
453: div128by128(&zExp, &zSig0, &zSig1, xExp, xSig0, xSig1);
454: }
455:
456: decimal_log(_T("BEFORE: zExp = %04x, zSig0 = %16llx, zSig1 = %16llx\n"),zExp,zSig0,zSig1);
457:
458: decSig = float128toint64(aSign, zExp, zSig0, zSig1, status);
459:
460: decimal_log(_T("AFTER: decSig = %llu\n"),decSig);
461:
462: if (ictr == 0) {
463:
464: decX = tentointdec(len - 1);
465:
466: if (decSig < decX) { // z < x
467: ilog -= 1;
468: ictr = 1;
469: goto try_again;
470: }
471:
472: decX *= 10;
473:
474: if (decSig > decX) { // z > x
475: ilog += 1;
476: ictr = 1;
477: goto try_again;
478: }
479: }
480:
481: decSign = aSign;
482: decExp = (ilog < 0) ? -ilog : ilog;
483: if (decExp > 999) {
484: float_raise(float_flag_invalid, status);
485: }
486: if (ilog < 0) decExp |= 0x4000;
487:
488: *k = len;
489:
490: return packFloatx80(decSign, decExp, decSig);
491: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.