CP Notebook

← all snippets

ContinuedFractions

Given N and a real number x ≥ 0, finds the closest rational approximation p/q with p, q ≤ N. It will obey |p/q - x| ≤ 1/qN. For consecutive convergents, p_k+1q_k - q_k+1p_k = (-1)^k. (p_k/q_k alternates between >x and <x.) If x is rational, y eventually becomes ∞; if x is the root of a degree 2 polynomial the a's eventually become cyclic.

Time: O(log N) 21 lines stress-tested for n <= 300

content/number-theory/ContinuedFractions.h — Simon Lindholm, source: Wikipedia

typedef double d; // for N ~ 1e7; long double for N ~ 1e9
pair<ll, ll> approximate(d x, ll N) {
	ll LP = 0, LQ = 1, P = 1, Q = 0, inf = LLONG_MAX; d y = x;
	for (;;) {
		ll lim = min(P ? (N-LP) / P : inf, Q ? (N-LQ) / Q : inf),
		   a = (ll)floor(y), b = min(a, lim),
		   NP = b*P + LP, NQ = b*Q + LQ;
		if (a > b) {
			// If b > a/2, we have a semi-convergent that gives us a
			// better approximation; if b = a/2, we *may* have one.
			// Return {P, Q} here for a more canonical approximation.
			return (abs(x - (d)NP / (d)NQ) < abs(x - (d)P / (d)Q)) ?
				make_pair(NP, NQ) : make_pair(P, Q);
		}
		if (abs(y = 1/(y - (d)a)) > 3*N) {
			return {NP, NQ};
		}
		LP = P; P = NP;
		LQ = Q; Q = NQ;
	}
}