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.122018454301963435591038946477905736722308507360552962445074448170103302686224355942322411 1.122018454301963435591038946788236577875789772716985136 1.1220184543019634355910389464779057367223085073605529649771222079485 1.122018454301963435591038946477905736722308507360552962445074448170103302686224355942322410693 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% --