// Magma program to compute terms of A341046 and A341047
// Jon E. Schoenfield
// 14 Mar 2021

// Running on the Magma Calculator at
//
//     http://magma.maths.usyd.edu.au/calc/
//
// the program below finds the first 22 terms (n = 0..21)
// of A341046 in about 0.01 seconds.


nMax:=21; // max index of terms to compute
DP:=50;   // digits of precision
R:=RealField(DP);
SetDefaultRealField(R);

c:=Pi(R); // A341047(n) = floor(c*A341046(n)) where c = Pi

cF:=c-Round(c); // the "signed fractional part" of c:
                //   -0.5 <= cF < 0.5

// Initialize the integer step sizes to use (for increasing x)
//  to yield small increases/decreases in the signed
//  fractional part of y=c*x:

if cF gt 0.0 then
   dxUp:=1; // "upstep" (change to x to increase y-Round(y))
   dxDn:=Ceiling(0.5/cF); // "downstep" (to decrease y-Round(y))
else
   dxDn:=1;
   dxUp:=Ceiling(-0.5/cF);
end if;

// Get sizes of changes to y and its signed fractional part
//  that will result from incrementing x by dxUp and dxDn:

dyUp:=c*dxUp;
dyDn:=c*dxDn;
dyUpF:=dyUp-Round(dyUp);
dyDnF:=dyDn-Round(dyDn);

// Initial pair of terms:
//  x=A341046(0)=1, yI=A341047(0)=floor(Pi).

x:=1;
y:=c*x;
yI:=Floor(y); // not Round(y) here; here, we must truncate
yF:=y-yI; // (unsigned) fractional part: 0 <= yF < 1

b:=true; // x & yI will be stored as terms in A340146 & A340147
n:=0; // offset of A340146 & A340147
X:=[];
YI:=[];

while true do

   if b then // store the value of x and its corresponding yI

      X[#X+1]:=x;
      YI[#YI+1]:=yI;
      if n ge nMax then // term nMax has been stored;
         break;         // exit the loop
      end if;
      n+:=1; // index of next term
      yFlim:=0.1^n; // update upper lim on signed fractional part
      b:=(yF lt yFlim); // true iff curr term is also next term

   else // current yF value is not in [0, 10^-n]

      D:=Ceiling((yFlim-yF)/dyDnF); // D = min # downsteps to get
                                    //  down to the upper limit.
      b:=(yF gt -D*dyDnF); // true iff D downsteps would get
                           //  below upper lim but not below 0
                           //  (i.e., would yield a term)
      if not b then        // D downsteps would overshoot
         D-:=1;            //  (i.e., go below 0), so go 1 fewer.
      end if;

      if D gt 0 then    // If there are any downsteps to take,
         x+:=D*dxDn;    //  then take them (i.e., D downsteps),
         y:=c*x;        //  and update y,
         yI:=Floor(y);  //  the integer part of y,
         yF:=y-yI;      //  and the fractional part of y.
      end if;

      if not b then // If the downsteps (if any) didn't yield a
                    //  term, refine the upstep:

         D:=Floor(dyUpF/-dyDnF); // (this is a different "D")
                    // D = max # of downsteps whose addition
                    //  to the upstep wouldn't turn it into a
                    //  downstep (i.e., wouldn't yield a net
                    //  decrease in the fractional part of y).

         if D gt 0 then    // If D is positive, then
            dxUp+:=D*dxDn; //  grow curr upstep by D downsteps,
            dyUp:=c*dxUp;  //  and update the resulting y-step
            dyUpF:=dyUp-Round(dyUp); // and its fractional part.
         end if;

         U:=Ceiling(-(yF+dyDnF)/dyUpF);
                    // U = min # of upsteps whose addition to
                    //  one downstep would prevent the "signed
                    //  fractional part" of y from going below 0.

         dxDn+:=U*dxUp;  // Grow current downstep by U upsteps,
         dyDn:=c*dxDn;   //  and update the resulting y-step
         dyDnF:=dyDn-Round(dyDn); // and its fractional part.

      end if;

   end if;

end while;

X; // this line outputs terms of A341046
// YI; // this line, if uncommented, would output A341047