5

Suppose I have two floating point numbers, x and y, with their values being very close.

There's a discrete number of floating point numbers representable on a computer, so we can enumerate them in increasing order: f_1, f_2, f_3, .... I wish to find the distance of x and y in this list (i.e. are they 1, 2, 3, ... or n discrete steps apart?)

Is it possible to do this by only using arithmetic operations (+-*/), and not looking at the binary representation? I'm primarily interested in how this works on x86.

Is the following approximation correct, assuming that y > x and that x and y are only a few steps (say, < 100) apart? (Probably not ...)

(y-x) / x / eps

Here eps denotes the machine epsilon. (The machine epsilon is the difference between 1.0 and the next smallest floating point number.)

Szabolcs
  • 24,728
  • 9
  • 85
  • 174
  • @ShreevatsaR, that's why I divided by `x` too (not only `eps`). `(y-x) / x / eps = (y-x) / (x*eps)`. I'm still not sure it's correct though. – Szabolcs May 31 '11 at 14:33
  • Ah, ok, sorry. Your formula may be right; let us think. :-) In the meantime you can read [What Every Computer Scientist Should Know About Floating Point Arithmetic](http://cr.yp.to/2005-590/goldberg.pdf). – ShreevatsaR May 31 '11 at 14:51

2 Answers2

4

Floats are lexicographically ordered, therefore:

int steps(float a, float b){

  int ai = *(int*)&a;  // reinterpret as integer
  int bi = *(int*)&b;  // reinterpret as integer
  return bi - ai;
}

steps(5.0e-1, 5.0000054e-1);  // returns 9

Such a technique is used when comparing floating point numbers.

Arlen
  • 6,641
  • 4
  • 29
  • 61
  • 1
    as long as the numbers are > 0. Some logic is still needed to cover negative values and 0 (including the "-0"). – Alexey Frunze Jan 16 '12 at 05:54
  • +1, Interesting observation about lexicographic ordering. What I'd really like to know though is still if this is possible using only floating point arithmetic operations. The language I was using when I asked this doesn't support reinterpreting bit patterns as a different type. – Szabolcs Jan 16 '12 at 07:19
  • @Szabolcs what language were you using? I wrote the code in C++, but it should work in C and D. – Arlen Jan 16 '12 at 07:24
  • My question was not really about how to do it in that language, I was merely curious if it can be done by using floating point operations only, without resorting to examining the bit-representation of the float in some way. I was using Mathematica. Even if not a direct answer to my question, you answer was interesting and I appreciate it (especially considering how old the question is). – Szabolcs Jan 16 '12 at 07:28
1

You do not have to examine the binary representation directly, but you do have to rely on it to get an exact answer, I think.

Start by using frexp() to break x into exponent exp and mantissa. I believe the next float bigger than x is x + eps * 2^(exp-1). (The "-1" is because frexp returns a mantissa in the range [1/2, 1) and not [1, 2).)

If x and y have the same exponent, you are basically done. Otherwise you need to count how many steps there are per power of 2, which is just 1.0/eps. In other words, the number of steps between 2^n and 2^(n+1) is 1.0/eps.

So, for y > x, count how many steps there are from x to the next power of two; then count how many more steps it takes to get to the largest power of 2 less than y; then count how many more steps it takes to get from there up to y. All of these are pretty easily expressible in terms of eps, I think.

Nemo
  • 70,042
  • 10
  • 116
  • 153