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