Square Root: Digit by Digit Algorithm - Printable Version +- HP Forums ( https://www.hpmuseum.org/forum)+-- Forum: HP Calculators (and very old HP Computers) ( /forum-3.html)+--- Forum: General Forum ( /forum-4.html)+--- Thread: Square Root: Digit by Digit Algorithm ( /thread-19201.html) |

Square Root: Digit by Digit Algorithm - Thomas Klemm - 11-25-2022 06:02 AM
The following program for the WP-34S is inspired by the this snippet of the HP-45 ROM source code: Code: ` 501 L01250: .1111...1. sqt15: c + 1 -> c[p]` The registers are mapped to the stack as follows: Z: p Y: c X: a While the program could easily be translated to HP calculators, the INC and SDL commands make it shorter: Code: `001 INC Y` Example 41 A 6.4031242374 ^{12}It is recommended to use a number in the range [1, 100), although it works outside of it as well. However, it is more efficient to apply a decimal shift first. So use 47.11 instead of 4711 or 0.4711. Make sure to shift an even number of digits. The program can be helpful to understand the algorithm if you go through an example step by step. References - HP-45 ROM source code from United States Patent 4,001,569
- HP35 square root algorithm
- Personal Calculator Algorithms I: Square Roots
- Patent US3576983
- Calculator Patents - Hewlett-Packard
By now you may have noticed that I didn't pay attention to the exponent. But that's not so interesting in my opinion. RE: Square Root: Digit by Digit Algorithm - brouhaha - 11-25-2022 06:33 AM
I'm off-topic here, but you cited William Egbert's HP Journal article "Personal Calculator Algorithms I: Square Roots", and that reminds me of a very optimistic statement Egbert made in the second of the four articles, "Personal Calculator Algorithms II: Trigonometric Functions", which I quote: "The problem with this scaling process is that in current computers numbers can be expressed only to a limited number of digits," I'm still looking forward to the future computers that can express numbers to an unlimited number of digits. RE: Square Root: Digit by Digit Algorithm - Thomas Klemm - 11-25-2022 07:51 AM
(11-25-2022 06:33 AM)brouhaha Wrote: I'm off-topic hereNo problem with that. From a comment in the linked HP-45 ROM: Quote:keyed in by Eric Smith on 9-Mar-1995Thanks a lot for this. (11-25-2022 06:33 AM)brouhaha Wrote: I'm still looking forward to the future computers that can express numbers to an unlimited number of digits. What exactly do you have in mind apart from unlimited integers in Python and real and complex floating-point arithmetic with arbitrary precision in mpmath. The calculations could also be carried out lazily and only at the end the required precision is specified. Similar to WolframAlpha where more digits can be requested. Or the precision could be changed during the calculation. RE: Square Root: Digit by Digit Algorithm - brouhaha - 11-25-2022 09:30 AM
(11-25-2022 07:51 AM)Thomas Klemm Wrote: What exactly do you have in mind apart from unlimited integers in Python and real and complex floating-point arithmetic with arbitrary precision in mpmath. Arbitrary precistion, but still limited, though there may not be a specific limit set by the language implementation. (I haven't actually looked at how the implementation defines the length, but if there's a fixed-length field that specifies the length of the variable-length number, then there is an implementation limit, though it might be bigger than the amount of DRAM in existence.) Egbert said "in current computers numbers can be expressed only to a limited number of digits", implying the possibility of computers in some other era (as opposed to "current") that can express numbers to an unlimited number of digits. I can't even begin to imagine how that could be done. RE: Square Root: Digit by Digit Algorithm - Thomas Klemm - 11-25-2022 10:59 AM
(11-25-2022 09:30 AM)brouhaha Wrote: if there's a fixed-length field that specifies the length of the variable-length number, then there is an implementation limit, though it might be bigger than the amount of DRAM in existence. As a playground I recommend to use an online Jupyter-Lab. There mpmath is preinstalled. Code: `from mpmath import *` The Naïve Way Code: `mp.dps = 10` (mpf('1.41421356237'), mpf('0.3465735902791')) Code: `mp.dps = 100` (mpf('1.414213562369695864617824554443359375'), mpf('0.3465735902791493572294712066650390625')) First I was surprised but then I noticed: nah, that aren't 100 digits. The Lazy Way Using functions as numbers we can postpone the evaluation. Code: `mp.dps = 10` (mpf('1.41421356237'), mpf('0.3465735902791')) Code: `mp.dps = 100` (mpf('1.414213562373095048801688724209698078569671875376948073176679737990732478462107038850387534327641572735'), mpf('0.3465735902799726547086160607290882840377500671801276270603400047466968109848473578029316634982093437698')) The usage is rather clumsy but could be improved by wrapping these "lazy numbers" in a class. Digits that have been calculated could be cached. Maybe someone already did this? (11-25-2022 09:30 AM)brouhaha Wrote: Egbert said "in current computers numbers can be expressed only to a limited number of digits", implying the possibility of computers in some other era (as opposed to "current") that can express numbers to an unlimited number of digits. I can't even begin to imagine how that could be done. You enter 2 and press the √ button. And then you die before reading all the digits. RE: Square Root: Digit by Digit Algorithm - Thomas Klemm - 11-25-2022 11:57 AM
(11-25-2022 10:59 AM)Thomas Klemm Wrote: Using functions as numbers we can postpone the evaluation. This reminded me of a recent thread: (11-17-2022 08:58 PM)halbitton Wrote: On the 50G, you can build an un-evaluated equation right on the stack using RPN keystrokes, with the following caveat... Using an algebraic expression is another way to achieve the same thing. RE: Square Root: Digit by Digit Algorithm - Albert Chan - 11-25-2022 02:28 PM
(11-25-2022 10:59 AM)Thomas Klemm Wrote: The usage is rather clumsy but could be improved by wrapping these "lazy numbers" in a class. Personally, I like the "clumsy" approach. It reminds user the delayed evaluation part. mpmath does use the "clean" approach with many constants, by wrapping it in a class. >>> import mpmath >>> mpmath.ln2 <ln(2): 0.693147~> >>> type(_) <class 'mpmath.ctx_mp_python.constant'> ctx_mp_python.py: Code: `class _constant(_mpf):` RE: Square Root: Digit by Digit Algorithm - ttw - 11-25-2022 02:40 PM
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. RE: Square Root: Digit by Digit Algorithm - Albert Chan - 11-25-2022 03:48 PM
(11-25-2022 02:40 PM)ttw Wrote: 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. If we let x0 = P0/Q0, x1 = P1/Q1, continued fraction coefficient, k ≥ 1: x2 = P2/Q2 = (k*P1+P0)/(k*Q1+Q0) = x1 - (x1-x0)* (Q0/(Q0+k*Q1)) Since Q1 > Q0, 0 < last term < 1/2, x2 between ((x0+x1)/2, x1) In other words, convergent is always better estimate than mean of previously 2 convergents. RE: Square Root: Digit by Digit Algorithm - EdS2 - 11-25-2022 05:31 PM
Infinite digits is asking a lot, but digits bounded only by available memory is a possibility. See perhaps this previous thread: Arbitrary precision - two or three approaches RE: Square Root: Digit by Digit Algorithm - Thomas Klemm - 11-25-2022 07:05 PM
(11-25-2022 02:40 PM)ttw Wrote: One can get (with a pretty large speed penalty) nice accuracy using continued fractions to represent integers. Not sure if this is rather "representing real numbers" or maybe "fractions"? Otherwise I don't understand the sentence and you could explain it in more detail. (11-25-2022 02:40 PM)ttw Wrote: Sqrt(2) is the continued fraction (1; 2,2,2,2,2,2...) and approximations have terms like 17/12 Can we digress any more? C.f. Convergents of a Continued Fraction The program can be used to get the approximations by entering the coefficients of a continued fraction one by one. (11-25-2022 02:40 PM)ttw Wrote: I haven't looked yet, but Gosper's algorithms for continued fractions may speed things up a bit. We should be able to apply arbitrary (or maybe just analytic?) functions. These can be approximated by polynomials. I think the basic arithmetic with sequence of digits is easy. Even the square root, as we have seen, can be taken digit by digit. We can also invert any polynomial using the same method. At least in a certain neighbourhood. Is this also possible with continued fractions? In the end, we might want a representation using digits rather than continued fractions, since their coefficients can be any integer, making interpretation difficult for us. RE: Square Root: Digit by Digit Algorithm - Thomas Klemm - 11-26-2022 01:15 PM
(11-25-2022 10:59 AM)Thomas Klemm Wrote: Digits that have been calculated could be cached. When I first tried to implement the Log-Arcsine Algorithm it became very slow after only a few iterations. The reason is the same as with this recursive implementation of the Fibonacci sequence: Code: `def fib(n):` The solution here is to use a cache decorator: Code: `from functools import cache` We can do the same: the numbers have an additional parameter dps (decimal places): Code: `from mpmath import *` Code: `mp.dps = 10` But we assume that it is consistent with the global setting mp.dps. The returned function f is cached: Code: `from functools import cache` It is only used in the print statement for debugging purposes. Now we can iterate the following cell: Code: `a, b = b, s(a, b, mp.dps, n)` Let's assume we want to do that until two consecutive values agree: n=0 dps=10 u=mpf('0.0') v=mpf('1.0') mpf('1.41421356237') n=1 dps=10 u=mpf('1.0') v=mpf('1.41421356237') mpf('1.530733729451') n=2 dps=10 u=mpf('1.41421356237') v=mpf('1.530733729451') mpf('1.56072257612') n=3 dps=10 u=mpf('1.530733729451') v=mpf('1.56072257612') mpf('1.568274245263') n=4 dps=10 u=mpf('1.56072257612') v=mpf('1.568274245263') mpf('1.570165578465') n=5 dps=10 u=mpf('1.568274245263') v=mpf('1.570165578465') mpf('1.570638625461') … n=15 dps=10 u=mpf('1.570796324319') v=mpf('1.570796326123') mpf('1.57079632656') n=16 dps=10 u=mpf('1.570796326123') v=mpf('1.57079632656') mpf('1.570796326676') n=17 dps=10 u=mpf('1.57079632656') v=mpf('1.570796326676') mpf('1.570796326705') n=18 dps=10 u=mpf('1.570796326676') v=mpf('1.570796326705') mpf('1.570796326705') n=19 dps=10 u=mpf('1.570796326705') v=mpf('1.570796326705') mpf('1.570796326705') We can increase the decimal places to 20 or more and calculate the numbers again: Code: `mp.dps = 20` n=0 dps=20 u=mpf('0.0') v=mpf('1.0') n=1 dps=20 u=mpf('1.0') v=mpf('1.4142135623730950488011') n=2 dps=20 u=mpf('1.4142135623730950488011') v=mpf('1.5307337294603590869117') n=3 dps=20 u=mpf('1.5307337294603590869117') v=mpf('1.5607225761290261427843') n=4 dps=20 u=mpf('1.5607225761290261427843') v=mpf('1.5682742452729696319058') n=5 dps=20 u=mpf('1.5682742452729696319058') v=mpf('1.570165578477376456156') n=6 dps=20 u=mpf('1.570165578477376456156') v=mpf('1.5706386254663864340288') n=7 dps=20 u=mpf('1.5706386254663864340288') v=mpf('1.5707569005721505381635') n=8 dps=20 u=mpf('1.5707569005721505381635') v=mpf('1.5707864701835456920665') n=9 dps=20 u=mpf('1.5707864701835456920665') v=mpf('1.5707938626385798503144') n=10 dps=20 u=mpf('1.5707938626385798503144') v=mpf('1.5707957107555999869997') n=11 dps=20 u=mpf('1.5707957107555999869997') v=mpf('1.5707961727850588711724') n=12 dps=20 u=mpf('1.5707961727850588711724') v=mpf('1.5707962882924363328452') n=13 dps=20 u=mpf('1.5707962882924363328452') v=mpf('1.5707963171692814945527') n=14 dps=20 u=mpf('1.5707963171692814945527') v=mpf('1.5707963243884928347474') n=15 dps=20 u=mpf('1.5707963243884928347474') v=mpf('1.5707963261932956729051') n=16 dps=20 u=mpf('1.5707963261932956729051') v=mpf('1.5707963266444963826394') n=17 dps=20 u=mpf('1.5707963266444963826394') v=mpf('1.5707963267572965600852') n=18 dps=20 u=mpf('1.5707963267572965600852') v=mpf('1.5707963267854966044475') n=19 dps=20 u=mpf('1.5707963267854966044475') v=mpf('1.5707963267925466155385') (mpf('1.5707963267925466155385'), mpf('1.5707963267943091183104')) These intermediate values are all cached. So, if we evaluate the values again, we don't see the debug statements: Code: `mp.dps = 20` (mpf('1.5707963267925466155385'), mpf('1.5707963267943091183104')) The calculations can be continued with higher precision: Code: `a, b = b, s(a, b, mp.dps, n)` n=20 dps=20 u=mpf('1.5707963267925466155385') v=mpf('1.5707963267943091183104') mpf('1.5707963267947497440042') Even with 1000 decimal places, the calculation is reasonably fast. Without caching, the kernel would just spin. RE: Square Root: Digit by Digit Algorithm - ttw - 11-27-2022 07:34 PM
I finally found the extension of continued fractions (like using Newton's method on square roots) to polynomials of arbitrary degree. It might be useful. https://hal.inria.fr/inria-00116990v4/document RE: Square Root: Digit by Digit Algorithm - BruceH - 11-27-2022 10:44 PM
(11-25-2022 09:30 AM)brouhaha Wrote: Arbitrary precistion, but still limited, though there may not be a specific limit set by the language implementation. (I haven't actually looked at how the implementation defines the length, but if there's a fixed-length field that specifies the length of the variable-length number, then there is an implementation limit, though it might be bigger than the amount of DRAM in existence.) Clearly everything will be limited by RAM, or the number of atoms in the Universe, but that aside it's been discussed before. https://www.hpmuseum.org/forum/thread-11865-post-108358.html#pid108358 |