2

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?

tachyon
  • 471
  • 3
  • 14
  • For `y` = 1 and `x` = 49, `s = y / x` will not work. It sets `s` to 0.0204081632653061208204636756136096664704382419586181640625, and then `floor(x * s)` is 0, because `x * s` is 0.99999999999999988897769753748434595763683319091796875. – Eric Postpischil Jun 24 '21 at 13:55
  • 2
    I suspect you could use `s = (y+.5) / x` to get a scale that works in most cases of interest, but further consideration is necessary, along with bounds on `y` and `x`. (For example, `s = (y+.5) / x` will fail when `y` is 2^64 and `x` is 49, because the `+.5` will have no effect, so the case is equivalent to the 1/49 failure above.) – Eric Postpischil Jun 24 '21 at 13:57
  • Link to formula now reports "AccessDenied" – chux - Reinstate Monica Jun 28 '21 at 12:57

2 Answers2

2

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
chux - Reinstate Monica
  • 143,097
  • 13
  • 135
  • 256
  • Interesting. So we can find out the range of `s` but analytically can not narrow it down further. I wonder if we can apply bisection method in this range. – tachyon Jun 25 '21 at 13:01
  • @tachyon I do not think bisection method would be sufficient.. – chux - Reinstate Monica Jun 28 '21 at 04:26
  • @chux-reinstate-monica Is it because of the limitation of the floating point precision or the fact that we can not correctly represent all real values (despite being in the range of precision)? What if we use arbitrary precision math library (such as [mpmath](https://mpmath.org/))? – tachyon Jun 28 '21 at 10:21
0

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
ibarrond
  • 6,617
  • 4
  • 26
  • 45
  • In you example `floor( x * y/x ) > y`. Is this behavior always the case or in some cases `floor( x * y/x ) < y` can also be a possibility? – tachyon Jun 24 '21 at 13:12
  • 1
    `y/x` fails for `y` = 1 and `x` = 49. Then `y/x` is computed as 0.0204081632653061208204636756136096664704382419586181640625, and multiplying that by `x` produces 0.99999999999999988897769753748434595763683319091796875, for which `floor` returns 0, not the desired 1. Also “Python's 64b floating points follow the IEEE-754 norm” is incorrect; the Python documentation leaves the floating-point format and some aspects of behavior up to the implementation. – Eric Postpischil Jun 24 '21 at 13:51
  • 1
    The precision of the IEEE-754 binary64 format is 53 bits, not 15-17. 15-17 is often cited as the number of decimal digits of precision, but this is [nonsensical and wrong](https://stackoverflow.com/questions/61609276/how-to-calculate-float-type-precision-and-does-it-make-sense/61614323#61614323). It is incorrect to say a precision is reduced by 1 bit; precision is a measure of the number of bits used to represent a value, not a measure of the accuracy of the value that is represented. – Eric Postpischil Jun 24 '21 at 13:53