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): 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. |
|||
« Next Oldest | Next Newest »
|
User(s) browsing this thread: 1 Guest(s)