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