Error in rational -> float conversion

[snunez1]

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

Interestingly, taking one digit off the rational produces a correct result:

CL-USER> (float 4110710000054127/10000000000)
411071.0

[informatimago]

Indeed.

If you want a double-float conversion, you can use (float x 0.0d0).

But it would be nice if coerce used float...

On macOS intel:

cl-user> (float 41107100000541273/100000000000)
411071.03
cl-user> (float 41107100000541273/100000000000 0.0d0)
411071.00000541273D0

It would be nice too if we had: (= (float (float 41107100000541273/100000000000 0.0d0) 0.0e0) (float 41107100000541273/100000000000 0.0e0)

cl-user> (float (float 41107100000541273/100000000000 0.0d0) 0.0e0)
411071.0

Otherwise note that the result obtained is not too bad. It's not the closest single-float, but it's the closest above:

cl-user> (let ((f 411071.03))
           (multiple-value-bind (signif expon sign) (integer-decode-float f)
             (* sign (scale-float (float signif f) expon))))
           
411071.03
cl-user> (let ((f 411071.03))
           (multiple-value-bind (signif expon sign) (integer-decode-float f)
             (* sign (scale-float (float (1+ signif) f) expon))))
411071.06
cl-user> (let ((f 411071.03))
           (multiple-value-bind (signif expon sign) (integer-decode-float f)
             (* sign (scale-float (float (1- signif) f) expon))))
411071.0
cl-user> 

[snunez1]

Are there any workarounds? Returning 411071.03 when the value is 411071.00000541273 is a rather significant error in many situations. Just confirmed that SBCL behaviour is 'correct' (i.e., returns the closest).

[informatimago]

Are there any workarounds? Returning 411071.03 when the value is 411071.00000541273 is a rather significant error in many situations. Just confirmed that SBCL behaviour is 'correct' (i.e., returns the closest).

Yes, use: (float (float 41107100000541273/100000000000 0.0d0) 0.0e0)
411071.0

[metayan]

Found a big difference:

(let ((a 411071000000000101) (b 100000000000)) (- (float (/ a b)) (float (/ a b) 0d0)))
0.24999999860301614D0

Lisp	Difference
SBCL	-0.0000000009313225746154785
ECL	-0.0000000009313225746154785
CLISP	-0.0000000009313226
ABCL	-0.0000000013969838619232178
CCL	0.24999999860301614

Formatted with "~f".

[kwccoin]

Found a big difference:

(let ((a 411071000000000101) (b 100000000000)) (- (float (/ a b)) (float (/ a b) 0d0)))
0.24999999860301614D0

Lisp Difference
SBCL -0.0000000009313225746154785
ECL -0.0000000009313225746154785
CLISP -0.0000000009313226
ABCL -0.0000000013969838619232178
CCL 0.24999999860301614
Formatted with "~f".

I noted ECL actually give the same answers. I used mainly SBCL, except ECL has support over IOS ... wonder I may try to side test also ECL. Let us see.

[mirkov]

I don't understand what the issue is - maybe someone can explain it:

Floating point numbers carry between 6-7 digits of precision (as per IBM page).

The single precision result of 411071.03 is correct to 7 digits. The last digit shown (3) is the eight digit and should not be taken into account.

To me this looks like a user error - not understanding floating point precision. But then again, maybe I am missing something.

Thanks,

[olopierpa]

On Sun, Jul 3, 2022 at 2:29 AM Mirko ***@***.***> wrote: I don't understand what the issue is - maybe someone can explain it: Floating point numbers carry between 6-7 digits of precision (as per IBM page <https://www.ibm.com/support/pages/single-precision-floating-point-accuracy#:~:text=A%20single%2Dprecision%20float%20only,left%20for%20floating%20part%20precision.> .

SINGLE floats. The problem is with DOUBLE floats.

…

Message ID: ***@***.***>

[mirkov]

Sorry, I still don't get it:

; SLIME 2.26.1
CL-USER> (lisp-implementation-type)
"Clozure Common Lisp"
CL-USER> (lisp-implementation-version)
"Version 1.12.1 (v1.12.1) LinuxX8664"
CL-USER> (float 41107100000541273/100000000000 0.d0)
411071.00000541273D0

[rtoy]

I think the issue is that 411071.0 is closer to the fraction than 411071.03. 411071.03 and 411071 differ by 1 bit.

Consider:

(abs (- 41107100000541273/100000000000 411071))
541273/100000000000
(abs (- 41107100000541273/100000000000 (rational 441071.03)))
3124458727/100000000000

It's pretty clear that the error from the first is less than the second. So you would want (float 41107100000541273/100000000000) to be 411071.0 and not 411071.03 since it's closer to the rational number than 411071.03.

[LdBeth]

This is demonstrated by:

;; Clozure Common Lisp Version 1.12.2  DarwinX8664
CL-USER> (/ 41107100000541273.0 100000000000.0 )
411071.03
CL-USER> (/ 41107100000541273.D0 100000000000.D0 )
411071.00000541273D0
CL-USER> (= (float 41107100000541273/100000000000) (/ 41107100000541273.0 100000000000.0 ))
T
CL-USER> (coerce 41107100000541273/100000000000 'double-float)
411071.00000541273D0

So, it appears for CCL the float converts the rational to single-float by first convert the denominator and numerator to single precision float, thus the loss of precision.

;; This is SBCL 2.3.0
* (float 41107100000541273/100000000000)
411071.0
* (/ 41107100000541273.0 100000000000.0) 
411071.03
* (coerce 41107100000541273/100000000000 'single-float)
411071.0

SBCL's float which called coerce function used more clever tricks. Unfortunately it is buried into too may layers of abstractions so I cannot get it, but for CCL, the actual job is done by %short-float-ratio defined in level-0/l0-float.lisp, and below is the 64 bit version which appears can prove by guess.

https://github.com/Clozure/ccl/blob/fa2bbe8bcb51eef65d4e37cea9d0a86c9a8120d5/level-0/l0-float.lisp#L491C18-L491C18

If anyone have more experience with SBCL, could you please give me some hints on where I could find the rational to float conversion?

[LdBeth]

I think this is SBCL's implementation that coverts rational to float

https://github.com/sbcl/sbcl/blob/ac267f21721663b196aefe4bfd998416e3cc4332/src/code/float.lisp#L444

See the comments above it on how they attempt to avoid inaccuracies.

It seems what SBCL did is in fact very simple: down scale the numerator and denominator by power of 2 if they exceeds the mantissa size.

When I have time I might write a patch to CCL to fix the issue.

[svspire]

There are [at least] two issues here:

1. #'%short-float-ratio makes assumptions that were true in 32-bit CCL but are not true in 64-bit CCL.
2. CCL's #'rationalize is not very good.

Let's take these one at a time:

There are several places in CCL that assume a fixnum has roughly the same number of significant bits as a short-float, which was roughly true in the days of CCL-32:

? (integer-length x8632::target-most-positive-fixnum)
29

? (float-digits 1.0s0)
24

Thus in CCL-32 if you convert a fixnum to a short-float you won't lose much precision. Unfortunately, this assumption is quite invalid in CCL-64:

? (integer-length x8664::target-most-positive-fixnum)
60

So in CCL-64, if you convert a fixnum to a short-float, you can lose a great deal of precision. #'%short-float-ratio makes this invalid assumption. Ouch.

In principle it's easy to fix #'%short-float-ratio as follows:

#+64-bit-target
(defun %short-float-ratio (number)
  (let* ((num (numerator number))
         (den (denominator number)))
    ;; dont error if result is floatable when either top or bottom is
    ;; not.  maybe do usual first, catching error
    (if (not (or (bignump num)(bignump den)))
        (%short-float (/ (the double-float (%double-float num)) ;;; <--- this is the fix
                         (the double-float (%double-float den))))
      (let* ((numlen (integer-length num))
             (denlen (integer-length den))
             (exp (- numlen denlen))
             (minusp (minusp num)))
        (if (and (<= numlen IEEE-single-float-bias)
                 (<= denlen IEEE-single-float-bias)
                 #|(not (minusp exp))|# 
                 (<= (abs exp) IEEE-single-float-mantissa-width))
          (/ (the short-float (%short-float num))
             (the short-float (%short-float den)))
          (if (> exp IEEE-single-float-mantissa-width)
            (progn  (%short-float (round num den)))
            (if (>= exp 0)
              ; exp between 0 and 23 and nums big
              (let* ((shift (- IEEE-single-float-digits exp))
                     (num (if minusp (- num) num))
                     (int (round (ash num shift) den)) ; gaak
                     (intlen (integer-length int))
                     (new-exp (+ intlen (- IEEE-single-float-bias shift))))
		(when (> intlen IEEE-single-float-digits)
                  (setq shift (1- shift))
                  (setq int (round (ash num shift) den))
                  (setq intlen (integer-length int))
                  (setq new-exp (+ intlen (- IEEE-single-float-bias shift))))
                (when (> new-exp IEEE-single-float-normal-exponent-max)
                  (error (make-condition 'floating-point-overflow
                                         :operation 'short-float
                                         :operands (list number))))
                (make-short-float-from-fixnums 
                   (ldb (byte IEEE-single-float-digits  (- intlen  IEEE-single-float-digits)) int)
                   new-exp
                   (if minusp -1 1)))
              ; den > num - exp negative
              (progn  
                (float-rat-neg-exp num den (if minusp -1 1) nil t))))))))

After this fix,
(float 41107100000541273/100000000000)
returns
411071.0

And now CCL-64 passes @metayan's test above too, returning the same value as ABCL.

But there's a problem, which brings us to the second issue. Now CCL's RATIONALIZE.1 test fails.
Here's the failure:

? (float (rationalize 8.131823E-18) 8.131823E-18)
8.131824E-18

Only the least significant digit is wrong, but it's enough to make (float (rationalize x) x) not be true, and the ANSI spec requires it to be true.

I've investigated this issue enough to know that nothing is wrong with #'/ in CCL; the problem is CCL's #'rationalize implementation. It may be that #'rationalize is making similar (i.e. wrong) assumptions about short-floats and fixnums. In any case, here's what I think needs to be done:

1. Investigate #'rationalize in CCL and figure out why it's losing. Look at Bruno Haible's #'rationalize in the SBCL source and compare it. It's pretty good, and maybe we should just adopt it.
2. Look for other places in the CCL-64 source that assume that short-floats have similar precision as fixnums. Correct these mistakes.
3. Adopt the above patch to #'%short-float-ratio.

[LdBeth]

@svspire I think the rational to double conversion needs also been patched and follows SBCL’s implementation that to single and to double conversion is treated in the same way (only the idea, not necessarily using same implementation). It is a bad idea assuming nothing would be wrong if double-precision is used.

I think the proper way to fix this is replace the bignump test by more proper conditions. Which as you observed is not a valid assumption for x86-64

also, I doubt the in consistencies for rationalize is caused by incorrect rounding in rational to float conversion instead of a flaw in rationalize implementation. I have only viewed the function definition but cannot notice any problems yet, but I’ll test them after I get back to the computer.

[LdBeth]

Adopting the solution from SBCL is easier than I think:

LdBeth@faaed9f

and this does not suffer from the rationalize problem that the ad hoc fix @svspire would produce

Clozure Common Lisp Version 1.12.2  DarwinX8664

For more information about CCL, please see http://ccl.clozure.com.

CCL is free software.  It is distributed under the terms of the Apache
Licence, Version 2.0.
? (float (rationalize 8.131823E-18) 8.131823E-18)
8.131823E-18
? (float 41107100000541273/100000000000)
411071.0
? (format t "~f" (let ((a 411071000000000101) (b 100000000000)) (- (float (/ a b)) (float (/ a b) 0d0))))
-0.0000000013969838619232178

Now the question would be should I do the same to %double-float function which shares almost same code with %short-float-ratio.