Fun with Unbounded Rational Numbers

Fun with Unbounded Rational Numbers

[Jeffrey R. Carter]

I decided to see what happens when you calculate square roots using Newton's 
method and the unbounded rational numbers in PragmARC.Rational_Numbers. Some 
interesting things happen.

For example, trying to calculate sqrt(2) to a small enough accuracy (0.000001 
will do it) results in the calculation of the next estimate

          M := Two * X;
          X := X - Y / M;

taking a Long Time (it does conclude eventually).

More alarming, though, is that for some values, even a fairly large accuracy 
will quickly result in Storage_Error. For example, sqrt(1.28) to an accuracy of 
0.01 will do it.

I realize that square roots are often irrational, but the accuracies seem so 
large that I didn't anticipate consuming GB of memory. Just naive, I guess.

-- 
Jeff Carter
"C++ is vast and dangerous, a sort of Mordor of
programming languages."
Jason R. Fruit
120

Re: Fun with Unbounded Rational Numbers

[Dmitry A. Kazakov]

On 2017-04-08 12:37, Jeffrey R. Carter wrote:
> I decided to see what happens when you calculate square roots using
> Newton's method and the unbounded rational numbers in
> PragmARC.Rational_Numbers. Some interesting things happen.
>
> For example, trying to calculate sqrt(2) to a small enough accuracy
> (0.000001 will do it) results in the calculation of the next estimate
>
>          M := Two * X;
>          X := X - Y / M;
>
> taking a Long Time (it does conclude eventually).
>
> More alarming, though, is that for some values, even a fairly large
> accuracy will quickly result in Storage_Error. For example, sqrt(1.28)
> to an accuracy of 0.01 will do it.
>
> I realize that square roots are often irrational, but the accuracies
> seem so large that I didn't anticipate consuming GB of memory. Just
> naive, I guess.

To use unbounded rational numbers should have the purpose of having 
exact computations <=> both absolutely accurate and absolutely precise. 
Otherwise, if accuracy is not absolute, you are allowed to choose from a 
set of rational numbers within the accuracy interval. Doing so you can 
take one with a less demanding representation, i.e. you are fiddling 
with the precision. And conversely for any algorithm involving 
irrational numbers you can always find an approximation requiring 
however much memory to waste.

In short unbounded rational numbers is *always* bad idea unless proven 
otherwise.

For transferring algorithms from real to unbounded rational numbers one 
could consider rounding the results within the accuracy constraint to 
some simpler (shorter) number. OK, this is what the floating-point 
arithmetic does already, right? (:-))

-- 
Regards,
Dmitry A. Kazakov
http://www.dmitry-kazakov.de

Re: Fun with Unbounded Rational Numbers

[Robert Eachus]

On Saturday, April 8, 2017 at 6:37:48 AM UTC-4, Jeffrey R. Carter wrote:
> I decided to see what happens when you calculate square roots using Newton's 
> method and the unbounded rational numbers in PragmARC.Rational_Numbers. Some 
> interesting things happen.
> 
> For example, trying to calculate sqrt(2) to a small enough accuracy (0.000001 
> will do it) results in the calculation of the next estimate
> 
>           M := Two * X;
>           X := X - Y / M;
> 
> taking a Long Time (it does conclude eventually).
> 
> More alarming, though, is that for some values, even a fairly large accuracy 
> will quickly result in Storage_Error. For example, sqrt(1.28) to an accuracy of 
> 0.01 will do it.
> 
> I realize that square roots are often irrational, but the accuracies seem so 
> large that I didn't anticipate consuming GB of memory. Just naive, I guess.

Sigh! You are running through a minefield.  Slowing to a walk will keep you alive a few minutes longer, but to survive you have to take each step very carefully.

First, you are trying to calculate an approximation to an irrational number.  The algorithm promises AT LEAST quadratic convergence which means that there will be at least twice as many bits/digits/words to store at each step.

But that is only half the story.  Your approximation will be a rational number. (Add, subtract, multiply, divide on rationals result in rationals.) The difference between your estimate and the actual irrational number is also an irrational number.  So when you test to a tolerance that subtraction should always raise Storage_Error.  But you did not say what you are using as an approximation of sqrt(2) for testing.  I suspect you are squaring your approximate value and testing that.  Sounds good, but...

For some reason you are using decimal limits.  These can not be represented accurately in binary, but using a rational package they will be represented by two numbers, with the divisor being a power of ten.  How do you compare two numbers stored as the ratio of integers? A/B ? C/D evaluate A*D ? B*C both integers, so now compare.  The rational arithmetic package will do those two multiplications, then throw the results away.  Unfortunately, the fragments left on the heap will be too small for any of the results in the next iteration.  This is probably what is running you out of memory.

You could try increasing the heap size.  Set your paging file to some large number, and try again with heap size set approximately to your main memory size.  You do not want to actually use that pagefile, set it to 3x you main memory size just to avoid system crashes.  On Windows, allowing the system to grow the pagefile, or putting it on C: are both bad news.  A cheap thumb drive will do nicely, and you don't have to worry about backing it up.

A much better idea is to throw away the limit, and have your program print the approximation at every step.  I don't know what you are looking for here, but a dozen iterations should give you all the accuracy you might need. (Even if you are measuring a diagonal brace for a house the size of the universe.)

Re: Fun with Unbounded Rational Numbers

[antispam]

Jeffrey R. Carter <spam.jrcarter.not@spam.not.acm.org> wrote:
> I decided to see what happens when you calculate square roots using Newton's 
> method and the unbounded rational numbers in PragmARC.Rational_Numbers. Some 
> interesting things happen.
> 
> For example, trying to calculate sqrt(2) to a small enough accuracy (0.000001 
> will do it) results in the calculation of the next estimate
> 
>          M := Two * X;
>          X := X - Y / M;
> 
> taking a Long Time (it does conclude eventually).
> 
> More alarming, though, is that for some values, even a fairly large accuracy 
> will quickly result in Storage_Error. For example, sqrt(1.28) to an accuracy of 
> 0.01 will do it.
> 
> I realize that square roots are often irrational, but the accuracies seem so 
> large that I didn't anticipate consuming GB of memory. Just naive, I guess.

There is somthing wrong with your code (or less likely with
PragmARC.Rational_Numbers).  Using computer algebra system
FriCAS and the following naive code:

newton(y : Fraction(Integer), eps : Fraction(Integer)) : Fraction(Integer) ==
    x : Fraction(Integer) := 1
    while abs(y - x^2) >= eps repeat
        x := (x + y/x)/2
    x

I can compute newton(2, 1/10^10000) with almost no time.  More
precisely time for actual computation is of order of timing resolution,
that is 10ms, while time needed to print result (that is binary to
decimal convertion) is of order 20ms.  The result has 12544 characters.

FriCAS after each rational operation cancels common factors between
numerator and denominator, this limits growth of numbers and is crucial
for bigger calculations.  But even without canceling common factors
your computation should work.  Namely size of numbers after
each operation is approximately sum of sizes of arguments.  Since
y above does not change size of y/x is bounded by size of x plus
a constant.  So in iteration n + 1 size of x is bounded by two
times size of x in iteration n plus a constant.  Each iteration
of Newton iteration approximately doubles accuracy.  The total
effect is that size of final result is proportional to accuracy,
and low accuracy (relatively large eps) result should be small.

-- 
                              Waldek Hebisch

Re: Fun with Unbounded Rational Numbers

[Jeffrey R. Carter]

On 04/09/2017 02:30 AM, antispam@math.uni.wroc.pl wrote:
>
> newton(y : Fraction(Integer), eps : Fraction(Integer)) : Fraction(Integer) ==
>     x : Fraction(Integer) := 1
>     while abs(y - x^2) >= eps repeat
>         x := (x + y/x)/2
>     x

Maybe I'm missing something, but this doesn't look like Newton's method to me.

-- 
Jeff Carter
"This trial is a travesty. It's a travesty of a mockery of a
sham of a mockery of a travesty of two mockeries of a sham. ...
Do you realize there's not a single homosexual on that jury?"
Bananas
27

Re: Fun with Unbounded Rational Numbers

[antispam]

Jeffrey R. Carter <spam.jrcarter.not@spam.not.acm.org> wrote:
> On 04/09/2017 02:30 AM, antispam@math.uni.wroc.pl wrote:
> >
> > newton(y : Fraction(Integer), eps : Fraction(Integer)) : Fraction(Integer) ==
> >     x : Fraction(Integer) := 1
> >     while abs(y - x^2) >= eps repeat
> >         x := (x + y/x)/2
> >     x
> 
> Maybe I'm missing something, but this doesn't look like Newton's method to me.
> 

We are sloving equation f(x) - y = 0 where f(x) = x^2.  Newton says:

x_{n+1} = x_{n} - (f(x_n) - y)/f'(x_n)

we have

f'(x) = 2x

so

x_{n+1} = x_{n} - (x_n^2 - y)/(2*(x_n)) =
          x_{n} - x_n^2/(2*(x_n)) + y/(2*(x_n)) =
          x_{n} - x_n/2 + y/(2*(x_n)) =
          x_{n}/2 + y/(2*(x_n)) = (x_{n} + y/x_n)/2
       
So indeed given the same initial approximation the formula produces
exactly the same sequence of numbers as direct use of Newton formula,
but is slightly cheaper to compute.

-- 
                              Waldek Hebisch

Re: Fun with Unbounded Rational Numbers

[Jeffrey R. Carter]

On 04/09/2017 09:25 PM, antispam@math.uni.wroc.pl wrote:
> Jeffrey R. Carter <spam.jrcarter.not@spam.not.acm.org> wrote:
>> On 04/09/2017 02:30 AM, antispam@math.uni.wroc.pl wrote:
>>>
>>> newton(y : Fraction(Integer), eps : Fraction(Integer)) : Fraction(Integer) ==
>>>     x : Fraction(Integer) := 1
>>>     while abs(y - x^2) >= eps repeat
>>>         x := (x + y/x)/2
>>>     x
>>
>> Maybe I'm missing something, but this doesn't look like Newton's method to me.
>>
>
> We are sloving equation f(x) - y = 0 where f(x) = x^2.  Newton says:
>
> x_{n+1} = x_{n} - (f(x_n) - y)/f'(x_n)
>
> we have
>
> f'(x) = 2x
>
> so
>
> x_{n+1} = x_{n} - (x_n^2 - y)/(2*(x_n)) =
>           x_{n} - x_n^2/(2*(x_n)) + y/(2*(x_n)) =
>           x_{n} - x_n/2 + y/(2*(x_n)) =
>           x_{n}/2 + y/(2*(x_n)) = (x_{n} + y/x_n)/2
>
> So indeed given the same initial approximation the formula produces
> exactly the same sequence of numbers as direct use of Newton formula,
> but is slightly cheaper to compute.

As I learned it, Newton's method is to iteratively improve an estimate of the 
(positive) root of

Y = X**2 - R

At each step, the next estimate is the root of the line through (X, Y) with slope

M = Y' = 2 * X

which works out to

X - Y/M

Substituting in the definitions of Y and M, it's easy to show this is the same 
as your version. I guess I'd never seen it simplified further.

As to being cheaper to calculate, is that true? For rational math, division is 
multiplication and subtraction is addition, so both have 2 multiplications and 
an addition.

-- 
Jeff Carter
"I didn't squawk about the steak, dear. I
merely said I didn't see that old horse
that used to be tethered outside here."
Never Give a Sucker an Even Break
103

Re: Fun with Unbounded Rational Numbers

[antispam]

Jeffrey R. Carter <spam.jrcarter.not@spam.not.acm.org> wrote:
> On 04/09/2017 09:25 PM, antispam@math.uni.wroc.pl wrote:
> > Jeffrey R. Carter <spam.jrcarter.not@spam.not.acm.org> wrote:
> >> On 04/09/2017 02:30 AM, antispam@math.uni.wroc.pl wrote:
> >>>
> >>> newton(y : Fraction(Integer), eps : Fraction(Integer)) : Fraction(Integer) ==
> >>>     x : Fraction(Integer) := 1
> >>>     while abs(y - x^2) >= eps repeat
> >>>         x := (x + y/x)/2
> >>>     x
> >>
> >> Maybe I'm missing something, but this doesn't look like Newton's method to me.
> >>
> >
> > We are sloving equation f(x) - y = 0 where f(x) = x^2.  Newton says:
> >
> > x_{n+1} = x_{n} - (f(x_n) - y)/f'(x_n)
> >
> > we have
> >
> > f'(x) = 2x
> >
> > so
> >
> > x_{n+1} = x_{n} - (x_n^2 - y)/(2*(x_n)) =
> >           x_{n} - x_n^2/(2*(x_n)) + y/(2*(x_n)) =
> >           x_{n} - x_n/2 + y/(2*(x_n)) =
> >           x_{n}/2 + y/(2*(x_n)) = (x_{n} + y/x_n)/2
> >
> > So indeed given the same initial approximation the formula produces
> > exactly the same sequence of numbers as direct use of Newton formula,
> > but is slightly cheaper to compute.
> 
> As I learned it, Newton's method is to iteratively improve an estimate of the 
> (positive) root of
> 
> Y = X**2 - R
> 
> At each step, the next estimate is the root of the line through (X, Y) with slope
> 
> M = Y' = 2 * X
> 
> which works out to
> 
> X - Y/M
> 
> Substituting in the definitions of Y and M, it's easy to show this is the same 
> as your version. I guess I'd never seen it simplified further.
> 
> As to being cheaper to calculate, is that true? For rational math, division is 
> multiplication and subtraction is addition, so both have 2 multiplications and 
> an addition.
> 

For rational math key factor is size of numbers.  Y has size approximatly
twice size of X plus size of R.  M is about the same as X, so Y/M
is approximatly 3 times size of X plus size of R before canceling
common factors.  Without canceling common factors X - Y/M
would be approximatly 4 times size of X plus size of R, that
is twice as large (assuming small R) as the version I gave.  And
without canceling common factors the difference in size will
accumulate leading to much worse behaviour.   With canceling common
factors the most expensive part is computing GCD-s and your version
leads to larger GCD-s.  But the difference will be moderate.

-- 
                              Waldek Hebisch

Re: Fun with Unbounded Rational Numbers

[Paul Rubin]

"Jeffrey R. Carter" <spam.jrcarter.not@spam.not.acm.org> writes:
> For example, trying to calculate sqrt(2)...
> taking a Long Time (it does conclude eventually).

Something seems wrong with that calculation or implementation.  I did
a simple Haskell version and it did 10 iterations (resulting in numbers
of around 750 digits) in almost no time:

    {-# LANGUAGE BangPatterns #-}
    import Data.Ratio
    type RI = Ratio Integer

    msqrt :: RI -> RI
    msqrt x = go (3/2) 10 where
      go g 0 = g
      go !g fuel =
        let y = g*g - x in
        go (g - y/(2*g)) (fuel-1)

    main = print . msqrt $ 2

("fuel" is the number of iterations).  With 5 iterations it gets the
more manageable ratio:

  1572584048032918633353217 % 1111984844349868137938112

which is approximately: x = 1.414213562373095
where |x*x - 2| is about 4.44e-16.

did you remember to recompute Y in your loop?  i.e.

         M := Two * X;
         Y = X * X - 2     -- <---- this might have been missing?
         X := X - Y / M;

Re: Fun with Unbounded Rational Numbers

[Jeffrey R. Carter]

On 04/09/2017 09:15 AM, Paul Rubin wrote:
>
> did you remember to recompute Y in your loop?  i.e.
>
>          M := Two * X;
>          Y = X * X - 2     -- <---- this might have been missing?
>          X := X - Y / M;

The implementation is

       X := R / Two;

       All_Iterations : for I in 1 .. 15 loop
          Y := X ** 2 - R;

          exit All_Iterations when abs Y < A;

          M := Two * X;
          X := X - Y / M;
       end loop All_Iterations;

       return X;

where R is the argument and A the accuracy.

-- 
Jeff Carter
"This trial is a travesty. It's a travesty of a mockery of a
sham of a mockery of a travesty of two mockeries of a sham. ...
Do you realize there's not a single homosexual on that jury?"
Bananas
27

Re: Fun with Unbounded Rational Numbers

[Paul Rubin]

"Jeffrey R. Carter" <spam.jrcarter.not@spam.not.acm.org> writes:
> The implementation is
>
>       X := R / Two;
>       All_Iterations : for I in 1 .. 15 loop
>          Y := X ** 2 - R;
>          exit All_Iterations when abs Y < A;   ...


Can you instrument this loop to check how many iterations are actually
being executed, and preferably how long each one takes?

The same calculation in compiled Haskell takes about 0.012 cpu seconds
including printing the answer, which is about 50k chars (i.e. 25k digits
for each of the numerator and denominator).

Doing 25 iterations takes about 2 minutes and uses around 250MB of
memory, giving an answer of around 50M chars.

You should get within 1e-8 or whatever it was in less than 5 iterations.

Re: Fun with Unbounded Rational Numbers

[Jeffrey R. Carter]

On 04/09/2017 11:18 PM, Paul Rubin wrote:
> "Jeffrey R. Carter" <spam.jrcarter.not@spam.not.acm.org> writes:
>> The implementation is
>>
>>       X := R / Two;
>>       All_Iterations : for I in 1 .. 15 loop
>>          Y := X ** 2 - R;
>>          exit All_Iterations when abs Y < A;   ...
>
>
> Can you instrument this loop to check how many iterations are actually
> being executed, and preferably how long each one takes?

As I said in another post, the thing that was taking the time was subtracting R 
(2) in the calculation of Y. I made some improvements to the addition algorithm 
and now it calculates sqrt(2) to within 0.000001 very quickly, in 4-5 
iterations. The result is 665857/470832.

-- 
Jeff Carter
"I didn't squawk about the steak, dear. I
merely said I didn't see that old horse
that used to be tethered outside here."
Never Give a Sucker an Even Break
103

Re: Fun with Unbounded Rational Numbers

[Paul Rubin]

"Jeffrey R. Carter" <spam.jrcarter.not@spam.not.acm.org> writes:
> As I said in another post, the thing that was taking the time was
> subtracting R (2) in the calculation of Y.

I'm trying to understand why that caused a bottleneck in the Ada
implementation when it didn't with other implementations.

> I made some improvements to the addition algorithm 

Oh, ok, you mean the rational number functions were something you were
working on rather than a released library.  That makes more sense--there
was some kind of algorithmic bug in it that's taken care of now.

I'd be interested to know if you get reasonable answers after 10 or 15
iterations.  Fwiw, for depth 7, I get:

  4946041176255201878775086487573351061418968498177    /
  3497379255757941172020851852070562919437964212608         

Higher depths get bigger numbers very quickly but I can supply the
values if you want to check them against your own implementation.

Re: Fun with Unbounded Rational Numbers

[Jeffrey R. Carter]

On 04/09/2017 09:15 AM, Paul Rubin wrote:
>          M := Two * X;
>          Y = X * X - 2     -- <---- this might have been missing?
>          X := X - Y / M;

Further investigation shows that the slow operation is the subtraction of R=2 in

Y := X ** 2 - R;

for X=665857/470832 (X**2=443365544449/221682772224≈2.00000000000451095044494277).

-- 
Jeff Carter
"This trial is a travesty. It's a travesty of a mockery of a
sham of a mockery of a travesty of two mockeries of a sham. ...
Do you realize there's not a single homosexual on that jury?"
Bananas
27