// Magma program for generating terms of 
//  OEIS A088306:  "Integers n with tan n > |n|, ordered by |n|."
//
// Jon E. Schoenfield, Aug 19 2014; updated Nov 07 2014 so that
//  the output will reflect the change of the sequence's Name from
//
//     Positive integers n with |tan n| > n.
// to
//     Integers n with tan n > |n|, ordered by |n|.
//
//
// Given a positive integer n, define A = n - Pi * floor(n/Pi); then A
//  is the angle in the interval (0, Pi) such that tan A = tan n.  The
//  inequality |tan n| > n is satisfied iff A - Pi/2 lies within the
//  success interval (arctan(-1/n), arctan(1/n)).  As n increases and
//  the success interval becomes narrower, an exhaustive search for
//  values of n satisfying |tan n| > n can be expedited by
//  successively incrementing n by only those increments that will
//  raise or lower A by small amounts; as n increases, the largest
//  changes to A that need to be considered become smaller, and the
//  corresponding increments to n become larger.
//
// (After each positive integer n for which |tan n| > n is found,
//  the value given as output is
//      n  if tan n > 0
//     -n  if tan n < 0
//  to agree with the changed Name of the sequence.)

OutputDP:=0; // total number of dgts of precision for output values;
             //  use 0 to output exact values
OutputExactMax:=10^OutputDP;

nMaxSearchDgts:=73; // (there are 100 terms < 10^73)
nMaxSearch:=10^nMaxSearchDgts;

DP:=2*nMaxSearchDgts+10; // would probably be okay w/o the +10
R:=RealField(DP);
SetDefaultRealField(R);
pi:=Pi(R);
halfPi:=0.5*pi;

dnUp:=1; // initialize step size that will increase angle A
dnDn:=2; // initialize step size that will decrease angle A
nTermsFound:=0;
n:=0;
s:=1;
bSuccessAtN:=false;

while n lt nMaxSearch do

   if bSuccessAtN then
      // Does a step AWAY FROM the A=Pi/2 centerline
      //  yield ANOTHER successful point?
      if s eq 1 then // sign(t)=1, so A < Pi/2
         nTest:=n+dnDn; // make angle A even lower
      else // A > Pi/2
         nTest:=n+dnUp; // make angle A even higher
      end if;
      if Abs(Tan(nTest)) gt nTest then // another success!
         // Is this even possible?  (Not seen for any nTest < 10^2000.)
         "Both n =", n, "and nTest =", nTest, "are in the sequence;";
         "terminating.";
         break;
      end if;
   end if;

   // Test the result of taking 1 step in the inward (i.e.,
   //  toward A=Pi/2) direction
   if s eq 1 then // A < Pi/2
      dnIn:=dnUp;
   else // A > Pi/2
      dnIn:=dnDn;
   end if;

   nTest:=n+dnIn;
   tTest:=Tan(nTest);
   sTest:=Sign(tTest);

   if Abs(tTest) gt nTest then // success!

      n:=nTest;
      t:=tTest;
      s:=sTest;
      bSuccessAtN:=true;
      nTermsFound+:=1;
      if (OutputDP eq 0) or (n le OutputExactMax) then
         nTermsFound, n*s; // "*s" (here and below) for signed output 
      else
         nTermsFound, ChangePrecision(n*s*1.0, OutputDP);
      end if;

   elif sTest eq s then // out of bounds on same side of centerline
                        // (i.e., undershot the success interval)

      // take the minimum number of steps in the inward direction
      //  that will at least reach (and possibly overshoot)
      //  the success interval

      A:=n-pi*Floor(n/pi);
      AThresh:=halfPi-s/n; // i.e., subtract s*(1/n); using 1/n instead
                           //  of Arctan(1/n) will make the k value
                           //  (found below) a lower bound; so will
                           //  the fact that |Arctan(1/n)| decreases
                           //  as n increases
      dAIn:=dnIn-pi*Round(dnIn/pi);
      k:=Ceiling((AThresh-A)/dAIn);

      n+:=k*dnIn;
      t:=Tan(n);
      s:=Sign(t);
      bSuccessAtN:=(Abs(t) gt n);

      if bSuccessAtN then
         nTermsFound+:=1;
         if (OutputDP eq 0) or (n le OutputExactMax) then
            nTermsFound, n*s;
         else
            nTermsFound, ChangePrecision(n*s*1.0, OutputDP);
         end if;
      end if;

   else // out of bounds on opposite side of centerline
        //  (i.e., overshot the success interval)

      if s eq 1 then // A < Pi/2
         dnOut:=dnDn;
      else // A > Pi/2
         dnOut:=dnUp;
      end if;

      dAIn:=dnIn-pi*Round(dnIn/pi);
      dAOut:=dnOut-pi*Round(dnOut/pi);

      if Abs(dAIn) lt Abs(dAOut) then // inward step is "narrower"

         // proceed with the (known overshot) step

         n:=nTest;
         t:=tTest;
         s:=sTest;
         bSuccessAtN:=false;

      else // inward step is "wider" than outward step

         A:=n-pi*Floor(n/pi);

         // 1 step dnIn would land at angle A+dAIn;
         //  compute a lower bound on the number k of subsequent
         //  steps dnOut to get back to success interval:

         AThresh:=halfPi+s/n; // far side of success interval (and, as
                              //  before, using 1/n instead of Arctan)
         k:=Ceiling((AThresh-(A+dAIn))/dAOut);

         if k gt Abs(dAIn/dAOut) then
            // augmenting dnIn by k*dnOut would cause the resulting
            //  dnIn to yield a dAIn in the wrong direction, i.e., out
            k-:=1;
         end if;

         if s eq 1 then // A < Pi/2; dnIn is dnUp;
            dnUp+:=k*dnDn; // augment it by k*dnOut, i.e., k*dnDn
            n+:=dnUp;
         else // A > Pi/2; dnIn is dnDn;;
            dnDn+:=k*dnUp; // augment it by k*dnOut, i.e., k*dnUp
            n+:=dnDn;
         end if;

         t:=Tan(n);
         s:=Sign(t);
         bSuccessAtN:=(Abs(t) gt n);

         if bSuccessAtN then
            nTermsFound+:=1;
            if (OutputDP eq 0) or (n le OutputExactMax) then
               nTermsFound, n*s;
            else
               nTermsFound, ChangePrecision(n*s*1.0, OutputDP);
            end if;
         end if;

      end if;

   end if;

end while;