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