Correctly determining scaling factor

Viewed 267

A function determines y(integer) from given x (integer) and s (float) as follows:

floor(x * s)

If x and y are known how to calculate s so that floor(x * s) is guaranteed to be exactly equal to y.

If I simply perform s = y / x is there any chance that floor(x * s) won't be equal to y due to floating point operations?

2 Answers

If I simply perform s = y / x is there any chance that floor(x * s) won't be equal to y due to floating point operations?

Yes, there is a chance it won't be equal. @Eric Postpischil offer a simple counter example: y = 1 and x = 49.


(For discussion, let us limit x,y > 0.)

To find a scale factor s for a given x,y, that often works, we need to reverse y = floor(x * s) mathematically. We need to account for the multiplication error (see ULP) and floor truncation.

 # Pseudo code

 e = ULP(x*s)
 y <  (x*s + 0.5*e) + 1
 y >= (x*s - 0.5*e)

 # Estimate e
 est = ULP((float)y)
 s_lower = ((float)y - 1 - 0.5*est)/(float)x
 s_upper = ((float)y     + 0.5*est)/(float)x 

A candidate s will lie s_lower < s <= s_upper.

Perform the above with higher precision routines. Then I recommend to use the float closest to the mid-point of s_lower, s_upper.

Alternatively, an initial stab at s could use:

s_first_attempt = ((float)y - 0.5)/(float)x

If we rephrase your question, you are wondering if the equation y = floor( x * y/x ) holds for x and y integers, where y/x translates in python into a 64-bit floating-point, and the subsequent multiplication also generates a 64b floating point value.

Python's 64b floating points follow the IEEE-754 norm, which gives them 15-17 bits of decimal precision. To perform the division and multiplication, both x and y are converted into floats, and these operations might reduce the minimum precision in up to 1 bit (really worst case), but they will for sure not increase the precision. As such, you can only expect up to 15-17 bits of precision in this operation. This means that y values above 10^15 might/will present rounding errors.

More practically, one example of this can be (and you can reuse this code for other examples):

import numpy as np
print("{:f}".format(np.floor(1.3 * (1.1e24 / 1.3))))
#> 1100000000000000008388608.000000
Related