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

use Test::More;
use Math::Prime::Util qw/bernfrac bernreal harmfrac harmreal stirling sumdigits/;
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/);
my @A001008 = (qw/1 3 11 25 137 49 363 761 7129 7381 83711 86021 1145993 1171733 1195757 2436559 42142223 14274301 275295799 55835135 18858053 19093197 444316699 1347822955 34052522467 34395742267 312536252003 315404588903 9227046511387/);
my @A002805 = (qw/1 2 6 12 60 20 140 280 2520 2520 27720 27720 360360 360360 360360 720720 12252240 4084080 77597520 15519504 5173168 5173168 118982864 356948592 8923714800 8923714800 80313433200 80313433200 2329089562800/);
if (!$extra) {
  $#A000367 = 20;
  $#A002445 = 20;
  $#A001008 = 20;
  $#A002805 = 20;
}

my @bern_reals = (1,1/2,1/6,0,-1/30,0,1/42,0,-1/30,0,5/66,0,-691/2730,0,7/6,0,-3617/510,0,43867/798,0,-174611/330,0,854513/138,0,-236364091/2730);
my @harm_reals = (0/1,1/1,3/2,11/6,50/24,274/120,1764/720,13068/5040,109584/40320,1026576/362880,10628640/3628800,120543840/39916800,1486442880/479001600,19802759040/6227020800,283465647360/87178291200,4339163001600/1307674368000,70734282393600/20922789888000,1223405590579200/355687428096000,22376988058521600/6402373705728000,431565146817638400/121645100408832000,8752948036761600000/2432902008176640000);

# 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/],
);

if (!$extra) {
  $#stirling1 = 12;
  $#stirling2 = 12;
}

plan tests => 2*$extra
            + 2 + scalar(@bern_reals)
            + 2 + scalar(@harm_reals)
            + 2 + scalar(@stirling2) + scalar(@stirling1) + 2*$extra
            + 2*$extra;


if ($extra) {
  like( bernreal(46), qr/2115074863808199160560.145/, "bernreal(46)" );
  like( harmreal(46), qr/4.416687245986104750714329/, "harmreal(46)" );
}

{
  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" );
}
for my $n (0 .. $#bern_reals) {
  cmp_closeto( bernreal($n), $bern_reals[$n], 1e-8, "bernreal($n)" );
}

{
  my @num = map { (harmfrac(1+$_))[0] }  0 .. $#A001008;
  my @den = map { (harmfrac(1+$_))[1] }  0 .. $#A002805;
  is_deeply( \@num, \@A001008, "H_n numerators 0 .. $#A001008" );
  is_deeply( \@den, \@A002805, "H_n denominators 0 .. $#A002805" );
}
for my $n (0 .. $#harm_reals) {
  cmp_closeto( harmreal($n), $harm_reals[$n], 1e-8, "harmreal($n)" );
}


sub cmp_closeto {
  my $got = shift;
  my $expect = shift;
  my $tolerance = shift;
  my $message = shift;
  cmp_ok( abs($got - $expect), '<=', $tolerance, $message );
}


# These are compiler errors:
#ok(!eval { stirling() }, "stirling with no args");
#ok(!eval { stirling(4) }, "stirling with one arg");
#ok(!eval { stirling(1,2,3,4) }, "stirling with four args");
ok(!eval { stirling(-4, -3) }, "Expected fail: stirling with negative args");
ok(!eval { stirling(4,3,4) }, "Expected fail: stirling type 4");
# These just generate warnings due to our prototype
#ok(!eval { stirling(undef, undef) }, "stirling given undefs");
#ok(!eval { stirling("x","y") }, "stirling x y");

{
  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++;
  }
}

if ($extra) {
  # Random large values
  is( stirling(114,85,2), '722095587897382907118640452680242028195738761915144254970925658656935934040', "S(114,85)" );
  is( stirling(132,67,1), '-6132458966070920781607687809239433538883836871765225500351514785120957322534135782514155513931693375104995311496306605620444680401484569675682191339176710', "s(132,67)" );
}

if ($extra) {
  is(sumdigits( (bernfrac(502))[0], 157), 27893, "sumdigits(bernfrac(502) numerator) base 157");
  is(sumdigits(stirling(234,95)), 1485, "sumdigits(stirling(234,95))");
}