Add support for compressed power-of-10 table

[vitaut]

Add support for compressed power-of-10 table similarly to Dragonbox. It is important for {fmt} which is often used on embedded systems.

[dougallj]

I did a few different prototypes of this a while ago. If you need exact values, which I think you do unless you re-verify the fast-path, the options I can come up with are >128-bit arithmetic, or a fix-up table.

These fix-up tables have the least values that'd need to be reverified if you removed the fix-up. Admittedly if you did the maths you might know some constraints I didn't and be able to generate something better (rounding in a different direction or whatever):
https://gist.github.com/dougallj/82284f0990473c43d2b42e233b6ad01a

These fix-up tables alter more outputs, but can be faster (some for 32-bit targets, some just by using powers of two for the mod): https://gist.github.com/dougallj/d73bd203bb6630233283f611d7390754

This is the >128-bit arithmetic version – more multiplies, but I think it might actually be faster since it loads one extra byte per compressed table entry rather than one extra bit per output value? Don't quite remember either way:
https://gist.github.com/dougallj/ce8230b8a327745453b134daa0e64099

[dougallj]

(And, while cataloguing old gists, these are slow table-free versions. As a drop-in for embedded devices where perf doesn't matter they might be interesting, but in a context where this was a good idea, I'd imagine you'd want to get rid of the zmij fast path and accept looser rounding, unless you also have a tiny result cache or something. I did do some better versions I can try to dig up sometime: https://gist.github.com/dougallj/9bd9834dcefe05f12bdc676c5447f1f9 )

[xjb714]

Regarding the conversion from f32 to string, I have a very ingenious algorithm that has been tested successfully. Only a 64-byte lookup table is required, and it is for little-endian CPUs. Code link: https://github.com/xjb714/xjb/blob/main/bench/xjb/xjb_comp.cpp

[xjb714]

@dougallj
In the dtoa algorithm, it is necessary to calculate the 128-bit value of pow10. How to balance the size of instructions and the size of lookup tables while maintaining a relatively good performance. How should this algorithm be designed? Are you interested in conducting research?

[xjb714]

@dougallj @vitaut
In the dtoa algorithm, the lookup table can be compressed to 368 bytes, with a compression ratio of approximately 27 times. It is necessary to set direct_via_lut_get_pow5 to 0. Code link: https://github.com/xjb714/xjb/blob/main/bench/xjb/xjb_comp.cpp

[vitaut]

Dragonbox-based compression has been implemented. Thanks @dougallj for the suggestions around fix-ups, I'll take a look. Alternative solution is to switch to plain Schubfach when optimizing for binary size. It only needs 126 bits of precision.

[vitaut]

For future reference: with the Dragonbox compression 5 cases are off by -1 and 305 cases are off by +1.

[jk-jeon]

@dougallj Is it correct understanding that, basically what you did is, instead of computing 5^offset * (5^(index/compression_ratio) aligned to the MSB) and then realigning it to the MSB, you compute (5^offset aligned to the MSB) * (5^(index/compression_ratio) aligned to the MSB) and then shift the result by 1-bit if needed? That's clever, why I did it in the dumb way 🤣

Regarding the table-free one, is there a reason to split the case of positive/negative exponents? Why wouldn't you just start with the largest negative exponent and multiply 5? Is 192-bit not enough precision if you do so?

By the way it's quite likely that the fixup table is not needed as in the case in Dragonbox. The fixup table was how I initially did it, but some error analysis revealed that it's okay to just unconditionally add 1 to the lower half of the resulting pow10 computation. Something similar is likely true here too.

To actually verify, we would need to do some math (as it's impossible to test all f64) but we could do random testing to know the answer beforehand we actually do the math. @vitaut I suppose you already did something like that, as you seem to had implemented the compression algorithm of Dragonbox before you replaced it by a better one?

[dougallj]

@jk-jeon

Is it correct understanding that, basically what you did is, instead of computing 5^offset * (5^(index/compression_ratio) aligned to the MSB) and then realigning it to the MSB, you compute (5^offset aligned to the MSB) * (5^(index/compression_ratio) aligned to the MSB) and then shift the result by 1-bit if needed?

Yep, that's the idea – it came pretty naturally from multiplying numbers straight from the table (or 192-significant-bit values) and dropping the zero-filled low bits.

Regarding the table-free one, is there a reason to split the case of positive/negative exponents? Why wouldn't you just start with the largest negative exponent and multiply 5? Is 192-bit not enough precision if you do so?

You absolutely can, there is enough precision (if you're a bit careful). It works and that's what my code to generate the table at compile time does (https://gist.github.com/dougallj/1ecb82ea57f7b8b0628229a0d2a603d4). So, why don't I?

1. My assumption is that exponents near zero are most common, numbers like 0.5 or 1.1 or whatever, and if we're using this at runtime rather than compile time that matters. We probably don't mind so much if subnormals and near-infinities are more expensive.
2. It's nice to avoid weird numbers: on ARM64 with 192-bit values it takes 9 instructions to load up that largest negative exponent, but only 6 to load up both 5 and 1/5.
3. We can get any exponent in O(log(n)) time rather than O(n) time, but the maths for that requires squaring to double the exponent, so you want to start with 5¹ or 5^-1 and square from there. I'm assuming the reason you'd want to multiply by 5 is that it is cheaper than multiplying by a number with more significant bits. Squaring isn't as cheap, but it is significantly cheaper than a full 192-bit multiplication (as, for example, you can reuse an x₀y₁ multiplication result for x₁y₀ if x=y). The version I need to dig up uses only squaring and multiplications by five to get to any exponent in O(log(n)) time, but still starts with either 5¹ and goes up, or 5^-1 and goes down.

By the way it's quite likely that the fixup table is not needed as in the case in Dragonbox.

I'd hope so, but I don't know how to prove it, and the fast-path is unnervingly close to not verifying (even as is, there's 1 case that the verification script thinks is wrong, but that turns out to give the correct result anyway).

[dougallj]

Ah, here it is, I think, my smallest/fastest O(log(n)) table-free version: https://gist.github.com/dougallj/2c02d12c6d8d5e5736e6c2f5a801ee19 (get_pow10_significand_tens)

~66 instruction or ~264 bytes on arm64: https://www.godbolt.org/z/aa7cY6xT1

[vitaut]

I suppose you already did something like that, as you seem to had implemented the compression algorithm of Dragonbox before you replaced it by a better one?

I didn't, which is why compression was initially disabled. My plan was to implement some form of fixup but Dougall beat me to it.