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