// Magma program to generate b-file for OEIS sequence A227777
//
// Program converts the n-th and (n+1)st partial sums described at
//      https://oeis.org/A227777
// i.e.,
//       n-th_partial_sum = 1/0! + 1/1! + 1/2! + ... + 1/n!
// and
//    (n+1)st_partial_sum = 1/0! + 1/1! + 1/2! + ... + 1/(n+1)!
//
// to continued fractions, compares them term-by-term to build the
// continued fraction representation of the fraction c/d (where c
// and d are both integers) that minimizes the denominator d
// subject to the constraint that
//
//      n-th_partial_sum <= c/d < (n+1)st_partial_sum
//
// and, for each index from n=0 to n=nMax outputs the index n and
// the corresponding denominator d, i.e., A227777(n).
//
// (The logic could probably be simplified considerably, but the
// program generates hundreds of terms very quickly as it is.)
//
// Jon E. Schoenfield, 27 Jun 2015

nMax:=200;

PartialSumNth:=1;
CFCurr:=ContinuedFraction(PartialSumNth);
NPlus1Factorial:=1;

for n in [0..nMax] do

   NPlus1Factorial*:=(n+1);
   PartialSumNthplus1st:=PartialSumNth+1/NPlus1Factorial;
   CFNext:=ContinuedFraction(PartialSumNthplus1st);
   CFSoln:=[];

   if #CFCurr eq #CFNext then // CFs have same number of terms

      for j in [1..#CFCurr] do
         if CFCurr[j] eq CFNext[j] then
            CFSoln[j]:=CFCurr[j];
         else // 1st difference found
            if j eq #CFCurr then // 1st diff is at last term
               CFSoln[j]:=CFCurr[j];
            else // 1st diff is before last term
               CFSoln[j]:=Min(CFCurr[j],CFNext[j])+1;
               break;
            end if;
         end if;
      end for;

   elif #CFCurr lt #CFNext then // CFCurr has fewer terms

      for j in [1..#CFCurr] do
         if CFCurr[j] eq CFNext[j] then
            CFSoln[j]:=CFCurr[j];
         else // 1st difference found
            if j eq #CFCurr then // 1st diff is at last term of CFCurr
               if CFCurr[j] lt CFNext[j] then
                  CFSoln[j]:=CFCurr[j];
               else // CFCurr[j] > CFNext[j]
                  CFSoln[j]:=CFNext[j]+1;
               end if;
            else // 1st diff is before last term of CFCurr
               CFSoln[j]:=Min(CFCurr[j],CFNext[j])+1;
               break;
            end if;
         end if;
      end for;

   else // CFNext has fewer terms

      j1stDiff:=0; // init

      for j in [1..#CFNext] do
         if CFCurr[j] eq CFNext[j] then
            CFSoln[j]:=CFCurr[j];
         else
            j1stDiff:=j;
            break;
         end if;
      end for;

      if j1stDiff eq 0 then // no diffs within the CFs' common length

         if #CFCurr eq #CFNext+1 then // CFCurr has just 1 more term
            CFSoln[#CFNext+1]:=CFCurr[#CFNext+1];
         else // CFCurr has more than 1 more term
            CFSoln[#CFNext+1]:=CFCurr[#CFNext+1]+1;
         end if;

      elif j1stDiff eq #CFNext then // 1st diff is at end of CFNext

         if CFCurr[j1stDiff] eq CFNext[j1stDiff]-1 then
            CFSoln[j1stDiff]:=CFCurr[j1stDiff];
            CFSoln[j1stDiff+1]:=2;
         else
            CFSoln[j1stDiff]:=Min(CFCurr[j1stDiff],CFNext[j1stDiff])+1;
         end if;

      else // 1st diff is at j1stDiff, somewhere before end of CFNext

         CFSoln[j1stDiff]:=Min(CFCurr[j1stDiff],CFNext[j1stDiff])+1;

      end if;

   end if;

   // convert the continued fraction to the rational fraction c/d

   frac:=CFSoln[#CFSoln];
   for j in [#CFSoln-1..1 by -1] do
      frac:=CFSoln[j]+1/frac;
   end for;

   n, Denominator(frac);

   PartialSumNth:=PartialSumNthplus1st;
   CFCurr:=CFNext;

end for;