http://qs1969.pair.com?node_id=134468


in reply to Re: Re: Help w/ Code Optimization
in thread Help w/ Code Optimization

There's an algorithm here to compute natural logarithms to any precision. It works for me, you might have to combine it with Math::BigFloat to get the precision you want (Oh, you're already using that...)
:-)

However, I don't know if going this route is any better than your current algorithm for getting the nth root of a number...

Update: By request, here's what I did to make the algorithm in the link above work (in the Read More link at the bottom of this post). This is not a perfect implementation, I'm pretty sure I'm losing precision somewhere, and you could probably make better use of Math::BigFloat, but it seems to be on the right track. And I don't know how this compares with the other ways of finding the root of a number...

Update: Works alot better now, not losing so much precision anymore (in fact, it looks like results are accurate to the precision you specify now) :-)

Updated code to reflect new version of Math::BigFloat
Also, did some rough benchmarks computing the 60th root of 1000. "Rough," because sifukurt's algorithm specifies 'iterations' to get precision, and I specify decimal places. For 50 decimal places on mine, and for 1 of his iterations (which is lt 40 decimal places), his wins, for 2 iterations (which is slightly gt 50 decimal places), mine wins slightly, and for 3 iterations mine wins by alot (but then his is accurate to much gt 50 places). And I'm still using operator overloading, which is a speed hit (though I'm not sure how much, and so far we're both using it).
For 90 decimal places, his takes 3 iterations, and the time mine takes is right in the middle of the time his takes to do 2 or 3 of his iterations.

use strict; use warnings; use Math::BigFloat; # number my $n = 10; # root my $r = 2; # precision my $p = 40; my $root = root($n, $r, $p); print "$root\n"; print $root*$root,"\n"; sub root { my ($n, $r, $p) = @_; $p++; my $log = LN($n, $p)/$r; EXP($log, $p)->bfround(1-$p); } BEGIN { # Cache the maximum precision so far calculated my $max_p=0; my $max_econst; sub ECONST { my ($P) = @_; $P = 9 unless defined $P; if ($P <= $max_p) { return $max_econst->copy->bfround(-$P); } my $Eps = 0.5 * Math::BigFloat->new("1"."0"x$P); my $N = Math::BigFloat->new("1")->bfround(-$P); my $D = Math::BigFloat->new("1")->bfround(-$P); my $J = Math::BigFloat->new("1")->bfround(-$P); { $N = $J * $N + 1; $D = $J * $D; if ($D >= $Eps) { $max_p = $P; return $max_econst = $N / $D; } $J++; redo; } } } sub EXP { my ($X, $P) = @_; $X = ref($X) ? $X->copy : Math::BigFloat->new($X); $P = 9 unless defined $P; $X->bfround(-$P); my $Y = $X->copy->bfround(0); $Y->bfround(-$P); $Y += ( 0 cmp $X ) if abs($X - $Y) > 0.5; $X = $X - $Y; my $Sum = Math::BigFloat->new("1")->bfround(-$P); my $Term = Math::BigFloat->new("1")->bfround(-$P); my $J = Math::BigFloat->new("1")->bfround(-$P); { $Term *= $X / $J; my $NewSum = $Sum + $Term; last unless $NewSum cmp $Sum; $Sum = $NewSum; $J++; redo; } return $Sum unless $Y cmp 0; my $E = ECONST($P); my $E_Y = 1; $E_Y *= $E for 1..$Y; $E_Y *= $Sum; return $E_Y; } sub LN { my ($X, $P) = @_; $X = ref($X) ? $X->copy : Math::BigFloat->new($X); $P = 9 unless defined $P; $X->bfround(-$P); return -LN(1 / $X, $P) if $X < 1; my $M = 0; ++$M until (2 ** $M) > $X; $M--; my $Z = $X / (2 ** $M); my $Zeta = (1 - $Z) / (1 + $Z); my $N = $Zeta; my $Ln = $Zeta; my $Zetasup2 = $Zeta * $Zeta; my $J = 1; { $N = $N * $Zetasup2; my $NewLn = $Ln + $N / (2 * $J + 1); return $M * LN2P($P) - 2 * $Ln unless $NewLn cmp $Ln; $Ln = $NewLn; $J++; redo; } } BEGIN { # Cache the maximum precision so far calculated my $max_p = 0; my $max_ln2p; sub LN2P { my ($P) = @_; if ($P <= $max_p) { return $max_ln2p->copy->bfround(-$P); } my $one = Math::BigFloat->new("1")->bfround(-$P); my $two = Math::BigFloat->new("2")->bfround(-$P); my $N = $one / 3; my $Ln = $N; my $Zetasup2 = $one / 9; my $J = 1; { $N = $N * $Zetasup2; my $NewLn = $Ln + $N / (2 * $J + 1); unless ($NewLn cmp $Ln) { $max_p = $P; return $max_ln2p = $Ln * 2; } $Ln = $NewLn; $J++; redo; } } } # Benchmark code and results: use Benchmark qw(cmpthese); # number my $n = 1000; # root my $r = 60; # precision my $p = 90; # A slight cheat - by calling root() before # the benchmark code I am precalculating # and caching the constants e and ln(2) my $root = root($n, $r, $p); print "$root\n"; my $Root1 = Root($n, $r, 1); print "$Root1\n"; my $Root2 = Root($n, $r, 2); print "$Root2\n"; my $Root3 = Root($n, $r, 3); print "$Root3\n"; cmpthese(-30, { MINE=>sub { root($n, $r, $p) }, NEWTON1=>sub { Root($n, $r, 1) }, NEWTON2=>sub { Root($n, $r, 2) }, NEWTON3=>sub { Root($n, $r, 3) }, }); >./tst 1.12201845430196343559103894647790573672230850736055296244507444817010 +3302686224355942322411 1.122018454301963435591038946788236577875789772716985136 1.1220184543019634355910389464779057367223085073605529649771222079485 1.12201845430196343559103894647790573672230850736055296244507444817010 +3302686224355942322410693 Benchmark: running MINE, NEWTON1, NEWTON2, NEWTON3, each for at least +120 CPU se conds... MINE: 128 wallclock secs (127.63 usr + 0.10 sys = 127.73 CPU) @ + 0.08/s ( n=10) NEWTON1: 125 wallclock secs (125.67 usr + 0.02 sys = 125.69 CPU) @ + 0.58/s ( n=73) NEWTON2: 120 wallclock secs (120.07 usr + 0.00 sys = 120.07 CPU) @ + 0.11/s ( n=13) NEWTON3: 128 wallclock secs (128.07 usr + 0.00 sys = 128.07 CPU) @ + 0.05/s ( n=6) s/iter NEWTON3 MINE NEWTON2 NEWTON1 NEWTON3 21.3 -- -40% -57% -92% MINE 12.8 67% -- -28% -87% NEWTON2 9.24 131% 38% -- -81% NEWTON1 1.72 1140% 642% 436% --