Square Root: Digit by Digit Algorithm

[ttw]

One can get (with a pretty large speed penalty) nice accuracy using continued fractions to represent integers. For symbolic computation, the symbols can be kept and numerically evaluated as needed. For example, Sqrt(2) is the continued fraction (1; 2,2,2,2,2,2...) and approximations have terms like 17/12; the accuracy is quadratic in the denominator at worst.

Some 50 or so years ago, David Matula used such a system; actually, 2 systems: one using a computer word each for the numerator and denominator; the other having a "floating slash" type of notation (needs extra space for the slash location.) I thought that on a system like the HP50g, three items per number might be good: the integer part, the numerator, and the denominator of the fractions. So far, I haven't tried using integers.

There's a hidden point that's not obvious (at least it wasn't to me). Any continued fraction P2/Q2 lies between the fractions P1/Q1 (the previous approximation) and (P1+P0)/(Q1+Q0), the mediant of the previous two approximations. It's sort of like doing interval arithmetic but with mediants rather than half-interval representatives.

Overflow in this type of system is handled by reducing a fraction Pk/Qk where one of these is too large, to the nearest P/Q which fits the implementation. It's expensive the one still has lots of precision. Some experiments with ill-conditioned matrices showed that even though the intermediate computations were approximate, the final answers were very good. One could get the inverse of high-order Hilbert matrices (the inverse is all integers) even though the intermediates don't fit with the system. The combination of good precision (1/Q^2 or better for P/Q) means that integer results popped out automatically.

I haven't looked yet, but Gosper's algorithms for continued fractions may speed things up a bit.

I've managed to get the cases I cared about to work by just storing the partial quotients in a list (HP50g) and converting these back and forth to fractions. Any system with large integers (Pari, HP50g, etc) can do well.