Post Reply 
Square Root: Digit by Digit Algorithm
11-26-2022, 01:15 PM
Post: #12
RE: Square Root: Digit by Digit Algorithm
(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.
Find all posts by this user
Quote this message in a reply
Post Reply 


Messages In This Thread
RE: Square Root: Digit by Digit Algorithm - Thomas Klemm - 11-26-2022 01:15 PM



User(s) browsing this thread: 1 Guest(s)