1 // Copyright 2011 The Go Authors. All rights reserved.
2 // Use of this source code is governed by a BSD-style
3 // license that can be found in the LICENSE file.
11 // An extFloat represents an extended floating-point number, with more
12 // precision than a float64. It does not try to save bits: the
13 // number represented by the structure is mant*(2^exp), with a negative
14 // sign if neg is true.
15 type extFloat
struct {
21 // Powers of ten taken from double-conversion library.
22 // https://code.google.com/p/double-conversion/
24 firstPowerOfTen
= -348
28 var smallPowersOfTen
= [...]extFloat
{
29 {1 << 63, -63, false}, // 1
30 {0xa << 60, -60, false}, // 1e1
31 {0x64 << 57, -57, false}, // 1e2
32 {0x3e8 << 54, -54, false}, // 1e3
33 {0x2710 << 50, -50, false}, // 1e4
34 {0x186a0 << 47, -47, false}, // 1e5
35 {0xf4240 << 44, -44, false}, // 1e6
36 {0x989680 << 40, -40, false}, // 1e7
39 var powersOfTen
= [...]extFloat
{
40 {0xfa8fd5a0081c0288, -1220, false}, // 10^-348
41 {0xbaaee17fa23ebf76, -1193, false}, // 10^-340
42 {0x8b16fb203055ac76, -1166, false}, // 10^-332
43 {0xcf42894a5dce35ea, -1140, false}, // 10^-324
44 {0x9a6bb0aa55653b2d, -1113, false}, // 10^-316
45 {0xe61acf033d1a45df, -1087, false}, // 10^-308
46 {0xab70fe17c79ac6ca, -1060, false}, // 10^-300
47 {0xff77b1fcbebcdc4f, -1034, false}, // 10^-292
48 {0xbe5691ef416bd60c, -1007, false}, // 10^-284
49 {0x8dd01fad907ffc3c, -980, false}, // 10^-276
50 {0xd3515c2831559a83, -954, false}, // 10^-268
51 {0x9d71ac8fada6c9b5, -927, false}, // 10^-260
52 {0xea9c227723ee8bcb, -901, false}, // 10^-252
53 {0xaecc49914078536d, -874, false}, // 10^-244
54 {0x823c12795db6ce57, -847, false}, // 10^-236
55 {0xc21094364dfb5637, -821, false}, // 10^-228
56 {0x9096ea6f3848984f, -794, false}, // 10^-220
57 {0xd77485cb25823ac7, -768, false}, // 10^-212
58 {0xa086cfcd97bf97f4, -741, false}, // 10^-204
59 {0xef340a98172aace5, -715, false}, // 10^-196
60 {0xb23867fb2a35b28e, -688, false}, // 10^-188
61 {0x84c8d4dfd2c63f3b, -661, false}, // 10^-180
62 {0xc5dd44271ad3cdba, -635, false}, // 10^-172
63 {0x936b9fcebb25c996, -608, false}, // 10^-164
64 {0xdbac6c247d62a584, -582, false}, // 10^-156
65 {0xa3ab66580d5fdaf6, -555, false}, // 10^-148
66 {0xf3e2f893dec3f126, -529, false}, // 10^-140
67 {0xb5b5ada8aaff80b8, -502, false}, // 10^-132
68 {0x87625f056c7c4a8b, -475, false}, // 10^-124
69 {0xc9bcff6034c13053, -449, false}, // 10^-116
70 {0x964e858c91ba2655, -422, false}, // 10^-108
71 {0xdff9772470297ebd, -396, false}, // 10^-100
72 {0xa6dfbd9fb8e5b88f, -369, false}, // 10^-92
73 {0xf8a95fcf88747d94, -343, false}, // 10^-84
74 {0xb94470938fa89bcf, -316, false}, // 10^-76
75 {0x8a08f0f8bf0f156b, -289, false}, // 10^-68
76 {0xcdb02555653131b6, -263, false}, // 10^-60
77 {0x993fe2c6d07b7fac, -236, false}, // 10^-52
78 {0xe45c10c42a2b3b06, -210, false}, // 10^-44
79 {0xaa242499697392d3, -183, false}, // 10^-36
80 {0xfd87b5f28300ca0e, -157, false}, // 10^-28
81 {0xbce5086492111aeb, -130, false}, // 10^-20
82 {0x8cbccc096f5088cc, -103, false}, // 10^-12
83 {0xd1b71758e219652c, -77, false}, // 10^-4
84 {0x9c40000000000000, -50, false}, // 10^4
85 {0xe8d4a51000000000, -24, false}, // 10^12
86 {0xad78ebc5ac620000, 3, false}, // 10^20
87 {0x813f3978f8940984, 30, false}, // 10^28
88 {0xc097ce7bc90715b3, 56, false}, // 10^36
89 {0x8f7e32ce7bea5c70, 83, false}, // 10^44
90 {0xd5d238a4abe98068, 109, false}, // 10^52
91 {0x9f4f2726179a2245, 136, false}, // 10^60
92 {0xed63a231d4c4fb27, 162, false}, // 10^68
93 {0xb0de65388cc8ada8, 189, false}, // 10^76
94 {0x83c7088e1aab65db, 216, false}, // 10^84
95 {0xc45d1df942711d9a, 242, false}, // 10^92
96 {0x924d692ca61be758, 269, false}, // 10^100
97 {0xda01ee641a708dea, 295, false}, // 10^108
98 {0xa26da3999aef774a, 322, false}, // 10^116
99 {0xf209787bb47d6b85, 348, false}, // 10^124
100 {0xb454e4a179dd1877, 375, false}, // 10^132
101 {0x865b86925b9bc5c2, 402, false}, // 10^140
102 {0xc83553c5c8965d3d, 428, false}, // 10^148
103 {0x952ab45cfa97a0b3, 455, false}, // 10^156
104 {0xde469fbd99a05fe3, 481, false}, // 10^164
105 {0xa59bc234db398c25, 508, false}, // 10^172
106 {0xf6c69a72a3989f5c, 534, false}, // 10^180
107 {0xb7dcbf5354e9bece, 561, false}, // 10^188
108 {0x88fcf317f22241e2, 588, false}, // 10^196
109 {0xcc20ce9bd35c78a5, 614, false}, // 10^204
110 {0x98165af37b2153df, 641, false}, // 10^212
111 {0xe2a0b5dc971f303a, 667, false}, // 10^220
112 {0xa8d9d1535ce3b396, 694, false}, // 10^228
113 {0xfb9b7cd9a4a7443c, 720, false}, // 10^236
114 {0xbb764c4ca7a44410, 747, false}, // 10^244
115 {0x8bab8eefb6409c1a, 774, false}, // 10^252
116 {0xd01fef10a657842c, 800, false}, // 10^260
117 {0x9b10a4e5e9913129, 827, false}, // 10^268
118 {0xe7109bfba19c0c9d, 853, false}, // 10^276
119 {0xac2820d9623bf429, 880, false}, // 10^284
120 {0x80444b5e7aa7cf85, 907, false}, // 10^292
121 {0xbf21e44003acdd2d, 933, false}, // 10^300
122 {0x8e679c2f5e44ff8f, 960, false}, // 10^308
123 {0xd433179d9c8cb841, 986, false}, // 10^316
124 {0x9e19db92b4e31ba9, 1013, false}, // 10^324
125 {0xeb96bf6ebadf77d9, 1039, false}, // 10^332
126 {0xaf87023b9bf0ee6b, 1066, false}, // 10^340
129 // floatBits returns the bits of the float64 that best approximates
130 // the extFloat passed as receiver. Overflow is set to true if
131 // the resulting float64 is ±Inf.
132 func (f
*extFloat
) floatBits(flt
*floatInfo
) (bits
uint64, overflow
bool) {
137 // Exponent too small.
138 if exp
< flt
.bias
+1 {
139 n
:= flt
.bias
+ 1 - exp
144 // Extract 1+flt.mantbits bits from the 64-bit mantissa.
145 mant
:= f
.mant
>> (63 - flt
.mantbits
)
146 if f
.mant
&(1<<(62-flt
.mantbits
)) != 0 {
151 // Rounding might have added a bit; shift down.
152 if mant
== 2<<flt
.mantbits
{
158 if exp
-flt
.bias
>= 1<<flt
.expbits
-1 {
161 exp
= 1<<flt
.expbits
- 1 + flt
.bias
163 } else if mant
&(1<<flt
.mantbits
) == 0 {
168 bits
= mant
& (uint64(1)<<flt
.mantbits
- 1)
169 bits |
= uint64((exp
-flt
.bias
)&(1<<flt
.expbits
-1)) << flt
.mantbits
171 bits |
= 1 << (flt
.mantbits
+ flt
.expbits
)
176 // AssignComputeBounds sets f to the floating point value
177 // defined by mant, exp and precision given by flt. It returns
178 // lower, upper such that any number in the closed interval
179 // [lower, upper] is converted back to the same floating point number.
180 func (f
*extFloat
) AssignComputeBounds(mant
uint64, exp
int, neg
bool, flt
*floatInfo
) (lower
, upper extFloat
) {
182 f
.exp
= exp
- int(flt
.mantbits
)
184 if f
.exp
<= 0 && mant
== (mant
>>uint(-f
.exp
))<<uint(-f
.exp
) {
186 f
.mant
>>= uint(-f
.exp
)
190 expBiased
:= exp
- flt
.bias
192 upper
= extFloat
{mant
: 2*f
.mant
+ 1, exp
: f
.exp
- 1, neg
: f
.neg
}
193 if mant
!= 1<<flt
.mantbits || expBiased
== 1 {
194 lower
= extFloat
{mant
: 2*f
.mant
- 1, exp
: f
.exp
- 1, neg
: f
.neg
}
196 lower
= extFloat
{mant
: 4*f
.mant
- 1, exp
: f
.exp
- 2, neg
: f
.neg
}
201 // Normalize normalizes f so that the highest bit of the mantissa is
202 // set, and returns the number by which the mantissa was left-shifted.
203 func (f
*extFloat
) Normalize() uint {
204 // bits.LeadingZeros64 would return 64
208 shift
:= bits
.LeadingZeros64(f
.mant
)
209 f
.mant
<<= uint(shift
)
214 // Multiply sets f to the product f*g: the result is correctly rounded,
215 // but not normalized.
216 func (f
*extFloat
) Multiply(g extFloat
) {
217 fhi
, flo
:= f
.mant
>>32, uint64(uint32(f
.mant
))
218 ghi
, glo
:= g
.mant
>>32, uint64(uint32(g
.mant
))
224 // f.mant*g.mant is fhi*ghi << 64 + (cross1+cross2) << 32 + flo*glo
225 f
.mant
= fhi
*ghi
+ (cross1
>> 32) + (cross2
>> 32)
226 rem
:= uint64(uint32(cross1
)) + uint64(uint32(cross2
)) + ((flo
* glo
) >> 32)
230 f
.mant
+= (rem
>> 32)
231 f
.exp
= f
.exp
+ g
.exp
+ 64
234 var uint64pow10
= [...]uint64{
235 1, 1e1
, 1e2
, 1e3
, 1e4
, 1e5
, 1e6
, 1e7
, 1e8
, 1e9
,
236 1e10
, 1e11
, 1e12
, 1e13
, 1e14
, 1e15
, 1e16
, 1e17
, 1e18
, 1e19
,
239 // AssignDecimal sets f to an approximate value mantissa*10^exp. It
240 // reports whether the value represented by f is guaranteed to be the
241 // best approximation of d after being rounded to a float64 or
242 // float32 depending on flt.
243 func (f
*extFloat
) AssignDecimal(mantissa
uint64, exp10
int, neg
bool, trunc
bool, flt
*floatInfo
) (ok
bool) {
244 const uint64digits
= 19
246 errors
:= 0 // An upper bound for error, computed in errorscale*ulp.
248 // the decimal number was truncated.
249 errors
+= errorscale
/ 2
256 // Multiply by powers of ten.
257 i
:= (exp10
- firstPowerOfTen
) / stepPowerOfTen
258 if exp10
< firstPowerOfTen || i
>= len(powersOfTen
) {
261 adjExp
:= (exp10
- firstPowerOfTen
) % stepPowerOfTen
263 // We multiply by exp%step
264 if adjExp
< uint64digits
&& mantissa
< uint64pow10
[uint64digits
-adjExp
] {
265 // We can multiply the mantissa exactly.
266 f
.mant
*= uint64pow10
[adjExp
]
270 f
.Multiply(smallPowersOfTen
[adjExp
])
271 errors
+= errorscale
/ 2
274 // We multiply by 10 to the exp - exp%step.
275 f
.Multiply(powersOfTen
[i
])
279 errors
+= errorscale
/ 2
282 shift
:= f
.Normalize()
285 // Now f is a good approximation of the decimal.
286 // Check whether the error is too large: that is, if the mantissa
287 // is perturbated by the error, the resulting float64 will change.
288 // The 64 bits mantissa is 1 + 52 bits for float64 + 11 extra bits.
290 // In many cases the approximation will be good enough.
291 denormalExp
:= flt
.bias
- 63
293 if f
.exp
<= denormalExp
{
294 // f.mant * 2^f.exp is smaller than 2^(flt.bias+1).
295 extrabits
= 63 - flt
.mantbits
+ 1 + uint(denormalExp
-f
.exp
)
297 extrabits
= 63 - flt
.mantbits
300 halfway
:= uint64(1) << (extrabits
- 1)
301 mant_extra
:= f
.mant
& (1<<extrabits
- 1)
303 // Do a signed comparison here! If the error estimate could make
304 // the mantissa round differently for the conversion to double,
305 // then we can't give a definite answer.
306 if int64(halfway
)-int64(errors
) < int64(mant_extra
) &&
307 int64(mant_extra
) < int64(halfway
)+int64(errors
) {
313 // Frexp10 is an analogue of math.Frexp for decimal powers. It scales
314 // f by an approximate power of ten 10^-exp, and returns exp10, so
315 // that f*10^exp10 has the same value as the old f, up to an ulp,
316 // as well as the index of 10^-exp in the powersOfTen table.
317 func (f
*extFloat
) frexp10() (exp10
, index
int) {
318 // The constants expMin and expMax constrain the final value of the
319 // binary exponent of f. We want a small integral part in the result
320 // because finding digits of an integer requires divisions, whereas
321 // digits of the fractional part can be found by repeatedly multiplying
325 // Find power of ten such that x * 10^n has a binary exponent
326 // between expMin and expMax.
327 approxExp10
:= ((expMin
+expMax
)/2 - f
.exp
) * 28 / 93 // log(10)/log(2) is close to 93/28.
328 i
:= (approxExp10
- firstPowerOfTen
) / stepPowerOfTen
331 exp
:= f
.exp
+ powersOfTen
[i
].exp
+ 64
341 // Apply the desired decimal shift on f. It will have exponent
342 // in the desired range. This is multiplication by 10^-exp10.
343 f
.Multiply(powersOfTen
[i
])
345 return -(firstPowerOfTen
+ i
*stepPowerOfTen
), i
348 // frexp10Many applies a common shift by a power of ten to a, b, c.
349 func frexp10Many(a
, b
, c
*extFloat
) (exp10
int) {
350 exp10
, i
:= c
.frexp10()
351 a
.Multiply(powersOfTen
[i
])
352 b
.Multiply(powersOfTen
[i
])
356 // FixedDecimal stores in d the first n significant digits
357 // of the decimal representation of f. It returns false
358 // if it cannot be sure of the answer.
359 func (f
*extFloat
) FixedDecimal(d
*decimalSlice
, n
int) bool {
367 panic("strconv: internal error: extFloat.FixedDecimal called with n == 0")
369 // Multiply by an appropriate power of ten to have a reasonable
370 // number to process.
372 exp10
, _
:= f
.frexp10()
374 shift
:= uint(-f
.exp
)
375 integer
:= uint32(f
.mant
>> shift
)
376 fraction
:= f
.mant
- (uint64(integer
) << shift
)
377 ε
:= uint64(1) // ε is the uncertainty we have on the mantissa of f.
379 // Write exactly n digits to d.
380 needed
:= n
// how many digits are left to write.
381 integerDigits
:= 0 // the number of decimal digits of integer.
382 pow10
:= uint64(1) // the power of ten by which f was scaled.
383 for i
, pow
:= 0, uint64(1); i
< 20; i
++ {
384 if pow
> uint64(integer
) {
391 if integerDigits
> needed
{
392 // the integral part is already large, trim the last digits.
393 pow10
= uint64pow10
[integerDigits
-needed
]
394 integer
/= uint32(pow10
)
395 rest
-= integer
* uint32(pow10
)
400 // Write the digits of integer: the digits of rest are omitted.
403 for v
:= integer
; v
> 0; {
407 buf
[pos
] = byte(v
+ '0')
410 for i
:= pos
; i
< len(buf
); i
++ {
415 d
.dp
= integerDigits
+ exp10
419 if rest
!= 0 || pow10
!= 1 {
420 panic("strconv: internal error, rest != 0 but needed > 0")
422 // Emit digits for the fractional part. Each time, 10*fraction
423 // fits in a uint64 without overflow.
426 ε
*= 10 // the uncertainty scales as we multiply by ten.
428 // the error is so large it could modify which digit to write, abort.
431 digit
:= fraction
>> shift
432 d
.d
[nd
] = byte(digit
+ '0')
433 fraction
-= digit
<< shift
440 // We have written a truncation of f (a numerator / 10^d.dp). The remaining part
441 // can be interpreted as a small number (< 1) to be added to the last digit of the
444 // If rest > 0, the amount is:
445 // (rest<<shift | fraction) / (pow10 << shift)
446 // fraction being known with a ±ε uncertainty.
447 // The fact that n > 0 guarantees that pow10 << shift does not overflow a uint64.
449 // If rest = 0, pow10 == 1 and the amount is
450 // fraction / (1 << shift)
451 // fraction being known with a ±ε uncertainty.
453 // We pass this information to the rounding routine for adjustment.
455 ok
:= adjustLastDigitFixed(d
, uint64(rest
)<<shift|fraction
, pow10
, shift
, ε
)
459 // Trim trailing zeros.
460 for i
:= d
.nd
- 1; i
>= 0; i
-- {
469 // adjustLastDigitFixed assumes d contains the representation of the integral part
470 // of some number, whose fractional part is num / (den << shift). The numerator
471 // num is only known up to an uncertainty of size ε, assumed to be less than
474 // It will increase the last digit by one to account for correct rounding, typically
475 // when the fractional part is greater than 1/2, and will return false if ε is such
476 // that no correct answer can be given.
477 func adjustLastDigitFixed(d
*decimalSlice
, num
, den
uint64, shift
uint, ε
uint64) bool {
478 if num
> den
<<shift
{
479 panic("strconv: num > den<<shift in adjustLastDigitFixed")
481 if 2*ε
> den
<<shift
{
482 panic("strconv: ε > (den<<shift)/2")
484 if 2*(num
+ε
) < den
<<shift
{
487 if 2*(num
-ε
) > den
<<shift
{
509 // ShortestDecimal stores in d the shortest decimal representation of f
510 // which belongs to the open interval (lower, upper), where f is supposed
511 // to lie. It returns false whenever the result is unsure. The implementation
512 // uses the Grisu3 algorithm.
513 func (f
*extFloat
) ShortestDecimal(d
*decimalSlice
, lower
, upper
*extFloat
) bool {
520 if f
.exp
== 0 && *lower
== *f
&& *lower
== *upper
{
524 for v
:= f
.mant
; v
> 0; {
527 buf
[n
] = byte(v
+ '0')
531 nd
:= len(buf
) - n
- 1
532 for i
:= 0; i
< nd
; i
++ {
536 for d
.nd
> 0 && d
.d
[d
.nd
-1] == '0' {
546 // Uniformize exponents.
547 if f
.exp
> upper
.exp
{
548 f
.mant
<<= uint(f
.exp
- upper
.exp
)
551 if lower
.exp
> upper
.exp
{
552 lower
.mant
<<= uint(lower
.exp
- upper
.exp
)
553 lower
.exp
= upper
.exp
556 exp10
:= frexp10Many(lower
, f
, upper
)
557 // Take a safety margin due to rounding in frexp10Many, but we lose precision.
561 // The shortest representation of f is either rounded up or down, but
562 // in any case, it is a truncation of upper.
563 shift
:= uint(-upper
.exp
)
564 integer
:= uint32(upper
.mant
>> shift
)
565 fraction
:= upper
.mant
- (uint64(integer
) << shift
)
567 // How far we can go down from upper until the result is wrong.
568 allowance
:= upper
.mant
- lower
.mant
569 // How far we should go to get a very precise result.
570 targetDiff
:= upper
.mant
- f
.mant
572 // Count integral digits: there are at most 10.
573 var integerDigits
int
574 for i
, pow
:= 0, uint64(1); i
< 20; i
++ {
575 if pow
> uint64(integer
) {
581 for i
:= 0; i
< integerDigits
; i
++ {
582 pow
:= uint64pow10
[integerDigits
-i
-1]
583 digit
:= integer
/ uint32(pow
)
584 d
.d
[i
] = byte(digit
+ '0')
585 integer
-= digit
* uint32(pow
)
586 // evaluate whether we should stop.
587 if currentDiff
:= uint64(integer
)<<shift
+ fraction
; currentDiff
< allowance
{
589 d
.dp
= integerDigits
+ exp10
591 // Sometimes allowance is so large the last digit might need to be
592 // decremented to get closer to f.
593 return adjustLastDigit(d
, currentDiff
, targetDiff
, allowance
, pow
<<shift
, 2)
600 // Compute digits of the fractional part. At each step fraction does not
601 // overflow. The choice of minExp implies that fraction is less than 2^60.
603 multiplier
:= uint64(1)
607 digit
= int(fraction
>> shift
)
608 d
.d
[d
.nd
] = byte(digit
+ '0')
610 fraction
-= uint64(digit
) << shift
611 if fraction
< allowance
*multiplier
{
612 // We are in the admissible range. Note that if allowance is about to
613 // overflow, that is, allowance > 2^64/10, the condition is automatically
614 // true due to the limited range of fraction.
615 return adjustLastDigit(d
,
616 fraction
, targetDiff
*multiplier
, allowance
*multiplier
,
617 1<<shift
, multiplier
*2)
622 // adjustLastDigit modifies d = x-currentDiff*ε, to get closest to
623 // d = x-targetDiff*ε, without becoming smaller than x-maxDiff*ε.
624 // It assumes that a decimal digit is worth ulpDecimal*ε, and that
625 // all data is known with an error estimate of ulpBinary*ε.
626 func adjustLastDigit(d
*decimalSlice
, currentDiff
, targetDiff
, maxDiff
, ulpDecimal
, ulpBinary
uint64) bool {
627 if ulpDecimal
< 2*ulpBinary
{
628 // Approximation is too wide.
631 for currentDiff
+ulpDecimal
/2+ulpBinary
< targetDiff
{
633 currentDiff
+= ulpDecimal
635 if currentDiff
+ulpDecimal
<= targetDiff
+ulpDecimal
/2+ulpBinary
{
636 // we have two choices, and don't know what to do.
639 if currentDiff
< ulpBinary || currentDiff
> maxDiff
-ulpBinary
{
643 if d
.nd
== 1 && d
.d
[0] == '0' {
644 // the number has actually reached zero.