[Openmcl-devel] Error in rational -> float conversion
Shannon Spires
svs at bearlanding.com
Mon Sep 20 11:44:01 PDT 2021
I wondered the same thing. This is exactly how #'%short-float-ratio
works and I assume it was done for speed reasons. The other big mistake
in #'%short-float-ratio is probably that it immediately checks whether
the numerator or denominator is a bignum. That check might have been
sensible in 32-bit implementations since in those implementations the
number of significant bits in a fixnum is roughly the same as those in a
short-float. But such a check doesn't make sense in 64-bit implementations:
Clozure Common Lisp Version 1.12 (v1.12-2-g5d13fc7d) DarwinX8664
? (integer-length most-positive-fixnum)
60
? (float-digits 1.0s0)
24
Thus by first converting the numerator and denominator to short-floats,
#'%short-float-ratio is losing a great deal of precision immediately,
even when both numerator and denominator are fixnums.
Assuming we fixed that, we are left with the monotonicity issue that Ron
points out:
; Here's Ron's original example that uses rational numbers (more
precisely, /ratios/):
(let ((a 4110710000054) (b 100000000000))
(list (float (/ a b)) (float (/ (+ a 1) b)) (float (/ (+ a 2) b))))
; (41.1071 41.107098 41.1071) ; weird nonmonotonicity
; Reframed to not use any ratios:
(let ((a 4110710000054) (b 100000000000))
(list (float (/ (float a 1.0s0)
(float b 1.0s0)))
(float (/ (float (+ a 1) 1.0s0)
(float b 1.0s0)))
(float (/ (float (+ a 2) 1.0s0)
(float b 1.0s0)))))
; (41.1071 41.1071 41.1071) ; [non-strict] monotonicity returns
This might be what Ron meant when he said "I'm an idiot" but it also
took me a while to figure it out: As Ron demonstrated, ratios are
[effectively] always stored in lowest terms. Because of this, they
bounce around a lot when you increment the numerator. Which in turn
causes unpredictable losses of precision: The variance of said losses
will be large if you always convert numerator and denominator to short
floats before you divide them. But if you force that lowest terms
mechanism not to happen, monotonicity returns (second example above). So
doing ratio conversions with double-floats will be a big help. It won't
make the nonmonotonicity go away entirely -- there's no way to do that
because the lowest terms mechanism will still be present -- but it will
decrease the nonmonotonicity dramatically.
The obvious next step is to remove the bignum check in
#'%short-float-ratio. I'll try to test that when I have time. The test
is not trivial because it requires a full rebuild of CCL.
-SS
On 9/20/21 11:28 AM, Ron Garret wrote:
> Is there actually a reason, when converting rationals to
> single-floats, that the computation is done by converting the
> numerator and denominator to single-floats and dividing, rather than
> converting them to double-floats, dividing, and then converting the
> result to a single-float? The former doesn’t seem to be appreciably
> faster on modern hardware:
>
> ? (time (dotimes (i 1000000)
> (let ((x (/ (random 1000000000000) (random 1000000000000))))
> (coerce x 'single-float))))
> (DOTIMES (I 1000000) (LET ((X (/ (RANDOM 1000000000000) (RANDOM
> 1000000000000)))) (COERCE X 'SINGLE-FLOAT)))
> took 525,951 microseconds (0.525951 seconds) to run.
> 3,967 microseconds (0.003967 seconds, 0.75%) of which was spent
> in GC.
> During that period, and with 4 available CPU cores,
> 534,009 microseconds (0.534009 seconds) were spent in user mode
> 6,094 microseconds (0.006094 seconds) were spent in system mode
> 32,000,000 bytes of memory allocated.
> 42 minor page faults, 0 major page faults, 0 swaps.
> NIL
> ? (time (dotimes (i 1000000)
> (let ((x (/ (random 1000000000000) (random 1000000000000))))
> (coerce (the double-float
> (/ (the doble-float (coerce (numerator x)
> 'double-float))
> (the double-float (coerce (denominator x)
> 'double-float))))
> 'single-float))))
> ;Compiler warnings :
> ; In an anonymous lambda form at position 158: Unknown type
> DOBLE-FLOAT, declaration ignored
> (DOTIMES (I 1000000) (LET ((X (/ (RANDOM 1000000000000) (RANDOM
> 1000000000000)))) (COERCE (THE DOUBLE-FLOAT (/ (THE DOBLE-FLOAT
> (COERCE (NUMERATOR X) 'DOUBLE-FLOAT)) (THE DOUBLE-FLOAT (COERCE
> (DENOMINATOR X) 'DOUBLE-FLOAT)))) 'SINGLE-FLOAT)))
> took 542,226 microseconds (0.542226 seconds) to run.
> 8,720 microseconds (0.008720 seconds, 1.61%) of which was spent
> in GC.
> During that period, and with 4 available CPU cores,
> 548,880 microseconds (0.548880 seconds) were spent in user mode
> 14,421 microseconds (0.014421 seconds) were spent in system mode
> 64,000,000 bytes of memory allocated.
> 2,926 minor page faults, 0 major page faults, 0 swaps.
>
>
> On Sep 20, 2021, at 10:16 AM, Ron Garret <ron at flownet.com> wrote:
>
>> Oh never mind, I’m an idiot.
>>
>> ? 4110710000054/100000000000
>> 2055355000027/50000000000
>> ? 4110710000055/100000000000
>> 822142000011/20000000000
>> ? 4110710000056/100000000000
>> 513838750007/12500000000
>>
>> On Sep 20, 2021, at 9:53 AM, Ron Garret <ron at flownet.com> wrote:
>>
>>> All this is true, but it does not account for the anomaly. This
>>> appears to be a legitimate bug:
>>>
>>> ? (let ((a 4110710000054) (b 100000000000))
>>> (list (float (/ a b)) (float (/ (+ a 1) b)) (float (/ (+ a 2) b))))
>>> (41.1071 41.107098 41.1071)
>>> ? (mapcar 'decode-float *)
>>> (0.64229846 0.6422984 0.64229846)
>>>
>>> It is very data dependent, only manifesting itself for very
>>> particular values of A. I have not yet been able to discern any
>>> pattern, but here are a few interesting edge cases:
>>>
>>> ? (let ((a 7110710000054) (b 100000000000))
>>> (list (float (/ a b)) (float (/ (+ a 1) b)) (float (/ (+ a 2) b))))
>>> (71.1071 71.10709 71.1071)
>>>
>>> ? (let ((a 7110720000054) (b 100000000000))
>>> (list (float (/ a b)) (float (/ (+ a 1) b)) (float (/ (+ a 2) b))))
>>> (71.1072 71.1072 71.1072)
>>>
>>> ? (let ((a 4110709000054) (b 100000000000))
>>> (list (float (/ a b)) (float (/ (+ a 1) b)) (float (/ (+ a 2) b))))
>>> (41.10709 41.10709 41.10709)
>>>
>>> AFAICT the bug is in ccl::%short-float-ratio. I’ll track it down
>>> when I have time, but with the above test cases it should not be too
>>> hard to figure it out.
>>>
>>> rg
>>>
>>> On Sep 19, 2021, at 10:56 PM, Raymond Wiker <rwiker at gmail.com> wrote:
>>>
>>>> For IEEE-754 single-precision floating-point, the mantissa is held
>>>> in 23 bits. log10 of 2^23 is about 6.9237, so you should expect
>>>> less than 7 decimal digits of accuracy. 411071.0 is 7 digits, so
>>>> anything after that can (and should) be ignored. The difference
>>>> between 411071.0 and 411071.03 expressed as a float32 is a single
>>>> bit: the least significant bit of the mantissa.
>>>>
>>>> The reason that coerce gives a slightly worse result is that you
>>>> compute the value as float32, then convert to float64, and when you
>>>> print the number, more digits will be printed.
>>>>
>>>>
>>>>> On 19 Sep 2021, at 15:05, Steven Nunez <steve_nunez at yahoo.com> wrote:
>>>>>
>>>>> Hi folks,
>>>>>
>>>>> I just posted an issue regarding what looks like a bug in
>>>>> rational->float conversion
>>>>> <https://github.com/Clozure/ccl/issues/390>, and I'm hoping
>>>>> someone here can try this on 1.12 (preferably on MS Windows), or
>>>>> perhaps provide some insight as to what might be going on here for
>>>>> a quick fix.
>>>>>
>>>>> CL-USER> (lisp-implementation-version)
>>>>> "Version 1.11.8 (v1.11.8-2-gd411e378) WindowsX8664"
>>>>> CL-USER> (float 41107100000541273/100000000000)
>>>>> 411071.03
>>>>>
>>>>> Hmm. If you work out the decimal points, the float is:
>>>>> 411071.00000541273, that's a long way from .03. Maybe coerce?
>>>>> Nope, that's even worse:
>>>>>
>>>>> CL-USER> (coerce 41107100000541273/100000000000 'double-float)
>>>>> 411071.03125D0
>>>>>
>>>>> Any ideas? That's not a small rounding error.
>>>>>
>>>>> Regards,
>>>>> Steve
>>>>> _______________________________________________
>>>>> Openmcl-devel mailing list
>>>>> Openmcl-devel at clozure.com
>>>>> https://lists.clozure.com/mailman/listinfo/openmcl-devel
>>>>
>>>> _______________________________________________
>>>> Openmcl-devel mailing list
>>>> Openmcl-devel at clozure.com
>>>> https://lists.clozure.com/mailman/listinfo/openmcl-devel
>>>
>>> _______________________________________________
>>> Openmcl-devel mailing list
>>> Openmcl-devel at clozure.com
>>> https://lists.clozure.com/mailman/listinfo/openmcl-devel
>>
>> _______________________________________________
>> Openmcl-devel mailing list
>> Openmcl-devel at clozure.com
>> https://lists.clozure.com/mailman/listinfo/openmcl-devel
>
>
> _______________________________________________
> Openmcl-devel mailing list
> Openmcl-devel at clozure.com
> https://lists.clozure.com/mailman/listinfo/openmcl-devel
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.clozure.com/pipermail/openmcl-devel/attachments/20210920/45cb7448/attachment.htm>
More information about the Openmcl-devel
mailing list