[libmath-prime-util-perl] 03/23: Increased precision for bignum zeta and R
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:45:54 UTC 2015
This is an automated email from the git hooks/post-receive script.
ppm-guest pushed a commit to annotated tag v0.14
in repository libmath-prime-util-perl.
commit d574abdb6851c2d2c2fbca722e81ded0d264a234
Author: Dana Jacobsen <dana at acm.org>
Date: Wed Nov 21 00:42:52 2012 -0800
Increased precision for bignum zeta and R
---
TODO | 5 +
lib/Math/Prime/Util.pm | 59 ++++---
lib/Math/Prime/Util/PP.pm | 385 ++++++++++++++++++++++++++++++++++++++--------
3 files changed, 357 insertions(+), 92 deletions(-)
diff --git a/TODO b/TODO
index b5a4c35..e60e8c2 100644
--- a/TODO
+++ b/TODO
@@ -34,3 +34,8 @@
- Add Lehmer in PP (I have a version, just needs some polishing).
- Move AKS tests into their own test, and add some provable prime tests.
+
+- For bignums, RiemannZeta and RiemmannR are slow and give questionable
+ precision. We should be able to do better. One problem is the accuracy
+ bug in Math::BigFloat. Perhaps check for Math::MPFR installed and use it?
+ Alternately write real-number pow function using GMP.
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index b4b4ebb..7e570a9 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -207,7 +207,7 @@ sub _validate_positive_integer {
sub _upgrade_to_float {
my($n) = @_;
return $n unless defined $Math::BigInt::VERSION || defined $Math::BigFloat::VERSION;
- do { require Math::BigFloat; Math::BigFloat->import(try=>'GMP,Pari') }
+ do { require Math::BigFloat; Math::BigFloat->import() }
if defined $Math::BigInt::VERSION && !defined $Math::BigFloat::VERSION;
return Math::BigFloat->new($n);
}
@@ -578,8 +578,8 @@ sub primes {
return random_nbit_prime($k) if $k <= $p0;
eval {
- require Math::BigInt; Math::BigInt->import( try=>'GMP,Pari' );
- require Math::BigFloat; Math::BigFloat->import( try=>'GMP,Pari' );
+ require Math::BigInt; Math::BigInt->import( try=>'GMP,Pari' );
+ require Math::BigFloat; Math::BigFloat->import();
1;
} or do {
croak "Cannot load Math::BigInt and Math::BigFloat";
@@ -659,8 +659,6 @@ sub primes {
return $n;
}
- no Math::BigFloat;
- no Math::BigInt;
}
}
@@ -1079,6 +1077,12 @@ sub prime_count_approx {
return $_prime_count_small[$x] if $x <= $#_prime_count_small;
+ # Below 2^58th or so, all differences between the high precision and C double
+ # precision result are less than 0.5.
+ if ($x <= $_XS_MAXVAL && $x <= 144115188075855872) {
+ return int(_XS_RiemannR($x) + 0.5);
+ }
+
# Turn on high precision FP if they gave us a big number.
$x = _upgrade_to_float($x) if ref($x) eq 'Math::BigInt';
@@ -1097,7 +1101,8 @@ sub prime_count_approx {
# my $result = int(LogarithmicIntegral($x) - LogarithmicIntegral(sqrt($x))/2);
- my $result = RiemannR($x) + 0.5;
+ my $tol = 10**-(length(int($x))+1);
+ my $result = RiemannR($x, $tol) + 0.5;
return Math::BigInt->new($result->bfloor->bstr()) if ref($result) eq 'Math::BigFloat';
return int($result);
@@ -1341,24 +1346,18 @@ sub RiemannZeta {
my($n) = @_;
croak("Invalid input to ReimannZeta: x must be > 0") if $n <= 0;
- #return Math::Prime::Util::PP::RiemannZeta($n) if defined $bignum::VERSION || ref($n) eq 'Math::BigFloat';
- return Math::Prime::Util::PP::RiemannZeta($n) if !$_Config{'xs'};
- return _XS_RiemannZeta($n);
+ return Math::Prime::Util::PP::RiemannZeta($n) if defined $bignum::VERSION || ref($n) eq 'Math::BigFloat';
+ return _XS_RiemannZeta($n) if $n <= $_XS_MAXVAL;
+ return Math::Prime::Util::PP::RiemannZeta($n);
}
sub RiemannR {
- my($n) = @_;
+ my($n, $tol) = @_;
croak("Invalid input to ReimannR: x must be > 0") if $n <= 0;
- return Math::Prime::Util::PP::RiemannR($n) if defined $bignum::VERSION || ref($n) eq 'Math::BigFloat';
- return Math::Prime::Util::PP::RiemannR($n) if !$_Config{'xs'};
- return _XS_RiemannR($n);
-
- # We could make a new object, like:
- # require Math::BigFloat;
- # my $bign = new Math::BigFloat "$n";
- # my $result = Math::Prime::Util::PP::RiemannR($bign);
- # return $result;
+ return Math::Prime::Util::PP::RiemannR($n, $tol) if defined $bignum::VERSION || ref($n) eq 'Math::BigFloat';
+ return _XS_RiemannR($n) if $n <= $_XS_MAXVAL;
+ return Math::Prime::Util::PP::RiemannR($n, $tol);
}
sub ExponentialIntegral {
@@ -2282,18 +2281,17 @@ Accuracy should be at least 14 digits.
my $z = RiemannZeta($s);
-Given a floating point input C<s> where C<s E<gt>= 0.5>, returns the floating
+Given a floating point input C<s> where C<s E<gt> 0>, returns the floating
point value of ζ(s)-1, where ζ(s) is the Riemann zeta function. One is
subtracted to ensure maximum precision for large values of C<s>. The zeta
-function is the sum from k=1 to infinity of C<1 / k^s>.
-
-Since the argument and the result are real, this is technically the Euler Zeta
-function.
+function is the sum from k=1 to infinity of C<1 / k^s>. This function only
+uses real arguments, so is basically the Euler Zeta function.
-Accuracy should be at least 14 digits, but currently does not increase
-accuracy with big floats. Small integer values are returned from a table,
-values between 0.5 and 5 use rational Chebyshev approximation, and larger
-values use a series.
+Accuracy should be at least 14 digits with native numbers and 35 digits with
+bignum or a BigInt/BigFloat argument. Small integer values are returned from
+a table. For native-precision numbers, a rational Chebyshev approximation is
+used between 0.5 and 5, while larger values use a series. Multiple precision
+numbers use Borwein (1991) algorithm 2 or the basic series.
=head2 RiemannR
@@ -2304,9 +2302,8 @@ Given a positive non-zero floating point input, returns the floating
point value of Riemann's R function. Riemann's R function gives a very close
approximation to the prime counting function.
-Accuracy should be at least 14 digits. The current implementation isn't
-correctly storing constants as big floats, so is not giving increased accuracy
-with big numbers like it should.
+Accuracy should be at least 14 digits for native numbers and 35 digits with
+bignum or a BigInt/BigFloat argument.
=head1 EXAMPLES
diff --git a/lib/Math/Prime/Util/PP.pm b/lib/Math/Prime/Util/PP.pm
index a1e5406..e3d154e 100644
--- a/lib/Math/Prime/Util/PP.pm
+++ b/lib/Math/Prime/Util/PP.pm
@@ -1288,6 +1288,9 @@ sub holf_factor {
@factors;
}
+
+
+
my $_const_euler = 0.57721566490153286060651209008240243104215933593992;
my $_const_li2 = 1.045163780117492784844588889194613136522615578151;
@@ -1378,96 +1381,356 @@ sub LogarithmicIntegral {
}
# Riemann Zeta function for integers, used for computing Riemann R
+# So many terms and digits are used so we can quickly do bignum R.
my @_Riemann_Zeta_Table = (
- 0.6449340668482264364724151666460251892, # zeta(2)
- 0.2020569031595942853997381615114499908,
- 0.0823232337111381915160036965411679028,
- 0.0369277551433699263313654864570341681,
- 0.0173430619844491397145179297909205279,
- 0.0083492773819228268397975498497967596,
- 0.0040773561979443393786852385086524653,
- 0.0020083928260822144178527692324120605,
- 0.0009945751278180853371459589003190170,
- 0.0004941886041194645587022825264699365,
- 0.0002460865533080482986379980477396710,
- 0.0001227133475784891467518365263573957,
- 0.0000612481350587048292585451051353337,
- 0.0000305882363070204935517285106450626,
- 0.0000152822594086518717325714876367220,
- 0.0000076371976378997622736002935630292,
- 0.0000038172932649998398564616446219397,
- 0.0000019082127165539389256569577951013,
- 0.0000009539620338727961131520386834493,
- 0.0000004769329867878064631167196043730,
- 0.0000002384505027277329900036481867530,
- 0.0000001192199259653110730677887188823,
- 0.0000000596081890512594796124402079358,
- 0.0000000298035035146522801860637050694,
- 0.0000000149015548283650412346585066307,
- 0.0000000074507117898354294919810041706,
- 0.0000000037253340247884570548192040184,
- 0.0000000018626597235130490064039099454,
- 0.0000000009313274324196681828717647350,
- 0.0000000004656629065033784072989233251,
- 0.0000000002328311833676505492001455976,
- 0.0000000001164155017270051977592973835,
- 0.0000000000582077208790270088924368599,
- 0.0000000000291038504449709968692942523,
- 0.0000000000145519218910419842359296322,
- 0.0000000000072759598350574810145208690,
- 0.0000000000036379795473786511902372363,
- 0.0000000000018189896503070659475848321,
- 0.0000000000009094947840263889282533118,
+ '0.64493406684822643647241516664602518921894990', # zeta(2) - 1
+ '0.20205690315959428539973816151144999076498629',
+ '0.082323233711138191516003696541167902774750952',
+ '0.036927755143369926331365486457034168057080920',
+ '0.017343061984449139714517929790920527901817490',
+ '0.0083492773819228268397975498497967595998635606',
+ '0.0040773561979443393786852385086524652589607906',
+ '0.0020083928260822144178527692324120604856058514',
+ '0.00099457512781808533714595890031901700601953156',
+ '0.00049418860411946455870228252646993646860643576',
+ '0.00024608655330804829863799804773967096041608846',
+ '0.00012271334757848914675183652635739571427510590',
+ '0.000061248135058704829258545105135333747481696169',
+ '0.000030588236307020493551728510645062587627948707',
+ '0.000015282259408651871732571487636722023237388990',
+ '0.0000076371976378997622736002935630292130882490903',
+ '0.0000038172932649998398564616446219397304546972190',
+ '0.0000019082127165539389256569577951013532585711448',
+ '0.00000095396203387279611315203868344934594379418741',
+ '0.00000047693298678780646311671960437304596644669478',
+ '0.00000023845050272773299000364818675299493504182178',
+ '0.00000011921992596531107306778871888232638725499778',
+ '0.000000059608189051259479612440207935801227503918837',
+ '0.000000029803503514652280186063705069366011844730920',
+ '0.000000014901554828365041234658506630698628864788168',
+ '0.0000000074507117898354294919810041706041194547190319',
+ '0.0000000037253340247884570548192040184024232328930593',
+ '0.0000000018626597235130490064039099454169480616653305',
+ '0.00000000093132743241966818287176473502121981356795514',
+ '0.00000000046566290650337840729892332512200710626918534',
+ '0.00000000023283118336765054920014559759404950248298228',
+ '0.00000000011641550172700519775929738354563095165224717',
+ '0.000000000058207720879027008892436859891063054173122605',
+ '0.000000000029103850444970996869294252278840464106981987',
+ '0.000000000014551921891041984235929632245318420983808894',
+ '0.0000000000072759598350574810145208690123380592648509256',
+ '0.0000000000036379795473786511902372363558732735126460284',
+ '0.0000000000018189896503070659475848321007300850305893096',
+ '0.00000000000090949478402638892825331183869490875386000099',
+ '0.00000000000045474737830421540267991120294885703390452991',
+ '0.00000000000022737368458246525152268215779786912138298220',
+ '0.00000000000011368684076802278493491048380259064374359028',
+ '0.000000000000056843419876275856092771829675240685530571589',
+ '0.000000000000028421709768893018554550737049426620743688265',
+ '0.000000000000014210854828031606769834307141739537678698606',
+ '0.0000000000000071054273952108527128773544799568000227420436',
+ '0.0000000000000035527136913371136732984695340593429921456555',
+ '0.0000000000000017763568435791203274733490144002795701555086',
+ '0.00000000000000088817842109308159030960913863913863256088715',
+ '0.00000000000000044408921031438133641977709402681213364596031',
+ '0.00000000000000022204460507980419839993200942046539642366543',
+ '0.00000000000000011102230251410661337205445699213827024832229',
+ '0.000000000000000055511151248454812437237365905094302816723551',
+ '0.000000000000000027755575621361241725816324538540697689848904',
+ '0.000000000000000013877787809725232762839094906500221907718625',
+ '0.0000000000000000069388939045441536974460853262498092748358742',
+ '0.0000000000000000034694469521659226247442714961093346219504706',
+ '0.0000000000000000017347234760475765720489729699375959074780545',
+ '0.00000000000000000086736173801199337283420550673429514879071415',
+ '0.00000000000000000043368086900206504874970235659062413612547801',
+ '0.00000000000000000021684043449972197850139101683209845761574010',
+ '0.00000000000000000010842021724942414063012711165461382589364744',
+ '0.000000000000000000054210108624566454109187004043886337150634224',
+ '0.000000000000000000027105054312234688319546213119497764318887282',
+ '0.000000000000000000013552527156101164581485233996826928328981877',
+ '0.0000000000000000000067762635780451890979952987415566862059812586',
+ '0.0000000000000000000033881317890207968180857031004508368340311585',
+ '0.0000000000000000000016940658945097991654064927471248619403036418',
+ '0.00000000000000000000084703294725469983482469926091821675222838642',
+ '0.00000000000000000000042351647362728333478622704833579344088109717',
+ '0.00000000000000000000021175823681361947318442094398180025869417612',
+ '0.00000000000000000000010587911840680233852265001539238398470699902',
+ '0.000000000000000000000052939559203398703238139123029185055866375629',
+ '0.000000000000000000000026469779601698529611341166842038715592556134',
+ '0.000000000000000000000013234889800848990803094510250944989684323826',
+ '0.0000000000000000000000066174449004244040673552453323082200147137975',
+ '0.0000000000000000000000033087224502121715889469563843144048092764894',
+ '0.0000000000000000000000016543612251060756462299236771810488297723589',
+ '0.00000000000000000000000082718061255303444036711056167440724040096811',
+ '0.00000000000000000000000041359030627651609260093824555081412852575873',
+ '0.00000000000000000000000020679515313825767043959679193468950443365312',
+ '0.00000000000000000000000010339757656912870993284095591745860911079606',
+ '0.000000000000000000000000051698788284564313204101332166355512893608164',
+ '0.000000000000000000000000025849394142282142681277617708450222269121159',
+ '0.000000000000000000000000012924697071141066700381126118331865309299779',
+ '0.0000000000000000000000000064623485355705318034380021611221670660356864',
+ '0.0000000000000000000000000032311742677852653861348141180266574173608296',
+ '0.0000000000000000000000000016155871338926325212060114057052272720509148',
+ '0.00000000000000000000000000080779356694631620331587381863408997398684847',
+ '0.00000000000000000000000000040389678347315808256222628129858130379479700',
+ '0.00000000000000000000000000020194839173657903491587626465673047518903728',
+ '0.00000000000000000000000000010097419586828951533619250700091044144538432',
+ '0.000000000000000000000000000050487097934144756960847711725486604360898735',
+ '0.000000000000000000000000000025243548967072378244674341937966175648398693',
+ '0.000000000000000000000000000012621774483536189043753999660777148710632765',
+ '0.0000000000000000000000000000063108872417680944956826093943332037500694712',
+ '0.0000000000000000000000000000031554436208840472391098412184847972814371270',
+ '0.0000000000000000000000000000015777218104420236166444327830159601782237092',
+);
+
+# Select n = 55, good for 46ish digits of accuracy.
+my $_Borwein_n = 55;
+my @_Borwein_dk = (
+ '1',
+ '6051',
+ '6104451',
+ '2462539971',
+ '531648934851',
+ '71301509476803',
+ '6504925195108803',
+ '429144511928164803',
+ '21392068013887742403',
+ '832780518854440804803',
+ '25977281563850106233283',
+ '662753606729324750201283',
+ '14062742362385399866745283',
+ '251634235316509414702211523',
+ '3841603462178827861104812483',
+ '50535961819850087101900022211',
+ '577730330374203014014104003011',
+ '5782012706584553297863989289411',
+ '50984922488525881477588707205571',
+ '398333597655022403279683908035011',
+ '2770992240330783259897072664469955',
+ '17238422988353715312442126057365955',
+ '96274027751337344115352100618133955',
+ '484350301573059857715727453968687555',
+ '2201794236784087151947175826243477955',
+ '9068765987529892610841571032285864387',
+ '33926582279822401059328069515697217987',
+ '115535262182820447663793177744255246787',
+ '358877507711760077538925500462137369027',
+ '1018683886695854101193095537014797787587',
+ '2646951832121008166346437186541363159491',
+ '6306464665572570713623910486640730071491',
+ '13799752848354341643763498672558481367491',
+ '27780237373991939435100856211039992177091',
+ '51543378762608611361377523633779417047491',
+ '88324588911945720951614452340280439890371',
+ '140129110249040241501243929391690331218371',
+ '206452706984942815385219764876242498642371',
+ '283527707823296964404071683165658912154051',
+ '364683602811933600833512164561308162744771',
+ '441935796522635816776473230396154031661507',
+ '508231717051242054487234759342047053767107',
+ '559351463001010719709990637083458540691907',
+ '594624787018881191308291683229515933311427',
+ '616297424973434835299724300924272199623107',
+ '628083443816135918099559567176252011864515',
+ '633714604276098212796088600263676671320515',
+ '636056734158553360761837806887547188568515',
+ '636894970116484676875895417679248215794115',
+ '637149280289288581322870186196318041432515',
+ '637213397278310656625865036925470191411651',
+ '637226467136294189739463288384528579584451',
+ '637228536449134002301138291602841035366851',
+ '637228775173095037281299181461988671775171',
+ '637228793021615488494769154535569803469251',
+ '637228793670652595811622608101881844621763',
);
+# "An Efficient Algorithm for the Riemann Zeta Function", Borwein, 1991.
+# About 1.3n terms are needed for n digits of accuracy.
+sub _Recompute_Dk {
+ my $nterms = shift;
+ $_Borwein_n = $nterms;
+ @_Borwein_dk = ();
+ foreach my $k (0 .. $nterms) {
+ my $dsum = Math::BigFloat->bzero;
+ $dsum->accuracy(2*$_Borwein_n);
+ my $n = Math::BigInt->new($nterms-1)->bfac;
+ my $d = Math::BigInt->new($nterms)->bfac;
+ foreach my $i (0 .. $k) {
+ my $term = Math::BigFloat->bone;
+ $term->accuracy(2*$_Borwein_n);
+ $term->bmul($n)->bdiv($d);
+ $dsum += $term;
+ $n->bmul($nterms+$i)->bmul(4);
+ $d->bdiv($nterms-$i)->bmul(2*$i+1)->bmul(2*$i+2);
+ }
+ my $dk = ($nterms * $dsum + 1e-20)->as_int;
+ $_Borwein_dk[$k] = $dk;
+ #print " '$dk',\n";
+ }
+}
-# Compute Riemann Zeta function minus 1.
sub RiemannZeta {
my($x, $tol) = @_;
- $tol = 1e-16 unless defined $tol;
- my $sum = 0.0;
- my($y, $t);
- my $c = 0.0;
- for my $k (2 .. 1000000) {
- my $term = (2*$k+1) ** -$x;
+ if (!defined $bignum::VERSION && ref($x) !~ /^Math::Big/) {
+ return 0.0 + $_Riemann_Zeta_Table[int($x)-2]
+ if $x == int($x) && defined $_Riemann_Zeta_Table[int($x)-2];
+ $tol = 1e-16 unless defined $tol;
+ my($y, $t);
+ my $sum = 0.0;
+ my $c = 0.0;
+
+ for my $k (2 .. 1000000) {
+ my $term = (2*$k+1) ** -$x;
+ $y = $term-$c; $t = $sum+$y; $c = ($t-$sum)-$y; $sum = $t;
+ last if $term < abs($tol*$sum);
+ }
+ my $term = 3 ** -$x;
$y = $term-$c; $t = $sum+$y; $c = ($t-$sum)-$y; $sum = $t;
- last if $term < abs($tol*$sum);
+ $t = 1.0 / (1.0 - (2 ** -$x));
+ $sum *= $t;
+ $sum += ($t - 1.0);
+ return $sum;
}
- my $term = 3 ** -$x;
- $y = $term-$c; $t = $sum+$y; $c = ($t-$sum)-$y; $sum = $t;
- $t = 1.0 / (1.0 - (2 ** -$x));
- $sum *= $t;
- $sum += ($t - 1.0);
- $sum;
+ return Math::BigFloat->new($_Riemann_Zeta_Table[int($x)-2])
+ if $x == int($x) && defined $_Riemann_Zeta_Table[int($x)-2];
+
+ $tol = 1e-40 unless defined $tol;
+
+ # Trying to work around Math::BigFloat bugs RT 43692 and RT 77105 which make
+ # a right mess of things. Watch this:
+ # my $n = Math::BigFloat->new(11); $n->accuracy(64); say $n**1.1; # 13.98
+ # my $n = Math::BigFloat->new(11); $n->accuracy(67); say $n**1.1; # 29.98
+ # We can fix some issues with large exponents (e.g. 6^-40.5) by turning it
+ # into (6^-(40.5/4))^4 (assuming the base is positive). Without that hack,
+ # none of this would work at all.
+
+ $x = Math::BigFloat->new($x);
+ my $superx = 1;
+ my $subx = Math::BigFloat->new($x);
+ while ($subx > 8) {
+ $superx *= 2;
+ $subx /= 2;
+ }
+
+ # Go with the basic formula for large x, as it best works around the mess,
+ # though is unfortunately much slower.
+ if ($x > 30) {
+ my $sum = 0.0;
+ for my $k (4 .. 1000) {
+ my $term = ( $k ** -$subx ) ** $superx;
+ $sum += $term;
+ last if $term < ($sum*$tol);
+ }
+ for my $k (3, 2) {
+ my $term = ( $k ** -$subx ) ** $superx;
+ $sum += $term;
+ }
+ return $sum;
+ }
+ #if ($x > 25) {
+ # my $sum = 0.0;
+ # my $divisor = 1.0 - ((2 ** -$subx) ** $superx);
+ # for my $k (2 .. 1000) {
+ # my $term = ( (2*$k+1) ** -$subx ) ** $superx;
+ # $sum += $term;
+ # last if $term < ($tol*$divisor);
+ # }
+ # $sum += (3 ** -$subx) ** $superx;
+ # my $t = 1.0 / $divisor;
+ # $sum *= $t;
+ # $sum += ($t - 1.0);
+ # return $sum;
+ #}
+
+ # If we wanted to change the Borwein series being used:
+ # _Recompute_Dk(55);
+
+ if (ref $_Borwein_dk[0] ne 'Math::BigInt') {
+ @_Borwein_dk = map { Math::BigInt->new("$_") } @_Borwein_dk;
+ }
+
+ my $n = $_Borwein_n;
+ my $intermediate_accuracy = undef;
+ my $one = Math::BigFloat->bone;
+ $one->accuracy($intermediate_accuracy) if defined $intermediate_accuracy;
+
+ my $d1 = Math::BigFloat->new(2);
+ $d1->accuracy($intermediate_accuracy) if defined $intermediate_accuracy;
+ # with bignum on, $d1->bpow($one-$x) doesn't change d1 !
+ $d1 = $d1 ** ($one - $x);
+ my $divisor = $one->copy->bsub($d1)->bmul(-$_Borwein_dk[$n]);
+
+ $tol = $divisor->copy->bmul($tol)->babs();
+
+ my $sum = Math::BigFloat->bzero;
+ $sum->accuracy($intermediate_accuracy) if defined $intermediate_accuracy;
+ foreach my $k (1 .. $n-1) {
+ my $term = Math::BigFloat->new( $_Borwein_dk[$k] - $_Borwein_dk[$n] );
+ $term *= -1 if $k % 2;
+ $term->accuracy($intermediate_accuracy) if defined $intermediate_accuracy;
+ my $den = Math::BigFloat->new($k+1);
+ $den->accuracy($intermediate_accuracy) if defined $intermediate_accuracy;
+ $den = ($den ** $subx) ** $superx;
+ $term /= $den;
+ $sum += $term;
+ last if $term->copy->babs() < $tol;
+ }
+ $sum += Math::BigFloat->new( $one - $_Borwein_dk[$n] ); # term k=0
+ $sum->bdiv( $divisor );
+ $sum->bsub(1);
+ return $sum;
}
# Riemann R function
sub RiemannR {
my($x, $tol) = @_;
- $tol = 1e-16 unless defined $tol;
- my $sum = 0.0;
my($y, $t);
- my $c = 0.0;
croak "Invalid input to ReimannR: x must be > 0" if $x <= 0;
- $x = new Math::BigFloat "$x" if defined $bignum::VERSION && ref($x) ne 'Math::BigFloat';
+ if (!defined $bignum::VERSION && ref($x) !~ /^Math::Big/) {
+ $tol = 1e-16 unless defined $tol;
+ my $sum = 0.0;
+ my $c = 0.0;
+
+ $y = 1.0-$c; $t = $sum+$y; $c = ($t-$sum)-$y; $sum = $t;
+ my $flogx = log($x);
+ my $part_term = 1.0;
+ for my $k (1 .. 10000) {
+ # Small k from table, larger k from function
+ my $zeta = ($k <= $#_Riemann_Zeta_Table)
+ ? $_Riemann_Zeta_Table[$k+1-2]
+ : RiemannZeta($k+1);
+ $part_term *= $flogx / $k;
+ my $term = $part_term / ($k + $k * $zeta);
+ $y = $term-$c; $t = $sum+$y; $c = ($t-$sum)-$y; $sum = $t;
+ last if abs($term/$sum) < $tol;
+ }
+ return $sum;
+ }
+
+ $x = new Math::BigFloat "$x" if ref($x) ne 'Math::BigFloat';
+
+ $tol = 1e-35 unless defined $tol;
+ my $sum = Math::BigFloat->bzero;
+ my $c = Math::BigFloat->bzero;
$y = 1.0-$c; $t = $sum+$y; $c = ($t-$sum)-$y; $sum = $t;
my $flogx = log($x);
- my $part_term = 1.0;
+ my $part_term = Math::BigFloat->bone;
for my $k (1 .. 10000) {
# Small k from table, larger k from function
- my $zeta = ($k <= $#_Riemann_Zeta_Table) ? $_Riemann_Zeta_Table[$k+1-2]
- : RiemannZeta($k+1);
+ my $zeta = ($k <= $#_Riemann_Zeta_Table)
+ ? Math::BigFloat->new($_Riemann_Zeta_Table[$k+1-2])
+ : RiemannZeta($k+1);
$part_term *= $flogx / $k;
my $term = $part_term / ($k + $k * $zeta);
$y = $term-$c; $t = $sum+$y; $c = ($t-$sum)-$y; $sum = $t;
last if abs($term/$sum) < $tol;
}
- $sum;
+ return $sum;
}
1;
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/pkg-perl/packages/libmath-prime-util-perl.git
More information about the Pkg-perl-cvs-commits
mailing list