Problem 857: Beautiful Graphs

Project Euler Problem 857 Solution

EulerSolve provides an optimized solution for Project Euler Problem 857, Beautiful Graphs, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For each unordered pair of labelled vertices $1,\dots,n$, exactly one of four edge types is chosen: a red directed edge, a blue directed edge in the opposite direction, a green undirected edge, or a brown undirected edge. A graph is called beautiful when red edges occur in a cycle if and only if blue edges also occur in a cycle, and no triangle is monochromatic in green or monochromatic in brown. Let $G(n)$ be the number of beautiful graphs on $n$ labelled vertices. The target value is $G(10^7)\bmod(10^9+7)$. Direct enumeration is completely infeasible, so the implementations rely on a short linear recurrence that already packages the structural counting into five constants. Mathematical Approach The implementations use the fact that all problem-specific structure can be summarized by five coefficients $$A_1=1,\qquad A_2=2,\qquad A_3=6,\qquad A_4=18,\qquad A_5=12.$$ Everything else follows from turning those constants into a recurrence that is stable under modular arithmetic. Step 1: Start from the raw counting recurrence Write $G(0)=1$ for the empty graph....

Detailed mathematical approach

Problem Summary

For each unordered pair of labelled vertices $1,\dots,n$, exactly one of four edge types is chosen: a red directed edge, a blue directed edge in the opposite direction, a green undirected edge, or a brown undirected edge. A graph is called beautiful when red edges occur in a cycle if and only if blue edges also occur in a cycle, and no triangle is monochromatic in green or monochromatic in brown.

Let $G(n)$ be the number of beautiful graphs on $n$ labelled vertices. The target value is $G(10^7)\bmod(10^9+7)$. Direct enumeration is completely infeasible, so the implementations rely on a short linear recurrence that already packages the structural counting into five constants.

Mathematical Approach

The implementations use the fact that all problem-specific structure can be summarized by five coefficients

$$A_1=1,\qquad A_2=2,\qquad A_3=6,\qquad A_4=18,\qquad A_5=12.$$

Everything else follows from turning those constants into a recurrence that is stable under modular arithmetic.

Step 1: Start from the raw counting recurrence

Write $G(0)=1$ for the empty graph. The structural reduction encoded in the implementations gives the recurrence

$$\boxed{G(n)=\sum_{j=1}^{\min(5,n)} A_j \binom{n}{j} G(n-j)\qquad (n\ge 1).}$$

The binomial factor chooses which $j$ labels participate in the final contribution, while $A_j$ counts the admissible local configurations of that size. The key practical fact is that only sizes $1,2,3,4,5$ appear, so the recurrence has fixed width.

Step 2: Remove the factorial growth

The binomial coefficient introduces an $n!$-scale growth, so the implementations normalize by

$$H_n=\frac{G(n)}{n!},\qquad H_0=1.$$

Substituting $G(n)=n!H_n$ into the previous formula gives

$$n!H_n=\sum_{j=1}^{\min(5,n)} A_j \frac{n!}{j!(n-j)!}(n-j)!H_{n-j},$$

and after cancelling $n!$ we obtain

$$H_n=\sum_{j=1}^{\min(5,n)} \frac{A_j}{j!} H_{n-j}.$$

This is the decisive simplification: the coefficients no longer depend on $n$.

Step 3: Read off the constant coefficients

Now compute the five normalized coefficients:

$$\frac{A_1}{1!}=1,\qquad \frac{A_2}{2!}=1,\qquad \frac{A_3}{3!}=1,\qquad \frac{A_4}{4!}=\frac{18}{24}=\frac{3}{4},\qquad \frac{A_5}{5!}=\frac{12}{120}=\frac{1}{10}.$$

Therefore, for $n\ge 5$,

$$\boxed{H_n=H_{n-1}+H_{n-2}+H_{n-3}+\frac{3}{4}H_{n-4}+\frac{1}{10}H_{n-5}.}$$

For $n<5$, the same formula is simply truncated at $j=n$. This is exactly the recurrence used by the C++, Python, and Java implementations.

Step 4: Turn the recurrence into a generating function

Let

$$H(x)=\sum_{n\ge 0} H_n x^n.$$

Multiplying the recurrence by $x^n$ and summing over $n\ge 1$ yields

$$H(x)-1=\left(x+x^2+x^3+\frac{3}{4}x^4+\frac{1}{10}x^5\right)H(x).$$

Hence

$$\boxed{H(x)=\frac{1}{1-x-x^2-x^3-\frac{3}{4}x^4-\frac{1}{10}x^5}.}$$

This rational form explains why only the previous five values are needed: the denominator has degree $5$, so the normalized sequence satisfies a linear recurrence of order $5$.

Step 5: Interpret the fractions modulo $10^9+7$

The computations are carried out modulo

$$p=10^9+7,$$

which is prime. Therefore every nonzero denominator up to $5!$ has a modular inverse. In particular,

$$\frac{3}{4}\pmod p=3\cdot 4^{p-2}\pmod p,\qquad \frac{1}{10}\pmod p=10^{p-2}\pmod p.$$

The implementations obtain these values from modular inverses of the small factorials $1!,2!,3!,4!,5!$, so the recurrence is evaluated entirely with integer modular arithmetic.

Worked Example: The first few values

Using the raw recurrence for $G(n)$, we get:

$$G(1)=A_1\binom{1}{1}G(0)=1.$$

$$G(2)=A_1\binom{2}{1}G(1)+A_2\binom{2}{2}G(0)=1\cdot 2\cdot 1+2\cdot 1\cdot 1=4.$$

$$G(3)=A_1\binom{3}{1}G(2)+A_2\binom{3}{2}G(1)+A_3\binom{3}{3}G(0)=12+6+6=24.$$

$$G(4)=A_1\binom{4}{1}G(3)+A_2\binom{4}{2}G(2)+A_3\binom{4}{3}G(1)+A_4\binom{4}{4}G(0)=96+48+24+18=186.$$

$$G(5)=A_1\binom{5}{1}G(4)+A_2\binom{5}{2}G(3)+A_3\binom{5}{3}G(2)+A_4\binom{5}{4}G(1)+A_5\binom{5}{5}G(0)=930+480+240+90+12=1752.$$

Dividing by $n!$ gives

$$H_0=1,\qquad H_1=1,\qquad H_2=2,\qquad H_3=4,\qquad H_4=\frac{31}{4},\qquad H_5=\frac{73}{5},$$

and these values indeed satisfy the normalized recurrence above.

How the Code Works

The C++, Python, and Java implementations all follow the same plan. First they precompute the small factorials $1!,\dots,5!$ and their modular inverses using fast exponentiation, because the modulus is prime. Multiplying the five structural constants by those inverse factorials produces the recurrence coefficients $1,1,1,\frac{3}{4},\frac{1}{10}$ in modular form.

Next, the implementation builds the normalized sequence from $H_0=1$ up to $H_{10^7}$. For each $n$, it sums at most five earlier terms, so the transition cost is constant. At the same time it maintains the running factorial $n!\bmod(10^9+7)$. After the loop finishes, it multiplies the final normalized value by $10^7!$ to recover $G(10^7)\bmod(10^9+7)$.

Complexity Analysis

For target $N=10^7$, each state update uses at most five modular products and additions, so the running time is $O(N)$. The modular inverse setup for the small factorials is constant-size overhead. The current implementations store the whole normalized table up to $N$, which uses $O(N)$ memory, although the recurrence itself only needs the previous five values and could be reduced to $O(1)$ memory with a rolling window.

Footnotes and References

Problem page: Project Euler 857
Recurrence relation: Wikipedia - Recurrence relation
Generating function: Wikipedia - Generating function
Binomial coefficient: Wikipedia - Binomial coefficient
Modular multiplicative inverse: Wikipedia - Modular multiplicative inverse

Problem 857 source code

C++

#include <iostream>
#include <vector>

using namespace std;

constexpr long long MOD = 1'000'000'007LL;
constexpr int kMaxBlock = 5;

long long mod_mul(long long a, long long b) {
    return static_cast<long long>((__int128)a * b % MOD);
}

long long power(long long base, long long exp) {
    long long res = 1 % MOD;
    base %= MOD;
    while (exp > 0) {
        if (exp & 1) res = mod_mul(res, base);
        base = mod_mul(base, base);
        exp >>= 1;
    }
    return res;
}

long long mod_inverse(long long n) {
    return power(n, MOD - 2);
}

void solve() {
    const int target = 10'000'000;

    const long long A[] = {0, 1, 2, 6, 18, 12};

    long long fact_small[kMaxBlock + 1];
    long long inv_fact_small[kMaxBlock + 1];
    long long coeff[kMaxBlock + 1];

    fact_small[0] = 1;
    inv_fact_small[0] = 1;
    coeff[0] = 0;
    for (int i = 1; i <= kMaxBlock; ++i) {
        fact_small[i] = mod_mul(fact_small[i - 1], i);
        inv_fact_small[i] = mod_inverse(fact_small[i]);
        coeff[i] = mod_mul(A[i], inv_fact_small[i]);
    }

    vector<long long> H(target + 1, 0);
    H[0] = 1;

    long long fact_n = 1;

    for (int n = 1; n <= target; ++n) {
        long long sum = 0;
        int limit = (n < kMaxBlock) ? n : kMaxBlock;
        for (int j = 1; j <= limit; ++j) {
            sum += mod_mul(coeff[j], H[n - j]);
            if (sum >= MOD) sum -= MOD;
        }
        H[n] = sum;
        fact_n = mod_mul(fact_n, n);
    }

    long long result = mod_mul(H[target], fact_n);
    cout << result << '\n';
}

int main() {
    ios_base::sync_with_stdio(false);
    cin.tie(nullptr);
    solve();
    return 0;
}

Python

kMod = 1000000007
kMaxBlock = 5

def solve():
    target = 10000000
    A = [0, 1, 2, 6, 18, 12]

    fact_small = [1] * (kMaxBlock + 1)
    inv_fact_small = [1] * (kMaxBlock + 1)
    coeff = [0] * (kMaxBlock + 1)

    for i in range(1, kMaxBlock + 1):
        fact_small[i] = (fact_small[i - 1] * i) % kMod
        inv_fact_small[i] = pow(fact_small[i], kMod - 2, kMod)
        coeff[i] = (A[i] * inv_fact_small[i]) % kMod

    H = [0] * (target + 1)
    H[0] = 1

    fact_n = 1

    for n in range(1, target + 1):
        total = 0
        limit = min(n, kMaxBlock)
        for j in range(1, limit + 1):
            total += coeff[j] * H[n - j]
            if total >= kMod:
                total %= kMod
        H[n] = total
        fact_n = (fact_n * n) % kMod

    result = (H[target] * fact_n) % kMod
    return str(result)

if __name__ == "__main__":
    print(solve())

Java

public class Euler857 {
    static final long MOD = 1000000007L;
    static final int kMaxBlock = 5;

    static long modPow(long base, long exp) {
        long res = 1;
        base %= MOD;
        while (exp > 0) {
            if ((exp & 1) == 1)
                res = (res * base) % MOD;
            base = (base * base) % MOD;
            exp >>= 1;
        }
        return res;
    }

    static long modInverse(long n) {
        return modPow(n, MOD - 2);
    }

    public static String solve() {
        int target = 10000000;
        long[] A = { 0, 1, 2, 6, 18, 12 };

        long[] factSmall = new long[kMaxBlock + 1];
        long[] invFactSmall = new long[kMaxBlock + 1];
        long[] coeff = new long[kMaxBlock + 1];

        factSmall[0] = 1;
        invFactSmall[0] = 1;
        coeff[0] = 0;

        for (int i = 1; i <= kMaxBlock; ++i) {
            factSmall[i] = (factSmall[i - 1] * i) % MOD;
            invFactSmall[i] = modInverse(factSmall[i]);
            coeff[i] = (A[i] * invFactSmall[i]) % MOD;
        }

        long[] H = new long[target + 1];
        H[0] = 1;
        long factN = 1;

        for (int n = 1; n <= target; ++n) {
            long sum = 0;
            int limit = Math.min(n, kMaxBlock);
            for (int j = 1; j <= limit; ++j) {
                sum += (coeff[j] * H[n - j]) % MOD;
            }
            H[n] = sum % MOD;
            factN = (factN * n) % MOD;
        }

        long result = (H[target] * factN) % MOD;
        return Long.toString(result);
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}