In article <firstname.lastname@example.org>, ddtl
> Hello everybody.
> I am trying to find a real number x so that x is the smallest
> number bigger than 1 which can be represented by the Matlab but
> which doesn't have an inverse which could be represented in Matlab.
> Now, the smallest number t which could be added to 1 so that 1 ~= (1+t)
> is 2^-52. So, I wrote the following program:
> x = 1;
> y = x;
> while x == y
> x = x + 2^-52;
> y = 1/(1/x);
> Obviously, that program should find the required number. The
> problem is that it take a *long* time to find that number (I tried
> to run the program for 6 hours on Athlon 1333, and still no success).
> Is it somehow possible to find the required number quicker?
I can't imagine what led you to such a problem as this, ddtl, but it
does have a definite answer, and that is the number x =
1.414213573163048743, or in matlab's format hex, 3ff6a09e6964b6ac. I
don't wonder that your while loop wasn't finished after a mere six hours,
since it would have required some 1.865e+15 (almost two quadrillion) trips
through the loop to reach that number, starting as far back as x = 1 and
only adding 2^(-52) each time.
What makes a solution possible in practical terms is the fact that no
value of x below sqrt(2) can possibly fail to achieve equality of x and
1/(1/x). To demonstrate this, start with the assumption that 1 < x < 2.
Then you can write the result of the first reciprocal as:
t = 1/x + e1
where e1 is the round off error, and where, since 1/2 < t < 1, we have
abs(e1) <= 2^(-54). This is a guaranteed property of division with double
precision IEEE 754 numbers - that is, the round off error cannot exceed
half the value of the least bit position, which for t is worth 2^(-53).
The next step is:
y = 1/t + e2
where, since 1 < y < 2, we have abs(e2) <= 2^(-53), again because of round
off precision requirements.
Given the meaning of e1 and e2, the previous two equations are
mathematically exact equations, and we can manipulate with them as such.
Combining them gives
y = 1/(1/x + e1) + e2 = x - (x^2/(1+e1*x))*e1 + e2,
and if we multiply both sides by 2^52 and do a transpose, we obtain
2^52*y - 2^52*x = -(x^2/(1+e1*x))*(2^52*e1) + (2^52*e2)
Because the fractional parts of x and y have only 52 bits, the left side
of this last equation must be an integer, and therefore so must the right
side. However, abs(2^52*e1) <= 1/4 and abs(2^52*e2) <= 1/2. This tells
us that the only way we will ever get a non zero integer value for the
right side is to have the factor x^2/(1+e1*x) be greater than or equal to
2. That in turn means that x must be at least the square root of 2 (or
within round off error of it.)
That fact being established, the only modification needed on your
while-loop method to speed things up is to start with x = sqrt(2), or to
play safe, a small amount below that. I started with x = 1.414213562,
and let it run until it finally stopped some ten or twenty minutes later
(I forgot to time it) and arrived at the number above.