Project Euler

Lissajous Curves

For coprime positive integers a and b, define the curve x=cos(at), y=cos (b(t-(pi)/(10))), 0 <= t <= 2pi. Let d(a,b) be the sum of x^2+y^2 over all self-intersection points of the curve, and let s(...

Source sync Apr 19, 2026

Problem #0777

Level Level 38

Solved By 153

Languages C++, Python

Answer 2.533018434e23

Length 670 words

number_theorymodular_arithmeticlinear_algebra

For coprime positive integers $a$ and $b$, let $C_{a,b}$ be the curve defined by: \[ \begin {cases} x = \cos \left (at\right ) \\ y = \cos \left (b\left (t-\frac {\pi }{10}\right )\right ) \end {cases} \] where $t$ varies between 0 and $2\pi $.

For example, the images below show $C_{2,5}$ (left) and $C_{7,4}$ (right):

Define $d(a,b) = \displaystyle \sum (x^2 + y^2)$, where the sum is over all points $(x, y)$ at which $C_{a,b}$ crosses itself.

For example, in the case of $C_{2,5}$ illustrated above, the curve crosses itself at two points: $(0.31, 0)$ and $(-0.81, 0)$, rounding coordinates to two decimal places, yielding $d(2, 5)=0.75$. Some other examples are $d(2,3)=4.5$, $d(7,4)=39.5$, $d(7,5)=52$, and $d(10,7)=23.25$.

Let $s(m) = \displaystyle \sum d(a,b)$, where this sum is over all pairs of coprime integers $a,b$ with $2\le a\le m$ and $2\le b\le m$.

You are given that $s(10) = 1602.5$ and $s(100) = 24256505$.

Find $s(10^6)$. Give your answer in scientific notation rounded to $10$ significant digits; for example $s(100)$ would be given as $2.425650500e7$.

Problem 777: Lissajous Curves

1. When Do Two Parameters Give the Same Point?

Suppose t_1 != t_2 and the curve takes the same value at both parameters. Since

\cos u = \cos v \iff u \equiv \pm v \pmod{2\pi},

we must have

a t_1 \equiv \pm a t_2 \pmod{2\pi}, \qquad b\left(t_1-\frac{\pi}{10}\right)\equiv \pm b\left(t_2-\frac{\pi}{10}\right)\pmod{2\pi}.

The (+,+) case gives t_1 \equiv t_2 (mod 2\pi) because gcd(a,b)=1, so it is trivial.

The two nontrivial cases are therefore:

Case B

a(t_1-t_2)=2\pi m,\qquad b\left(t_1+t_2-\frac{\pi}{5}\right)=2\pi n.

Case C

a(t_1+t_2)=2\pi m,\qquad b(t_1-t_2)=2\pi n.

The (-,-) case will matter only when the whole curve is retraced; that happens exactly when 10 | ab.

2. Coordinates of the Candidate Intersections

Let m,n be the integers from Cases B and C.

For Case B, after solving for t_1,t_2 and simplifying, the point satisfies

x^2 = \cos^2\!\left(\pi\left(\frac{a}{10} + \frac{an}{b}\right)\right), \qquad y^2 = \cos^2\!\left(\pi\frac{bm}{a}\right).

Because gcd(a,b)=1, multiplication by a permutes residues modulo b, and multiplication by b permutes residues modulo a. Hence the Case-B contributions are exactly

x^2 \in \left\{ \cos^2\!\left(\pi\left(\frac{a}{10} + \frac{k}{b}\right)\right) : 0 \le k < b \right\},

and

y^2 \in \left\{ \cos^2\!\left(\pi\frac{j}{a}\right) : 1 \le j < a \right\}.

Similarly, Case C gives

x^2 \in \left\{ \cos^2\!\left(\pi\frac{k}{b}\right) : 1 \le k < b \right\},

and

y^2 \in \left\{ \cos^2\!\left(\pi\left(\frac{j}{a} - \frac{b}{10}\right)\right) : 0 \le j < a \right\}.

3. Which Repeated Points Are Genuine Crossings?

Different parameter values can also describe the same local branch, which does not count as a crossing.

In Case B the two tangent vectors are

(x'(t_1),y'(t_1)) \quad\text{and}\quad (x'(t_1),-y'(t_1)),

while in Case C they are

(x'(t_1),y'(t_1)) \quad\text{and}\quad (-x'(t_1),y'(t_1)).

So a repeated point is a genuine crossing exactly when neither derivative component vanishes, i.e.

x^2 \ne 1 \qquad\text{and}\qquad y^2 \ne 1.

This is the only exclusion we need.

4. Formula for `d(a,b)`

4.1 The Generic Case `10 ∤ ab`

If 10 ∤ ab, then the curve is not retraced, and every valid Case-B or Case-C point is counted exactly once.

Using

\sum_{j=1}^{a-1}\cos^2\!\left(\pi\frac{j}{a}\right)=\frac{a}{2}-1, \qquad \sum_{k=0}^{b-1}\cos^2\!\left(\alpha+\pi\frac{k}{b}\right)=\frac{b}{2},

we get

\text{Case B} = b\left(\frac{a}{2}-1\right) + (a-1)\frac{b}{2} = ab - \frac{3b}{2},

and

\text{Case C} = a\left(\frac{b}{2}-1\right) + (b-1)\frac{a}{2} = ab - \frac{3a}{2}.

Therefore

d(a,b)=2ab-\frac{3(a+b)}{2} \qquad\text{if } 10 \nmid ab.

4.2 The Retraced Case `10 | ab`

If 10 | ab, the omitted (-,-) congruence has a solution, so the whole curve is traced twice. Then every non-tangential crossing coming from Cases B and C is counted four times in the raw parameter count.

The only tangential raw points are the ones with x^2=1 in Case B and the ones with y^2=1 in Case C. Their raw total is

\left((a-1)+\left(\frac{a}{2}-1\right)\right) + \left((b-1)+\left(\frac{b}{2}-1\right)\right) = \frac{3(a+b)}{2} - 4.

d(a,b) = \frac{1}{4} \left( 2ab-\frac{3(a+b)}{2} -\left(\frac{3(a+b)}{2}-4\right) \right) = \frac{2ab-3a-3b+4}{4}.

Hence:

d(a,b)= \begin{cases} 2ab-\dfrac{3(a+b)}{2}, & 10 \nmid ab,\\[1.2ex] \dfrac{2ab-3a-3b+4}{4}, & 10 \mid ab. \end{cases}

This matches all examples in the statement.

5. Reducing `s(m)` to Arithmetic Sums

Let, over coprime pairs 2 <= a,b <= m,

C = \sum 1,\qquad A = \sum a,\qquad M = \sum ab.

Let C_{10}, A_{10}, M_{10} be the same sums restricted to 10 | ab.

Because the summation is symmetric in a and b, we also have \sum b = A.

From the formula for d(a,b) above,

4s(m)=8M-12A+4C_{10}-6M_{10}+6A_{10}.

So it remains to compute these six arithmetic totals.

6. Möbius Inversion

For any filter P(n), define:

c_P(d) = number of n <= m such that d | n and P(n) holds,
s_P(d) = sum of those n.

Then over 1 <= a,b <= m,

\sum_{\substack{P(a),\,P(b)\\ \gcd(a,b)=1}} 1 = \sum_{d=1}^{m} \mu(d)\, c_P(d)^2,

\sum_{\substack{P(a),\,P(b)\\ \gcd(a,b)=1}} a = \sum_{d=1}^{m} \mu(d)\, s_P(d)\, c_P(d),

\sum_{\substack{P(a),\,P(b)\\ \gcd(a,b)=1}} ab = \sum_{d=1}^{m} \mu(d)\, s_P(d)^2.

After computing these on 1 <= a,b <= m, we subtract the boundary pairs with a=1 or b=1, because the problem starts at 2.

We use four filters:

all
odd
non5 = numbers not divisible by 5
coprime10 = numbers coprime to 10

Their c_P(d) and s_P(d) are closed-form:

all: plain counts and arithmetic sums of multiples of d
odd: only odd multiples of d
non5: multiples of d excluding those divisible by 5
coprime10: odd multiples of d excluding odd multiples of 5

Finally,

1_{10 \mid ab} = 1 - 1_{\text{both odd}} - 1_{\text{both non5}} + 1_{\text{both coprime10}},

X_{10}=X_{\text{all}}-X_{\text{odd}}-X_{\text{non5}}+X_{\text{coprime10}}

for each of X in {C,A,M}.

This yields s(m) in O(m) time after a linear sieve for \mu.

Correctness

Theorem. The algorithm returns the correct value of s(m).

Proof. Sections 1 through 4 prove the exact closed form for d(a,b), separating the cases 10 \nmid ab and 10 \mid ab. Section 5 rewrites s(m) in terms of the six arithmetic totals C, A, M, C_{10}, A_{10}, M_{10}. Section 6 computes each of these exactly by Möbius inversion together with inclusion-exclusion for the condition 10 \mid ab. Therefore every summand is exact, so the final value of s(m) is exact as well. $\square$

7. Complexity Analysis

Möbius sieve: O(m)
The four filtered divisor sums: O(m)
Memory: O(m)

For m = 10^6, this is easily feasible.

Answer

The computation gives

\boxed{2.533018434e23}

C++ project_euler/problem_777/solution.cpp

#include <cmath>
#include <cassert>
#include <iomanip>
#include <iostream>
#include <sstream>
#include <string>
#include <vector>

namespace {

constexpr int TARGET = 1'000'000;
const std::string TARGET_ANSWER = "2.533018434e23";

using i128 = __int128_t;

struct Stats {
    i128 count;
    i128 sum_first;
    i128 sum_product;
};

i128 triangular(i128 n) {
    return n * (n + 1) / 2;
}

std::vector<int> mobius_sieve(int limit) {
    std::vector<int> mu(limit + 1);
    std::vector<int> primes;
    std::vector<bool> composite(limit + 1, false);
    mu[1] = 1;

    for (int value = 2; value <= limit; ++value) {
        if (!composite[value]) {
            primes.push_back(value);
            mu[value] = -1;
        }
        for (int prime : primes) {
            long long next = 1LL * value * prime;
            if (next > limit) {
                break;
            }
            composite[static_cast<int>(next)] = true;
            if (value % prime == 0) {
                mu[static_cast<int>(next)] = 0;
                break;
            }
            mu[static_cast<int>(next)] = -mu[value];
        }
    }

    return mu;
}

std::pair<i128, i128> count_and_sum(const std::string& filter_name, int limit, int divisor) {
    const i128 n = limit / divisor;

    if (filter_name == "all") {
        return {n, static_cast<i128>(divisor) * triangular(n)};
    }

    if (filter_name == "odd") {
        if (divisor % 2 == 0) {
            return {0, 0};
        }
        const i128 odd_count = (n + 1) / 2;
        return {odd_count, static_cast<i128>(divisor) * odd_count * odd_count};
    }

    if (filter_name == "non5") {
        if (divisor % 5 == 0) {
            return {0, 0};
        }
        return {
            n - n / 5,
            static_cast<i128>(divisor) * (triangular(n) - 5 * triangular(n / 5)),
        };
    }

    if (filter_name == "coprime10") {
        if (divisor % 2 == 0 || divisor % 5 == 0) {
            return {0, 0};
        }
        const i128 odd_count = (n + 1) / 2;
        const i128 odd_multiple_5_count = (n / 5 + 1) / 2;
        const i128 count = odd_count - odd_multiple_5_count;
        const i128 total =
            static_cast<i128>(divisor)
            * (odd_count * odd_count - 5 * odd_multiple_5_count * odd_multiple_5_count);
        return {count, total};
    }

    assert(false && "unknown filter");
    return {0, 0};
}

Stats pair_stats(int limit, const std::string& filter_name, const std::vector<int>& mu) {
    i128 count_full = 0;
    i128 sum_first_full = 0;
    i128 sum_product_full = 0;

    for (int divisor = 1; divisor <= limit; ++divisor) {
        const int mobius = mu[divisor];
        if (mobius == 0) {
            continue;
        }
        const auto [count, total] = count_and_sum(filter_name, limit, divisor);
        if (count == 0) {
            continue;
        }
        count_full += static_cast<i128>(mobius) * count * count;
        sum_first_full += static_cast<i128>(mobius) * total * count;
        sum_product_full += static_cast<i128>(mobius) * total * total;
    }

    const auto [total_count, total_sum] = count_and_sum(filter_name, limit, 1);
    return {
        count_full - (2 * total_count - 1),
        sum_first_full - (total_count + total_sum - 1),
        sum_product_full - (2 * total_sum - 1),
    };
}

i128 solve_quarters(int limit = TARGET) {
    const std::vector<int> mu = mobius_sieve(limit);

    const Stats all_stats = pair_stats(limit, "all", mu);
    const Stats odd_stats = pair_stats(limit, "odd", mu);
    const Stats non5_stats = pair_stats(limit, "non5", mu);
    const Stats coprime10_stats = pair_stats(limit, "coprime10", mu);

    const i128 count_10 = all_stats.count - odd_stats.count - non5_stats.count
                        + coprime10_stats.count;
    const i128 sum_first_10 = all_stats.sum_first - odd_stats.sum_first - non5_stats.sum_first
                            + coprime10_stats.sum_first;
    const i128 sum_product_10 =
        all_stats.sum_product - odd_stats.sum_product - non5_stats.sum_product
        + coprime10_stats.sum_product;

    return 8 * all_stats.sum_product
         - 12 * all_stats.sum_first
         + 4 * count_10
         - 6 * sum_product_10
         + 6 * sum_first_10;
}

long double to_long_double(i128 value) {
    return static_cast<long double>(value);
}

long double d_formula(int a, int b) {
    if ((1LL * a * b) % 10 == 0) {
        return static_cast<long double>(2LL * a * b - 3LL * a - 3LL * b + 4) / 4.0L;
    }
    return static_cast<long double>(4LL * a * b - 3LL * a - 3LL * b) / 2.0L;
}

long double brute_sum(int limit) {
    long double total = 0.0L;
    for (int a = 2; a <= limit; ++a) {
        for (int b = 2; b <= limit; ++b) {
            int x = a;
            int y = b;
            while (y != 0) {
                const int temp = x % y;
                x = y;
                y = temp;
            }
            if (x == 1) {
                total += d_formula(a, b);
            }
        }
    }
    return total;
}

std::string format_quarters(i128 quarters, int digits_after_decimal = 9) {
    const long double value = to_long_double(quarters) / 4.0L;
    std::ostringstream out;
    out << std::scientific << std::setprecision(digits_after_decimal) << value;
    const std::string text = out.str();
    const std::size_t split = text.find('e');
    const std::string mantissa = text.substr(0, split);
    const int exponent = std::stoi(text.substr(split + 1));
    return mantissa + "e" + std::to_string(exponent);
}

bool nearly_equal(long double a, long double b, long double eps = 1e-12L) {
    return fabsl(a - b) <= eps;
}

void self_test() {
    assert(nearly_equal(d_formula(2, 3), 4.5L));
    assert(nearly_equal(d_formula(2, 5), 0.75L));
    assert(nearly_equal(d_formula(7, 4), 39.5L));
    assert(nearly_equal(d_formula(7, 5), 52.0L));
    assert(nearly_equal(d_formula(10, 7), 23.25L));

    for (int limit : {3, 4, 5, 10, 20}) {
        assert(nearly_equal(to_long_double(solve_quarters(limit)) / 4.0L, brute_sum(limit)));
    }

    assert(solve_quarters(10) == static_cast<i128>(6410));
    assert(solve_quarters(100) == static_cast<i128>(97026020));
    assert(format_quarters(solve_quarters()) == TARGET_ANSWER);
}

}  // namespace

int main() {
    self_test();
    std::cout << format_quarters(solve_quarters()) << '\n';
    return 0;
}

Python project_euler/problem_777/solution.py

"""Project Euler 777: Lissajous Curves."""

from __future__ import annotations

from dataclasses import dataclass
from decimal import Decimal, localcontext
from math import gcd, pi

TARGET = 10**6
TARGET_ANSWER = "2.533018434e23"


def mobius_sieve(limit: int) -> list[int]:
    mu = [0] * (limit + 1)
    mu[1] = 1
    primes: list[int] = []
    composite = [False] * (limit + 1)
    for value in range(2, limit + 1):
        if not composite[value]:
            primes.append(value)
            mu[value] = -1
        for prime in primes:
            nxt = value * prime
            if nxt > limit:
                break
            composite[nxt] = True
            if value % prime == 0:
                mu[nxt] = 0
                break
            mu[nxt] = -mu[value]
    return mu


def triangular(n: int) -> int:
    return n * (n + 1) // 2


@dataclass(frozen=True)
class Stats:
    count: int
    sum_first: int
    sum_product: int


def count_and_sum(filter_name: str, limit: int, divisor: int) -> tuple[int, int]:
    n = limit // divisor

    if filter_name == "all":
        return n, divisor * triangular(n)

    if filter_name == "odd":
        if divisor % 2 == 0:
            return 0, 0
        odd_count = (n + 1) // 2
        return odd_count, divisor * odd_count * odd_count

    if filter_name == "non5":
        if divisor % 5 == 0:
            return 0, 0
        return n - n // 5, divisor * (triangular(n) - 5 * triangular(n // 5))

    if filter_name == "coprime10":
        if divisor % 2 == 0 or divisor % 5 == 0:
            return 0, 0
        odd_count = (n + 1) // 2
        odd_multiple_5_count = ((n // 5) + 1) // 2
        count = odd_count - odd_multiple_5_count
        total = divisor * (
            odd_count * odd_count - 5 * odd_multiple_5_count * odd_multiple_5_count
        )
        return count, total

    raise ValueError(f"Unknown filter: {filter_name}")


def pair_stats(limit: int, filter_name: str, mu: list[int]) -> Stats:
    count_full = 0
    sum_first_full = 0
    sum_product_full = 0

    for divisor in range(1, limit + 1):
        mobius = mu[divisor]
        if mobius == 0:
            continue
        count, total = count_and_sum(filter_name, limit, divisor)
        if count == 0:
            continue
        count_full += mobius * count * count
        sum_first_full += mobius * total * count
        sum_product_full += mobius * total * total

    total_count, total_sum = count_and_sum(filter_name, limit, 1)
    return Stats(
        count=count_full - (2 * total_count - 1),
        sum_first=sum_first_full - (total_count + total_sum - 1),
        sum_product=sum_product_full - (2 * total_sum - 1),
    )


def solve_quarters(limit: int = TARGET) -> int:
    """Return 4 * s(limit) exactly."""
    mu = mobius_sieve(limit)
    stats = {
        name: pair_stats(limit, name, mu)
        for name in ("all", "odd", "non5", "coprime10")
    }

    all_stats = stats["all"]
    odd_stats = stats["odd"]
    non5_stats = stats["non5"]
    coprime10_stats = stats["coprime10"]

    count_10 = (
        all_stats.count
        - odd_stats.count
        - non5_stats.count
        + coprime10_stats.count
    )
    sum_first_10 = (
        all_stats.sum_first
        - odd_stats.sum_first
        - non5_stats.sum_first
        + coprime10_stats.sum_first
    )
    sum_product_10 = (
        all_stats.sum_product
        - odd_stats.sum_product
        - non5_stats.sum_product
        + coprime10_stats.sum_product
    )

    return (
        8 * all_stats.sum_product
        - 12 * all_stats.sum_first
        + 4 * count_10
        - 6 * sum_product_10
        + 6 * sum_first_10
    )


def format_quarters(quarters: int, digits_after_decimal: int = 9) -> str:
    with localcontext() as ctx:
        ctx.prec = 50
        value = Decimal(quarters) / Decimal(4)
    mantissa, exponent = f"{value:.{digits_after_decimal}e}".split("e")
    return f"{mantissa}e{int(exponent)}"


def solve(limit: int = TARGET) -> str:
    return format_quarters(solve_quarters(limit))


def d_formula(a: int, b: int) -> Decimal:
    if (a * b) % 10 == 0:
        return Decimal(2 * a * b - 3 * a - 3 * b + 4) / Decimal(4)
    return Decimal(4 * a * b - 3 * a - 3 * b) / Decimal(2)


def brute_sum(limit: int) -> Decimal:
    total = Decimal(0)
    for a in range(2, limit + 1):
        for b in range(2, limit + 1):
            if gcd(a, b) == 1:
                total += d_formula(a, b)
    return total


def self_test() -> None:
    assert d_formula(2, 3) == Decimal("4.5")
    assert d_formula(2, 5) == Decimal("0.75")
    assert d_formula(7, 4) == Decimal("39.5")
    assert d_formula(7, 5) == Decimal("52")
    assert d_formula(10, 7) == Decimal("23.25")

    for limit in (3, 4, 5, 10, 20):
        assert Decimal(solve_quarters(limit)) / Decimal(4) == brute_sum(limit)

    assert Decimal(solve_quarters(10)) / Decimal(4) == Decimal("1602.5")
    assert Decimal(solve_quarters(100)) / Decimal(4) == Decimal("24256505")
    assert solve() == TARGET_ANSWER


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

Lissajous Curves

Problem Statement

Problem 777: Lissajous Curves

1. When Do Two Parameters Give the Same Point?

Case B

Case C

2. Coordinates of the Candidate Intersections

3. Which Repeated Points Are Genuine Crossings?

4. Formula for `d(a,b)`

4.1 The Generic Case `10 ∤ ab`

4.2 The Retraced Case `10 | ab`

5. Reducing `s(m)` to Arithmetic Sums

6. Möbius Inversion

Correctness

7. Complexity Analysis

Answer

Code

Problem Statement

Problem 777: Lissajous Curves

1. When Do Two Parameters Give the Same Point?

Case B

Case C

2. Coordinates of the Candidate Intersections

3. Which Repeated Points Are Genuine Crossings?

4. Formula for d(a,b)

4.1 The Generic Case 10 ∤ ab

4.2 The Retraced Case 10 | ab

5. Reducing s(m) to Arithmetic Sums

6. Möbius Inversion

Correctness

7. Complexity Analysis

Answer

Code

4. Formula for `d(a,b)`

4.1 The Generic Case `10 ∤ ab`

4.2 The Retraced Case `10 | ab`

5. Reducing `s(m)` to Arithmetic Sums