f

#### Smallest real number which has no inverse in Matlab.

```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);
end
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?

```
 0
this (338)
11/15/2006 7:33:26 PM
comp.soft-sys.matlab 211266 articles. 19 followers. lunamoonmoon (257) is leader.

2 Replies
551 Views

Similar Articles

[PageSpeed] 27

```ddtl wrote:
> 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 ~=3D (1+t)
> is 2^-52. So, I wrote the following program:
>
> x =3D 1;
> y =3D x;
> while x =3D=3D y
>   x =3D x + 2^-52;
>   y =3D 1/(1/x);
> end
> 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?

Generally, you will always get a round off error when representing
floating point numbers on a computer. This especially the case in your
example. The inverse of x will almost allways suffer from this, and is
also depending on the standard that is used for floating point
calculations. You can't even be sure that 1/1 equals 1 exactly (even if
it does when I test it in Matlab). When comparing x with y, you should
see if the difference is below a precision limit of your choice. You
should also take a look at the eps function if you haven't done so

Regards,=20
J=F8ger

```
 0
jogerh2032 (17)
11/15/2006 7:58:11 PM
```In article <obqml21jgi6shm35bkg4mub9n93vbn3edv@4ax.com>, ddtl
<this@is.invalid> wrote:

> 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);
> end
> 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

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.

Roger Stafford
```
 0
11/16/2006 6:24:28 AM