#!/usr/bin/env perl
use strict;
use warnings;

use Test::More;
use Math::Prime::Util::GMP qw/bernfrac bernreal stirling harmfrac harmreal/;
my $extra = defined $ENV{EXTENDED_TESTING} && $ENV{EXTENDED_TESTING};

my @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/);
my @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/);

# 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"))
my @stirling2 = (
[qw/1 0/],
[qw/0 1 0/],
[qw/0 1 1 0/],
[qw/0 1 3 1 0/],
[qw/0 1 7 6 1 0/],
[qw/0 1 15 25 10 1 0/],
[qw/0 1 31 90 65 15 1 0/],
[qw/0 1 63 301 350 140 21 1 0/],
[qw/0 1 127 966 1701 1050 266 28 1 0/],
[qw/0 1 255 3025 7770 6951 2646 462 36 1 0/],
[qw/0 1 511 9330 34105 42525 22827 5880 750 45 1 0/],
[qw/0 1 1023 28501 145750 246730 179487 63987 11880 1155 55 1 0/],
[qw/0 1 2047 86526 611501 1379400 1323652 627396 159027 22275 1705 66 1 0/],
[qw/0 1 4095 261625 2532530 7508501 9321312 5715424 1899612 359502 39325 2431 78 1 0/],
[qw/0 1 8191 788970 10391745 40075035 63436373 49329280 20912320 5135130 752752 66066 3367 91 1 0/],
[qw/0 1 16383 2375101 42355950 210766920 420693273 408741333 216627840 67128490 12662650 1479478 106470 4550 105 1 0/],
[qw/0 1 32767 7141686 171798901 1096190550 2734926558 3281882604 2141764053 820784250 193754990 28936908 2757118 165620 6020 120 1 0/],
[qw/0 1 65535 21457825 694337290 5652751651 17505749898 25708104786 20415995028 9528822303 2758334150 512060978 62022324 4910178 249900 7820 136 1 0/],
[qw/0 1 131071 64439010 2798806985 28958095545 110687251039 197462483400 189036065010 106175395755 37112163803 8391004908 1256328866 125854638 8408778 367200 9996 153 1 0/],
[qw/0 1 262143 193448101 11259666950 147589284710 693081601779 1492924634839 1709751003480 1144614626805 477297033785 129413217791 23466951300 2892439160 243577530 13916778 527136 12597 171 1 0/],
[qw/0 1 524287 580606446 45232115901 749206090500 4306078895384 11143554045652 15170932662679 12011282644725 5917584964655 1900842429486 411016633391 61068660380 6302524580 452329200 22350954 741285 15675 190 1 0/],
);

my @stirling1 = (
[qw/1 0/],
[qw/0 1 0/],
[qw/0 -1 1 0/],
[qw/0 2 -3 1 0/],
[qw/0 -6 11 -6 1 0/],
[qw/0 24 -50 35 -10 1 0/],
[qw/0 -120 274 -225 85 -15 1 0/],
[qw/0 720 -1764 1624 -735 175 -21 1 0/],
[qw/0 -5040 13068 -13132 6769 -1960 322 -28 1 0/],
[qw/0 40320 -109584 118124 -67284 22449 -4536 546 -36 1 0/],
[qw/0 -362880 1026576 -1172700 723680 -269325 63273 -9450 870 -45 1 0/],
[qw/0 3628800 -10628640 12753576 -8409500 3416930 -902055 157773 -18150 1320 -55 1 0/],
[qw/0 -39916800 120543840 -150917976 105258076 -45995730 13339535 -2637558 357423 -32670 1925 -66 1 0/],
[qw/0 479001600 -1486442880 1931559552 -1414014888 657206836 -206070150 44990231 -6926634 749463 -55770 2717 -78 1 0/],
[qw/0 -6227020800 19802759040 -26596717056 20313753096 -9957703756 3336118786 -790943153 135036473 -16669653 1474473 -91091 3731 -91 1 0/],
[qw/0 87178291200 -283465647360 392156797824 -310989260400 159721605680 -56663366760 14409322928 -2681453775 368411615 -37312275 2749747 -143325 5005 -105 1 0/],
[qw/0 -1307674368000 4339163001600 -6165817614720 5056995703824 -2706813345600 1009672107080 -272803210680 54631129553 -8207628000 928095740 -78558480 4899622 -218400 6580 -120 1 0/],
[qw/0 20922789888000 -70734282393600 102992244837120 -87077748875904 48366009233424 -18861567058880 5374523477960 -1146901283528 185953177553 -23057159840 2185031420 -156952432 8394022 -323680 8500 -136 1 0/],
[qw/0 -355687428096000 1223405590579200 -1821602444624640 1583313975727488 -909299905844112 369012649234384 -110228466184200 24871845297936 -4308105301929 577924894833 -60202693980 4853222764 -299650806 13896582 -468180 10812 -153 1 0/],
[qw/0 6402373705728000 -22376988058521600 34012249593822720 -30321254007719424 17950712280921504 -7551527592063024 2353125040549984 -557921681547048 102417740732658 -14710753408923 1661573386473 -147560703732 10246937272 -549789282 22323822 -662796 13566 -171 1 0/],
[qw/0 -121645100408832000 431565146817638400 -668609730341153280 610116075740491776 -371384787345228000 161429736530118960 -52260903362512720 12953636989943896 -2503858755467550 381922055502195 -46280647751910 4465226757381 -342252511900 20692933630 -973941900 34916946 -920550 16815 -190 1 0/],
);

plan tests => 2 + scalar(@stirling2) + scalar(@stirling1) + 2 + 2+6 + 3;

{
  my @num = map { (bernfrac(2*$_))[0] }  0 .. $#A000367;
  my @den = map { (bernfrac(2*$_))[1] }  0 .. $#A002445;
  is_deeply( \@num, \@A000367, "B_2n numerators 0 .. $#A000367" );
  is_deeply( \@den, \@A002445, "B_2n denominators 0 .. $#A002445" );
}

{
  my $n = 0;
  foreach my $narr (@stirling2) {
    my @s2 = map { stirling($n,$_,2) } 0..$n+1;
    is_deeply( \@s2, $narr, "Stirling 2: S($n,0..". ($n+1) .")" );
    $n++;
  }
}
{
  my $n = 0;
  foreach my $narr (@stirling1) {
    my @s1 = map { stirling($n,$_,1) } 0..$n+1;
    is_deeply( \@s1, $narr, "Stirling 1: s($n,0..". ($n+1) .")" );
    $n++;
  }
}

# Random large values
is( stirling(137,14,2), '119921091552849030298712472784655311636455680166733181212686613915193327160253320582242126677455942135033082921631146726468348015037515052954267400', "S(137,14)" );
is( stirling(99,14,1), '-76185801962487294910690331431878395972434237854033124053130281967496159110075633618163409236750708320666097319441407130017141175404355750702612480000000', "s(99,14)" );

###########
is_deeply( [harmfrac(27)], [qw/312536252003 80313433200/], "harmfrac(27)");
is_deeply( [harmfrac(172)], [qw/79501058066082280981403896527589637839225193778798494740482914683623969349 13880309301456766718169645534485312116208863916210197754275864148562288000/], "harmfrac(172)");

is( harmreal(5,5), "2.28333", "harmreal(5,5)" );
is( harmreal(15,2), "3.32", "harmreal(15,2)" );
is( harmreal(15,24), "3.318228993228993228993229", "harmreal(15,24)" );
is( harmreal(1500,84), "7.890769348288132237024982466533383516587334672793034866299830645277849644355893942483", "harmreal(1500,84)" );
is( harmreal(2502,763), "8.4022611827442616317889358234638600819223549876775756009040579007028304622396704717826054849631311959827336370652350363210214348075460454127806820991887787301726524318112301788039510819911216154119126624804602690267244260972035688366573000214774520753765479437246871327307674220381391918202939677188395567028236867891536047084719881902222764381394522556877103822410093096790602691889769666625414585123756378881628875119610391914967976757016187471943182230789317973147668458937306028722119815966094143929327181951971889166127652438346064969183176053197768023014774945977838767778096364258279792884478027384432419591251420003615953125404288714656235208556151078286913808883501164776686214636296222092373622320260770103224594844688053478029558802715467037435761134100", "harmreal(2502,763)" );
is( harmreal(2502,764), "8.40226118274426163178893582346386008192235498767757560090405790070283046223967047178260548496313119598273363706523503632102143480754604541278068209918877873017265243181123017880395108199112161541191266248046026902672442609720356883665730002147745207537654794372468713273076742203813919182029396771883955670282368678915360470847198819022227643813945225568771038224100930967906026918897696666254145851237563788816288751196103919149679767570161874719431822307893179731476684589373060287221198159660941439293271819519718891661276524383460649691831760531977680230147749459778387677780963642582797928844780273844324195912514200036159531254042887146562352085561510782869138088835011647766862146362962220923736223202607701032245948446880534780295588027154670374357611340997", "harmreal(2502,764)" );

###########
is( bernreal(24), "-86580.2531135531135531135531135531135531135531", "bern(24)" );
is( bernreal(16,5), "-7.09216", "bern(16,5)" );
is( bernreal(200,7), "-364707726451913543621383088655499449048682346861910587376827309036218502061077098613909538979906038749680879817942730412840299537621084437107198005516180654930958621012134647126652378010626301508404512970954087034967.1585337", "bern(200,7)" );