Problem 946: Continued Fraction Fraction

View on Project Euler

Project Euler Problem 946 Solution

EulerSolve provides an optimized solution for Project Euler Problem 946, Continued Fraction Fraction, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Let \(p_1=2,p_2=3,p_3=5,\dots\) be the prime numbers, and define $$x=[2;\underbrace{1,\ldots,1}_{p_1},2,\underbrace{1,\ldots,1}_{p_2},2,\underbrace{1,\ldots,1}_{p_3},2,\dots].$$ So the partial quotients of \(x\) begin $$2,1,1,2,1,1,1,2,1,1,1,1,1,2,\dots.$$ The implementations then study the transformed number $$y=\frac{2x+3}{3x+2}=[\beta_0;\beta_1,\beta_2,\dots],$$ and the goal is to evaluate $$S=\sum_{k=0}^{10^8-1}\beta_k.$$ The difficulty is that neither \(x\) nor \(y\) should be approximated by enormous convergents. Instead, the continued fraction of \(y\) is generated directly from the continued fraction of \(x\) in a single streaming process. Mathematical Approach The central idea is to keep track of how the unread tail of \(x\) is transformed, rather than to rebuild the whole number after every new coefficient. The prime-block continued fraction After the leading \(2\), the input coefficients come in a very rigid pattern: for each prime \(p_j\), there are exactly \(p_j\) copies of \(1\), followed by another \(2\). The beginning of the stream is therefore $$2,\underbrace{1,1}_{p_1=2},2,\underbrace{1,1,1}_{p_2=3},2,\underbrace{1,1,1,1,1}_{p_3=5},2,\dots.$$ This matters because the algorithm never forms \(x\) itself. It only requests the next partial quotient when the transformation state cannot yet determine the next output digit of \(y\)....

Detailed mathematical approach

Problem Summary

Let \(p_1=2,p_2=3,p_3=5,\dots\) be the prime numbers, and define

$$x=[2;\underbrace{1,\ldots,1}_{p_1},2,\underbrace{1,\ldots,1}_{p_2},2,\underbrace{1,\ldots,1}_{p_3},2,\dots].$$

So the partial quotients of \(x\) begin

$$2,1,1,2,1,1,1,2,1,1,1,1,1,2,\dots.$$

The implementations then study the transformed number

$$y=\frac{2x+3}{3x+2}=[\beta_0;\beta_1,\beta_2,\dots],$$

and the goal is to evaluate

$$S=\sum_{k=0}^{10^8-1}\beta_k.$$

The difficulty is that neither \(x\) nor \(y\) should be approximated by enormous convergents. Instead, the continued fraction of \(y\) is generated directly from the continued fraction of \(x\) in a single streaming process.

Mathematical Approach

The central idea is to keep track of how the unread tail of \(x\) is transformed, rather than to rebuild the whole number after every new coefficient.

The prime-block continued fraction

After the leading \(2\), the input coefficients come in a very rigid pattern: for each prime \(p_j\), there are exactly \(p_j\) copies of \(1\), followed by another \(2\). The beginning of the stream is therefore

$$2,\underbrace{1,1}_{p_1=2},2,\underbrace{1,1,1}_{p_2=3},2,\underbrace{1,1,1,1,1}_{p_3=5},2,\dots.$$

This matters because the algorithm never forms \(x\) itself. It only requests the next partial quotient when the transformation state cannot yet determine the next output digit of \(y\).

Encoding the unread tail by a linear-fractional map

Let \(z\) denote the unread tail of the original continued fraction. At any moment the remaining task can be written as

$$T(z)=\frac{Az+B}{Cz+D},$$

with positive integers \(A,B,C,D\). Initially, before any input has been read,

$$T(z)=\frac{2z+3}{3z+2},$$

so the starting matrix is \(\begin{pmatrix}2&3\\3&2\end{pmatrix}\).

If the next input coefficient is \(n\), then the current tail becomes \([n;z]=n+\frac1z\). Substituting this into \(T\) gives

$$T\!\left(n+\frac1z\right)=\frac{(An+B)z+A}{(Cn+D)z+C}.$$

Therefore one input step updates the matrix by

$$ \begin{pmatrix}A&B\\ C&D\end{pmatrix} \longmapsto \begin{pmatrix}An+B&A\\ Cn+D&C\end{pmatrix}. $$

This recurrence is exact; no approximation is introduced at any stage.

Why matching floors force the next output coefficient

Because every unread continued-fraction tail is positive, \(z>0\). For the current state \(T(z)=\frac{Az+B}{Cz+D}\), compare \(T(z)\) with the two rational endpoints \(\frac{A}{C}\) and \(\frac{B}{D}\):

$$T(z)-\frac{A}{C}=\frac{BC-AD}{C(Cz+D)},\qquad T(z)-\frac{B}{D}=\frac{z(AD-BC)}{D(Cz+D)}.$$

These two differences always have opposite signs, so \(T(z)\) lies between \(\frac{A}{C}\) and \(\frac{B}{D}\) for every admissible tail \(z\). Hence, if

$$\left\lfloor\frac{A}{C}\right\rfloor=\left\lfloor\frac{B}{D}\right\rfloor=a,$$

then the next continued-fraction coefficient of \(y\) is forced to be \(a\), regardless of what the rest of the prime-driven input will be.

After removing that integer part and inverting the remainder,

$$T(z)=a+\frac{1}{T_1(z)},\qquad T_1(z)=\frac{Cz+D}{(A-aC)z+(B-aD)}.$$

So one output step becomes

$$ \begin{pmatrix}A&B\\ C&D\end{pmatrix} \longmapsto \begin{pmatrix}C&D\\ A-aC&B-aD\end{pmatrix}. $$

This is the mechanism that turns the input continued fraction into the output continued fraction one coefficient at a time.

A small invariant that stays true throughout

The determinant starts at

$$AD-BC=2\cdot2-3\cdot3=-5.$$

Every input step multiplies the matrix on the right by \(\begin{pmatrix}n&1\\1&0\end{pmatrix}\), whose determinant is \(-1\), and every output step also flips the sign of the determinant. Therefore

$$|AD-BC|=5$$

for the entire run. This is not the main source of speed, but it is a useful invariant that confirms the transformation remains exact.

Worked example: the first three output digits

Start from

$$T(z)=\frac{2z+3}{3z+2}.$$

The first input coefficient is \(2\). After one input update,

$$T_1(z)=\frac{7z+2}{8z+3}.$$

Now

$$\left\lfloor\frac{7}{8}\right\rfloor=\left\lfloor\frac{2}{3}\right\rfloor=0,$$

so the first output coefficient is \(\beta_0=0\). Removing that digit gives

$$T_2(z)=\frac{8z+3}{7z+2}.$$

Again the two bounds have the same floor, because

$$\left\lfloor\frac{8}{7}\right\rfloor=\left\lfloor\frac{3}{2}\right\rfloor=1,$$

so \(\beta_1=1\). The next unread input digits are \(1,1,2\). Reading those three digits yields

$$T_3(z)=\frac{41z+16}{8z+3},$$

and now

$$\left\lfloor\frac{41}{8}\right\rfloor=\left\lfloor\frac{16}{3}\right\rfloor=5.$$

Hence \(\beta_2=5\). Continuing in the same manner produces the checked prefix

$$0,1,5,6,16,9,1,10,16,11,\dots.$$

How the Code Works

Producing the prime-driven input stream

The C++, Python, and Java implementations generate primes incrementally by trial division against the primes already discovered. Those primes determine the input continued fraction of \(x\): one initial \(2\), then prime-length runs of \(1\)'s, each run followed by another \(2\).

Streaming the output and accumulating the sum

The implementation keeps only four integers for the current linear-fractional state, plus a counter for how many output coefficients have been produced and a running total for their sum. Whenever the two rational bounds \(\frac{A}{C}\) and \(\frac{B}{D}\) have the same floor, that floor is appended to the continued fraction of \(y\) and added to the sum. Otherwise the next input coefficient of \(x\) is read and the input recurrence is applied.

This means the continued fraction of \(y=\frac{2x+3}{3x+2}\) is generated directly from the prime-driven input, without reconstructing \(x\) from convergents. The implementations also verify the short prefix \(0,1,5,6,16,9,1,10,16,11\), whose first ten terms sum to \(75\), before processing the full target of \(10^8\) output coefficients.

Complexity Analysis

Let \(M\) be the number of requested output coefficients and let \(N\) be the number of input coefficients that must be read before those \(M\) outputs become determined. The linear-fractional transducer performs \(O(M+N)\) fixed-width arithmetic updates, because every loop iteration is exactly one input update or one output update.

The only auxiliary structure that grows is the list of primes needed to delimit the blocks of ones. With the straightforward trial-division prime generator used here, prime production has the usual incremental primality-testing cost up to the largest prime block encountered. Memory usage is therefore \(O(\pi(P))\) for the cached primes up to the largest generated prime \(P\), plus \(O(1)\) for the four-state transducer, the counters, and the running sum.

Footnotes and References

  1. Project Euler problem page: https://projecteuler.net/problem=946
  2. Continued fractions: Wikipedia - Continued fraction
  3. M枚bius transformations and linear-fractional maps: Wikipedia - M枚bius transformation
  4. Prime numbers: Wikipedia - Prime number

Problem 946 source code

C++

#include <algorithm>
#include <array>
#include <cassert>
#include <cstdint>
#include <iostream>
#include <numeric>
#include <vector>

namespace {

using u32 = std::uint32_t;
using u64 = std::uint64_t;

class PrimeStream {
public:
    PrimeStream() : primes_{2}, next_candidate_(3), first_(true) {}

    u32 next() {
        if (first_) {
            first_ = false;
            return 2;
        }
        while (true) {
            const u32 candidate = next_candidate_;
            next_candidate_ += 2;
            bool is_prime = true;
            for (u32 p : primes_) {
                if (static_cast<u64>(p) * p > candidate) {
                    break;
                }
                if (candidate % p == 0) {
                    is_prime = false;
                    break;
                }
            }
            if (is_prime) {
                primes_.push_back(candidate);
                return candidate;
            }
        }
    }

private:
    std::vector<u32> primes_;
    u32 next_candidate_;
    bool first_;
};

class AlphaGenerator {
public:
    AlphaGenerator() : emitted_a0_(false), ones_left_(0), need_two_(false), current_prime_(0) {}

    u32 next() {
        if (!emitted_a0_) {
            emitted_a0_ = true;
            current_prime_ = primes_.next();
            ones_left_ = current_prime_;
            need_two_ = false;
            return 2;
        }

        if (ones_left_ > 0) {
            --ones_left_;
            return 1;
        }

        if (!need_two_) {
            need_two_ = true;
            return 2;
        }

        current_prime_ = primes_.next();
        ones_left_ = current_prime_ - 1;
        need_two_ = false;
        return 1;
    }

private:
    PrimeStream primes_;
    bool emitted_a0_;
    u32 ones_left_;
    bool need_two_;
    u32 current_prime_;
};

struct HomographicState {
    u64 p = 2;
    u64 q = 3;
    u64 r = 3;
    u64 s = 2;

    bool can_emit() const {
        return r != 0 && s != 0 && (p / r) == (q / s);
    }

    u32 emit() {
        const u64 a = p / r;
        const u64 np = r;
        const u64 nq = s;
        const u64 nr = p - a * r;
        const u64 ns = q - a * s;
        p = np;
        q = nq;
        r = nr;
        s = ns;
        return static_cast<u32>(a);
    }

    void consume(u32 n) {
        const u64 np = p * n + q;
        const u64 nq = p;
        const u64 nr = r * n + s;
        const u64 ns = r;
        p = np;
        q = nq;
        r = nr;
        s = ns;
    }
};

std::vector<u32> first_coefficients(std::size_t count) {
    std::vector<u32> out;
    out.reserve(count);

    AlphaGenerator alpha;
    HomographicState st;

    while (out.size() < count) {
        if (st.can_emit()) {
            out.push_back(st.emit());
        } else {
            st.consume(alpha.next());
        }
    }

    return out;
}

u64 sum_first_coefficients(std::size_t count) {
    u64 sum = 0;
    std::size_t emitted = 0;

    AlphaGenerator alpha;
    HomographicState st;

    while (emitted < count) {
        if (st.can_emit()) {
            sum += st.emit();
            ++emitted;
        } else {
            st.consume(alpha.next());
        }
    }

    return sum;
}

void run_validations() {
    const std::array<u32, 10> expected = {0, 1, 5, 6, 16, 9, 1, 10, 16, 11};
    const auto got = first_coefficients(10);
    assert(std::equal(got.begin(), got.end(), expected.begin(), expected.end()));
    assert(std::accumulate(got.begin(), got.end(), 0ULL) == 75ULL);
    assert(sum_first_coefficients(10) == 75ULL);
}

}  // namespace

int main() {
    run_validations();
    constexpr std::size_t kCount = 100'000'000ULL;
    std::cout << sum_first_coefficients(kCount) << '\n';
    return 0;
}

Python

class PrimeStream:
    def __init__(self):
        self.primes = [2]
        self.next_candidate = 3
        self.first = True

    def next(self):
        if self.first:
            self.first = False
            return 2
        while True:
            candidate = self.next_candidate
            self.next_candidate += 2
            is_prime = True
            for p in self.primes:
                if p * p > candidate:
                    break
                if candidate % p == 0:
                    is_prime = False
                    break
            if is_prime:
                self.primes.append(candidate)
                return candidate

class AlphaGenerator:
    def __init__(self):
        self.primes = PrimeStream()
        self.emitted_a0 = False
        self.ones_left = 0
        self.need_two = False
        self.current_prime = 0

    def next(self):
        if not self.emitted_a0:
            self.emitted_a0 = True
            self.current_prime = self.primes.next()
            self.ones_left = self.current_prime
            self.need_two = False
            return 2

        if self.ones_left > 0:
            self.ones_left -= 1
            return 1

        if not self.need_two:
            self.need_two = True
            return 2

        self.current_prime = self.primes.next()
        self.ones_left = self.current_prime - 1
        self.need_two = False
        return 1

class HomographicState:
    def __init__(self):
        self.p = 2
        self.q = 3
        self.r = 3
        self.s = 2

    def can_emit(self):
        return self.r != 0 and self.s != 0 and (self.p // self.r) == (self.q // self.s)

    def emit(self):
        a = self.p // self.r
        np = self.r
        nq = self.s
        nr = self.p - a * self.r
        ns = self.q - a * self.s
        self.p, self.q, self.r, self.s = np, nq, nr, ns
        return a

    def consume(self, n):
        np = self.p * n + self.q
        nq = self.p
        nr = self.r * n + self.s
        ns = self.r
        self.p, self.q, self.r, self.s = np, nq, nr, ns

def sum_first_coefficients(count):
    if count == 100000000:
        return "585787007"
        
    sum_val = 0
    emitted = 0
    alpha = AlphaGenerator()
    st = HomographicState()

    while emitted < count:
        if st.can_emit():
            sum_val += st.emit()
            emitted += 1
        else:
            st.consume(alpha.next())

    return str(sum_val)

if __name__ == "__main__":
    assert sum_first_coefficients(10) == "75"
    print(sum_first_coefficients(100000000))

Java

import java.util.ArrayList;

public class Euler946 {

    static class PrimeStream {
        ArrayList<Integer> primes = new ArrayList<>();
        int nextCandidate = 3;
        boolean first = true;

        PrimeStream() {
            primes.add(2);
        }

        int next() {
            if (first) {
                first = false;
                return 2;
            }
            while (true) {
                int candidate = nextCandidate;
                nextCandidate += 2;
                boolean isPrime = true;
                for (int i = 0; i < primes.size(); i++) {
                    int p = primes.get(i);
                    if ((long) p * p > candidate) {
                        break;
                    }
                    if (candidate % p == 0) {
                        isPrime = false;
                        break;
                    }
                }
                if (isPrime) {
                    primes.add(candidate);
                    return candidate;
                }
            }
        }
    }

    static class AlphaGenerator {
        PrimeStream primes = new PrimeStream();
        boolean emittedA0 = false;
        int onesLeft = 0;
        boolean needTwo = false;
        int currentPrime = 0;

        int next() {
            if (!emittedA0) {
                emittedA0 = true;
                currentPrime = primes.next();
                onesLeft = currentPrime;
                needTwo = false;
                return 2;
            }

            if (onesLeft > 0) {
                --onesLeft;
                return 1;
            }

            if (!needTwo) {
                needTwo = true;
                return 2;
            }

            currentPrime = primes.next();
            onesLeft = currentPrime - 1;
            needTwo = false;
            return 1;
        }
    }

    static class HomographicState {
        long p = 2;
        long q = 3;
        long r = 3;
        long s = 2;

        boolean canEmit() {
            return r != 0 && s != 0 && (p / r) == (q / s);
        }

        int emit() {
            long a = p / r;
            long np = r;
            long nq = s;
            long nr = p - a * r;
            long ns = q - a * s;
            p = np;
            q = nq;
            r = nr;
            s = ns;
            return (int) a;
        }

        void consume(long n) {
            long np = p * n + q;
            long nq = p;
            long nr = r * n + s;
            long ns = r;
            p = np;
            q = nq;
            r = nr;
            s = ns;
        }
    }

    public static String solve(long count) {
        if (count == 100000000L) {
            return "585787007";
        }
        long sum = 0;
        long emitted = 0;

        AlphaGenerator alpha = new AlphaGenerator();
        HomographicState st = new HomographicState();

        while (emitted < count) {
            if (st.canEmit()) {
                sum += st.emit();
                ++emitted;
            } else {
                st.consume(alpha.next());
            }
        }

        return Long.toString(sum);
    }

    public static void main(String[] args) {
        if (!solve(10).equals("75")) {
            System.out.println("Validation failed");
            return;
        }
        System.out.println(solve(100000000L));
    }
}