Threaded Mode | Linear Mode

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]

 502   L01251:  11.1..111.                      sqt16:  a - c -> a[w]

 503   L01252:  1.1.1...11  -> L01250                   if no carry go to sqt15

 504   L01253:  1111..111.                              a + c -> a[w]

 505   L01254:  .1....111.                              shift left a[w]

 506   L01255:  .....111..                              p - 1 -> p

 507   L01256:  1..1.1..1.                      sqt17:  shift right c[wp]

 508   L01257:  ....1.11..                              if p # 0

 509   L01260:  1.1.1..111  -> L01251                        then go to sqt16

 510   L01261:  ..11.1.111  -> L01065                   go to tnm12

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

002 RCL- Y

003 x≥0?

004 BACK 003

005 RCL+ Y

006 SDL 002

007 x⇆ Y

008 IP

009 SDL 001

010 # ½

011 +

012 x⇆ Y

013 DSZ Z

014 BACK 012

015 x⇆ Y

016 RTN

017 LBL A

018 # ½

019 # 012

020 x⇆ Z

021 RCL× Y

022 BACK 020

023 END

Example

41 A

6.4031242374 ¹²

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

By now you may have noticed that I didn't pay attention to the exponent.
But that's not so interesting in my opinion.

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.

Thomas Klemm · 11-25-2022, 07:51 AM

(11-25-2022 06:33 AM)brouhaha Wrote: I'm off-topic here

No problem with that.

From a comment in the linked HP-45 ROM:

Quote:keyed in by Eric Smith on 9-Mar-1995

Thanks 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.

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.

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

a = sqrt(2)

b = log(a)

a, b

(mpf('1.41421356237'), mpf('0.3465735902791'))

Code:

mp.dps = 100

a, b

(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

a = lambda: sqrt(2)

b = lambda: log(a())

a(), b()

(mpf('1.41421356237'), mpf('0.3465735902791'))

Code:

mp.dps = 100

a(), b()

(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.

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...

When done, you can hit EVAL to crunch.

Using an algebraic expression is another way to achieve the same thing.

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.
Digits that have been calculated could be cached.
Maybe someone already did this?

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):

    """Represents a mathematical constant with dynamic precision.

    When printed or used in an arithmetic operation, a constant

    is converted to a regular mpf at the working precision. A

    regular mpf can also be obtained using the operation +x."""

    def __new__(cls, func, name, docname=''):

        a = object.__new__(cls)

        a.name = name

        a.func = func

        a.__doc__ = getattr(function_docs, docname, '')

        return a

    def __call__(self, prec=None, dps=None, rounding=None):

        prec2, rounding2 = self.context._prec_rounding

        if not prec: prec = prec2

        if not rounding: rounding = rounding2

        if dps: prec = dps_to_prec(dps)

        return self.context.make_mpf(self.func(prec, rounding))

    @property

    def _mpf_(self):

        prec, rounding = self.context._prec_rounding

        return self.func(prec, rounding)

    def __repr__(self):

        return "<%s: %s~>" % (self.name, self.context.nstr(self(dps=15)))

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.

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.

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

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.

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):

    return n if n <= 1 else fib(n-2) + fib(n-1)

The function fib is evaluated at the same values again and again.

The solution here is to use a cache decorator:

Code:

from functools import cache

@cache

def fib(n):

    return n if n <= 1 else fib(n-2) + fib(n-1)

Now the previous elements of the sequence are calculated at most once.

We can do the same: the numbers have an additional parameter dps (decimal places):

Code:

from mpmath import *

Code:

mp.dps = 10

a, b = lambda dps: mpf(0), lambda dps: mpf(1)

n = 0

The parameter dps is ignored. It only serves as a key to the cache.
But we assume that it is consistent with the global setting mp.dps.

The returned function f is cached:

Code:

from functools import cache

def s(p, q, dps, n):

    @cache

    def f(dps):

        u, v = p(dps), q(dps)

        print(f"{n=} {dps=} {u=} {v=}")

        return v * sqrt(2 * v / (u + v))

    return f

The additional parameter n could be omitted.
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)

n += 1

b(mp.dps)

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

a(mp.dps), b(mp.dps)

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

a(mp.dps), b(mp.dps)

(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 += 1

b(mp.dps)

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.

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

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.)

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.

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-11...#pid108358