// Magma program to compute terms of OEIS A322557
//   using the (seemingly nearly certain to give all correct terms) formula
//
//        a(n) = ceiling(z - 1/2 - 1/(12*z))
//               where z = 6/(Pi^2 - (floor(Pi*10^n)/10^n)^2) 
//
//   and confirm their correctness using
//      series approximations for all terms computed
//      and also using exact sums for the first few terms.
//
// Jon E. Schoenfield, Aug 31 2019

DP:=200; R:=RealField(DP); SetDefaultRealField(R);
pi:=Pi(R); L:=pi*pi/6;

numer:=[1, -1,  1, -1,  5, -691, 7, -3617, 43867, -174611, 854513, -236364091]; // from A000367
denom:=[6, 30, 42, 30, 66, 2730, 6,   510,   798,     330,    138,       2730]; // from A002445

a:=[];
nMaxComputeExactSum:=4;
for n in [0..19] do
   "n =", n;
   x := Floor(pi*10^n)/10^n;
   z := 6/(pi^2 - x^2);
   aTest := Ceiling(z - 1/2 - 1/(12*z));
   if n le nMaxComputeExactSum then
      exactSum1:=0.0;
      for k in [1..aTest-1] do
         exactSum1+:=6.0/k^2;
      end for;
      exactSum2:=exactSum1+6.0/aTest^2;
   end if;
   for k in [aTest-1, aTest] do
      approxSum := L - 1/k + (1/2)/k^2;
      for j in [1..#numer] do
         approxSum -:= (numer[j]/denom[j])/k^(2*j+1);
      end for;
      k, ChangePrecision(Sqrt(6*approxSum), 42), "(series approx)";
      if n le nMaxComputeExactSum then
         if k eq aTest-1 then
            k, ChangePrecision(Sqrt(exactSum1),42), "(exact sum)";
         else
            k, ChangePrecision(Sqrt(exactSum2),42), "(exact sum)";
         end if;
      end if;
   end for;
   "";
end for;