/base/ryu/shortest.jl

https://github.com/JuliaLang/julia · Julia · 494 lines · 453 code · 18 blank · 23 comment · 114 complexity · 8dc653bbafe32439ff79648398cb709f MD5 · raw file

  1. """
  2. b, e10 = reduce_shortest(f[, maxsignif])
  3. Reduce to shortest decimal representation where `abs(f) == b * 10^e10` and `b` is an
  4. integer. If a `maxsignif` argument is provided, then `b < maxsignif`.
  5. """
  6. @inline function reduce_shortest(f::T, maxsignif=nothing) where {T}
  7. U = uinttype(T)
  8. uf = reinterpret(U, f)
  9. m = uf & significand_mask(T)
  10. e = ((uf & exponent_mask(T)) >> significand_bits(T)) % Int
  11. ## Step 1
  12. # mf * 2^ef == f
  13. mf = (one(U) << significand_bits(T)) | m
  14. ef = e - exponent_bias(T) - significand_bits(T)
  15. f_isinteger = mf & ((one(U) << -ef) - one(U)) == 0
  16. if ef > 0 || ef < -Base.significand_bits(T) || !f_isinteger
  17. # fixup subnormals
  18. if e == 0
  19. ef = 1 - exponent_bias(T) - significand_bits(T)
  20. mf = m
  21. end
  22. ## Step 2
  23. # u * 2^e2 == (f + prevfloat(f))/2
  24. # v * 2^e2 == f
  25. # w * 2^e2 == (f + nextfloat(f))/2
  26. e2 = ef - 2
  27. mf_iseven = iseven(mf) # trailing bit of significand is zero
  28. v = U(4) * mf
  29. w = v + U(2)
  30. u_shift_half = m == 0 && e > 1 # if first element of binade, other than first normal one
  31. u = v - U(2) + u_shift_half
  32. ## Step 3
  33. # a == floor(u * 2^e2 / 10^e10), exact if a_allzero
  34. # b == floor(v * 2^e2 / 10^e10), exact if b_allzero
  35. # c == floor(w * 2^e2 / 10^e10)
  36. a_allzero = false
  37. b_allzero = false
  38. b_lastdigit = 0x00
  39. if e2 >= 0
  40. q = log10pow2(e2) - (T == Float64 ? (e2 > 3) : 0)
  41. e10 = q
  42. k = pow5_inv_bitcount(T) + pow5bits(q) - 1
  43. i = -e2 + q + k
  44. a, b, c = mulshiftinvsplit(T, u, v, w, q, i)
  45. if T == Float32 || T == Float16
  46. if q != 0 && div(c - 1, 10) <= div(a, 10)
  47. l = pow5_inv_bitcount(T) + pow5bits(q - 1) - 1
  48. mul = pow5invsplit_lookup(T, q-1)
  49. b_lastdigit = (mulshift(v, mul, -e2 + q - 1 + l) % 10) % UInt8
  50. end
  51. end
  52. if q <= qinvbound(T)
  53. if ((v % UInt32) - 5 * div(v, 5)) == 0
  54. b_allzero = pow5(v, q)
  55. elseif mf_iseven
  56. a_allzero = pow5(u, q)
  57. else
  58. c -= pow5(w, q)
  59. end
  60. end
  61. else
  62. q = log10pow5(-e2) - (T == Float64 ? (-e2 > 1) : 0)
  63. e10 = q + e2
  64. i = -e2 - q
  65. k = pow5bits(i) - pow5_bitcount(T)
  66. j = q - k
  67. a, b, c = mulshiftsplit(T, u, v, w, i, j)
  68. if T == Float32 || T == Float16
  69. if q != 0 && div(c - 1, 10) <= div(a, 10)
  70. j = q - 1 - (pow5bits(i + 1) - pow5_bitcount(T))
  71. mul = pow5split_lookup(T, i+1)
  72. b_lastdigit = (mulshift(v, mul, j) % 10) % UInt8
  73. end
  74. end
  75. if q <= 1
  76. b_allzero = true
  77. if mf_iseven
  78. a_allzero = !u_shift_half
  79. else
  80. c -= 1
  81. end
  82. elseif q < qbound(T)
  83. b_allzero = pow2(v, q - (T != Float64))
  84. end
  85. end
  86. ## Step 4: reduction
  87. if a_allzero || b_allzero
  88. # a) slow loop
  89. while true
  90. c_div10 = div(c, 10)
  91. a_div10 = div(a, 10)
  92. if c_div10 <= a_div10
  93. break
  94. end
  95. a_mod10 = (a % UInt32) - UInt32(10) * (a_div10 % UInt32)
  96. b_div10 = div(b, 10)
  97. b_mod10 = (b % UInt32) - UInt32(10) * (b_div10 % UInt32)
  98. a_allzero &= a_mod10 == 0
  99. b_allzero &= b_lastdigit == 0
  100. b_lastdigit = b_mod10 % UInt8
  101. b = b_div10
  102. c = c_div10
  103. a = a_div10
  104. e10 += 1
  105. end
  106. if a_allzero
  107. while true
  108. a_div10 = div(a, 10)
  109. a_mod10 = (a % UInt32) - UInt32(10) * (a_div10 % UInt32)
  110. if a_mod10 != 0 && (maxsignif === nothing || b < maxsignif)
  111. break
  112. end
  113. c_div10 = div(c, 10)
  114. b_div10 = div(b, 10)
  115. b_mod10 = (b % UInt32) - UInt32(10) * (b_div10 % UInt32)
  116. b_allzero &= b_lastdigit == 0
  117. b_lastdigit = b_mod10 % UInt8
  118. b = b_div10
  119. c = c_div10
  120. a = a_div10
  121. e10 += 1
  122. end
  123. end
  124. if b_allzero && b_lastdigit == 5 && iseven(b)
  125. b_lastdigit = UInt8(4)
  126. end
  127. roundup = (b == a && (!mf_iseven || !a_allzero)) || b_lastdigit >= 5
  128. else
  129. # b) specialized for common case (99% Float64, 96% Float32)
  130. roundup = b_lastdigit >= 5
  131. c_div100 = div(c, 100)
  132. a_div100 = div(a, 100)
  133. if c_div100 > a_div100
  134. b_div100 = div(b, 100)
  135. b_mod100 = (b % UInt32) - UInt32(100) * (b_div100 % UInt32)
  136. roundup = b_mod100 >= 50
  137. b = b_div100
  138. c = c_div100
  139. a = a_div100
  140. e10 += 2
  141. end
  142. while true
  143. c_div10 = div(c, 10)
  144. a_div10 = div(a, 10)
  145. if c_div10 <= a_div10
  146. break
  147. end
  148. b_div10 = div(b, 10)
  149. b_mod10 = (b % UInt32) - UInt32(10) * (b_div10 % UInt32)
  150. roundup = b_mod10 >= 5
  151. b = b_div10
  152. c = c_div10
  153. a = a_div10
  154. e10 += 1
  155. end
  156. roundup = (b == a || roundup)
  157. end
  158. if maxsignif !== nothing && b > maxsignif
  159. # reduce to max significant digits
  160. while true
  161. b_div10 = div(b, 10)
  162. b_mod10 = (b % UInt32) - UInt32(10) * (b_div10 % UInt32)
  163. if b <= maxsignif
  164. break
  165. end
  166. b = b_div10
  167. roundup = (b_allzero && iseven(b)) ? b_mod10 > 5 : b_mod10 >= 5
  168. b_allzero &= b_mod10 == 0
  169. e10 += 1
  170. end
  171. b = b + roundup
  172. # remove trailing zeros
  173. while true
  174. b_div10 = div(b, 10)
  175. b_mod10 = (b % UInt32) - UInt32(10) * (b_div10 % UInt32)
  176. if b_mod10 != 0
  177. break
  178. end
  179. b = b_div10
  180. e10 += 1
  181. end
  182. else
  183. b = b + roundup
  184. end
  185. else
  186. # c) specialized f an integer < 2^53
  187. b = mf >> -ef
  188. e10 = 0
  189. if maxsignif !== nothing && b > maxsignif
  190. b_allzero = true
  191. # reduce to max significant digits
  192. while true
  193. b_div10 = div(b, 10)
  194. b_mod10 = (b % UInt32) - UInt32(10) * (b_div10 % UInt32)
  195. if b <= maxsignif
  196. break
  197. end
  198. b = b_div10
  199. roundup = (b_allzero && iseven(b)) ? b_mod10 > 5 : b_mod10 >= 5
  200. b_allzero &= b_mod10 == 0
  201. e10 += 1
  202. end
  203. b = b + roundup
  204. end
  205. while true
  206. b_div10 = div(b, 10)
  207. b_mod10 = (b % UInt32) - UInt32(10) * (b_div10 % UInt32)
  208. if b_mod10 != 0
  209. break
  210. end
  211. b = b_div10
  212. e10 += 1
  213. end
  214. end
  215. return b, e10
  216. end
  217. function writeshortest(buf::Vector{UInt8}, pos, x::T,
  218. plus=false, space=false, hash=true,
  219. precision=-1, expchar=UInt8('e'), padexp=false, decchar=UInt8('.'),
  220. typed=false, compact=false) where {T}
  221. @assert 0 < pos <= length(buf)
  222. # special cases
  223. if x == 0
  224. if typed && x isa Float16
  225. buf[pos] = UInt8('F')
  226. buf[pos + 1] = UInt8('l')
  227. buf[pos + 2] = UInt8('o')
  228. buf[pos + 3] = UInt8('a')
  229. buf[pos + 4] = UInt8('t')
  230. buf[pos + 5] = UInt8('1')
  231. buf[pos + 6] = UInt8('6')
  232. buf[pos + 7] = UInt8('(')
  233. pos += 8
  234. end
  235. pos = append_sign(x, plus, space, buf, pos)
  236. buf[pos] = UInt8('0')
  237. pos += 1
  238. if hash
  239. buf[pos] = decchar
  240. pos += 1
  241. end
  242. if precision == -1
  243. buf[pos] = UInt8('0')
  244. pos += 1
  245. if typed && x isa Float32
  246. buf[pos] = UInt8('f')
  247. buf[pos + 1] = UInt8('0')
  248. pos += 2
  249. end
  250. if typed && x isa Float16
  251. buf[pos] = UInt8(')')
  252. pos += 1
  253. end
  254. return pos
  255. end
  256. while hash && precision > 1
  257. buf[pos] = UInt8('0')
  258. pos += 1
  259. precision -= 1
  260. end
  261. if typed && x isa Float32
  262. buf[pos] = UInt8('f')
  263. buf[pos + 1] = UInt8('0')
  264. pos += 2
  265. end
  266. if typed && x isa Float16
  267. buf[pos] = UInt8(')')
  268. pos += 1
  269. end
  270. return pos
  271. elseif isnan(x)
  272. pos = append_sign(x, plus, space, buf, pos)
  273. buf[pos] = UInt8('N')
  274. buf[pos + 1] = UInt8('a')
  275. buf[pos + 2] = UInt8('N')
  276. if typed
  277. if x isa Float32
  278. buf[pos + 3] = UInt8('3')
  279. buf[pos + 4] = UInt8('2')
  280. elseif x isa Float16
  281. buf[pos + 3] = UInt8('1')
  282. buf[pos + 4] = UInt8('6')
  283. end
  284. end
  285. return pos + 3 + (typed && x isa Union{Float32, Float16} ? 2 : 0)
  286. elseif !isfinite(x)
  287. pos = append_sign(x, plus, space, buf, pos)
  288. buf[pos] = UInt8('I')
  289. buf[pos + 1] = UInt8('n')
  290. buf[pos + 2] = UInt8('f')
  291. if typed
  292. if x isa Float32
  293. buf[pos + 3] = UInt8('3')
  294. buf[pos + 4] = UInt8('2')
  295. elseif x isa Float16
  296. buf[pos + 3] = UInt8('1')
  297. buf[pos + 4] = UInt8('6')
  298. end
  299. end
  300. return pos + 3 + (typed && x isa Union{Float32, Float16} ? 2 : 0)
  301. end
  302. output, nexp = reduce_shortest(x, compact ? 999_999 : nothing)
  303. if typed && x isa Float16
  304. buf[pos] = UInt8('F')
  305. buf[pos + 1] = UInt8('l')
  306. buf[pos + 2] = UInt8('o')
  307. buf[pos + 3] = UInt8('a')
  308. buf[pos + 4] = UInt8('t')
  309. buf[pos + 5] = UInt8('1')
  310. buf[pos + 6] = UInt8('6')
  311. buf[pos + 7] = UInt8('(')
  312. pos += 8
  313. end
  314. pos = append_sign(x, plus, space, buf, pos)
  315. olength = decimallength(output)
  316. exp_form = true
  317. pt = nexp + olength
  318. if -4 < pt <= (precision == -1 ? (T == Float16 ? 3 : 6) : precision) &&
  319. !(pt >= olength && abs(mod(x + 0.05, 10^(pt - olength)) - 0.05) > 0.05)
  320. exp_form = false
  321. if pt <= 0
  322. buf[pos] = UInt8('0')
  323. pos += 1
  324. buf[pos] = decchar
  325. pos += 1
  326. for _ = 1:abs(pt)
  327. buf[pos] = UInt8('0')
  328. pos += 1
  329. end
  330. # elseif pt >= olength
  331. # nothing to do at this point
  332. # else
  333. # nothing to do at this point
  334. end
  335. else
  336. pos += 1
  337. end
  338. i = 0
  339. ptr = pointer(buf)
  340. ptr2 = pointer(DIGIT_TABLE)
  341. if (output >> 32) != 0
  342. q = output ÷ 100000000
  343. output2 = (output % UInt32) - UInt32(100000000) * (q % UInt32)
  344. output = q
  345. c = output2 % UInt32(10000)
  346. output2 = div(output2, UInt32(10000))
  347. d = output2 % UInt32(10000)
  348. c0 = (c % 100) << 1
  349. c1 = (c ÷ 100) << 1
  350. d0 = (d % 100) << 1
  351. d1 = (d ÷ 100) << 1
  352. memcpy(ptr, pos + olength - 2, ptr2, c0 + 1, 2)
  353. memcpy(ptr, pos + olength - 4, ptr2, c1 + 1, 2)
  354. memcpy(ptr, pos + olength - 6, ptr2, d0 + 1, 2)
  355. memcpy(ptr, pos + olength - 8, ptr2, d1 + 1, 2)
  356. i += 8
  357. end
  358. output2 = output % UInt32
  359. while output2 >= 10000
  360. c = output2 % UInt32(10000)
  361. output2 = div(output2, UInt32(10000))
  362. c0 = (c % 100) << 1
  363. c1 = (c ÷ 100) << 1
  364. memcpy(ptr, pos + olength - i - 2, ptr2, c0 + 1, 2)
  365. memcpy(ptr, pos + olength - i - 4, ptr2, c1 + 1, 2)
  366. i += 4
  367. end
  368. if output2 >= 100
  369. c = (output2 % UInt32(100)) << 1
  370. output2 = div(output2, UInt32(100))
  371. memcpy(ptr, pos + olength - i - 2, ptr2, c + 1, 2)
  372. i += 2
  373. end
  374. if output2 >= 10
  375. c = output2 << 1
  376. buf[pos + 1] = DIGIT_TABLE[c + 2]
  377. buf[pos - exp_form] = DIGIT_TABLE[c + 1]
  378. else
  379. buf[pos - exp_form] = UInt8('0') + (output2 % UInt8)
  380. end
  381. if !exp_form
  382. if pt <= 0
  383. pos += olength
  384. precision -= olength
  385. while hash && precision > 0
  386. buf[pos] = UInt8('0')
  387. pos += 1
  388. precision -= 1
  389. end
  390. elseif pt >= olength
  391. pos += olength
  392. precision -= olength
  393. for _ = 1:nexp
  394. buf[pos] = UInt8('0')
  395. pos += 1
  396. precision -= 1
  397. end
  398. if hash
  399. buf[pos] = decchar
  400. pos += 1
  401. if precision < 0
  402. buf[pos] = UInt8('0')
  403. pos += 1
  404. end
  405. while precision > 0
  406. buf[pos] = UInt8('0')
  407. pos += 1
  408. precision -= 1
  409. end
  410. end
  411. else
  412. pointoff = olength - abs(nexp)
  413. memmove(ptr, pos + pointoff + 1, ptr, pos + pointoff, olength - pointoff + 1)
  414. buf[pos + pointoff] = decchar
  415. pos += olength + 1
  416. precision -= olength
  417. while hash && precision > 0
  418. buf[pos] = UInt8('0')
  419. pos += 1
  420. precision -= 1
  421. end
  422. end
  423. if typed && x isa Float32
  424. buf[pos] = UInt8('f')
  425. buf[pos + 1] = UInt8('0')
  426. pos += 2
  427. end
  428. else
  429. if olength > 1 || hash
  430. buf[pos] = decchar
  431. pos += olength
  432. precision -= olength
  433. end
  434. if hash && olength == 1
  435. buf[pos] = UInt8('0')
  436. pos += 1
  437. end
  438. while hash && precision > 0
  439. buf[pos] = UInt8('0')
  440. pos += 1
  441. precision -= 1
  442. end
  443. buf[pos] = expchar
  444. pos += 1
  445. exp2 = nexp + olength - 1
  446. if exp2 < 0
  447. buf[pos] = UInt8('-')
  448. pos += 1
  449. exp2 = -exp2
  450. elseif padexp
  451. buf[pos] = UInt8('+')
  452. pos += 1
  453. end
  454. if exp2 >= 100
  455. c = exp2 % 10
  456. memcpy(ptr, pos, ptr2, 2 * div(exp2, 10) + 1, 2)
  457. buf[pos + 2] = UInt8('0') + (c % UInt8)
  458. pos += 3
  459. elseif exp2 >= 10
  460. memcpy(ptr, pos, ptr2, 2 * exp2 + 1, 2)
  461. pos += 2
  462. else
  463. if padexp
  464. buf[pos] = UInt8('0')
  465. pos += 1
  466. end
  467. buf[pos] = UInt8('0') + (exp2 % UInt8)
  468. pos += 1
  469. end
  470. end
  471. if typed && x isa Float16
  472. buf[pos] = UInt8(')')
  473. pos += 1
  474. end
  475. return pos
  476. end