1#!/usr/bin/env perl
2use strict;
3use warnings;
4
5use Test::More;
6use Math::Prime::Util::GMP qw/bernfrac bernreal stirling harmfrac harmreal/;
7my $extra = defined $ENV{EXTENDED_TESTING} && $ENV{EXTENDED_TESTING};
8
9my @A000367 = (qw/1 1 -1 1 -1 5 -691 7 -3617 43867 -174611 854513 -236364091 8553103 -23749461029 8615841276005 -7709321041217 2577687858367 -26315271553053477373 2929993913841559 -261082718496449122051 1520097643918070802691 -27833269579301024235023 596451111593912163277961 -5609403368997817686249127547 495057205241079648212477525 -801165718135489957347924991853 29149963634884862421418123812691 -2479392929313226753685415739663229 84483613348880041862046775994036021 -1215233140483755572040304994079820246041491/);
10my @A002445 = (qw/1 6 30 42 30 66 2730 6 510 798 330 138 2730 6 870 14322 510 6 1919190 6 13530 1806 690 282 46410 66 1590 798 870 354 56786730 6 510 64722 30 4686 140100870 6 30 3318 230010 498 3404310 6 61410 272118 1410 6 4501770 6 33330 4326 1590 642 209191710 1518 1671270 42/);
11
12# Generated by gp 2.8.0: for(n=0,20,printf("[qw/");for(m=0,n+1,printf("%d ",stirling(n,m,2)));printf("/],\n"))
13my @stirling2 = (
14[qw/1 0/],
15[qw/0 1 0/],
16[qw/0 1 1 0/],
17[qw/0 1 3 1 0/],
18[qw/0 1 7 6 1 0/],
19[qw/0 1 15 25 10 1 0/],
20[qw/0 1 31 90 65 15 1 0/],
21[qw/0 1 63 301 350 140 21 1 0/],
22[qw/0 1 127 966 1701 1050 266 28 1 0/],
23[qw/0 1 255 3025 7770 6951 2646 462 36 1 0/],
24[qw/0 1 511 9330 34105 42525 22827 5880 750 45 1 0/],
25[qw/0 1 1023 28501 145750 246730 179487 63987 11880 1155 55 1 0/],
26[qw/0 1 2047 86526 611501 1379400 1323652 627396 159027 22275 1705 66 1 0/],
27[qw/0 1 4095 261625 2532530 7508501 9321312 5715424 1899612 359502 39325 2431 78 1 0/],
28[qw/0 1 8191 788970 10391745 40075035 63436373 49329280 20912320 5135130 752752 66066 3367 91 1 0/],
29[qw/0 1 16383 2375101 42355950 210766920 420693273 408741333 216627840 67128490 12662650 1479478 106470 4550 105 1 0/],
30[qw/0 1 32767 7141686 171798901 1096190550 2734926558 3281882604 2141764053 820784250 193754990 28936908 2757118 165620 6020 120 1 0/],
31[qw/0 1 65535 21457825 694337290 5652751651 17505749898 25708104786 20415995028 9528822303 2758334150 512060978 62022324 4910178 249900 7820 136 1 0/],
32[qw/0 1 131071 64439010 2798806985 28958095545 110687251039 197462483400 189036065010 106175395755 37112163803 8391004908 1256328866 125854638 8408778 367200 9996 153 1 0/],
33[qw/0 1 262143 193448101 11259666950 147589284710 693081601779 1492924634839 1709751003480 1144614626805 477297033785 129413217791 23466951300 2892439160 243577530 13916778 527136 12597 171 1 0/],
34[qw/0 1 524287 580606446 45232115901 749206090500 4306078895384 11143554045652 15170932662679 12011282644725 5917584964655 1900842429486 411016633391 61068660380 6302524580 452329200 22350954 741285 15675 190 1 0/],
35);
36# Generated by gp 2.8.0: lah(n,k)={if(k>n,0,if(k==n,1,binomial(n-1,k-1)*n!/k!))}; for(n=0,20,printf("[qw/");for(m=0,n+1,printf("%d ",lah(n,m)));printf("/],\n"))
37my @stirling3 = (
38[qw/1 0 /],
39[qw/0 1 0 /],
40[qw/0 2 1 0 /],
41[qw/0 6 6 1 0 /],
42[qw/0 24 36 12 1 0 /],
43[qw/0 120 240 120 20 1 0 /],
44[qw/0 720 1800 1200 300 30 1 0 /],
45[qw/0 5040 15120 12600 4200 630 42 1 0 /],
46[qw/0 40320 141120 141120 58800 11760 1176 56 1 0 /],
47[qw/0 362880 1451520 1693440 846720 211680 28224 2016 72 1 0 /],
48[qw/0 3628800 16329600 21772800 12700800 3810240 635040 60480 3240 90 1 0 /],
49[qw/0 39916800 199584000 299376000 199584000 69854400 13970880 1663200 118800 4950 110 1 0 /],
50[qw/0 479001600 2634508800 4390848000 3293136000 1317254400 307359360 43908480 3920400 217800 7260 132 1 0 /],
51[qw/0 6227020800 37362124800 68497228800 57081024000 25686460800 6849722880 1141620480 122316480 8494200 377520 10296 156 1 0 /],
52[qw/0 87178291200 566658892800 1133317785600 1038874636800 519437318400 155831195520 29682132480 3710266560 309188880 17177160 624624 14196 182 1 0 /],
53[qw/0 1307674368000 9153720576000 19833061248000 19833061248000 10908183686400 3636061228800 779155977600 111307996800 10821610800 721440720 32792760 993720 19110 210 1 0 /],
54[qw/0 20922789888000 156920924160000 366148823040000 396661224960000 237996734976000 87265469491200 20777492736000 3339239904000 371026656000 28857628800 1574052480 59623200 1528800 25200 240 1 0 /],
55[qw/0 355687428096000 2845499424768000 7113748561920000 8299373322240000 5394592659456000 2157837063782400 565147802419200 100919250432000 12614906304000 1121325004800 71357045760 3243502080 103958400 2284800 32640 272 1 0 /],
56[qw/0 6402373705728000 54420176498688000 145120470663168000 181400588328960000 126980411830272000 55024845126451200 15721384321843200 3088129063219200 428906814336000 42890681433600 3119322286080 165418606080 6362254080 174787200 3329280 41616 306 1 0 /],
57);
58
59my @stirling1 = (
60[qw/1 0/],
61[qw/0 1 0/],
62[qw/0 -1 1 0/],
63[qw/0 2 -3 1 0/],
64[qw/0 -6 11 -6 1 0/],
65[qw/0 24 -50 35 -10 1 0/],
66[qw/0 -120 274 -225 85 -15 1 0/],
67[qw/0 720 -1764 1624 -735 175 -21 1 0/],
68[qw/0 -5040 13068 -13132 6769 -1960 322 -28 1 0/],
69[qw/0 40320 -109584 118124 -67284 22449 -4536 546 -36 1 0/],
70[qw/0 -362880 1026576 -1172700 723680 -269325 63273 -9450 870 -45 1 0/],
71[qw/0 3628800 -10628640 12753576 -8409500 3416930 -902055 157773 -18150 1320 -55 1 0/],
72[qw/0 -39916800 120543840 -150917976 105258076 -45995730 13339535 -2637558 357423 -32670 1925 -66 1 0/],
73[qw/0 479001600 -1486442880 1931559552 -1414014888 657206836 -206070150 44990231 -6926634 749463 -55770 2717 -78 1 0/],
74[qw/0 -6227020800 19802759040 -26596717056 20313753096 -9957703756 3336118786 -790943153 135036473 -16669653 1474473 -91091 3731 -91 1 0/],
75[qw/0 87178291200 -283465647360 392156797824 -310989260400 159721605680 -56663366760 14409322928 -2681453775 368411615 -37312275 2749747 -143325 5005 -105 1 0/],
76[qw/0 -1307674368000 4339163001600 -6165817614720 5056995703824 -2706813345600 1009672107080 -272803210680 54631129553 -8207628000 928095740 -78558480 4899622 -218400 6580 -120 1 0/],
77[qw/0 20922789888000 -70734282393600 102992244837120 -87077748875904 48366009233424 -18861567058880 5374523477960 -1146901283528 185953177553 -23057159840 2185031420 -156952432 8394022 -323680 8500 -136 1 0/],
78[qw/0 -355687428096000 1223405590579200 -1821602444624640 1583313975727488 -909299905844112 369012649234384 -110228466184200 24871845297936 -4308105301929 577924894833 -60202693980 4853222764 -299650806 13896582 -468180 10812 -153 1 0/],
79[qw/0 6402373705728000 -22376988058521600 34012249593822720 -30321254007719424 17950712280921504 -7551527592063024 2353125040549984 -557921681547048 102417740732658 -14710753408923 1661573386473 -147560703732 10246937272 -549789282 22323822 -662796 13566 -171 1 0/],
80[qw/0 -121645100408832000 431565146817638400 -668609730341153280 610116075740491776 -371384787345228000 161429736530118960 -52260903362512720 12953636989943896 -2503858755467550 381922055502195 -46280647751910 4465226757381 -342252511900 20692933630 -973941900 34916946 -920550 16815 -190 1 0/],
81);
82
83plan tests => 2 + scalar(@stirling3) + scalar(@stirling2) + scalar(@stirling1) + 3 + 2+6 + 4;
84
85{
86  my @num = map { (bernfrac(2*$_))[0] }  0 .. $#A000367;
87  my @den = map { (bernfrac(2*$_))[1] }  0 .. $#A002445;
88  is_deeply( \@num, \@A000367, "B_2n numerators 0 .. $#A000367" );
89  is_deeply( \@den, \@A002445, "B_2n denominators 0 .. $#A002445" );
90}
91
92{
93  my $n = 0;
94  foreach my $narr (@stirling3) {
95    my @s2 = map { stirling($n,$_,3) } 0..$n+1;
96    is_deeply( \@s2, $narr, "Stirling 3: L($n,0..". ($n+1) .")" );
97    $n++;
98  }
99}
100{
101  my $n = 0;
102  foreach my $narr (@stirling2) {
103    my @s2 = map { stirling($n,$_,2) } 0..$n+1;
104    is_deeply( \@s2, $narr, "Stirling 2: S($n,0..". ($n+1) .")" );
105    $n++;
106  }
107}
108{
109  my $n = 0;
110  foreach my $narr (@stirling1) {
111    my @s1 = map { stirling($n,$_,1) } 0..$n+1;
112    is_deeply( \@s1, $narr, "Stirling 1: s($n,0..". ($n+1) .")" );
113    $n++;
114  }
115}
116
117# Random large values
118is( stirling(246,61,3), '16781089031289908648739894658187968503831038435579384473036844176021715959960063309748494851984040210112281911986359570651366793437675293554444919512994972542928372289858529733210971633662476717598808875335617639975100691881555975405754262383248739607900482185116908949278639778972088005321758840374848760552543345287066862821699035121991640194236911712299712310508321119592475284255246533432645238910745625755648000000000000000000000000000000000000000000000', "L(246,61)" );
119is( stirling(137,14,2), '119921091552849030298712472784655311636455680166733181212686613915193327160253320582242126677455942135033082921631146726468348015037515052954267400', "S(137,14)" );
120is( stirling(99,14,1), '-76185801962487294910690331431878395972434237854033124053130281967496159110075633618163409236750708320666097319441407130017141175404355750702612480000000', "s(99,14)" );
121
122###########
123is_deeply( [harmfrac(27)], [qw/312536252003 80313433200/], "harmfrac(27)");
124is_deeply( [harmfrac(172)], [qw/79501058066082280981403896527589637839225193778798494740482914683623969349 13880309301456766718169645534485312116208863916210197754275864148562288000/], "harmfrac(172)");
125
126is( harmreal(5,6), "2.28333", "harmreal(5,6)" );
127is( harmreal(15,3), "3.32", "harmreal(15,3)" );
128is( harmreal(15,25), "3.318228993228993228993229", "harmreal(15,25)" );
129is( harmreal(1500,85), "7.890769348288132237024982466533383516587334672793034866299830645277849644355893942483", "harmreal(1500,85)" );
130is( harmreal(2502,764), "8.4022611827442616317889358234638600819223549876775756009040579007028304622396704717826054849631311959827336370652350363210214348075460454127806820991887787301726524318112301788039510819911216154119126624804602690267244260972035688366573000214774520753765479437246871327307674220381391918202939677188395567028236867891536047084719881902222764381394522556877103822410093096790602691889769666625414585123756378881628875119610391914967976757016187471943182230789317973147668458937306028722119815966094143929327181951971889166127652438346064969183176053197768023014774945977838767778096364258279792884478027384432419591251420003615953125404288714656235208556151078286913808883501164776686214636296222092373622320260770103224594844688053478029558802715467037435761134100", "harmreal(2502,764)" );
131is( harmreal(2502,765), "8.40226118274426163178893582346386008192235498767757560090405790070283046223967047178260548496313119598273363706523503632102143480754604541278068209918877873017265243181123017880395108199112161541191266248046026902672442609720356883665730002147745207537654794372468713273076742203813919182029396771883955670282368678915360470847198819022227643813945225568771038224100930967906026918897696666254145851237563788816288751196103919149679767570161874719431822307893179731476684589373060287221198159660941439293271819519718891661276524383460649691831760531977680230147749459778387677780963642582797928844780273844324195912514200036159531254042887146562352085561510782869138088835011647766862146362962220923736223202607701032245948446880534780295588027154670374357611340997", "harmreal(2502,765)" );
132
133###########
134is( bernreal(24), "-86580.25311355311355311355311355311355311", "bern(24)" );
135is( bernreal(16,5), "-7.0922", "bern(16,5)" );
136is( bernreal(200,7), "-3647077E209", "bern(200,7)" );
137is( bernreal(222,260), "1427295874284878567714163200871224998971799130039664528732851045408568226351052228503180915783619081891275551756930326329197981112038254961419750276336559804985081989565349639529368205790276538674423444433770367136781907003735697448421366567512564135.0193252189", "bern(222,260)" );
138