File Coverage

blib/lib/Math/Prime/Util/ZetaBigFloat.pm
Criterion Covered Total %
statement 108 159 67.9
branch 26 56 46.4
condition 11 30 36.6
subroutine 8 8 100.0
pod 2 2 100.0
total 155 255 60.7


line stmt bran cond sub pod time code
1             package Math::Prime::Util::ZetaBigFloat;
2 2     2   12 use strict;
  2         4  
  2         68  
3 2     2   10 use warnings;
  2         4  
  2         101  
4              
5             BEGIN {
6 2     2   5 $Math::Prime::Util::ZetaBigFloat::AUTHORITY = 'cpan:DANAJ';
7 2         98 $Math::Prime::Util::ZetaBigFloat::VERSION = '0.69';
8             }
9              
10 0         0 BEGIN {
11 2 50   2   4121 do { require Math::BigInt; Math::BigInt->import(try=>"GMP,Pari"); }
  0         0  
  0         0  
12             unless defined $Math::BigInt::VERSION;
13 2     2   11 use Math::BigFloat;
  2         4  
  2         13  
14             }
15             #my $_oldacc = Math::BigFloat->accuracy();
16             #Math::BigFloat->accuracy(undef);
17              
18              
19             # Riemann Zeta($k) for integer $k.
20             # So many terms and digits are used so we can quickly do bignum R.
21             my @_Riemann_Zeta_Table = (
22             '0.64493406684822643647241516664602518921894990', # zeta(2) - 1
23             '0.20205690315959428539973816151144999076498629',
24             '0.082323233711138191516003696541167902774750952',
25             '0.036927755143369926331365486457034168057080920',
26             '0.017343061984449139714517929790920527901817490',
27             '0.0083492773819228268397975498497967595998635606',
28             '0.0040773561979443393786852385086524652589607906',
29             '0.0020083928260822144178527692324120604856058514',
30             '0.00099457512781808533714595890031901700601953156',
31             '0.00049418860411946455870228252646993646860643576',
32             '0.00024608655330804829863799804773967096041608846',
33             '0.00012271334757848914675183652635739571427510590',
34             '0.000061248135058704829258545105135333747481696169',
35             '0.000030588236307020493551728510645062587627948707',
36             '0.000015282259408651871732571487636722023237388990',
37             '0.0000076371976378997622736002935630292130882490903',
38             '0.0000038172932649998398564616446219397304546972190',
39             '0.0000019082127165539389256569577951013532585711448',
40             '0.00000095396203387279611315203868344934594379418741',
41             '0.00000047693298678780646311671960437304596644669478',
42             '0.00000023845050272773299000364818675299493504182178',
43             '0.00000011921992596531107306778871888232638725499778',
44             '0.000000059608189051259479612440207935801227503918837',
45             '0.000000029803503514652280186063705069366011844730920',
46             '0.000000014901554828365041234658506630698628864788168',
47             '0.0000000074507117898354294919810041706041194547190319',
48             '0.0000000037253340247884570548192040184024232328930593',
49             '0.0000000018626597235130490064039099454169480616653305',
50             '0.00000000093132743241966818287176473502121981356795514',
51             '0.00000000046566290650337840729892332512200710626918534',
52             '0.00000000023283118336765054920014559759404950248298228',
53             '0.00000000011641550172700519775929738354563095165224717',
54             '0.000000000058207720879027008892436859891063054173122605',
55             '0.000000000029103850444970996869294252278840464106981987',
56             '0.000000000014551921891041984235929632245318420983808894',
57             '0.0000000000072759598350574810145208690123380592648509256',
58             '0.0000000000036379795473786511902372363558732735126460284',
59             '0.0000000000018189896503070659475848321007300850305893096',
60             '0.00000000000090949478402638892825331183869490875386000099',
61             '0.00000000000045474737830421540267991120294885703390452991',
62             '0.00000000000022737368458246525152268215779786912138298220',
63             '0.00000000000011368684076802278493491048380259064374359028',
64             '0.000000000000056843419876275856092771829675240685530571589',
65             '0.000000000000028421709768893018554550737049426620743688265',
66             '0.000000000000014210854828031606769834307141739537678698606',
67             '0.0000000000000071054273952108527128773544799568000227420436',
68             '0.0000000000000035527136913371136732984695340593429921456555',
69             '0.0000000000000017763568435791203274733490144002795701555086',
70             '0.00000000000000088817842109308159030960913863913863256088715',
71             '0.00000000000000044408921031438133641977709402681213364596031',
72             '0.00000000000000022204460507980419839993200942046539642366543',
73             '0.00000000000000011102230251410661337205445699213827024832229',
74             '0.000000000000000055511151248454812437237365905094302816723551',
75             '0.000000000000000027755575621361241725816324538540697689848904',
76             '0.000000000000000013877787809725232762839094906500221907718625',
77             '0.0000000000000000069388939045441536974460853262498092748358742',
78             '0.0000000000000000034694469521659226247442714961093346219504706',
79             '0.0000000000000000017347234760475765720489729699375959074780545',
80             '0.00000000000000000086736173801199337283420550673429514879071415',
81             '0.00000000000000000043368086900206504874970235659062413612547801',
82             '0.00000000000000000021684043449972197850139101683209845761574010',
83             '0.00000000000000000010842021724942414063012711165461382589364744',
84             '0.000000000000000000054210108624566454109187004043886337150634224',
85             '0.000000000000000000027105054312234688319546213119497764318887282',
86             '0.000000000000000000013552527156101164581485233996826928328981877',
87             '0.0000000000000000000067762635780451890979952987415566862059812586',
88             '0.0000000000000000000033881317890207968180857031004508368340311585',
89             '0.0000000000000000000016940658945097991654064927471248619403036418',
90             '0.00000000000000000000084703294725469983482469926091821675222838642',
91             '0.00000000000000000000042351647362728333478622704833579344088109717',
92             '0.00000000000000000000021175823681361947318442094398180025869417612',
93             '0.00000000000000000000010587911840680233852265001539238398470699902',
94             '0.000000000000000000000052939559203398703238139123029185055866375629',
95             '0.000000000000000000000026469779601698529611341166842038715592556134',
96             '0.000000000000000000000013234889800848990803094510250944989684323826',
97             '0.0000000000000000000000066174449004244040673552453323082200147137975',
98             '0.0000000000000000000000033087224502121715889469563843144048092764894',
99             '0.0000000000000000000000016543612251060756462299236771810488297723589',
100             '0.00000000000000000000000082718061255303444036711056167440724040096811',
101             '0.00000000000000000000000041359030627651609260093824555081412852575873',
102             '0.00000000000000000000000020679515313825767043959679193468950443365312',
103             '0.00000000000000000000000010339757656912870993284095591745860911079606',
104             '0.000000000000000000000000051698788284564313204101332166355512893608164',
105             '0.000000000000000000000000025849394142282142681277617708450222269121159',
106             '0.000000000000000000000000012924697071141066700381126118331865309299779',
107             '0.0000000000000000000000000064623485355705318034380021611221670660356864',
108             '0.0000000000000000000000000032311742677852653861348141180266574173608296',
109             '0.0000000000000000000000000016155871338926325212060114057052272720509148',
110             '0.00000000000000000000000000080779356694631620331587381863408997398684847',
111             '0.00000000000000000000000000040389678347315808256222628129858130379479700',
112             '0.00000000000000000000000000020194839173657903491587626465673047518903728',
113             '0.00000000000000000000000000010097419586828951533619250700091044144538432',
114             '0.000000000000000000000000000050487097934144756960847711725486604360898735',
115             '0.000000000000000000000000000025243548967072378244674341937966175648398693',
116             '0.000000000000000000000000000012621774483536189043753999660777148710632765',
117             '0.0000000000000000000000000000063108872417680944956826093943332037500694712',
118             '0.0000000000000000000000000000031554436208840472391098412184847972814371270',
119             '0.0000000000000000000000000000015777218104420236166444327830159601782237092',
120             '0.00000000000000000000000000000078886090522101180735205378276604136878962534',
121             '0.00000000000000000000000000000039443045261050590335263935513575963608141044',
122             '0.00000000000000000000000000000019721522630525295156852383215213909988473843',
123             '0.000000000000000000000000000000098607613152626475748329967604159218377505181',
124             '0.000000000000000000000000000000049303806576313237862187667644776975622245754',
125             '0.000000000000000000000000000000024651903288156618927101395103287812527732549',
126             '0.000000000000000000000000000000012325951644078309462219884645277065145764150',
127             '0.0000000000000000000000000000000061629758220391547306663380205162648609383631',
128             '0.0000000000000000000000000000000030814879110195773651853009095507130250105264',
129             '0.0000000000000000000000000000000015407439555097886825433610878728841686496904',
130             '0.00000000000000000000000000000000077037197775489434125525075496895150086398231',
131             '0.00000000000000000000000000000000038518598887744717062214878116197893873445220',
132             '0.00000000000000000000000000000000019259299443872358530924885847349054449873362',
133             '0.000000000000000000000000000000000096296497219361792654015918534245633717541108',
134             '0.000000000000000000000000000000000048148248609680896326805122366289604787579935',
135             '0.000000000000000000000000000000000024074124304840448163334948882867065229914248',
136             '0.000000000000000000000000000000000012037062152420224081644937008007620275295506',
137             '0.0000000000000000000000000000000000060185310762101120408149560261951727031681191',
138             '0.0000000000000000000000000000000000030092655381050560204049738538280405431094080',
139             '0.0000000000000000000000000000000000015046327690525280102016522071575050028177934',
140             '0.00000000000000000000000000000000000075231638452626400510054786365991407868525313',
141             '0.00000000000000000000000000000000000037615819226313200255018118519034423181524371',
142             '0.00000000000000000000000000000000000018807909613156600127505967704863451341028548',
143             '0.000000000000000000000000000000000000094039548065783000637519533342138055875645097',
144             '0.000000000000000000000000000000000000047019774032891500318756331610342627662060287',
145             '0.000000000000000000000000000000000000023509887016445750159377020784929180405960294',
146             '0.000000000000000000000000000000000000011754943508222875079688128719050545728002924',
147             '0.0000000000000000000000000000000000000058774717541114375398439371350539247056872356',
148             '0.0000000000000000000000000000000000000029387358770557187699219261593698463000750878',
149             '0.0000000000000000000000000000000000000014693679385278593849609489436325511324487536',
150             '0.00000000000000000000000000000000000000073468396926392969248046975979881822702829326',
151             '0.00000000000000000000000000000000000000036734198463196484624023330922692333378216377',
152             '0.00000000000000000000000000000000000000018367099231598242312011613105596640698043218',
153             '0.000000000000000000000000000000000000000091835496157991211560057891008818116853335663',
154             '0.000000000000000000000000000000000000000045917748078995605780028887331354029547708393',
155             '0.000000000000000000000000000000000000000022958874039497802890014424274658671814201226',
156             '0.000000000000000000000000000000000000000011479437019748901445007205673656554920549667',
157             '0.0000000000000000000000000000000000000000057397185098744507225036006822706837980911955',
158             '0.0000000000000000000000000000000000000000028698592549372253612517996229494773449843879',
159             '0.0000000000000000000000000000000000000000014349296274686126806258995720794504878051247',
160             '0.00000000000000000000000000000000000000000071746481373430634031294970624129584900687276',
161             '0.00000000000000000000000000000000000000000035873240686715317015647482652117145953820656',
162             '0.00000000000000000000000000000000000000000017936620343357658507823740439409357478069335',
163             '0.000000000000000000000000000000000000000000089683101716788292539118699241549402394210037',
164             '0.000000000000000000000000000000000000000000044841550858394146269559348635608906198392806',
165             '0.000000000000000000000000000000000000000000022420775429197073134779673989415854766292332',
166             '0.000000000000000000000000000000000000000000011210387714598536567389836885245061272178142',
167             '0.0000000000000000000000000000000000000000000056051938572992682836949184061349085990997301',
168             '0.0000000000000000000000000000000000000000000028025969286496341418474591909049136205534180',
169             '0.0000000000000000000000000000000000000000000014012984643248170709237295913982765839445600',
170             '0.00000000000000000000000000000000000000000000070064923216240853546186479434774488319489698',
171             '0.00000000000000000000000000000000000000000000035032461608120426773093239672340797200498749',
172             '0.00000000000000000000000000000000000000000000017516230804060213386546619821154916280500674',
173             '0.000000000000000000000000000000000000000000000087581154020301066932733099055722973670007705',
174             '0.000000000000000000000000000000000000000000000043790577010150533466366549511177617590838630',
175             '0.000000000000000000000000000000000000000000000021895288505075266733183274750027519047364241',
176             '0.000000000000000000000000000000000000000000000010947644252537633366591637373159996274330429',
177             '0.0000000000000000000000000000000000000000000000054738221262688166832958186859620770540479841',
178             '0.0000000000000000000000000000000000000000000000027369110631344083416479093427750648326515819',
179             '0.0000000000000000000000000000000000000000000000013684555315672041708239546713188745182016542',
180             '0.00000000000000000000000000000000000000000000000068422776578360208541197733563655129305944821',
181             '0.00000000000000000000000000000000000000000000000034211388289180104270598866781064699118259780',
182             '0.00000000000000000000000000000000000000000000000017105694144590052135299433390278061047559013',
183             '0.000000000000000000000000000000000000000000000000085528470722950260676497166950542676865892145',
184             '0.000000000000000000000000000000000000000000000000042764235361475130338248583474988795642311765',
185             '0.000000000000000000000000000000000000000000000000021382117680737565169124291737400216890944447',
186             '0.000000000000000000000000000000000000000000000000010691058840368782584562145868668714802068411',
187             '0.0000000000000000000000000000000000000000000000000053455294201843912922810729343238928532329351',
188             '0.0000000000000000000000000000000000000000000000000026727647100921956461405364671584582440160440',
189             '0.0000000000000000000000000000000000000000000000000013363823550460978230702682335780663944745475',
190             '0.00000000000000000000000000000000000000000000000000066819117752304891153513411678864562139278223',
191             '0.00000000000000000000000000000000000000000000000000033409558876152445576756705839419361874822728',
192             '0.00000000000000000000000000000000000000000000000000016704779438076222788378352919705374539139236',
193             '0.000000000000000000000000000000000000000000000000000083523897190381113941891764598512518034789088',
194             '0.000000000000000000000000000000000000000000000000000041761948595190556970945882299251474130425513',
195             '0.000000000000000000000000000000000000000000000000000020880974297595278485472941149624142102889746',
196             '0.000000000000000000000000000000000000000000000000000010440487148797639242736470574811539397337203',
197             '0.0000000000000000000000000000000000000000000000000000052202435743988196213682352874055924806327115',
198             '0.0000000000000000000000000000000000000000000000000000026101217871994098106841176437027371676377257',
199             '0.0000000000000000000000000000000000000000000000000000013050608935997049053420588218513488929259862',
200             '0.00000000000000000000000000000000000000000000000000000065253044679985245267102941092566788283203421',
201             );
202             # Convert to BigFloat objects.
203             @_Riemann_Zeta_Table = map { Math::BigFloat->new($_) } @_Riemann_Zeta_Table;
204             # for k = 1 .. n : (1 / (zeta(k+1) * k + k)
205             # Makes RiemannR run about twice as fast.
206             my @_Riemann_Zeta_Premult;
207             my $_Riemann_Zeta_premult_accuracy = 0;
208              
209             # Select n = 55, good for 46ish digits of accuracy.
210             my $_Borwein_n = 55;
211             my @_Borwein_dk = (
212             '1',
213             '6051',
214             '6104451',
215             '2462539971',
216             '531648934851',
217             '71301509476803',
218             '6504925195108803',
219             '429144511928164803',
220             '21392068013887742403',
221             '832780518854440804803',
222             '25977281563850106233283',
223             '662753606729324750201283',
224             '14062742362385399866745283',
225             '251634235316509414702211523',
226             '3841603462178827861104812483',
227             '50535961819850087101900022211',
228             '577730330374203014014104003011',
229             '5782012706584553297863989289411',
230             '50984922488525881477588707205571',
231             '398333597655022403279683908035011',
232             '2770992240330783259897072664469955',
233             '17238422988353715312442126057365955',
234             '96274027751337344115352100618133955',
235             '484350301573059857715727453968687555',
236             '2201794236784087151947175826243477955',
237             '9068765987529892610841571032285864387',
238             '33926582279822401059328069515697217987',
239             '115535262182820447663793177744255246787',
240             '358877507711760077538925500462137369027',
241             '1018683886695854101193095537014797787587',
242             '2646951832121008166346437186541363159491',
243             '6306464665572570713623910486640730071491',
244             '13799752848354341643763498672558481367491',
245             '27780237373991939435100856211039992177091',
246             '51543378762608611361377523633779417047491',
247             '88324588911945720951614452340280439890371',
248             '140129110249040241501243929391690331218371',
249             '206452706984942815385219764876242498642371',
250             '283527707823296964404071683165658912154051',
251             '364683602811933600833512164561308162744771',
252             '441935796522635816776473230396154031661507',
253             '508231717051242054487234759342047053767107',
254             '559351463001010719709990637083458540691907',
255             '594624787018881191308291683229515933311427',
256             '616297424973434835299724300924272199623107',
257             '628083443816135918099559567176252011864515',
258             '633714604276098212796088600263676671320515',
259             '636056734158553360761837806887547188568515',
260             '636894970116484676875895417679248215794115',
261             '637149280289288581322870186196318041432515',
262             '637213397278310656625865036925470191411651',
263             '637226467136294189739463288384528579584451',
264             '637228536449134002301138291602841035366851',
265             '637228775173095037281299181461988671775171',
266             '637228793021615488494769154535569803469251',
267             '637228793670652595811622608101881844621763',
268             );
269             # "An Efficient Algorithm for the Riemann Zeta Function", Borwein, 1991.
270             # About 1.3n terms are needed for n digits of accuracy.
271             sub _Recompute_Dk {
272 2     2   4 my $nterms = shift;
273 2         3 $_Borwein_n = $nterms;
274 2         43 @_Borwein_dk = ();
275 2         8 my $orig_acc = Math::BigFloat->accuracy();
276 2         26 Math::BigFloat->accuracy($nterms);
277 2         35 foreach my $k (0 .. $nterms) {
278 139         1426 my $sum = Math::BigInt->bzero;
279 139         2733 my $num = Math::BigInt->new($nterms-1)->bfac();
280 139         83750 foreach my $i (0 .. $k) {
281 4942         736035 my $den = Math::BigInt->new($nterms - $i)->bfac * Math::BigInt->new(2*$i)->bfac;
282 4942         5007283 $sum += $num->copy->bdiv($den);
283 4942         2864595 $num->bmul(4 * ($nterms+$i));
284             }
285 139         22686 $sum->bmul($nterms);
286 139         16958 $_Borwein_dk[$k] = $sum;
287             }
288 2         15 Math::BigFloat->accuracy($orig_acc);
289             }
290              
291             sub RiemannZeta {
292 4     4 1 12 my($ix) = @_;
293              
294 4 50       25 my $x = (ref($ix) eq 'Math::BigFloat') ? $ix->copy : Math::BigFloat->new("$ix");
295 4 50       142 $x->accuracy($ix->accuracy) if $ix->accuracy;
296 4   33     65 my $xdigits = $ix->accuracy() || Math::BigFloat->accuracy() || Math::BigFloat->div_scale();
297              
298 4 50 66     189 if ($x == int($x) && $xdigits <= 44 && (int($x)-2) <= $#_Riemann_Zeta_Table) {
      66        
299 2         1872 my $izeta = $_Riemann_Zeta_Table[int($x)-2]->copy;
300 2         1161 $izeta->bround($xdigits);
301 2         585 return $izeta;
302             }
303              
304             # Note, this code likely will not work correctly without fixes for RTs:
305             #
306             # 43692 : blog and others broken
307             # 43460 : exp and powers broken
308             #
309             # E.g:
310             # my $n = Math::BigFloat->new(11); $n->accuracy(64); say $n**1.1; # 13.98
311             # my $n = Math::BigFloat->new(11); $n->accuracy(67); say $n**1.1; # 29.98
312             #
313             # There is a hack that tries to work around some of the problem, but it
314             # can't cover everything and it slows things down a lot. There just isn't
315             # any way to do this if the basic math operations don't work right.
316              
317 2         727 my $orig_acc = Math::BigFloat->accuracy();
318 2         25 my $extra_acc = 5;
319 2 100 66     6 if ($x > 15 && $x <= 50) { $extra_acc = 15; }
  1         737  
320              
321 2         219 $xdigits += $extra_acc;
322 2         13 Math::BigFloat->accuracy($xdigits);
323 2         64 $x->accuracy($xdigits);
324 2         142 my $zero= $x->copy->bzero;
325 2         131 my $one = $zero->copy->binc;
326 2         234 my $two = $one->copy->binc;
327              
328 2         198 my $tol = ref($x)->new('0.' . '0' x ($xdigits-1) . '1');
329              
330             # Note: with bignum on, $d1->bpow($one-$x) doesn't change d1 !
331              
332             # This is a hack to turn 6^-40.5 into (6^-(40.5/4))^4. It helps work around
333             # the two RTs listed earlier, though does not completely fix their bugs.
334             # It has the downside of making integer arguments very slow.
335              
336 2         484 my $superx = Math::BigInt->bone;
337 2         91 my $subx = $x->copy;
338 2         49 my $intx = int("$x");
339 2 50 33     133 if ($Math::BigFloat::VERSION < 1.9996 || $x != $intx) {
340 2         882 while ($subx > 1) {
341 8         5638 $superx->blsft(1);
342 8         1200 $subx /= $two;
343             }
344             }
345              
346 2 0 33     1508 if (1 && $x == $intx && $x >= 2 && !($intx & 1) && $intx < 100) {
      33        
      0        
347             # Mathworld equation 63. How fast this is relative to the others is
348             # dependent on the backend library and if we have MPUGMP.
349 0         0 $x = int("$x");
350 0         0 my $den = Math::Prime::Util::factorial($x);
351 0         0 $xdigits -= $extra_acc;
352 0         0 $extra_acc += length($den);
353 0         0 $xdigits += $extra_acc;
354 0         0 $one->accuracy($xdigits); $two->accuracy($xdigits);
  0         0  
355 0         0 Math::BigFloat->accuracy($xdigits);
356 0         0 $subx->accuracy($xdigits); $superx->accuracy($xdigits);
  0         0  
357 0         0 my $Pix = Math::Prime::Util::Pi($xdigits)->bpow($subx)->bpow($superx);
358 0 0       0 my $Bn = Math::Prime::Util::bernreal($x,$xdigits); $Bn = -$Bn if $Bn < 0;
  0         0  
359 0         0 my $twox1 = $two->copy->bpow($x-1);
360             #my $num = $Pix * $Bn * $twox1;
361             #my $res = $num->bdiv($den)->bdec->bround($xdigits - $extra_acc);
362 0         0 my $res = $Bn->bdiv($den)->bmul($Pix)->bmul($twox1)->bdec
363             ->bround($xdigits-$extra_acc);
364 0         0 Math::BigFloat->accuracy($orig_acc);
365 0         0 return $res;
366             }
367              
368             # Go with the basic formula for large x.
369 2 50       798 if (1 && $x >= 50) {
370 0         0 my $negsubx = $subx->copy->bneg;
371 0         0 my $sum = $zero->copy;
372 0         0 my $k = $two->copy->binc;
373 0         0 while ($k->binc <= 1000) {
374 0         0 my $term = $k->copy->bpow($negsubx)->bpow($superx);
375 0         0 $sum += $term;
376 0 0       0 last if $term < ($sum*$tol);
377             }
378 0         0 $k = $two+$two;
379 0         0 $k->bdec(); $sum += $k->copy->bpow($negsubx)->bpow($superx);
  0         0  
380 0         0 $k->bdec(); $sum += $k->copy->bpow($negsubx)->bpow($superx);
  0         0  
381 0         0 $sum->bround($xdigits-$extra_acc);
382 0         0 Math::BigFloat->accuracy($orig_acc);
383 0         0 return $sum;
384             }
385              
386             {
387 2         837 my $dig = int($_Borwein_n / 1.3)+1;
  2         7  
388 2 50       12 _Recompute_Dk( int($xdigits * 1.3) + 4 ) if $dig < $xdigits;
389             }
390              
391 2 50       63 if (ref $_Borwein_dk[0] ne 'Math::BigInt') {
392 0         0 @_Borwein_dk = map { Math::BigInt->new("$_") } @_Borwein_dk;
  0         0  
393             }
394              
395 2         5 my $n = $_Borwein_n;
396              
397 2         14 my $d1 = $two ** ($one - $x);
398 2         419415 my $divisor = ($one - $d1) * $_Borwein_dk[$n];
399 2         3004 $divisor->bneg;
400 2         32 $tol = ($divisor * $tol)->babs();
401              
402 2         399 my ($sum, $bigk) = ($zero->copy, $one->copy);
403 2         89 my $negsubx = $subx->copy->bneg;
404 2         66 foreach my $k (1 .. $n-1) {
405 135         25639 my $den = $bigk->binc()->copy->bpow($negsubx)->bpow($superx);
406 135 100       32696007 my $term = ($k % 2) ? ($_Borwein_dk[$n] - $_Borwein_dk[$k])
407             : ($_Borwein_dk[$k] - $_Borwein_dk[$n]);
408 135 50       23855 $term = Math::BigFloat->new($term) unless ref($term) eq 'Math::BigFloat';
409 135         16866 $sum += $term * $den;
410 135 50       172891 last if $term->copy->babs() < $tol;
411             }
412 2         375 $sum += $_Borwein_dk[0] - $_Borwein_dk[$n];
413 2         1910 $sum = $sum->bdiv($divisor);
414 2         1984 $sum->bdec->bround($xdigits-$extra_acc);
415 2         1305 Math::BigFloat->accuracy($orig_acc);
416 2         120 return $sum;
417             }
418              
419             # Riemann R function
420             sub RiemannR {
421 4     4 1 13 my($x) = @_;
422              
423 4 50       20 if (ref($x) eq 'Math::BigInt') {
424 0         0 my $xacc = $x->accuracy();
425 0         0 $x = Math::BigFloat->new($x);
426 0 0       0 $x->accuracy($xacc) if $xacc;
427             }
428 4 50       107 $x = Math::BigFloat->new("$x") if ref($x) ne 'Math::BigFloat';
429 4   33     23 my $xdigits = $x->accuracy || Math::BigFloat->accuracy() || Math::BigFloat->div_scale();
430 4         60 my $extra_acc = 2;
431 4         18 $xdigits += $extra_acc;
432 4         11 my $orig_acc = Math::BigFloat->accuracy();
433 4         91 Math::BigFloat->accuracy($xdigits);
434 4         81 $x->accuracy($xdigits);
435 4         182 my $tol = $x->copy->bone->brsft($xdigits-1, 10);
436 4         5547 my $sum = $x->copy->bone;
437              
438 4 50       532 if ($xdigits <= length($x->copy->as_int->bstr())) {
439              
440 0         0 for my $k (1 .. 1000) {
441 0         0 my $mob = Math::Prime::Util::moebius($k);
442 0 0       0 next if $mob == 0;
443 0         0 $mob = Math::BigFloat->new($mob);
444 0         0 my $term = $mob->bdiv($k) *
445             Math::Prime::Util::LogarithmicIntegral($x->copy->broot($k));
446 0         0 $sum += $term;
447             #warn "k = $k term = $term sum = $sum\n";
448 0 0       0 last if abs($term) < ($tol * abs($sum));
449             }
450              
451             } else {
452              
453 4         614 my ($flogx, $part_term, $fone, $bigk)
454             = (log($x), Math::BigFloat->bone, Math::BigFloat->bone, Math::BigInt->bone);
455              
456 4 100       210293 if ($_Riemann_Zeta_premult_accuracy < $xdigits) {
457 3         349 @_Riemann_Zeta_Premult = ();
458 3         10 $_Riemann_Zeta_premult_accuracy = $xdigits;
459             }
460              
461 4         17 for my $k (1 .. 10000) {
462 571         22440 my $zeta_term = $_Riemann_Zeta_Premult[$k-1];
463 571 100       1298 if (!defined $zeta_term) {
464 422 50       1204 my $zeta = ($xdigits > 44) ? undef : $_Riemann_Zeta_Table[$k-1];
465 422 50       1030 if (!defined $zeta) {
466 0         0 my $kz = $fone->copy->badd($bigk); # kz is k+1
467 0 0 0     0 if (($k+1) >= 100 && $xdigits <= 40) {
468             # For this accuracy level, two terms are more than enough. Also,
469             # we should be able to miss the Math::BigFloat accuracy bug. If we
470             # try to do this for higher accuracy, things will go very bad.
471 0         0 $zeta = Math::BigFloat->new(3)->bpow(-$kz)
472             + Math::BigFloat->new(2)->bpow(-$kz);
473             } else {
474 0         0 $zeta = Math::Prime::Util::ZetaBigFloat::RiemannZeta( $kz );
475             }
476             }
477 422         1107 $zeta_term = $fone / ($zeta * $bigk + $bigk);
478 422 50       829348 $_Riemann_Zeta_Premult[$k-1] = $zeta_term if defined $_Riemann_Zeta_Table[$k-1];
479             }
480 571         1678 $part_term *= $flogx / $bigk;
481 571         718983 my $term = $part_term * $zeta_term;
482 571         285055 $sum += $term;
483             #warn "k = $k term = $term sum = $sum\n";
484 571 100       277048 last if $term < ($tol*$sum);
485 567         167857 $bigk->binc;
486             }
487              
488             }
489 4         1542 $sum->bround($xdigits-$extra_acc);
490 4         1012 Math::BigFloat->accuracy($orig_acc);
491 4         110 return $sum;
492             }
493              
494             #Math::BigFloat->accuracy($_oldacc);
495             #undef $_oldacc;
496              
497             1;
498              
499             __END__