package Math::Prime::Util::ZetaBigFloat; use strict; use warnings; BEGIN { $Math::Prime::Util::ZetaBigFloat::AUTHORITY = 'cpan:DANAJ'; $Math::Prime::Util::ZetaBigFloat::VERSION = '0.14'; } use Math::BigFloat; # 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; # 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"; } } sub RiemannZeta { my($x) = @_; $x = new Math::BigFloat "$x" if ref($x) ne 'Math::BigFloat'; my $xdigits = $x->accuracy || Math::BigFloat->accuracy() || Math::BigFloat->div_scale(); return $_Riemann_Zeta_Table[int($x)-2] if $x == int($x) && defined $_Riemann_Zeta_Table[int($x)-2] && $xdigits <= 44; my $tol = 0.0 + "1e-$xdigits"; # 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. my $superx = 1; my $subx = $x->copy; 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 = Math::BigFloat->bzero; $sum->accuracy($xdigits); for my $k (4 .. 1000) { my $term = ( $k ** -$subx ) ** $superx; $sum += $term; last if $term < ($sum*$tol); } $sum += ( 3 ** -$subx ) ** $superx; $sum += ( 2 ** -$subx ) ** $superx; 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; #} { 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 $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); #$sum->fround($xdigits); return $sum; } # Riemann R function sub RiemannR { my($x) = @_; $x = new Math::BigFloat "$x" if ref($x) ne 'Math::BigFloat'; my $xdigits = $x->accuracy || Math::BigFloat->accuracy() || Math::BigFloat->div_scale(); my $tol = 0.0 + "1e-$xdigits"; # TODO: The default table is only 44 digits. if ( (scalar @_Riemann_Zeta_Premult == 0) || ($_Riemann_Zeta_Premult_accuracy < $xdigits) ) { $_Riemann_Zeta_Premult_accuracy = $xdigits; @_Riemann_Zeta_Premult = map { my $v = Math::BigFloat->bone; $v->accuracy($xdigits); $v / ($_Riemann_Zeta_Table[$_-1] * $_ + $_) } (1 .. @_Riemann_Zeta_Table); } my $sum = Math::BigFloat->bone; my $flogx = log($x); my $part_term = Math::BigFloat->bone; for my $k (1 .. 10000) { my $zeta_term = $_Riemann_Zeta_Premult[$k-1]; if (!defined $zeta_term) { my $zeta = $_Riemann_Zeta_Table[$k-1]; if (!defined $zeta) { my $kz = Math::BigFloat->new($k+1); $kz->accuracy($xdigits); if ($kz >= 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 = Math::BigFloat->bone / ($zeta * $k + $k); } $part_term *= $flogx / $k; my $term = $part_term * $zeta_term; #warn "k = $k term = $term sum = $sum\n"; $sum += $term; last if $term < ($tol*$sum); } return $sum; } 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.14 =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 Math::MPFR module is not available. =back If you use these functions a lot, I B recommend you install L, which the main L functions will find. These give B better performance, and better accuracy. You can also use L for the Riemann Zeta function. =head1 FUNCTIONS =head2 RiemannZeta my $z = RiemannZeta($s); Given a floating point input C where C= 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. 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 apply workarounds, 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 L L =head1 AUTHORS Dana Jacobsen Edana@acm.orgE =head1 COPYRIGHT Copyright 2012 by Dana Jacobsen Edana@acm.orgE This program is free software; you can redistribute it and/or modify it under the same terms as Perl itself. =cut