blob: d31f44142c2b99636c21694ff02cb1a957663fc6 [file] [log] [blame]
Russ Cox079c00a2008-11-17 12:34:03 -08001// Copyright 2009 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.
4
5// Binary to decimal floating point conversion.
6// Algorithm:
7// 1) store mantissa in multiprecision decimal
8// 2) shift decimal by exponent
9// 3) read digits out & format
10
11package strconv
12
Russ Cox3b864e42009-08-12 13:18:37 -070013import "math"
Russ Cox079c00a2008-11-17 12:34:03 -080014
15// TODO: move elsewhere?
Russ Cox8a7cbad2009-01-15 17:22:17 -080016type floatInfo struct {
Russ Cox094f1d52009-10-08 15:14:54 -070017 mantbits uint;
18 expbits uint;
19 bias int;
Russ Cox079c00a2008-11-17 12:34:03 -080020}
Russ Cox094f1d52009-10-08 15:14:54 -070021
22var float32info = floatInfo{23, 8, -127}
23var float64info = floatInfo{52, 11, -1023}
Russ Cox079c00a2008-11-17 12:34:03 -080024
Russ Coxb54133d2009-01-15 16:16:42 -080025func floatsize() int {
Russ Cox079c00a2008-11-17 12:34:03 -080026 // Figure out whether float is float32 or float64.
27 // 1e-35 is representable in both, but 1e-70
28 // is too small for a float32.
29 var f float = 1e-35;
30 if f*f == 0 {
Robert Griesemer40621d52009-11-09 12:07:39 -080031 return 32
Russ Cox079c00a2008-11-17 12:34:03 -080032 }
33 return 64;
34}
Russ Cox5bf0fbe2009-03-05 15:29:04 -080035
36// Floatsize gives the size of the float type, either 32 or 64.
Russ Cox839a6842009-01-20 14:40:40 -080037var FloatSize = floatsize()
Russ Cox079c00a2008-11-17 12:34:03 -080038
Russ Cox5bf0fbe2009-03-05 15:29:04 -080039// Ftoa32 converts the 32-bit floating-point number f to a string,
40// according to the format fmt and precision prec.
41//
42// The format fmt is one of
43// 'b' (-ddddp±ddd, a binary exponent),
44// 'e' (-d.dddde±dd, a decimal exponent),
45// 'f' (-ddd.dddd, no exponent), or
46// 'g' ('e' for large exponents, 'f' otherwise).
47//
48// The precision prec controls the number of digits
49// (excluding the exponent) printed by the 'e', 'f', and 'g' formats.
50// For 'e' and 'f' it is the number of digits after the decimal point.
51// For 'g' it is the total number of digits.
52// The special precision -1 uses the smallest number of digits
53// necessary such that Atof32 will return f exactly.
54//
55// Ftoa32(f) is not the same as Ftoa64(float32(f)),
56// because correct rounding and the number of digits
57// needed to identify f depend on the precision of the representation.
Russ Cox839a6842009-01-20 14:40:40 -080058func Ftoa32(f float32, fmt byte, prec int) string {
Robert Griesemer40621d52009-11-09 12:07:39 -080059 return genericFtoa(uint64(math.Float32bits(f)), fmt, prec, &float32info)
Russ Cox079c00a2008-11-17 12:34:03 -080060}
61
Russ Cox5bf0fbe2009-03-05 15:29:04 -080062// Ftoa64 is like Ftoa32 but converts a 64-bit floating-point number.
Russ Cox839a6842009-01-20 14:40:40 -080063func Ftoa64(f float64, fmt byte, prec int) string {
Robert Griesemer40621d52009-11-09 12:07:39 -080064 return genericFtoa(math.Float64bits(f), fmt, prec, &float64info)
Russ Cox079c00a2008-11-17 12:34:03 -080065}
66
Russ Cox5bf0fbe2009-03-05 15:29:04 -080067// Ftoa behaves as Ftoa32 or Ftoa64, depending on the size of the float type.
Russ Cox839a6842009-01-20 14:40:40 -080068func Ftoa(f float, fmt byte, prec int) string {
Russ Coxb54133d2009-01-15 16:16:42 -080069 if FloatSize == 32 {
Robert Griesemer40621d52009-11-09 12:07:39 -080070 return Ftoa32(float32(f), fmt, prec)
Russ Cox079c00a2008-11-17 12:34:03 -080071 }
Russ Cox8a7cbad2009-01-15 17:22:17 -080072 return Ftoa64(float64(f), fmt, prec);
Russ Cox079c00a2008-11-17 12:34:03 -080073}
74
Russ Cox8a7cbad2009-01-15 17:22:17 -080075func genericFtoa(bits uint64, fmt byte, prec int, flt *floatInfo) string {
Russ Cox094f1d52009-10-08 15:14:54 -070076 neg := bits >> flt.expbits >> flt.mantbits != 0;
77 exp := int(bits >> flt.mantbits)&(1 << flt.expbits - 1);
78 mant := bits&(uint64(1) << flt.mantbits - 1);
Russ Cox079c00a2008-11-17 12:34:03 -080079
80 switch exp {
Russ Cox094f1d52009-10-08 15:14:54 -070081 case 1 << flt.expbits - 1:
Russ Cox079c00a2008-11-17 12:34:03 -080082 // Inf, NaN
83 if mant != 0 {
Robert Griesemer40621d52009-11-09 12:07:39 -080084 return "NaN"
Russ Cox079c00a2008-11-17 12:34:03 -080085 }
86 if neg {
Robert Griesemer40621d52009-11-09 12:07:39 -080087 return "-Inf"
Russ Cox079c00a2008-11-17 12:34:03 -080088 }
89 return "+Inf";
90
91 case 0:
92 // denormalized
Robert Griesemer40621d52009-11-09 12:07:39 -080093 exp++
Russ Cox079c00a2008-11-17 12:34:03 -080094
95 default:
96 // add implicit top bit
Robert Griesemer40621d52009-11-09 12:07:39 -080097 mant |= uint64(1) << flt.mantbits
Russ Cox079c00a2008-11-17 12:34:03 -080098 }
99 exp += flt.bias;
100
101 // Pick off easy binary format.
102 if fmt == 'b' {
Robert Griesemer40621d52009-11-09 12:07:39 -0800103 return fmtB(neg, mant, exp, flt)
Russ Cox079c00a2008-11-17 12:34:03 -0800104 }
105
106 // Create exact decimal representation.
107 // The shift is exp - flt.mantbits because mant is a 1-bit integer
108 // followed by a flt.mantbits fraction, and we are treating it as
109 // a 1+flt.mantbits-bit integer.
Russ Cox094f1d52009-10-08 15:14:54 -0700110 d := newDecimal(mant).Shift(exp-int(flt.mantbits));
Russ Cox079c00a2008-11-17 12:34:03 -0800111
112 // Round appropriately.
113 // Negative precision means "only as much as needed to be exact."
Russ Cox0e198da2008-11-23 17:27:44 -0800114 shortest := false;
Russ Cox079c00a2008-11-17 12:34:03 -0800115 if prec < 0 {
Russ Cox0e198da2008-11-23 17:27:44 -0800116 shortest = true;
Russ Cox8a7cbad2009-01-15 17:22:17 -0800117 roundShortest(d, mant, exp, flt);
Russ Cox079c00a2008-11-17 12:34:03 -0800118 switch fmt {
Russ Coxa843b452009-08-31 16:38:30 -0700119 case 'e', 'E':
Robert Griesemer40621d52009-11-09 12:07:39 -0800120 prec = d.nd - 1
Russ Cox079c00a2008-11-17 12:34:03 -0800121 case 'f':
Robert Griesemer40621d52009-11-09 12:07:39 -0800122 prec = max(d.nd - d.dp, 0)
Russ Coxa843b452009-08-31 16:38:30 -0700123 case 'g', 'G':
Robert Griesemer40621d52009-11-09 12:07:39 -0800124 prec = d.nd
Russ Cox079c00a2008-11-17 12:34:03 -0800125 }
126 } else {
127 switch fmt {
Russ Coxa843b452009-08-31 16:38:30 -0700128 case 'e', 'E':
Robert Griesemer40621d52009-11-09 12:07:39 -0800129 d.Round(prec+1)
Russ Cox079c00a2008-11-17 12:34:03 -0800130 case 'f':
Robert Griesemer40621d52009-11-09 12:07:39 -0800131 d.Round(d.dp + prec)
Russ Coxa843b452009-08-31 16:38:30 -0700132 case 'g', 'G':
Russ Cox079c00a2008-11-17 12:34:03 -0800133 if prec == 0 {
Robert Griesemer40621d52009-11-09 12:07:39 -0800134 prec = 1
Russ Cox079c00a2008-11-17 12:34:03 -0800135 }
136 d.Round(prec);
137 }
138 }
139
140 switch fmt {
Russ Coxa843b452009-08-31 16:38:30 -0700141 case 'e', 'E':
Robert Griesemer40621d52009-11-09 12:07:39 -0800142 return fmtE(neg, d, prec, fmt)
Russ Cox079c00a2008-11-17 12:34:03 -0800143 case 'f':
Robert Griesemer40621d52009-11-09 12:07:39 -0800144 return fmtF(neg, d, prec)
Russ Coxa843b452009-08-31 16:38:30 -0700145 case 'g', 'G':
Russ Cox079c00a2008-11-17 12:34:03 -0800146 // trailing zeros are removed.
147 if prec > d.nd {
Robert Griesemer40621d52009-11-09 12:07:39 -0800148 prec = d.nd
Russ Cox079c00a2008-11-17 12:34:03 -0800149 }
150 // %e is used if the exponent from the conversion
151 // is less than -4 or greater than or equal to the precision.
Russ Cox0e198da2008-11-23 17:27:44 -0800152 // if precision was the shortest possible, use precision 6 for this decision.
153 eprec := prec;
154 if shortest {
Robert Griesemer40621d52009-11-09 12:07:39 -0800155 eprec = 6
Russ Cox0e198da2008-11-23 17:27:44 -0800156 }
Russ Cox079c00a2008-11-17 12:34:03 -0800157 exp := d.dp - 1;
Russ Cox0e198da2008-11-23 17:27:44 -0800158 if exp < -4 || exp >= eprec {
Robert Griesemer40621d52009-11-09 12:07:39 -0800159 return fmtE(neg, d, prec-1, fmt+'e'-'g')
Russ Cox079c00a2008-11-17 12:34:03 -0800160 }
Russ Cox8a7cbad2009-01-15 17:22:17 -0800161 return fmtF(neg, d, max(prec - d.dp, 0));
Russ Cox079c00a2008-11-17 12:34:03 -0800162 }
163
Russ Cox094f1d52009-10-08 15:14:54 -0700164 return "%"+string(fmt);
Russ Cox079c00a2008-11-17 12:34:03 -0800165}
166
167// Round d (= mant * 2^exp) to the shortest number of digits
168// that will let the original floating point value be precisely
169// reconstructed. Size is original floating point size (64 or 32).
Russ Cox8a7cbad2009-01-15 17:22:17 -0800170func roundShortest(d *decimal, mant uint64, exp int, flt *floatInfo) {
Russ Coxcf9b7f72008-11-19 12:50:34 -0800171 // If mantissa is zero, the number is zero; stop now.
172 if mant == 0 {
173 d.nd = 0;
174 return;
175 }
176
Russ Cox7732d802009-11-01 09:25:55 -0800177 // TODO(rsc): Unless exp == minexp, if the number of digits in d
178 // is less than 17, it seems likely that it would be
Russ Cox079c00a2008-11-17 12:34:03 -0800179 // the shortest possible number already. So maybe we can
180 // bail out without doing the extra multiprecision math here.
181
182 // Compute upper and lower such that any decimal number
183 // between upper and lower (possibly inclusive)
184 // will round to the original floating point number.
185
186 // d = mant << (exp - mantbits)
187 // Next highest floating point number is mant+1 << exp-mantbits.
188 // Our upper bound is halfway inbetween, mant*2+1 << exp-mantbits-1.
Russ Cox094f1d52009-10-08 15:14:54 -0700189 upper := newDecimal(mant*2 + 1).Shift(exp-int(flt.mantbits)-1);
Russ Cox079c00a2008-11-17 12:34:03 -0800190
191 // d = mant << (exp - mantbits)
192 // Next lowest floating point number is mant-1 << exp-mantbits,
193 // unless mant-1 drops the significant bit and exp is not the minimum exp,
194 // in which case the next lowest is mant*2-1 << exp-mantbits-1.
195 // Either way, call it mantlo << explo-mantbits.
196 // Our lower bound is halfway inbetween, mantlo*2+1 << explo-mantbits-1.
197 minexp := flt.bias + 1; // minimum possible exponent
198 var mantlo uint64;
199 var explo int;
Russ Cox094f1d52009-10-08 15:14:54 -0700200 if mant > 1 << flt.mantbits || exp == minexp {
201 mantlo = mant-1;
Russ Cox079c00a2008-11-17 12:34:03 -0800202 explo = exp;
203 } else {
Russ Cox094f1d52009-10-08 15:14:54 -0700204 mantlo = mant*2 - 1;
Russ Cox079c00a2008-11-17 12:34:03 -0800205 explo = exp-1;
206 }
Russ Cox094f1d52009-10-08 15:14:54 -0700207 lower := newDecimal(mantlo*2 + 1).Shift(explo-int(flt.mantbits)-1);
Russ Cox079c00a2008-11-17 12:34:03 -0800208
209 // The upper and lower bounds are possible outputs only if
210 // the original mantissa is even, so that IEEE round-to-even
211 // would round to the original mantissa and not the neighbors.
212 inclusive := mant%2 == 0;
213
214 // Now we can figure out the minimum number of digits required.
215 // Walk along until d has distinguished itself from upper and lower.
216 for i := 0; i < d.nd; i++ {
217 var l, m, u byte; // lower, middle, upper digits
218 if i < lower.nd {
Robert Griesemer40621d52009-11-09 12:07:39 -0800219 l = lower.d[i]
Russ Cox079c00a2008-11-17 12:34:03 -0800220 } else {
Robert Griesemer40621d52009-11-09 12:07:39 -0800221 l = '0'
Russ Cox079c00a2008-11-17 12:34:03 -0800222 }
223 m = d.d[i];
224 if i < upper.nd {
Robert Griesemer40621d52009-11-09 12:07:39 -0800225 u = upper.d[i]
Russ Cox079c00a2008-11-17 12:34:03 -0800226 } else {
Robert Griesemer40621d52009-11-09 12:07:39 -0800227 u = '0'
Russ Cox079c00a2008-11-17 12:34:03 -0800228 }
229
230 // Okay to round down (truncate) if lower has a different digit
231 // or if lower is inclusive and is exactly the result of rounding down.
232 okdown := l != m || (inclusive && l == m && i+1 == lower.nd);
233
234 // Okay to round up if upper has a different digit and
235 // either upper is inclusive or upper is bigger than the result of rounding up.
236 okup := m != u && (inclusive || i+1 < upper.nd);
237
238 // If it's okay to do either, then round to the nearest one.
239 // If it's okay to do only one, do it.
240 switch {
241 case okdown && okup:
242 d.Round(i+1);
243 return;
244 case okdown:
245 d.RoundDown(i+1);
246 return;
247 case okup:
248 d.RoundUp(i+1);
249 return;
250 }
251 }
252}
253
254// %e: -d.ddddde±dd
Russ Coxa843b452009-08-31 16:38:30 -0700255func fmtE(neg bool, d *decimal, prec int, fmt byte) string {
Russ Cox094f1d52009-10-08 15:14:54 -0700256 buf := make([]byte, 3 + max(prec, 0) + 30); // "-0." + prec digits + exp
257 w := 0; // write index
Russ Cox079c00a2008-11-17 12:34:03 -0800258
259 // sign
260 if neg {
261 buf[w] = '-';
262 w++;
263 }
264
265 // first digit
266 if d.nd == 0 {
Robert Griesemer40621d52009-11-09 12:07:39 -0800267 buf[w] = '0'
Russ Cox079c00a2008-11-17 12:34:03 -0800268 } else {
Robert Griesemer40621d52009-11-09 12:07:39 -0800269 buf[w] = d.d[0]
Russ Cox079c00a2008-11-17 12:34:03 -0800270 }
271 w++;
272
273 // .moredigits
274 if prec > 0 {
275 buf[w] = '.';
276 w++;
277 for i := 0; i < prec; i++ {
278 if 1+i < d.nd {
Robert Griesemer40621d52009-11-09 12:07:39 -0800279 buf[w] = d.d[1+i]
Russ Cox079c00a2008-11-17 12:34:03 -0800280 } else {
Robert Griesemer40621d52009-11-09 12:07:39 -0800281 buf[w] = '0'
Russ Cox079c00a2008-11-17 12:34:03 -0800282 }
283 w++;
284 }
285 }
286
287 // e±
Russ Coxa843b452009-08-31 16:38:30 -0700288 buf[w] = fmt;
Russ Cox079c00a2008-11-17 12:34:03 -0800289 w++;
290 exp := d.dp - 1;
291 if d.nd == 0 { // special case: 0 has exponent 0
Robert Griesemer40621d52009-11-09 12:07:39 -0800292 exp = 0
Russ Cox079c00a2008-11-17 12:34:03 -0800293 }
294 if exp < 0 {
295 buf[w] = '-';
296 exp = -exp;
297 } else {
Robert Griesemer40621d52009-11-09 12:07:39 -0800298 buf[w] = '+'
Russ Cox079c00a2008-11-17 12:34:03 -0800299 }
300 w++;
301
302 // dddd
303 // count digits
304 n := 0;
305 for e := exp; e > 0; e /= 10 {
Robert Griesemer40621d52009-11-09 12:07:39 -0800306 n++
Russ Cox079c00a2008-11-17 12:34:03 -0800307 }
308 // leading zeros
309 for i := n; i < 2; i++ {
310 buf[w] = '0';
311 w++;
312 }
313 // digits
314 w += n;
315 n = 0;
316 for e := exp; e > 0; e /= 10 {
317 n++;
318 buf[w-n] = byte(e%10 + '0');
319 }
320
321 return string(buf[0:w]);
322}
323
324// %f: -ddddddd.ddddd
Russ Cox8a7cbad2009-01-15 17:22:17 -0800325func fmtF(neg bool, d *decimal, prec int) string {
Russ Cox094f1d52009-10-08 15:14:54 -0700326 buf := make([]byte, 1 + max(d.dp, 1) + 1 + max(prec, 0));
Russ Cox079c00a2008-11-17 12:34:03 -0800327 w := 0;
328
329 // sign
330 if neg {
331 buf[w] = '-';
332 w++;
333 }
334
335 // integer, padded with zeros as needed.
336 if d.dp > 0 {
337 var i int;
338 for i = 0; i < d.dp && i < d.nd; i++ {
339 buf[w] = d.d[i];
340 w++;
341 }
342 for ; i < d.dp; i++ {
343 buf[w] = '0';
344 w++;
345 }
346 } else {
347 buf[w] = '0';
348 w++;
349 }
350
351 // fraction
352 if prec > 0 {
353 buf[w] = '.';
354 w++;
355 for i := 0; i < prec; i++ {
Russ Cox094f1d52009-10-08 15:14:54 -0700356 if d.dp + i < 0 || d.dp + i >= d.nd {
Robert Griesemer40621d52009-11-09 12:07:39 -0800357 buf[w] = '0'
Russ Cox079c00a2008-11-17 12:34:03 -0800358 } else {
Robert Griesemer40621d52009-11-09 12:07:39 -0800359 buf[w] = d.d[d.dp + i]
Russ Cox079c00a2008-11-17 12:34:03 -0800360 }
361 w++;
362 }
363 }
364
365 return string(buf[0:w]);
366}
367
368// %b: -ddddddddp+ddd
Russ Cox8a7cbad2009-01-15 17:22:17 -0800369func fmtB(neg bool, mant uint64, exp int, flt *floatInfo) string {
Russ Cox079c00a2008-11-17 12:34:03 -0800370 var buf [50]byte;
371 w := len(buf);
372 exp -= int(flt.mantbits);
373 esign := byte('+');
374 if exp < 0 {
375 esign = '-';
376 exp = -exp;
377 }
378 n := 0;
379 for exp > 0 || n < 1 {
380 n++;
381 w--;
382 buf[w] = byte(exp%10 + '0');
Russ Cox094f1d52009-10-08 15:14:54 -0700383 exp /= 10;
Russ Cox079c00a2008-11-17 12:34:03 -0800384 }
385 w--;
386 buf[w] = esign;
387 w--;
388 buf[w] = 'p';
389 n = 0;
390 for mant > 0 || n < 1 {
391 n++;
392 w--;
393 buf[w] = byte(mant%10 + '0');
394 mant /= 10;
395 }
396 if neg {
397 w--;
398 buf[w] = '-';
399 }
Russ Coxd47d8882008-12-18 22:37:22 -0800400 return string(buf[w:len(buf)]);
Russ Cox079c00a2008-11-17 12:34:03 -0800401}
402
Russ Cox8a7cbad2009-01-15 17:22:17 -0800403func max(a, b int) int {
Russ Cox079c00a2008-11-17 12:34:03 -0800404 if a > b {
Robert Griesemer40621d52009-11-09 12:07:39 -0800405 return a
Russ Cox079c00a2008-11-17 12:34:03 -0800406 }
407 return b;
408}