P is for Practical PerlMonks

### Re: Re: Re: Help w/ Code Optimization

by runrig (Abbot)
 on Dec 26, 2001 at 23:21 UTC ( #134468=note: print w/replies, xml ) Need Help??

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

Create A New User
Domain Nodelet?
Node Status?
node history
Node Type: note [id://134468]
help
Chatterbox?
and the web crawler heard nothing...

How do I use this? | Other CB clients
Other Users?
Others taking refuge in the Monastery: (2)
As of 2023-03-29 00:54 GMT
Sections?
Information?
Find Nodes?
Leftovers?
Voting Booth?
Which type of climate do you prefer to live in?

Results (70 votes). Check out past polls.

Notices?