// 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;