Take the 2-minute tour ×
Programmers Stack Exchange is a question and answer site for professional programmers interested in conceptual questions about software development. It's 100% free, no registration required.

Is the following guaranteed to return true for all numerical and non-zero values of x?

bool IsRoundTrip(double x)
{
    double y = 1 / (1 / x);
    return x == y;
}

What conditions would cause a discrepency?

share|improve this question
3  
Depends how good your compiler/interpreter is at optimizing. It could work out that the assignment to y should be x and thus this function always returns true. Thus performing no actual evaluation at runtime. –  Loki Astari Apr 1 '13 at 4:01
2  
Please, try it out. –  zxcdw Apr 1 '13 at 4:16
3  
What result do you expect if you pass in x of 0.0 ? –  Carson63000 Apr 1 '13 at 4:22
2  
I would expect it to fail for a sufficiently small value of x. For instance in gcc it fails if I set x=1e-310, due to an overflow of 1/x. –  JSchlather Apr 1 '13 at 7:06
1  
@Carson63000, I would expect 0.0 to fail, which is why I stipulated "non-zero". –  Hand-E-Food Apr 1 '13 at 21:27
show 3 more comments

7 Answers

up vote 24 down vote accepted

To simplify things by defining a concrete implementation, I will assume (as other answers do) that we're talking about IEEE 754 64-bit floating point.

Each floating point number has three parts: a sign, an exponent, and a mantissa. (Technical details about hidden bits are irrelevant to this discussion).

Reciprocation doesn't affect the sign

1 / (2**e * m) = (1 / 2**e) * (1 / m) = 2**-e * (1 / m) , so there are two ways in which the double-reciprocation can fail to provide a fixpoint. The easy one is that the exponent can be an extreme value such that we move from a denormalised number to one which overflows. The second is that the mantissa can be a non-fixpoint of the double-reciprocation.

I wrote a simple program to test random mantissas:

import java.util.Random;

strictfp class RoundTrip
{
    public static void main(String[] args)
    {
        long one = Double.doubleToLongBits(1.0);
        Random rnd = new Random();
        for (int i = 1; i < 1<<30; i++) {
            long mantissa = rnd.nextLong() & 0xfffffffffffffL;
            double x = Double.longBitsToDouble(one + mantissa);
            double y = 1 / (1 / x);
            if (x != y) {
                System.out.println(Long.toHexString(one + mantissa));
                System.out.println(x);
                System.out.println(y);
                break;
            }
        }
    }
}

It quickly gave some output:

3ffeca41c09ebb2b
1.9243791126461456
1.9243791126461458

The program can be expected to find an answer if as few as 1 in 2**30 mantissas fail. With a slight modification, I found that about 17.15% of mantissas fail.


Slightly handwavy analysis:

There are 2**52-1 mantissas covering the open range (1, 2), and they're uniformly spaced. The same uniformly spaced mantissas cover the open range (0.5, 1), which contains the reciprocals. Note that in this range one unit in the last place (1ulp), i.e. the difference between consecutive values, has an absolute value half that of the ulp in the range (1, 2). But reciprocation isn't a linear operation, so in some parts of the range the density of values required is higher than in others. Therefore we expect that the reciprocation will not be injective.

Suppose values x and x+dx, both in (1, 2), differ by 1ulp. If they map to the same reciprocal mantissa, at most one of them can round-trip. What is the probability of this collision?

x^-1 differentiates to -x^-2, so the difference between 1/x and 1/(x+dx) is approximately -dx/x^2, or -2dx/x^2 ulps, so a difference of one ulp before reciprocation gives a difference of -2/x^2 ulps after reciprocation. Given that the separation between two exactly representable values is 1ulp (by definition), and assuming (for simplification) no particular alignment between mantissas and reciprocal mantissas, we can estimate the probability of a collision as max(0, 1 - 2/x^2), and we can approximate the proportion of collisions as \int_1^2 max(0, 1 - 2/x^2) dx = \int_{\sqrt 2}^2 (1 - 2/x^2) dx = 3 - 2\sqrt 2 is approximately 0.1716. This is in very good agreement with my empirical results for the proportion of mantissas that don't round-trip, so it seems reasonable to hypothesise that a mantissa will round-trip unless its reciprocal collides with that of another mantissa, in which case only one of the two will round-trip.

share|improve this answer
3  
So many terrible answers here, but this one actually addresses the specific question that was asked. I can't believe the sheer number of people who just said 'no, because of floating point precision' and thought that was good enough, without any proof or examples to show that the rounding errors from the two reciprocations wouldn't always cancel out (or as you would put it, that some non-fixed points of double-reciprocation of mantissas exist :P). A big +1, although this feels a little terse and intimidating at present. –  Mark Amery Apr 1 '13 at 16:01
1  
@MarkAmery: Why do you think that an example of two rounding errors not cancelling out is necessary? There are billions of single precision floating point numbers, why would we expect that every single time they were doubly-inverted we would get an exactly-cancelling rounding error? –  Michael Shaw Apr 1 '13 at 21:21
1  
@MichaelShaw, there's no good reason to expect that the reciprocal should even be injective, which would be a precondition for double-reciprocal to round-trip. However, there's a difference between a) stating without proof that there's probably at least one number that doesn't round-trip; b) calculating and justifying an estimate for the number of mantissas which don't round-trip; c) proving that at least one mantissa doesn't round-trip. I would be quite interested to see a careful analysis which shows how good or bad my empirical measurement is. –  Peter Taylor Apr 1 '13 at 21:41
    
@PeterTaylor: Given the multiple different rounding schemes that exist and the large number of possibilities it would have to work for, I find the probabilistic argument pretty convincing, although as you pointed out, it isn't proof. –  Michael Shaw Apr 1 '13 at 22:23
add comment

Besides rounding errors with either the first or the second division, IEEE floating point has + and - infinity and NaN (Not a Number), as well as zero, which doesn't even work with your formula for real numbers (meaning the mathematical abstraction, rather than a computer storage format).

Without getting into the odd special cases (zero, infinity, NaN), it seems likely that the only way you're going to get a true result is if there are two rounding errors that happen to cancel out or if the number is exactly representable in IEEE floating point. Since there are several different rounding schemes for IEEE floating point, and hardware could store more bits than the representation to try to keep things accurate, predicting when rounding would cancel out seems both difficult and hardware-dependent.

This is not the sort of thing that floating point was invented to deal with.

share|improve this answer
add comment

No. As a counterexample, 1/(1/49) works out to 49.00000000000000710542735760100185871124267578125 on my machine.


For a more abstract argument, let N be the number of representable floating-point numbers in the interval [1, 2). (In IEEE 754 double-precision, N happens to be 2^52) There is a trivial one-to-one mapping between this set and the floating-point numbers in the interval [1/2, 1): Just subtract 1 from the exponent. Thus, there are also N floating-point numbers [1/2, 1).

Within each interval [2^k, 2^(k+1)), floating-point numbers are equally-spaced. So:

  • With N floating point numbers in the interval [1/2, 1):
    • ~N/3 are in the interval [1/2, 2/3)
    • ~N/3 are in the interval [2/3, 5/6)
    • ~N/3 are in the interval [5/6, 1)
  • With N floating-point numbers in the interval [1, 2):
    • N/2 are in the interval [1, 3/2)
    • N/2 are in the interval [3/2, 2)

The N/2 floating-point numbers in the interval [3/2, 2) have reciprocals in the interval (1/2, 2/3]. But there are only ~N/3 floating-point numbers available in this range. Therefore, by the pigeonhole principle, there exists a pair of distinct floating-point numbers that have the same floating-point reciprocal.

share|improve this answer
1  
Nice use of the pigeonhole principle. It can be extended to get a bound equal to my estimate for the number of collisions. For x in (1, 2), the range (x, 2) maps onto the range (0.5, x^-1), so given the difference in ulp the proportion of collisions is at least (2 - x) - 2(x^-1 - 0.5) = 3 - x - 2x^-1. Differentiate: this has extrema when -1 + 2x^-2 = 0 i.e. x^2 = 2. Only the positive root is in range, so we get that the proportion of collisions is at least 3 - sqrt(2) - 2/sqrt(2) = 3 - 2sqrt(2). –  Peter Taylor Apr 3 '13 at 8:18
add comment

Rule of thumb: If the language has a type named "double", this will probably be false for most values of x.

"Double" is usually short-hand to "Double Precision Floating Point Number". By saying "Double precision", it implies "finite precision", i.e. rounding will occur at some point.

Some languages support "Arbitrary Precision Numbers", which allows storing any value that has a finite representation, and can usually do some math with these numbers without loss of precision. This still doesn't cover all options, as some numbers just don't have a finite representation - for example, 1/3 is close to 1.3333, but it can't be written down without some

To handle this, some languages also support "Symbolic numbers", which is basically means they keep the number as "1/3", and can do math with this form.

Using arbitrary precision or symbolic numbers has a high cost in performance, which is why it's not common in "mainstream" languages. Languages that support these forms also have a tendency to be untyped, and will only have join all number types under the name "Number" (But this is by no means related to supporting high-precision).

I most languages that have a type "double", it's identical (Or at lease very similar) to IEEE 754 Double Precision type.

Note: Your code might not fail for x=3, because although the program can't hold 1/3 precisely, it may also make a similar rounding error when dividing 1 by what it can hold, and come out with 3 again.

Note 2: You will find libraries for arbitrary/symbolic math for most languages that don't have them built in.

share|improve this answer
3  
Of course, 1/3 has a finite representation. In base 3, for example, its representation is 0.1. Perfectly finite. –  Jörg W Mittag Apr 1 '13 at 10:55
    
1/3 is also a bad example for symbolic algebra, as simple fractions suffice to get by with 1/3. In fact, with fractions OP's function is no problem at all, and very cheap (at worst, it copies the numerator and denominator twice, switching them each time). –  delnan Apr 1 '13 at 11:17
2  
-1, because you've just explained how the limited precision of floating points can lead to rounding errors without showing, or even trying to show, that it can cause the specific problem the question asked about. Without either a little testing or a little mathematics, it's pretty non-obvious to me whether the rounding error on the second division will cancel out the rounding error on the first division sometimes, always, or never. –  Mark Amery Apr 1 '13 at 15:47
    
Empirically your rule of thumb is wrong: it returns true for slightly over 80% of values of x. –  Peter Taylor Apr 1 '13 at 16:42
add comment

There exists the possibility that they will not be the same due to how the floating point numbers are represented.

Information on their numerical representation can be found at ... http://en.wikipedia.org/wiki/Floating-point

Most of the time, floating point numbers are normalized. Should the number be a denormalized float, your algorithm can be expected to return false. Whether it will or not, will depend upon whether your compiler chooses to recognize the mathematical equivalency and optimize based on that, OR to implement exactly what you have requested.

There may be other cases where your algorithm could return false, but that is the first that comes to my mind.

share|improve this answer
add comment

I'd assume this will produce true whenever the compiler optimizes the calculation away. Apart from that checking for equality is considered a programming error in case of doubles (unless your compiler/hardware combination warrants the opposite, which I would not dare to believe in). This is not limited to your example, but may also occur in other situations. Like writing a double to the DB and reading it back and comparing it with the original value you may still have in memory.

share|improve this answer
add comment

No, it is not always a round trip. In mathematics (i.e., the kind your write on paper with a pencil), it will be, because mathematics doesn't have a problem with limited precision. But in any finite-precision system (e.g., any computer), there will be quantities of the form x != 1/(1/x). It's downright trivial to prove - when x is 3, finite-precision decimal arithmetic cannot make 1/(1/3) come out to 3. 1/3 => 0.33333..., and no matter how many digits you extend it to, you'll still be incorrect. Finite-precision binary arithmetic has the same problem with other values.

share|improve this answer
3  
You cannot conclude from the fact that 1/x can't be represented exactly that 1/(1/x) != x. To take your example, in 3.s.f. decimal floating point, 1.00 / 3.00 = 0.333, and 1.00 / 0.333 is 3.00. That particular mantissa is a fixed point of the double-reciprocal in that floating point scheme. –  Peter Taylor Apr 1 '13 at 13:39
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.