float.go 7.3 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278
  1. // MIT License
  2. //
  3. // Copyright (c) 2017 José Santos <henrique_1609@me.com>
  4. //
  5. // Permission is hereby granted, free of charge, to any person obtaining a copy
  6. // of this software and associated documentation files (the "Software"), to deal
  7. // in the Software without restriction, including without limitation the rights
  8. // to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
  9. // copies of the Software, and to permit persons to whom the Software is
  10. // furnished to do so, subject to the following conditions:
  11. //
  12. // The above copyright notice and this permission notice shall be included in all
  13. // copies or substantial portions of the Software.
  14. //
  15. // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  16. // IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  17. // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
  18. // AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  19. // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
  20. // OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
  21. // SOFTWARE.
  22. package fastprinter
  23. import (
  24. "io"
  25. "math"
  26. )
  27. type floatInfo struct {
  28. mantbits uint
  29. expbits uint
  30. bias int
  31. }
  32. var (
  33. float64info = floatInfo{52, 11, -1023}
  34. floatNaN = []byte("Nan")
  35. floatNinf = []byte("-Inf")
  36. floatPinf = []byte("+Inf")
  37. pool_floatBuffer = newByteSliceBufferPool(800)
  38. )
  39. func PrintFloat(w io.Writer, f float64) (int, error) {
  40. return PrintFloatPrecision(w, f, -1)
  41. }
  42. func PrintFloatPrecision(dst io.Writer, val float64, prec int) (int, error) {
  43. var bits uint64
  44. var flt *floatInfo
  45. bits = math.Float64bits(val)
  46. flt = &float64info
  47. neg := bits>>(flt.expbits+flt.mantbits) != 0
  48. exp := int(bits>>flt.mantbits) & (1<<flt.expbits - 1)
  49. mant := bits & (uint64(1)<<flt.mantbits - 1)
  50. switch exp {
  51. case 1<<flt.expbits - 1:
  52. switch {
  53. case mant != 0:
  54. return dst.Write(floatNaN)
  55. case neg:
  56. return dst.Write(floatNinf)
  57. default:
  58. return dst.Write(floatPinf)
  59. }
  60. case 0:
  61. // denormalized
  62. exp++
  63. default:
  64. // add implicit top bit
  65. mant |= uint64(1) << flt.mantbits
  66. }
  67. exp += flt.bias
  68. var digs decimalSlice
  69. ok := false
  70. // Negative precision means "only as much as needed to be exact."
  71. shortest := prec < 0
  72. if shortest {
  73. // Try Grisu3 algorithm.
  74. f := new(extFloat)
  75. lower, upper := f.AssignComputeBounds(mant, exp, neg, flt)
  76. var buf [32]byte
  77. digs.d = buf[:]
  78. ok = f.ShortestDecimal(&digs, &lower, &upper)
  79. if !ok {
  80. return bigFtoa(dst, prec, neg, mant, exp, flt)
  81. }
  82. // Precision for shortest representation mode.
  83. prec = max(digs.nd-digs.dp, 0)
  84. }
  85. if !ok {
  86. return bigFtoa(dst, prec, neg, mant, exp, flt)
  87. }
  88. return fmtF(dst, neg, digs, prec)
  89. }
  90. // bigFtoa uses multiprecision computations to format a float.
  91. func bigFtoa(dst io.Writer, prec int, neg bool, mant uint64, exp int, flt *floatInfo) (int, error) {
  92. d := new(decimal)
  93. d.Assign(mant)
  94. d.Shift(exp - int(flt.mantbits))
  95. var digs decimalSlice
  96. shortest := prec < 0
  97. if shortest {
  98. roundShortest(d, mant, exp, flt)
  99. digs = decimalSlice{d: d.d[:], nd: d.nd, dp: d.dp}
  100. prec = max(digs.nd-digs.dp, 0)
  101. } else {
  102. d.Round(d.dp + prec)
  103. digs = decimalSlice{d: d.d[:], nd: d.nd, dp: d.dp}
  104. }
  105. return fmtF(dst, neg, digs, prec)
  106. }
  107. // roundShortest rounds d (= mant * 2^exp) to the shortest number of digits
  108. // that will let the original floating point value be precisely reconstructed.
  109. func roundShortest(d *decimal, mant uint64, exp int, flt *floatInfo) {
  110. // If mantissa is zero, the number is zero; stop now.
  111. if mant == 0 {
  112. d.nd = 0
  113. return
  114. }
  115. // Compute upper and lower such that any decimal number
  116. // between upper and lower (possibly inclusive)
  117. // will round to the original floating point number.
  118. // We may see at once that the number is already shortest.
  119. //
  120. // Suppose d is not denormal, so that 2^exp <= d < 10^dp.
  121. // The closest shorter number is at least 10^(dp-nd) away.
  122. // The lower/upper bounds computed below are at distance
  123. // at most 2^(exp-mantbits).
  124. //
  125. // So the number is already shortest if 10^(dp-nd) > 2^(exp-mantbits),
  126. // or equivalently log2(10)*(dp-nd) > exp-mantbits.
  127. // It is true if 332/100*(dp-nd) >= exp-mantbits (log2(10) > 3.32).
  128. minexp := flt.bias + 1 // minimum possible exponent
  129. if exp > minexp && 332*(d.dp-d.nd) >= 100*(exp-int(flt.mantbits)) {
  130. // The number is already shortest.
  131. return
  132. }
  133. // d = mant << (exp - mantbits)
  134. // Next highest floating point number is mant+1 << exp-mantbits.
  135. // Our upper bound is halfway between, mant*2+1 << exp-mantbits-1.
  136. upper := new(decimal)
  137. upper.Assign(mant*2 + 1)
  138. upper.Shift(exp - int(flt.mantbits) - 1)
  139. // d = mant << (exp - mantbits)
  140. // Next lowest floating point number is mant-1 << exp-mantbits,
  141. // unless mant-1 drops the significant bit and exp is not the minimum exp,
  142. // in which case the next lowest is mant*2-1 << exp-mantbits-1.
  143. // Either way, call it mantlo << explo-mantbits.
  144. // Our lower bound is halfway between, mantlo*2+1 << explo-mantbits-1.
  145. var mantlo uint64
  146. var explo int
  147. if mant > 1<<flt.mantbits || exp == minexp {
  148. mantlo = mant - 1
  149. explo = exp
  150. } else {
  151. mantlo = mant*2 - 1
  152. explo = exp - 1
  153. }
  154. lower := new(decimal)
  155. lower.Assign(mantlo*2 + 1)
  156. lower.Shift(explo - int(flt.mantbits) - 1)
  157. // The upper and lower bounds are possible outputs only if
  158. // the original mantissa is even, so that IEEE round-to-even
  159. // would round to the original mantissa and not the neighbors.
  160. inclusive := mant%2 == 0
  161. // Now we can figure out the minimum number of digits required.
  162. // Walk along until d has distinguished itself from upper and lower.
  163. for i := 0; i < d.nd; i++ {
  164. l := byte('0') // lower digit
  165. if i < lower.nd {
  166. l = lower.d[i]
  167. }
  168. m := d.d[i] // middle digit
  169. u := byte('0') // upper digit
  170. if i < upper.nd {
  171. u = upper.d[i]
  172. }
  173. // Okay to round down (truncate) if lower has a different digit
  174. // or if lower is inclusive and is exactly the result of rounding
  175. // down (i.e., and we have reached the final digit of lower).
  176. okdown := l != m || inclusive && i+1 == lower.nd
  177. // Okay to round up if upper has a different digit and either upper
  178. // is inclusive or upper is bigger than the result of rounding up.
  179. okup := m != u && (inclusive || m+1 < u || i+1 < upper.nd)
  180. // If it's okay to do either, then round to the nearest one.
  181. // If it's okay to do only one, do it.
  182. switch {
  183. case okdown && okup:
  184. d.Round(i + 1)
  185. return
  186. case okdown:
  187. d.RoundDown(i + 1)
  188. return
  189. case okup:
  190. d.RoundUp(i + 1)
  191. return
  192. }
  193. }
  194. }
  195. type decimalSlice struct {
  196. d []byte
  197. nd, dp int
  198. neg bool
  199. }
  200. // %f: -ddddddd.ddddd
  201. func fmtF(dst io.Writer, neg bool, d decimalSlice, prec int) (n int, err error) {
  202. a := pool_floatBuffer.Get().(*byteSliceBuffer)
  203. i := 0
  204. // sign
  205. if neg {
  206. a.bytes[i] = '-'
  207. i++
  208. }
  209. // integer, padded with zeros as needed.
  210. if d.dp > 0 {
  211. m := min(d.nd, d.dp)
  212. copy(a.bytes[i:], d.d[:m])
  213. i += m
  214. for ; m < d.dp; m++ {
  215. a.bytes[i] = '0'
  216. i++
  217. }
  218. } else {
  219. a.bytes[i] = '0'
  220. i++
  221. }
  222. // fraction
  223. if prec > 0 {
  224. a.bytes[i] = '.'
  225. i++
  226. for j := 0; j < prec; j++ {
  227. ch := byte('0')
  228. if j := d.dp + j; 0 <= j && j < d.nd {
  229. ch = d.d[j]
  230. }
  231. a.bytes[i] = ch
  232. i++
  233. }
  234. }
  235. n, err = dst.Write(a.bytes[0:i])
  236. pool_floatBuffer.Put(a)
  237. return
  238. }
  239. func min(a, b int) int {
  240. if a < b {
  241. return a
  242. }
  243. return b
  244. }
  245. func max(a, b int) int {
  246. if a > b {
  247. return a
  248. }
  249. return b
  250. }