# When a design develops complications that are the result of a mental habit
# that is no longer adaptive, the developer has succumbed to Whorfian mind-lock.
# -- Jargon File
# This will probably take 20 minutes or so the first time you run it.
use strict;
use warnings;
use Math::GMP;

my $digits = 1_000_000;

my $A = Math::GMP->new(13591409);
my $B = Math::GMP->new(545140134);
my $C = Math::GMP->new(640320);
my $D = $C ** 3 / 24;

my $total_terms;
my $done_terms = 0;

my $pifile = 'bigpi';
my $pi;
if (!-f $pifile || -s _ < $digits) {
   open my $F, '>', $pifile or die "Can't write $pifile: $!\n";
   $| = 1;
   print "Calculating pi... this may take a while\n";
   # Calculate pi using Chudnovsky algorithm
   $total_terms = int($digits / 14.1816474627254776555);
   my $c = ($C * Math::GMP->new(10) ** ($digits * 2 - 2))->bsqrt();
   my ($p, $q) = partial_sum(0, $total_terms);
   print "\rAlmost...";
   $pi = $p * ($C / 12) * $c / ($q + $A * $p);
   $pi = "$pi";
   print "\rFinished!\n";
   print $F $pi;
   close $F;
}
else {
   open my $F, '<', $pifile or die "Can't read $pifile: $!\n";
   read $F, $pi, $digits;
   close $F;
}

my $yr = 2015;
my @month_days = ( 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 );
my %days;
foreach my $mon (1 .. 12) {
   foreach my $day (1 .. $month_days[$mon-1]) {
      $days{sprintf('%02d%02d', $mon, $day)} = 1;
   }
}

my $i = index($pi, $yr);
while ($i >= 0) {
   delete $days{substr($pi, $i + 4, 4)};
   $i = index($pi, $yr, $i + 4);
}

print $yr, $_, "\n" for sort keys %days;

sub partial_sum {
   my ($a, $b) = @_;
   if ($b <= $a) {
      die;
   }
   elsif ($b == $a + 1) {
      if (++$done_terms % 10000 == 0) {
         printf "\r%.2f%%", 100 * $done_terms / $total_terms;
      }
      my $p = Math::GMP->new($b) * $b * $b * $D;
      my $g = Math::GMP->new(2 * $b - 1) * (6 * $b - 5) * (6 * $b - 1);
      my $q = $g * ($A + $B * $b);
      $q = -$q if $b & 1;
      return ($p, $q, $g);
   }
   else {
      my $mid = ($a + $b) >> 1;
      my ($p1, $q1, $g1) = partial_sum($a, $mid);
      my ($p2, $q2, $g2) = partial_sum($mid, $b);
      my $gcd = $g1->bgcd($p2);
      if ($gcd > 1) { $g1 /= $gcd; $p2 /= $gcd; }
      return ($p1 * $p2, $q1 * $p2 + $q2 * $g1, $g1 * $g2);
   }
}
