// OEIS A346534: "Denominators of approximations j/k for Pi such that // abs(j/k - Pi)*sqrt(5)*k^2 < 1." // /////////////////////////////////////////////////////////////////////////////// // // The Magma program below runs on the Magma Calculator at // // http://magma.maths.usyd.edu.au/calc/ // // and computes the first 31 terms of A346534 in about 0.01 seconds. // // The algorithm seems very fast, even for large terms. Changing the value of // kStop to 10^1000 and the value of DP (which needs to be at least a little // larger than 2*log_10(kStop) to avoid errors due to rounding) to 2025 and // replacing the line // // a; // // at the end of the program (which would print out the entire array of // computed terms, but that would exceed the output limit for the Magma // Calculator) with the line // // "a(" * IntegerToString(#a) * ") =", ChangePrecision(1.0*a[#a], 10); // // and rerunning the program gives the result (in less than 1 second) that the // last term found is a(2031) = 5.529760409E999. // // A question: Is there a name for this basic approach to searching for rational // approximations to a given irrational number? (I developed it myself, but my // guess is that some similar approach, and possibly a significantly better one, // is widely known.) // // Jon E. Schoenfield (with apologies that I'm still not very good at Magma) // Aug 09 2021 // /////////////////////////////////////////////////////////////////////////////// // // Algorithm description: // // On the j-k plane, define the "solution zone" as the zone in which // |j/k - Pi| * sqrt(5) * k^2 < 1. (This zone lies on both sides of the // centerline j = Pi*k, and becomes increasingly narrow as k increases.) // Equivalently, the solution zone is the zone in which E/M < 1, where // E = |j/k - Pi| and M = 1/(sqrt(5)*k^2). // // An integer k > 0 is a term of A346534 if and only if the point (j,k), where // j = round(Pi*k), lies inside the solution zone. // // The search implemented below progresses in such a way that, at any point in // the search, whatever the current values of k, dkUp, and dkDn are, all // remaining terms can be reached from the current value of k by incrementing k // by some nonnegative number of "up-steps" of size dkUp and some nonnegative // number of "down-steps" dkDn using the following simple rules: // // 1. If incrementing k by either dkUp or dkDn yields a point (j,k) inside // the solution zone, then increment k by dkUp or dkDn (whichever one // yields a point inside the solution zone). // // 2. Otherwise, increment k by dkUp or dkDn if the current point (j,k) is // on the "low" or "high" side of the centerline, i.e., j/k < Pi or // j/k > Pi, respectively. // // E.g., at a point early in the search where dkUp = 113 and dkDn = 106, // each time k is incremented by dkUp or dkDn, j will be incremented by // round(pi*dkUp) = round(354.99996985...) = 355 or round(pi*dkDn) = // 333.00882128... = 333, respectively, with the result that the values // of the ratio j/k remain very close to Pi. // // But incrementing k by, e.g., 113 or 106 at a time would not make it practical // to reach large values of k. A key part of the algorithm implemented below is // the way it grows the increments dkUp and dkDn. The full set of rules used for // incrementing k and growing dkUp and dkDn is as follows: // // 1. If taking a single up-step (i.e., increasing k by dkUp) from the current // value of k would yield a point inside the solution zone AND taking a // single down-step (i.e., increasing k by dkDn) from the current value of // k would ALSO yield a point inside the solution zone, then report an // error and terminate. (I have seen this happen only if the number of // digits of precision specified was insufficient; I don't think it can // happen if adequate precision is used.) Otherwise: // // 2. If taking a single up-step (i.e., increasing k by dkUp) would yield a // point inside the solution zone, take it. Otherwise: // // 3. If taking a single down-step (i.e., increasing k by dkDn) would yield a // point inside the solution zone, take it. Otherwise: // // 4. If taking a single step in the direction of the solution zone's // centerline (j/k = Pi) -- i.e., a single up-step taken if j/k < Pi at // the current value of k, or a single down-step taken if j/k > Pi at the // current value of k -- yields a point that is on the same side of the // centerline (in other words, if the step taken would move j/k toward Pi, // but would not change the sign of j/k - Pi), then take that single step. // Otherwise: // // 5. The only remaining possibility is that a single step in the direction of // the solution zone's centerline would result in crossing the centerline // and continuing in that direction so far as to completely overshoot the // solution zone. // // Suppose the current value of k is on the "low" side of the centerline, // i.e., that j/k < Pi. Then a step toward the centerline is an up-step, // dkUp. If the net effect of a single up-step plus a single down-step // is to increase j/k, then, from the current value of k, the path to // every remaining term of the sequence (i.e., every remaining point (j,k) // that lies inside the solution zone) must include at least one down-step // for each up-step. We could then speed up the search by combining one // up-step with one down-step and making the result (dkUp + dkDn) the new // up-step. // // Similarly, if the current value of k is on the "high" side of the // centerline, and if the net effect of a single down-step plus a single // up-step is to decrease j/k, then, from the current value of k, the path // to every remaining term of the sequence must include at least one up-step // for each down-step, so we could speed up the search by combining one // down-step with one up-step and making the result (dkDn + dkUp) the new // down-step. // // But there are situations (these arise in conjunction with large terms in // the continued fraction expansion of Pi) in which the effect (on j/k) of a // single step toward the centerline is far greater than the effect of a // single step in the opposite direction. In such a case, combining a single // step toward the centerline with a single step in the opposite direction, // and using the resulting combined step size as the new step toward the // centerline, results in a new step toward the centerline whose effect on // j/k is only slightly smaller than before, and continues to be much greater // than the effect of a single step in the opposite direction. So, to speed // things up, the algorithm as implemented below does the following whenever // none of the conditions for paragraphs 1 through 4 above apply: // // Let m be the minimum (positive integer) number of steps in the direction // away from the centerline such that the net effect (on j/k) of a single // step toward the centerline plus m steps away from the centerline is a // move that does NOT completely overshoot the solution zone. If that net // effect is an increase in j/k, then use the combined step (i.e., one step // toward the centerline plus m steps away) as the new up-step; otherwise, // use it as the new down-step. // // The result is that, as the search progresses, the up-steps and down-steps // are repeatedly refined so that dkUp and dkDn become larger and larger (and // take values that are denominators of fractions that more and more closely // approximate Pi), and no points that lie inside the solution zone are missed. // /////////////////////////////////////////////////////////////////////////////// // Set default precision (note: need 10^DP larger than kStop^2) DP := 50; R := RealField(DP); SetDefaultRealField(R); kStop := 10^11; // stop the search when k exceeds this value pi := Pi(R); sqrt5 := Sqrt(5); // for computing parameter M = sqrt(5)*k^2 k := 1; // k=1 is the first term of the sequence // and the starting point for the search a := [1]; // store k=1 as the first term of the sequence // Set the initial values of the step sizes delta(k) such that the resulting // changes to j=round(k*Pi) (i.e., delta(j)) will give // // delta(j)/delta(k) < Pi (this delta(k) value is dkDn) // // and // // delta(j)/delta(k) > Pi (this delta(k) value is dkUp) dkDn := 1; // delta(k)=1 -> delta(j)=3; delta(j)/delta(k) = 3 < Pi dkUp := 6; // delta(k)=6 -> delta(j)=19; delta(j)/delta(k) = 19/6 > Pi // Note: initializing dkUp to 6 makes the search work properly; from k=1 // (= a(1)), increasing k by dkUp=6 moves to k=7 (= a(2)), and from there, the // logic implemented below advances the value of k through all values that need // to be tested. bUpStepTested := false; // an up-step (from the current k) has not been tested bDnStepTested := false; // a down-step (from the current k) has not been tested while k lt kStop do j := Round(pi*k); Esgnd := j/k - pi; // signed value of parameter "E = abs(j/k - Pi)" if not bUpStepTested then // Test the value of k at one up-step from the current value of k; // does that value yield a term of the sequence? kTestUp := k + dkUp; jTestUp := Round(pi * kTestUp); EsgndTestUp := jTestUp / kTestUp - pi; ratioTestUp := EsgndTestUp * sqrt5 * kTestUp^2; // ratio is (signed) E/M bUpTerm := (Abs(ratioTestUp) lt 1.0); // (Boolean variable) bUpStepTested := true; end if; if not bDnStepTested then // Test the value of k at one down-step from the current value of k; // does that value yield a term of the sequence? kTestDn := k + dkDn; jTestDn := Round(pi * kTestDn); EsgndTestDn := jTestDn / kTestDn - pi; ratioTestDn := EsgndTestDn * sqrt5 * kTestDn^2; // ratio is (signed) E/M bDnTerm := (Abs(ratioTestDn) lt 1.0); // (Boolean variable) bDnStepTested := true; end if; if bUpTerm and bDnTerm then "Error: both k + dkUp and k + dkDn were", "found to yield a term when k =", k, "(this may happen as a result of rounding error if DP", "is not sufficiently larger than 2*log_10(kStop))."; "Terminating."; // This has never happened in my testing, aside from situations resulting // from rounding error, caused by using a value of DP that was too low, // given the value of kStop. // *Can* this ever happen, aside from rounding error? quit; elif bUpTerm or bDnTerm then // an up-step or down-step from k yields a term if bUpTerm then // up-step yields a term k +:= dkUp; // take an up-step else // down-step yields a term k +:= dkDn; // take a down-step end if; a[#a+1] := k; // store the new value of k as a term bUpStepTested := false; // an up-step from (new) k has not been tested bDnStepTested := false; // a down-step from (new) k has not been tested else // neither step yields a term // Would a single step toward the centerline (i.e., E/M = 0) of the // solution zone overshoot the zone entirely? if Esgnd lt 0.0 then // j/k < Pi: on low side, so test an up-step bOvershoot := (ratioTestUp gt 1.0); // (Boolean variable) else // j/k > Pi: on high side, so test a down-step bOvershoot := (ratioTestDn lt -1.0); // (Boolean variable) end if; if not bOvershoot then // A step toward the centerline would fall short of reaching the // solution zone (but would be closer to it). if Esgnd lt 0.0 then // j/k < Pi k +:= dkUp; // take an up-step else // j/k > Pi k +:= dkDn; // take a down-step end if; bUpStepTested := false; // an up-step from (new) k has not been tested bDnStepTested := false; // a down-step from (new) k has not been tested else // A step toward the centerline would overshoot the entire solution // zone. Determine the minimum positive number m of steps in the // opposite direction that would have to be combined with a single step // toward the centerline to yield a point that is not beyond the // solution zone. // As a simple way to do this, solve first for the value of m that // would yield E = 0 (i.e., j/k = Pi); then adjust m as necessary. djUp := Round(pi*dkUp); djDn := Round(pi*dkDn); if Esgnd lt 0.0 then // j/k < Pi: an up-step is toward the centerline // need (j + djUp + m*djDn)/(k + dkUp + m*dkDn) = Pi; // solving for m gives // // m = (Pi*(k + dkUp) - j - djUp) / (djDn - Pi*dkDn) m := (pi * (k + dkUp) - j - djUp) / (djDn - pi*dkDn); m := Floor(m); // round down to integer if m lt 1 then // but not less than 1 m := 1; end if; // while (m > 1) and (1 up-step + (m-1) down-steps would not // overshoot the zone), decrement m by 1 while (m gt 1) do // test the result of 1 up-step plus (m-1) down-steps kTest := k + dkUp + (m-1)*dkDn; jTest := j + djUp + (m-1)*djDn; ratioTest := (jTest/kTest - pi) * sqrt5 * kTest^2; if ratioTest lt 1.0 then // 1 up-step + (m-1) down-steps would not overshoot the zone m -:= 1; else break; end if; end while; dkTest := dkUp + m*dkDn; else // j/k > Pi: a down-step is toward the centerline // need (j + djDn + m*djUp)/(k + dkDn + m*dkUp) = Pi; // solving for m gives // // m = (Pi*(k + dkDn) - j - djDn) / (djUp - Pi*dkUp) m := (pi * (k + dkDn) - j - djDn) / (djUp - pi*dkUp); m := Floor(m); // round down to integer if m lt 1 then // but not less than 1 m := 1; end if; // while (m > 1) and (1 down-step + (m-1) up-steps would not // overshoot the zone), decrement m by 1 while (m gt 1) do // test the result of 1 down-step plus (m-1) up-steps kTest := k + dkDn + (m-1)*dkUp; jTest := j + djDn + (m-1)*djUp; ratioTest := (jTest/kTest - pi) * sqrt5 * kTest^2; if ratioTest gt -1.0 then // 1 down-step + (m-1) up-steps would not overshoot the zone m -:= 1; else break; end if; end while; dkTest := dkDn + m*dkUp; end if; dkTestPi := dkTest * pi; if Round(dkTestPi) gt dkTestPi then // The combined step is in the direction of increasing j/k, so use // it as the new up-step. dkUp := dkTest; bUpStepTested := false; // (new) up-step from k has not been tested // (but the existing down-step has) else // The combined step is in the direction of decreasing j/k, so use // it as the new down-step. dkDn := dkTest; bDnStepTested := false; // (new) down-step from k has not been // tested (but the existing up-step has) end if; end if; end if; end while; a; // output the array of terms