do not ignore explicitly given mantissa width

[mnieper]

This fixes the issue raised in #866, making mantissa widths meaningful.

[mflatt]

I haven't yet understand the function completely, but it looks rounder function needs to be a little different to defend against very large bit widths. With this patch, 1|100000000000000000000000000000 runs out of memory, but that reads as 1.0 in current Chez Scheme.

[mnieper]

I haven't yet understand the function completely, but it looks rounder function needs to be a little different to defend against very large bit widths. With this patch, 1|100000000000000000000000000000 runs out of memory, but that reads as 1.0 in current Chez Scheme.

Chez also runs out of memory with #e0.1|10000000000000. That is somewhat unavoidable by the definition of ...|.... That number represents an exact rational number that is the best binary floating-point approximation to 0.1 with 100...000 significant binary digits, and that exact number has a huge numerator and denominator. The only way out here would be to raise an &implementation-restriction for large mantissa widths. Maybe that is the best course of action as mantissa widths in practice are much smaller than these huge numbers.

Your example, on the other hand, is an inexact number, so the final result won't take much space. The question is whether the calculation has to take that much space, especially in the less trivial case 0.1|1000000000000000000000. I would like to get the result of (inexact #e0.1|10000000000000000), which is the best 53-bit approximation to the best 10000000000000-bit approximation to 1/10.

[mflatt]

Chez also runs out of memory with #e0.1|10000000000000.

Here's what I'm seeing:

Chez Scheme Version 10.1.0-pre-release.2
Copyright 1984-2024 Cisco Systems, Inc.

> #e0.1|10000000000000
1/10

I don't see why there's inherently a problem here. As long as the number written before the | has a tractable number of digits, any request requested additional precision is just zeros, right?

[mnieper]

(Please excuse the delayed answer; it was night here.)

Have you tried #e0.1|100...000 with the (most recent version of the) patch I provided here? For example, I get

> #e0.1|100
1014120480182583521197362564301/10141204801825835211973625643008

and much larger denominators for higher mantissa widths.

It would be wrong to truncate the precision. The reason is that 1/10 has a period of length 4 in binary representation. In fact,

1/10 = 0.00011001100110011001100...

in binary representation. From this, we can deduce that

#e0.1|1 = 0.001
#e0.1|2 = 0.00011
#e0.1|3 = 0.00011
#e0.1|4 = 0.0001101
#e0.1|5 = 0.00011011
...

in binary representation, where I rounded to even. This means that larger and larger denominators (all powers of two) are needed the larger the mantissa width is.

[mflatt]

Sorry, I misunderstood what you meant by "Chez also", and i was confused about fractions and binary representations. Thank you for the tutorial! It makes sense that #e0.1|10000000000000 runs out of memory.

It still seems like 0.1|10000000000000 should not run out of memory. Is it a matter of setting a ceiling on precision when working toward for an inexact result, or is it more complex than that?

[mflatt]

Also, the results of #e1|10000000000000 and #e0.25|10000000000000 fit comfortably into memory, so it seems like they should be allowed, too. Is that a matter of detecting a power-of-two denominator, or are the cases when the result is representable more complicated to characterize?

[mnieper]

Sorry, I misunderstood what you meant by "Chez also"

Oh, I see. Indeed, what I wrote wasn't very clear.

It still seems like 0.1|10000000000000 should not run out of memory. Is it a matter of setting a ceiling on precision when working toward for an inexact result, or is it more complex than that?

Consider the following number in binary notation (where N* means to repeat the following binary digit N times):

0.1 52*0 1 N*0 1

Let x be its decimal representation.

If N is sufficiently larger than p and q sufficiently larger than N, we have that #x|p is 0.1 in binary and #x|q is 0.1 51*0 1 in binary. In other words, we cannot simply truncate a huge mantissa width like q to some smaller one p without examining the value of x.

[mnieper]

Also, the results of #e1|10000000000000 and #e0.25|10000000000000 fit comfortably into memory, so it seems like they should be allowed, too. Is that a matter of detecting a power-of-two denominator, or are the cases when the result is representable more complicated to characterize?

I do not have a full characterisation. But a denominator N means that the quotient has a period length of at most N - 1, so one should be able to reduce the case of an arbitrary mantissa width roughly to the case of a mantissa width <= 2*N.

But I wonder whether it makes sense to spend the time getting the details right and writing extra code for huge mantissa widths. In practice, the largest mantissa widths may come from when using libraries like GNU MPFR. Do you think someone would use floats that use megabytes of memory?

[mnieper]

Maybe it is not that complicated to actually implement mantissa truncation.

1. For exact numbers, this is only possible when the denominator is a power of two because otherwise, we have an infinite period in the binary representation. When the denominator is a power of two, we can truncate at the log2 of the denominator.
2. For inexact numbers, we calculate the period length P according to the denominator. We then add 53 and the log2 of the numerator, getting some number N. If the mantissa width p is then > N + P, we can replace it with p - P.

The estimates can possibly be off by 1 or 2, but one can code safely.

[mflatt]

That sounds really great. As we've established, I'm not clear on the math, but it certainly sounds plausible.

My experience with Scheme numbers is that these details end up being worthwhile, even though it means extra code, and even through the happy spaces often end up being complex (e.g., only power-of-two denominators). Unfortunately, my experience with Scheme numbers is also that I have to learn a lot of new things, and then I forget them soon afterward!

[mflatt]

Hi @mnieper — I pushed a commit to add precision bounding in (I think) the way you describe. Does it look right? Do I understand correctly that this captures all of the cases where the end result number can be represented with about the same amount of memory as the number without a precision adjustment?

[mnieper]

Thanks a lot, Matthew. I am going to take a look at it within the next few days. (I wanted to have come up with some code as well but haven't found the time.)

[mflatt]

@mnieper Have you had a chance to take a look? (I have been thinking that it would make sense to create a v10.1.0 release soon, but waiting for this change.)

[mnieper]

@mnieper Have you had a chance to take a look? (I have been thinking that it would make sense to create a v10.1.0 release soon, but waiting for this change.)

Sorry, this somehow fell off my radar. I have some comments (and, hopefully, an improvement of my earlier analysis), which I will post as soon as possible.

[mnieper]

My original analysis of the inexact case (which needed the period length and which you didn't implement) was too complicated; the simpler analysis here should suffice (and leads to a simple formula).

[mnieper]

If b is a power of two, the number of significant bits of n is just the number of significant bits of a. (In the issue, I wrote erroneously "denominator" instead of "numerator" at one point.) Thus, (min p a-bits) should also work.

[mnieper]

We are in the case where a binary fraction x is first rounded to p significant bits and then to 53 significant bits (round to even in case of a tie). This is, in general, not the same as rounding directly to 53 significant bits. The case where it may differ looks like X.1Y. Here, X represents 53 significant binary digits, the dot is the point of rounding and Y are all further binary digits. The rounding of X.1Y to 53 binary digits depends on whether Y vanishes or not. If we first round X.1Y to p significant digits, it may happen that a nonzero Y becomes zero after rounding. Because of this, we must make sure that if we truncate p, the answer to the question of whether Y becomes zero or not won't change.

Every binary fraction is a transient (non-repeating digits, possibly trivial) followed by a repetend (repeating digits, possibly just zeros). We are allowed to truncate p when p points into the repetend and remains there after truncation (because this won't change whether Y becomes zero after rounding). So, for this, we have to estimate the length of the transient. This has at most (+ a-bits b-bits) bits. Together with the 53 bits for X and 2 bits for a safety measure to exclude possible edge cases, the formula should therefore be (min p (+ a-bits b-bits 53 2)).

[mnieper]

PS I didn't mention that the repetend may start with zero digits. The analysis remains correct when we count the zero digits as part of the transient.

[mflatt]

@mnieper Thanks for the corrections! I've made adjustments and rebased. Can you check one last time to make sure I didn't mangle the change?

I wasn't able to find any inputs that produce different results before an after the last change, trying millions of samples. That doesn't mean much, since my grasp of the arithmetic is weak. But my fuzzing script was able to find counterexamples within 10000 samples when I mangled the exact case to use (- a-bits 1) instead of a-bits or mangled the inexact case to use (+ 53 2), (+ a-bits 53 2), or (+ b-bits 53 2) instead of (+ a-bits b-bits 53 2). So, I'm 75|1 percent confident that my script fuzzing could find a mismatch if there were one.

[mnieper]

My analysis was very conservative, so it may very well be that there is no counterexample for your earlier code. I just want to make sure that we are on the safe side. We should improve the estimate once we have a formal proof that a more tight estimate also works. But I think this does not need to be done in the coming Chez version.

I think I am 100|7 percent sure that you didn't mangle the change.