// Magma program to compute terms of OEIS A342171,
//  "Nonnegative integers k such that k < sec(k)*csc(k)."
//
// The program below outputs the first 100 terms of A342171
//  (all terms up to 10^77) in about 0.1 seconds when run on
//  the Magma Calculator at http://magma.maths.usyd.edu.au/calc/
//
// J. Schoenfield, Mar 14 2021


// Basic approach:
//
// Set initial values of two small increments to k: a "down-step"
//  dkDn=1 and an "up-step" dkUp=2 such that, when 2*k/Pi is
//  close to an integer, incrementing k by dkDn or dkUp will
//  cause the (signed) difference
//
//        2*k/Pi - round(2*k/Pi)
//
//  to decrease or increase, respectively.
//
// Initialize kCurr to 1; then, for every positive integer kTest
//  up to dkUp+dkDn, set kCurr to kTest if kTest is a term of
//  this sequence -- i.e., if
//
//        t(kTest) = sec(kTest)*csc(kTest) = 2/sin(2*kTest)
//
//  exceeds kTest -- or if
//
//        |t(kTest)| > |t(kCurr)|.
//
// Then, starting from the last k value to which kCurr was set as
//  described above, do the following until kCurr exceeds the
//  value kEnd (set below):
//
//    For each of two larger values of k,
//
//        kTestDn=kCurr+dkDn
//    and
//        kTestUp=kCurr+dkUp
//
//    compute
//
//        t(k) = sec(k)*csc(k) = 2/sin(2*k)
//
//    as above. If either value of k yields |t(k)| > k, then set
//    kCurr to that value of k (and if that t(k) is positive,
//    store that value of k as a term of this sequence).
//    Otherwise, if either value of k yields |t(k)| > |t(kCurr)|,
//    then set kCurr to that value of k. Otherwise, combine a
//    down-step and an up-step as dkTest=dkDn+dkUp and use the
//    result as the new (refined) down-step or up-step (depending
//    on whether 2*dkTest/Pi is smaller or larger than its
//    nearest integer, respectively).


a:=[];
kEnd:=10^77; // to get all terms up 10^77
DP:=Max(2*Ceiling(Log(10,kEnd))+20, 30);
                    // DP digits of precision should be adequate

R:=RealField(DP);
SetDefaultRealField(R);
pi:=Pi(R);

dkDnInit:=1; // given an integer k such that 2*k/Pi is very close
             //  to an integer, incrementing k by 1 will make
             //  2*k/Pi smaller than its nearest integer

dkUpInit:=2; // given an integer k such that 2*k/Pi is very close
             //  to an integer, incrementing k by 2 will make
             //  2*k/Pi larger than its nearest integer

kInitTestMax:=dkDnInit+dkUpInit;


// Do exhaustive search through very small values of k

kCurr:=1;
tCurr:=2/Sin(2*kCurr);

for kTest in [1..kInitTestMax] do

   tTest:=2/Sin(2*kTest);

   if kTest lt tTest then
      a[#a+1]:=kTest;
      bUpdateKCurr:=true;
   elif Abs(tTest) gt Abs(tCurr) then
      bUpdateKCurr:=true;
   else
      bUpdateKCurr:=false;
   end if;

   if bUpdateKCurr then
      kCurr:=kTest;
      tCurr:=tTest;
   end if;

end for;


dkDn:=dkDnInit;
dkUp:=dkUpInit;

while kCurr lt kEnd do

   // test a down-step and an up-step from kCurr

   kTestDn:=kCurr+dkDn;
   kTestUp:=kCurr+dkUp;

   tTestDn:=2/Sin(2*kTestDn);
   tTestUp:=2/Sin(2*kTestUp);

   bDnStepYieldsGoodAbsT:=(kTestDn lt Abs(tTestDn));
   bUpStepYieldsGoodAbsT:=(kTestUp lt Abs(tTestUp));

   if bDnStepYieldsGoodAbsT and bUpStepYieldsGoodAbsT then

      // Both types of steps yield k values such that |t(k)| > k.
      //  Can this ever happen?
      //  (It doesn't happen for any term < 10^1000.)

      "Error: unexpected result at";
      "       kCurr =", kCurr, "  kTestDn =", kTestDn,
            "  kTestUp =", kTestUp;
      "Both k=kTestDn and k=kTestUp yield |t(k)| > k.";
      "Terminating.";

      quit;

   elif bDnStepYieldsGoodAbsT then

      if tTestDn gt 0.0 then
         a[#a+1]:=kTestDn;
      end if;

      kCurr:=kTestDn;
      tCurr:=tTestDn;

   elif bUpStepYieldsGoodAbsT then

      if tTestUp gt 0.0 then
         a[#a+1]:=kTestUp;
      end if;

      kCurr:=kTestUp;
      tCurr:=tTestUp;

   else // neither type of step yields |t(k)| > k

      // From kCurr, reaching any larger k that is a term will
      //  require at least 1 up-step AND at least 1 down-step;
      //  test the step that would result from using 1 of each.

      // (Note: on some passes through this loop,
      //
      //       |2*dkDn/pi - round(2*dkDn/pi)|
      //  and
      //       |2*dkUp/pi - round(2*dkUp/pi)|
      //
      //  are very different in size; one can be thousands of
      //  times larger than the other, which can result in
      //  thousands of passes through the loop before the next
      //  term of the sequence is found. (E.g., if kEnd is set
      //  to a large enough value [such as 10^221] so that the
      //  program finds 321 or more terms of the sequence, there
      //  are more than 10000 passes through the loop between
      //  finding a(320) and finding a(321).) Code could be added
      //  here to combine multiple up-steps with a single
      //  down-step or vice versa to speed up the process.

      dkTest:=dkDn+dkUp;
      fracTest:=2*dkTest/pi;

      if fracTest lt Round(fracTest) then

         // combining 1 down-step with 1 up-step would result
         //  in a refined down-step

         dkDn:=dkTest;

      else

         dkUp:=dkTest;

      end if;

   end if;

end while;

a; // output all terms found