Tags: bigger, butwhich, everybody, inverse, matlab, number, programming, real, represented, smallest, smallestnumber

Smallest real number which has no inverse in Matlab.

On Programmer » Matlab

5,635 words with 2 Comments; publish: Mon, 12 May 2008 10:03:00 GMT; (200109.38, « »)

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?

All Comments

Leave a comment...

  • 2 Comments
    • 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

      already.=20

      Regards,=20

      J=F8ger

      #1; Mon, 12 May 2008 10:04:00 GMT
    • In article <obqml21jgi6shm35bkg4mub9n93vbn3edv.matlab.todaysummary.com.4ax.com>, ddtl

      <this.matlab.todaysummary.com.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

      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.

      Roger Stafford

      #2; Mon, 12 May 2008 10:05:00 GMT