610 lines
27 KiB
Perl
610 lines
27 KiB
Perl
package Math::Prime::Util::ZetaBigFloat;
|
|
use strict;
|
|
use warnings;
|
|
|
|
BEGIN {
|
|
$Math::Prime::Util::ZetaBigFloat::AUTHORITY = 'cpan:DANAJ';
|
|
$Math::Prime::Util::ZetaBigFloat::VERSION = '0.73';
|
|
}
|
|
|
|
BEGIN {
|
|
do { require Math::BigInt; Math::BigInt->import(try=>"GMP,Pari"); }
|
|
unless defined $Math::BigInt::VERSION;
|
|
use Math::BigFloat;
|
|
}
|
|
#my $_oldacc = Math::BigFloat->accuracy();
|
|
#Math::BigFloat->accuracy(undef);
|
|
|
|
|
|
# Riemann Zeta($k) for integer $k.
|
|
# So many terms and digits are used so we can quickly do bignum R.
|
|
my @_Riemann_Zeta_Table = (
|
|
'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',
|
|
'0.00000000000000000000000000000078886090522101180735205378276604136878962534',
|
|
'0.00000000000000000000000000000039443045261050590335263935513575963608141044',
|
|
'0.00000000000000000000000000000019721522630525295156852383215213909988473843',
|
|
'0.000000000000000000000000000000098607613152626475748329967604159218377505181',
|
|
'0.000000000000000000000000000000049303806576313237862187667644776975622245754',
|
|
'0.000000000000000000000000000000024651903288156618927101395103287812527732549',
|
|
'0.000000000000000000000000000000012325951644078309462219884645277065145764150',
|
|
'0.0000000000000000000000000000000061629758220391547306663380205162648609383631',
|
|
'0.0000000000000000000000000000000030814879110195773651853009095507130250105264',
|
|
'0.0000000000000000000000000000000015407439555097886825433610878728841686496904',
|
|
'0.00000000000000000000000000000000077037197775489434125525075496895150086398231',
|
|
'0.00000000000000000000000000000000038518598887744717062214878116197893873445220',
|
|
'0.00000000000000000000000000000000019259299443872358530924885847349054449873362',
|
|
'0.000000000000000000000000000000000096296497219361792654015918534245633717541108',
|
|
'0.000000000000000000000000000000000048148248609680896326805122366289604787579935',
|
|
'0.000000000000000000000000000000000024074124304840448163334948882867065229914248',
|
|
'0.000000000000000000000000000000000012037062152420224081644937008007620275295506',
|
|
'0.0000000000000000000000000000000000060185310762101120408149560261951727031681191',
|
|
'0.0000000000000000000000000000000000030092655381050560204049738538280405431094080',
|
|
'0.0000000000000000000000000000000000015046327690525280102016522071575050028177934',
|
|
'0.00000000000000000000000000000000000075231638452626400510054786365991407868525313',
|
|
'0.00000000000000000000000000000000000037615819226313200255018118519034423181524371',
|
|
'0.00000000000000000000000000000000000018807909613156600127505967704863451341028548',
|
|
'0.000000000000000000000000000000000000094039548065783000637519533342138055875645097',
|
|
'0.000000000000000000000000000000000000047019774032891500318756331610342627662060287',
|
|
'0.000000000000000000000000000000000000023509887016445750159377020784929180405960294',
|
|
'0.000000000000000000000000000000000000011754943508222875079688128719050545728002924',
|
|
'0.0000000000000000000000000000000000000058774717541114375398439371350539247056872356',
|
|
'0.0000000000000000000000000000000000000029387358770557187699219261593698463000750878',
|
|
'0.0000000000000000000000000000000000000014693679385278593849609489436325511324487536',
|
|
'0.00000000000000000000000000000000000000073468396926392969248046975979881822702829326',
|
|
'0.00000000000000000000000000000000000000036734198463196484624023330922692333378216377',
|
|
'0.00000000000000000000000000000000000000018367099231598242312011613105596640698043218',
|
|
'0.000000000000000000000000000000000000000091835496157991211560057891008818116853335663',
|
|
'0.000000000000000000000000000000000000000045917748078995605780028887331354029547708393',
|
|
'0.000000000000000000000000000000000000000022958874039497802890014424274658671814201226',
|
|
'0.000000000000000000000000000000000000000011479437019748901445007205673656554920549667',
|
|
'0.0000000000000000000000000000000000000000057397185098744507225036006822706837980911955',
|
|
'0.0000000000000000000000000000000000000000028698592549372253612517996229494773449843879',
|
|
'0.0000000000000000000000000000000000000000014349296274686126806258995720794504878051247',
|
|
'0.00000000000000000000000000000000000000000071746481373430634031294970624129584900687276',
|
|
'0.00000000000000000000000000000000000000000035873240686715317015647482652117145953820656',
|
|
'0.00000000000000000000000000000000000000000017936620343357658507823740439409357478069335',
|
|
'0.000000000000000000000000000000000000000000089683101716788292539118699241549402394210037',
|
|
'0.000000000000000000000000000000000000000000044841550858394146269559348635608906198392806',
|
|
'0.000000000000000000000000000000000000000000022420775429197073134779673989415854766292332',
|
|
'0.000000000000000000000000000000000000000000011210387714598536567389836885245061272178142',
|
|
'0.0000000000000000000000000000000000000000000056051938572992682836949184061349085990997301',
|
|
'0.0000000000000000000000000000000000000000000028025969286496341418474591909049136205534180',
|
|
'0.0000000000000000000000000000000000000000000014012984643248170709237295913982765839445600',
|
|
'0.00000000000000000000000000000000000000000000070064923216240853546186479434774488319489698',
|
|
'0.00000000000000000000000000000000000000000000035032461608120426773093239672340797200498749',
|
|
'0.00000000000000000000000000000000000000000000017516230804060213386546619821154916280500674',
|
|
'0.000000000000000000000000000000000000000000000087581154020301066932733099055722973670007705',
|
|
'0.000000000000000000000000000000000000000000000043790577010150533466366549511177617590838630',
|
|
'0.000000000000000000000000000000000000000000000021895288505075266733183274750027519047364241',
|
|
'0.000000000000000000000000000000000000000000000010947644252537633366591637373159996274330429',
|
|
'0.0000000000000000000000000000000000000000000000054738221262688166832958186859620770540479841',
|
|
'0.0000000000000000000000000000000000000000000000027369110631344083416479093427750648326515819',
|
|
'0.0000000000000000000000000000000000000000000000013684555315672041708239546713188745182016542',
|
|
'0.00000000000000000000000000000000000000000000000068422776578360208541197733563655129305944821',
|
|
'0.00000000000000000000000000000000000000000000000034211388289180104270598866781064699118259780',
|
|
'0.00000000000000000000000000000000000000000000000017105694144590052135299433390278061047559013',
|
|
'0.000000000000000000000000000000000000000000000000085528470722950260676497166950542676865892145',
|
|
'0.000000000000000000000000000000000000000000000000042764235361475130338248583474988795642311765',
|
|
'0.000000000000000000000000000000000000000000000000021382117680737565169124291737400216890944447',
|
|
'0.000000000000000000000000000000000000000000000000010691058840368782584562145868668714802068411',
|
|
'0.0000000000000000000000000000000000000000000000000053455294201843912922810729343238928532329351',
|
|
'0.0000000000000000000000000000000000000000000000000026727647100921956461405364671584582440160440',
|
|
'0.0000000000000000000000000000000000000000000000000013363823550460978230702682335780663944745475',
|
|
'0.00000000000000000000000000000000000000000000000000066819117752304891153513411678864562139278223',
|
|
'0.00000000000000000000000000000000000000000000000000033409558876152445576756705839419361874822728',
|
|
'0.00000000000000000000000000000000000000000000000000016704779438076222788378352919705374539139236',
|
|
'0.000000000000000000000000000000000000000000000000000083523897190381113941891764598512518034789088',
|
|
'0.000000000000000000000000000000000000000000000000000041761948595190556970945882299251474130425513',
|
|
'0.000000000000000000000000000000000000000000000000000020880974297595278485472941149624142102889746',
|
|
'0.000000000000000000000000000000000000000000000000000010440487148797639242736470574811539397337203',
|
|
'0.0000000000000000000000000000000000000000000000000000052202435743988196213682352874055924806327115',
|
|
'0.0000000000000000000000000000000000000000000000000000026101217871994098106841176437027371676377257',
|
|
'0.0000000000000000000000000000000000000000000000000000013050608935997049053420588218513488929259862',
|
|
'0.00000000000000000000000000000000000000000000000000000065253044679985245267102941092566788283203421',
|
|
);
|
|
# Convert to BigFloat objects.
|
|
@_Riemann_Zeta_Table = map { Math::BigFloat->new($_) } @_Riemann_Zeta_Table;
|
|
# for k = 1 .. n : (1 / (zeta(k+1) * k + k)
|
|
# Makes RiemannR run about twice as fast.
|
|
my @_Riemann_Zeta_Premult;
|
|
my $_Riemann_Zeta_premult_accuracy = 0;
|
|
|
|
# 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 = ();
|
|
my $orig_acc = Math::BigFloat->accuracy();
|
|
Math::BigFloat->accuracy($nterms);
|
|
foreach my $k (0 .. $nterms) {
|
|
my $sum = Math::BigInt->bzero;
|
|
my $num = Math::BigInt->new($nterms-1)->bfac();
|
|
foreach my $i (0 .. $k) {
|
|
my $den = Math::BigInt->new($nterms - $i)->bfac * Math::BigInt->new(2*$i)->bfac;
|
|
$sum += $num->copy->bdiv($den);
|
|
$num->bmul(4 * ($nterms+$i));
|
|
}
|
|
$sum->bmul($nterms);
|
|
$_Borwein_dk[$k] = $sum;
|
|
}
|
|
Math::BigFloat->accuracy($orig_acc);
|
|
}
|
|
|
|
sub RiemannZeta {
|
|
my($ix) = @_;
|
|
|
|
my $x = (ref($ix) eq 'Math::BigFloat') ? $ix->copy : Math::BigFloat->new("$ix");
|
|
$x->accuracy($ix->accuracy) if $ix->accuracy;
|
|
my $xdigits = $ix->accuracy() || Math::BigFloat->accuracy() || Math::BigFloat->div_scale();
|
|
|
|
if ($x == int($x) && $xdigits <= 44 && (int($x)-2) <= $#_Riemann_Zeta_Table) {
|
|
my $izeta = $_Riemann_Zeta_Table[int($x)-2]->copy;
|
|
$izeta->bround($xdigits);
|
|
return $izeta;
|
|
}
|
|
|
|
# Note, this code likely will not work correctly without fixes for RTs:
|
|
#
|
|
# 43692 : blog and others broken
|
|
# 43460 : exp and powers broken
|
|
#
|
|
# E.g:
|
|
# 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
|
|
#
|
|
# There is a hack that tries to work around some of the problem, but it
|
|
# can't cover everything and it slows things down a lot. There just isn't
|
|
# any way to do this if the basic math operations don't work right.
|
|
|
|
my $orig_acc = Math::BigFloat->accuracy();
|
|
my $extra_acc = 5;
|
|
if ($x > 15 && $x <= 50) { $extra_acc = 15; }
|
|
|
|
$xdigits += $extra_acc;
|
|
Math::BigFloat->accuracy($xdigits);
|
|
$x->accuracy($xdigits);
|
|
my $zero= $x->copy->bzero;
|
|
my $one = $zero->copy->binc;
|
|
my $two = $one->copy->binc;
|
|
|
|
my $tol = ref($x)->new('0.' . '0' x ($xdigits-1) . '1');
|
|
|
|
# Note: with bignum on, $d1->bpow($one-$x) doesn't change d1 !
|
|
|
|
# This is a hack to turn 6^-40.5 into (6^-(40.5/4))^4. It helps work around
|
|
# the two RTs listed earlier, though does not completely fix their bugs.
|
|
# It has the downside of making integer arguments very slow.
|
|
|
|
my $superx = Math::BigInt->bone;
|
|
my $subx = $x->copy;
|
|
my $intx = int("$x");
|
|
if ($Math::BigFloat::VERSION < 1.9996 || $x != $intx) {
|
|
while ($subx > 1) {
|
|
$superx->blsft(1);
|
|
$subx /= $two;
|
|
}
|
|
}
|
|
|
|
if (1 && $x == $intx && $x >= 2 && !($intx & 1) && $intx < 100) {
|
|
# Mathworld equation 63. How fast this is relative to the others is
|
|
# dependent on the backend library and if we have MPUGMP.
|
|
$x = int("$x");
|
|
my $den = Math::Prime::Util::factorial($x);
|
|
$xdigits -= $extra_acc;
|
|
$extra_acc += length($den);
|
|
$xdigits += $extra_acc;
|
|
$one->accuracy($xdigits); $two->accuracy($xdigits);
|
|
Math::BigFloat->accuracy($xdigits);
|
|
$subx->accuracy($xdigits); $superx->accuracy($xdigits);
|
|
my $Pix = Math::Prime::Util::Pi($xdigits)->bpow($subx)->bpow($superx);
|
|
my $Bn = Math::Prime::Util::bernreal($x,$xdigits); $Bn = -$Bn if $Bn < 0;
|
|
my $twox1 = $two->copy->bpow($x-1);
|
|
#my $num = $Pix * $Bn * $twox1;
|
|
#my $res = $num->bdiv($den)->bdec->bround($xdigits - $extra_acc);
|
|
my $res = $Bn->bdiv($den)->bmul($Pix)->bmul($twox1)->bdec
|
|
->bround($xdigits-$extra_acc);
|
|
Math::BigFloat->accuracy($orig_acc);
|
|
return $res;
|
|
}
|
|
|
|
# Go with the basic formula for large x.
|
|
if (1 && $x >= 50) {
|
|
my $negsubx = $subx->copy->bneg;
|
|
my $sum = $zero->copy;
|
|
my $k = $two->copy->binc;
|
|
while ($k->binc <= 1000) {
|
|
my $term = $k->copy->bpow($negsubx)->bpow($superx);
|
|
$sum += $term;
|
|
last if $term < ($sum*$tol);
|
|
}
|
|
$k = $two+$two;
|
|
$k->bdec(); $sum += $k->copy->bpow($negsubx)->bpow($superx);
|
|
$k->bdec(); $sum += $k->copy->bpow($negsubx)->bpow($superx);
|
|
$sum->bround($xdigits-$extra_acc);
|
|
Math::BigFloat->accuracy($orig_acc);
|
|
return $sum;
|
|
}
|
|
|
|
{
|
|
my $dig = int($_Borwein_n / 1.3)+1;
|
|
_Recompute_Dk( int($xdigits * 1.3) + 4 ) if $dig < $xdigits;
|
|
}
|
|
|
|
if (ref $_Borwein_dk[0] ne 'Math::BigInt') {
|
|
@_Borwein_dk = map { Math::BigInt->new("$_") } @_Borwein_dk;
|
|
}
|
|
|
|
my $n = $_Borwein_n;
|
|
|
|
my $d1 = $two ** ($one - $x);
|
|
my $divisor = ($one - $d1) * $_Borwein_dk[$n];
|
|
$divisor->bneg;
|
|
$tol = ($divisor * $tol)->babs();
|
|
|
|
my ($sum, $bigk) = ($zero->copy, $one->copy);
|
|
my $negsubx = $subx->copy->bneg;
|
|
foreach my $k (1 .. $n-1) {
|
|
my $den = $bigk->binc()->copy->bpow($negsubx)->bpow($superx);
|
|
my $term = ($k % 2) ? ($_Borwein_dk[$n] - $_Borwein_dk[$k])
|
|
: ($_Borwein_dk[$k] - $_Borwein_dk[$n]);
|
|
$term = Math::BigFloat->new($term) unless ref($term) eq 'Math::BigFloat';
|
|
$sum += $term * $den;
|
|
last if $term->copy->babs() < $tol;
|
|
}
|
|
$sum += $_Borwein_dk[0] - $_Borwein_dk[$n];
|
|
$sum = $sum->bdiv($divisor);
|
|
$sum->bdec->bround($xdigits-$extra_acc);
|
|
Math::BigFloat->accuracy($orig_acc);
|
|
return $sum;
|
|
}
|
|
|
|
# Riemann R function
|
|
sub RiemannR {
|
|
my($x) = @_;
|
|
|
|
if (ref($x) eq 'Math::BigInt') {
|
|
my $xacc = $x->accuracy();
|
|
$x = Math::BigFloat->new($x);
|
|
$x->accuracy($xacc) if $xacc;
|
|
}
|
|
$x = Math::BigFloat->new("$x") if ref($x) ne 'Math::BigFloat';
|
|
my $xdigits = $x->accuracy || Math::BigFloat->accuracy() || Math::BigFloat->div_scale();
|
|
my $extra_acc = 2;
|
|
$xdigits += $extra_acc;
|
|
my $orig_acc = Math::BigFloat->accuracy();
|
|
Math::BigFloat->accuracy($xdigits);
|
|
$x->accuracy($xdigits);
|
|
my $tol = $x->copy->bone->brsft($xdigits-1, 10);
|
|
my $sum = $x->copy->bone;
|
|
|
|
if ($xdigits <= length($x->copy->as_int->bstr())) {
|
|
|
|
for my $k (1 .. 1000) {
|
|
my $mob = Math::Prime::Util::moebius($k);
|
|
next if $mob == 0;
|
|
$mob = Math::BigFloat->new($mob);
|
|
my $term = $mob->bdiv($k) *
|
|
Math::Prime::Util::LogarithmicIntegral($x->copy->broot($k));
|
|
$sum += $term;
|
|
#warn "k = $k term = $term sum = $sum\n";
|
|
last if abs($term) < ($tol * abs($sum));
|
|
}
|
|
|
|
} else {
|
|
|
|
my ($flogx, $part_term, $fone, $bigk)
|
|
= (log($x), Math::BigFloat->bone, Math::BigFloat->bone, Math::BigInt->bone);
|
|
|
|
if ($_Riemann_Zeta_premult_accuracy < $xdigits) {
|
|
@_Riemann_Zeta_Premult = ();
|
|
$_Riemann_Zeta_premult_accuracy = $xdigits;
|
|
}
|
|
|
|
for my $k (1 .. 10000) {
|
|
my $zeta_term = $_Riemann_Zeta_Premult[$k-1];
|
|
if (!defined $zeta_term) {
|
|
my $zeta = ($xdigits > 44) ? undef : $_Riemann_Zeta_Table[$k-1];
|
|
if (!defined $zeta) {
|
|
my $kz = $fone->copy->badd($bigk); # kz is k+1
|
|
if (($k+1) >= 100 && $xdigits <= 40) {
|
|
# For this accuracy level, two terms are more than enough. Also,
|
|
# we should be able to miss the Math::BigFloat accuracy bug. If we
|
|
# try to do this for higher accuracy, things will go very bad.
|
|
$zeta = Math::BigFloat->new(3)->bpow(-$kz)
|
|
+ Math::BigFloat->new(2)->bpow(-$kz);
|
|
} else {
|
|
$zeta = Math::Prime::Util::ZetaBigFloat::RiemannZeta( $kz );
|
|
}
|
|
}
|
|
$zeta_term = $fone / ($zeta * $bigk + $bigk);
|
|
$_Riemann_Zeta_Premult[$k-1] = $zeta_term if defined $_Riemann_Zeta_Table[$k-1];
|
|
}
|
|
$part_term *= $flogx / $bigk;
|
|
my $term = $part_term * $zeta_term;
|
|
$sum += $term;
|
|
#warn "k = $k term = $term sum = $sum\n";
|
|
last if $term < ($tol*$sum);
|
|
$bigk->binc;
|
|
}
|
|
|
|
}
|
|
$sum->bround($xdigits-$extra_acc);
|
|
Math::BigFloat->accuracy($orig_acc);
|
|
return $sum;
|
|
}
|
|
|
|
#Math::BigFloat->accuracy($_oldacc);
|
|
#undef $_oldacc;
|
|
|
|
1;
|
|
|
|
__END__
|
|
|
|
|
|
# ABSTRACT: Perl Big Float versions of Riemann Zeta and R functions
|
|
|
|
=pod
|
|
|
|
=encoding utf8
|
|
|
|
|
|
=head1 NAME
|
|
|
|
Math::Prime::Util::ZetaBigFloat - Perl Big Float versions of Riemann Zeta and R functions
|
|
|
|
|
|
=head1 VERSION
|
|
|
|
Version 0.73
|
|
|
|
|
|
=head1 SYNOPSIS
|
|
|
|
Math::BigFloat versions`of the Riemann Zeta and Riemann R functions. These
|
|
are kept in a separate module because they use a lot of big tables that we'd
|
|
prefer to only load if needed.
|
|
|
|
|
|
=head1 DESCRIPTION
|
|
|
|
Pure Perl implementations of Riemann Zeta and Riemann R using Math::BigFloat.
|
|
These functions are used if:
|
|
|
|
=over 4
|
|
|
|
=item The input is a BigInt, a BigFloat, or the bignum module has been loaded.
|
|
|
|
=item The L<Math::Prime::Util::GMP> module is not available or old.
|
|
|
|
=back
|
|
|
|
If you use these functions a lot, I B<highly> recommend you install
|
|
L<Math::Prime::Util::GMP>, which the main L<Math::Prime::Util> functions
|
|
will find.
|
|
These give B<much> better performance, and better accuracy. You can also
|
|
use L<Math::Pari> and L<Math::MPFR> for the Riemann Zeta function.
|
|
|
|
|
|
=head1 FUNCTIONS
|
|
|
|
=head2 RiemannZeta
|
|
|
|
my $z = RiemannZeta($s);
|
|
|
|
Given a floating point input C<s> where C<s E<gt>= 0.5>, 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>
|
|
|
|
Results are calculated using either Borwein (1991) algorithm 2, or the basic
|
|
series. Full input accuracy is attempted, but there are defects in
|
|
Math::BigFloat with high accuracy computations that make this difficult.
|
|
|
|
|
|
=head2 RiemannR
|
|
|
|
my $r = RiemannR($x);
|
|
|
|
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 about 35 digits.
|
|
|
|
|
|
=head1 LIMITATIONS
|
|
|
|
Bugs in Math::BigFloat (RT 43692, RT 77105) cause many problems with this code.
|
|
I've attempted to work around them, but it is possible there are cases they
|
|
miss.
|
|
|
|
The accuracy goals (35 digits) are sometimes missed by a digit or two.
|
|
|
|
|
|
=head1 PERFORMANCE
|
|
|
|
Performance is quite bad.
|
|
|
|
|
|
=head1 SEE ALSO
|
|
|
|
L<Math::Prime::Util>
|
|
|
|
L<Math::Prime::Util::GMP>
|
|
|
|
L<Math::MPFR>
|
|
|
|
L<Math::Pari>
|
|
|
|
|
|
=head1 AUTHORS
|
|
|
|
Dana Jacobsen E<lt>dana@acm.orgE<gt>
|
|
|
|
|
|
=head1 COPYRIGHT
|
|
|
|
Copyright 2012 by Dana Jacobsen E<lt>dana@acm.orgE<gt>
|
|
|
|
This program is free software; you can redistribute it and/or modify it under the same terms as Perl itself.
|
|
|
|
=cut
|