Project Euler

Lattice Points Enclosed by Parabola and Line

For integers a and b, define D(a,b) = {(x,y) in R^2: x^2 <= y <= ax + b}. Let L(a,b) count the lattice points in D(a,b). Define S(N) = sum L(a,b) over all pairs (a,b) with |a|, |b| <= N such that A...

Source sync Apr 19, 2026

Problem #0403

Level Level 25

Solved By 422

Languages C++, Python

Answer 18224771

Length 320 words

modular_arithmeticgeometryalgebra

For integers $a$ and $a$, we define $D(a, b)$ as the domain enclosed enclosed by the parabola $y = x^2$ and the line $y = a \cdot x + b$: $D(a, b) = \{(x, y) | x^2 \leq y \leq a \cdot x + b\}$.

$L(a, b)$ is defined as the number of lattice points contained in $D(a, b)$.

For example, $L(1, 2) = 8$ and $L(2, -1) = 1$.

We also define $S(N)$ as the sum of $L(a, b)$ for all the pairs $(a, b)$ such that the area of $D(a, b)$ is a rational number and $|a|,|b| \leq N$.

We can verify that $S(5) = 344$ and $S(100) = 26709528$.

Find $S(10^{12})$. Give your answer mod $10^8$.

Problem 403: Lattice Points Enclosed by Parabola and Line

Mathematical Foundation

Theorem 1 (Rational area condition). The area of $D(a,b)$ is rational if and only if the discriminant $\Delta = a^2 + 4b$ is a non-negative perfect square.

Proof. The parabola $y = x^2$ and line $y = ax + b$ intersect at roots of $x^2 - ax - b = 0$ , namely $r_{1,2} = (a \mp \sqrt{\Delta})/2$ where $\Delta = a^2 + 4b$ . The enclosed area is

\operatorname{Area}(D(a,b)) = \int_{r_1}^{r_2}(ax + b - x^2)\,dx = \frac{(r_2 - r_1)^3}{6} = \frac{\Delta^{3/2}}{6}.

Write $\Delta^{3/2} = \Delta \cdot \Delta^{1/2}$ . If $\Delta = s^2$ for non-negative integer $s$ , then $\Delta^{3/2} = s^3 \in \mathbb{Q}$ . Conversely, if $\Delta > 0$ is not a perfect square, then $\sqrt{\Delta}$ is irrational, so $\Delta^{3/2}/6$ is irrational. The case $\Delta = 0$ gives area $0 \in \mathbb{Q}$ . $\square$

Lemma 1 (Parity constraint). If $\Delta = s^2$ and $b = (\Delta - a^2)/4 = (s^2 - a^2)/4$ is an integer, then $a \equiv s \pmod{2}$ .

Proof. We need $4 \mid (s^2 - a^2) = (s-a)(s+a)$ . If $a$ and $s$ have the same parity, then $(s-a)$ and $(s+a)$ are both even, so $4 \mid (s-a)(s+a)$ . If they have different parity, $(s-a)(s+a)$ is odd, so $4 \nmid (s^2 - a^2)$ . $\square$

Theorem 2 (Closed form for $L$ ). When $\Delta = s^2$ with $s \ge 0$ , the lattice point count depends only on $s$ :

L(s) = \frac{(s+1)(s^2 - s + 6)}{6}.

Proof. The roots $r_1 = (a-s)/2$ and $r_2 = (a+s)/2$ are integers by Lemma 1. For each integer $x \in [r_1, r_2]$ , the count of integers $y$ with $x^2 \le y \le ax + b$ is $ax + b - x^2 + 1$ . Substituting $t = x - r_1$ (so $x = r_1 + t$ and $t$ ranges from $0$ to $s$ ):

ax + b - x^2 = -(x - r_1)(x - r_2) = t(s - t).

Therefore

L(a,b) = \sum_{t=0}^{s}\bigl[t(s-t) + 1\bigr] = \sum_{t=0}^{s} t(s-t) + (s+1) = \frac{s^2(s+1)}{2} - \frac{s(s+1)(2s+1)}{6} + (s+1).

Simplifying: $\frac{s(s+1)}{6}[3s - (2s+1)] + (s+1) = \frac{s(s+1)(s-1)}{6} + (s+1) = \frac{(s+1)(s^2 - s + 6)}{6}$ . $\square$

Lemma 2 (Reparametrization). Setting $a = s + 2u$ for integer $u$ , we get $b = -u(s+u)$ . The constraints become $|s + 2u| \le N$ , $|u(s+u)| \le N$ , $s \ge 0$ .

Proof. From $s^2 = a^2 + 4b$ with $a = s + 2u$ : $b = (s^2 - a^2)/4 = (s^2 - (s+2u)^2)/4 = (-4su - 4u^2)/4 = -u(s+u)$ . $\square$

Theorem 3 (Prefix sum formula). Define $P(K) = \sum_{s=0}^{K} L(s)$ . Then

P(K) = \frac{(K+1)(K+2)(K^2 - K + 12)}{24}.

Proof. By direct summation of $L(s) = (s+1)(s^2 - s + 6)/6$ using the identities $\sum_{s=0}^{K}(s+1) = (K+1)(K+2)/2$ , $\sum_{s=0}^{K}s(s+1) = K(K+1)(K+2)/3$ , and $\sum_{s=0}^{K}s^2(s+1) = K(K+1)(K+2)(3K+1)/12$ . Combining: $P(K) = \frac{1}{6}[\frac{K(K+1)(K+2)(3K+1)}{12} - \frac{K(K+1)(K+2)}{3} + (K+1)(K+2) \cdot 3]$ . Factoring out $(K+1)(K+2)/6$ and simplifying the bracket yields $(K^2 - K + 12)/4$ . $\square$

Theorem 4 (Sub-linear summation). $S(N)$ can be computed in $O(\sqrt{N})$ time by splitting the sum over $u$ into three cases:

$u = 0$ : contribution $P(N)$ .
$u \ge 1$ : direct summation for $u \le \sqrt{N}$ .
$u \le -1$ (set $v = -u$ ): for $v \le \sqrt{N}$ , direct summation; for $v > \sqrt{N}$ , group by $q = \lfloor N/v \rfloor$ , yielding $O(\sqrt{N})$ groups each with $O(1)$ work via Faulhaber’s formulas on degree-4 polynomial sums.

Proof. For $u \ge 1$ : $|u(s+u)| \le N$ with $s \ge 0$ forces $u^2 \le N$ , so $u \le \sqrt{N}$ . For $v = -u > \sqrt{N}$ : the valid $s$ -range has width $O(N/v)$ , and $\lfloor N/v \rfloor$ takes $O(\sqrt{N})$ distinct values. Within each group, $P(v+q) - P(f(v))$ is a degree-4 polynomial in $v$ , summable in $O(1)$ via $\sum v^k$ for $k = 0,\ldots,4$ . $\square$

Editorial

We case u >= 1: direct iteration. We then case v = -u >= 1, v <= sqrt(N): direct iteration. Finally, clamp to v > sqrt(N) range.

Pseudocode

Case u >= 1: direct iteration
Case v = -u >= 1, v <= sqrt(N): direct iteration
Case v > sqrt(N): group by q = floor(N/v)
Clamp to v > sqrt(N) range
Sum P(v + q) - P(max(v - q, 2v - N) - 1) over v in [v_lo, v_hi]
Each term is degree-4 polynomial in v; use Faulhaber sums

Complexity Analysis

Time: $O(\sqrt{N})$ . Cases 1—2 require $O(\sqrt{N})$ iterations. Case 3 has $O(\sqrt{N})$ groups, each $O(1)$ via Faulhaber’s formulas.
Space: $O(1)$ .

Answer

\boxed{18224771}

C++ project_euler/problem_403/solution.cpp

#include <bits/stdc++.h>
using namespace std;

/*
 * Project Euler Problem 403: Lattice Points Enclosed by Parabola and Line
 *
 * Area of D(a,b) = (a^2+4b)^{3/2}/6, rational iff a^2+4b = s^2 (perfect square).
 * L(a,b) depends only on s: L(s) = (s+1)(s^2-s+6)/6.
 * Reparametrize: a = s+2u, b = u(s+u). Sum over (u,s) pairs.
 * Use prefix sum P(K) = (K+1)(K+2)(K^2-K+12)/24 and floor(N/v) grouping
 * for O(sqrt(N)) complexity.
 *
 * Answer: S(10^12) mod 10^8 = 18224771
 */

typedef long long ll;
typedef __int128 lll;

const ll MOD = 100000000LL; // 10^8
const ll BIGMOD = 24LL * MOD; // 2,400,000,000

// Reduce __int128 to [0, m)
ll mmod(lll x, ll m) {
    x %= m;
    if (x < 0) x += m;
    return (ll)x;
}

// P(K) mod MOD using the numerator-mod-BIGMOD trick
ll P_mod(ll K) {
    if (K < 0) return 0;
    lll a = mmod((lll)(K + 1), BIGMOD);
    lll b = mmod((lll)(K + 2), BIGMOD);
    lll c = mmod((lll)K * K - K + 12, BIGMOD);
    ll num = (ll)(a * b % BIGMOD * c % BIGMOD);
    return (num / 24) % MOD;
}

// Prefix power sums mod BIGMOD.
// We use the exact-division trick: compute numerator mod (denom * BIGMOD),
// then integer-divide by denom to get result mod BIGMOD.
ll prefix_pow(ll n, int k) {
    if (n < 0) return 0;
    lll nn = n;
    if (k == 0) {
        return mmod(nn + 1, BIGMOD);
    }
    if (k == 1) {
        // n(n+1)/2
        lll M = (lll)2 * BIGMOD;
        lll val = (lll)mmod(nn, M) * mmod(nn + 1, M) % M;
        return (ll)((val / 2) % BIGMOD);
    }
    if (k == 2) {
        // n(n+1)(2n+1)/6
        lll M = (lll)6 * BIGMOD;
        lll val = (lll)mmod(nn, M) * mmod(nn + 1, M) % M;
        val = val * mmod(2 * nn + 1, M) % M;
        return (ll)((val / 6) % BIGMOD);
    }
    if (k == 3) {
        // [n(n+1)/2]^2
        lll M2 = (lll)2 * BIGMOD;
        lll half = (lll)mmod(nn, M2) * mmod(nn + 1, M2) % M2;
        half = half / 2;
        half %= BIGMOD;
        return (ll)(half * half % BIGMOD);
    }
    if (k == 4) {
        // n(n+1)(2n+1)(3n^2+3n-1)/30
        lll M = (lll)30 * BIGMOD;
        lll a = mmod(nn, M);
        lll b = mmod(nn + 1, M);
        lll c = mmod(2 * nn + 1, M);
        lll d = mmod((lll)3 * nn * nn + 3 * nn - 1, M);
        lll val = a * b % M;
        val = val * c % M;
        val = val * d % M;
        return (ll)((val / 30) % BIGMOD);
    }
    return 0;
}

ll sum_pow_range(ll a, ll b, int k) {
    if (a > b) return 0;
    return (prefix_pow(b, k) - prefix_pow(a - 1, k) + BIGMOD) % BIGMOD;
}

// sum_{v=v1}^{v2} P(coeff*v + offset) mod MOD
//
// P(w) = (w^4 + 2w^3 + 11w^2 + 34w + 24) / 24
// w = c*v + d, where c = coeff, d = offset.
//
// The coefficient of v^j in the numerator of P(cv+d) is:
//   A_j = sum_{deg=j..4} weight[deg] * C(deg,j) * c^j * d^(deg-j)
// where weight = {24, 34, 11, 2, 1} for deg = {0,1,2,3,4}.
//
// We compute A_j mod BIGMOD, then:
//   sum_v P(cv+d) = sum_j A_j * (sum_{v=v1}^{v2} v^j) / 24
// Everything mod BIGMOD, divide by 24 at end.

ll sum_P_linear_mod(ll v1, ll v2, ll coeff, ll offset) {
    if (v1 > v2) return 0;

    // c^j mod BIGMOD for j = 0..4
    ll cpow[5];
    cpow[0] = 1;
    for (int j = 1; j <= 4; j++)
        cpow[j] = mmod((lll)cpow[j-1] * coeff, BIGMOD);

    // d^k mod BIGMOD for k = 0..4
    ll dpow[5];
    dpow[0] = 1;
    for (int j = 1; j <= 4; j++)
        dpow[j] = mmod((lll)dpow[j-1] * offset, BIGMOD);

    // Binomial coefficients (small, precomputed)
    int binom[5][5] = {
        {1, 0, 0, 0, 0},
        {1, 1, 0, 0, 0},
        {1, 2, 1, 0, 0},
        {1, 3, 3, 1, 0},
        {1, 4, 6, 4, 1}
    };

    int weight[5] = {24, 34, 11, 2, 1};

    lll total = 0;
    for (int j = 0; j <= 4; j++) {
        // A_j = sum_{deg=j}^{4} weight[deg] * binom[deg][j] * c^j * d^{deg-j}
        lll Aj = 0;
        for (int deg = j; deg <= 4; deg++) {
            lll term = (lll)weight[deg] * binom[deg][j];
            term = term % BIGMOD * cpow[j] % BIGMOD;
            term = term * dpow[deg - j] % BIGMOD;
            Aj = (Aj + term) % BIGMOD;
        }

        ll sp = sum_pow_range(v1, v2, j);
        total = (total + Aj * sp) % BIGMOD;
    }

    return (ll)((total / 24) % MOD);
}

int main() {
    ios_base::sync_with_stdio(false);
    cin.tie(NULL);

    auto S_compute = [&](ll N) -> ll {
        ll total = 0;
        ll sq = (ll)sqrtl((long double)N);
        while ((sq + 1) * (sq + 1) <= N) sq++;
        while (sq * sq > N) sq--;

        // Part 1: u = 0
        total = (total + P_mod(N)) % MOD;

        // Part 2: u = 1 .. sq
        for (ll u = 1; u <= sq; u++) {
            ll s_hi_b = (N - u * u) / u;
            ll s_hi_a = N - 2 * u;
            ll s_hi = min(s_hi_b, s_hi_a);
            if (s_hi >= 0) {
                total = (total + P_mod(s_hi)) % MOD;
            }
        }

        // Part 3: v = 1 .. sq (u = -v)
        for (ll v = 1; v <= sq; v++) {
            ll s_hi = v + N / v;
            ll s_lo = max(0LL, 2 * v - N);
            if (s_hi >= s_lo) {
                total = (total + P_mod(s_hi) - P_mod(s_lo - 1) + MOD) % MOD;
            }
        }

        // Part 4: v = sq+1 .. N, grouped by q = N/v
        ll v = sq + 1;
        while (v <= N) {
            ll q = N / v;
            if (q == 0) break;
            ll v_end = min(N / q, N);
            ll v_cut = N - q;

            ll v_start = v;
            ll v_end_A = min(v_cut, v_end);
            ll v_start_B = max(v_start, v_cut + 1);
            ll v_end_B = v_end;

            if (v_start <= v_end_A) {
                ll add = (sum_P_linear_mod(v_start, v_end_A, 1, q)
                        - sum_P_linear_mod(v_start, v_end_A, 1, -q - 1) + MOD) % MOD;
                total = (total + add) % MOD;
            }

            if (v_start_B <= v_end_B) {
                ll add = (sum_P_linear_mod(v_start_B, v_end_B, 1, q)
                        - sum_P_linear_mod(v_start_B, v_end_B, 2, -N - 1) + MOD) % MOD;
                total = (total + add) % MOD;
            }

            v = v_end + 1;
        }

        return total;
    };

    // Verify
    printf("S(5) mod 10^8 = %lld (expected 344)\n", S_compute(5));
    printf("S(100) mod 10^8 = %lld (expected 26709528)\n", S_compute(100));
    printf("S(1000) mod 10^8 = %lld (expected %lld)\n", S_compute(1000), 263967605900LL % MOD);

    // Main computation
    ll N = 1000000000000LL; // 10^12
    auto t0 = chrono::high_resolution_clock::now();
    ll ans = S_compute(N);
    auto t1 = chrono::high_resolution_clock::now();
    double elapsed = chrono::duration<double>(t1 - t0).count();

    printf("\nS(10^12) mod 10^8 = %lld\n", ans);
    printf("Time: %.3f s\n", elapsed);

    return 0;
}

Python project_euler/problem_403/solution.py

#!/usr/bin/env python3
"""
Project Euler Problem 403: Lattice points enclosed by parabola and line

For integers a and b, the domain D(a,b) = {(x,y) | x^2 <= y <= ax + b}.
L(a,b) counts the lattice points in D(a,b).
S(N) = sum of L(a,b) for all (a,b) with |a|,|b| <= N such that
       the area of D(a,b) is rational.

Find S(10^12) mod 10^8.

Key mathematical insights:
1. Area of D(a,b) = (a^2+4b)^(3/2) / 6, which is rational iff a^2+4b is a perfect square.
2. Let a^2+4b = s^2 (s >= 0). Then L(a,b) depends only on s:
   L(s) = (s+1)(s^2 - s + 6) / 6
3. Substituting a = s + 2u gives b = u(s+u), reducing the problem to counting
   valid (u,s) pairs and summing L(s).
4. Using prefix sums P(K) = sum_{s=0}^{K} L(s) = (K+1)(K+2)(K^2-K+12)/24,
   each group of u values contributes a closed-form expression.
5. The floor(N/v) grouping trick reduces the O(N) loop to O(sqrt(N)).
"""

import math
from math import comb
import time

# ---------------------------------------------------------------------------
# Core mathematical functions
# ---------------------------------------------------------------------------

def L(s):
    """Number of lattice points for discriminant root s."""
    return (s + 1) * (s * s - s + 6) // 6

def P(K):
    """Prefix sum: sum_{s=0}^{K} L(s)."""
    if K < 0:
        return 0
    return (K + 1) * (K + 2) * (K * K - K + 12) // 24

def sum_pow(a, b, k):
    """Exact sum of v^k for v = a, a+1, ..., b."""
    if a > b:
        return 0

    def prefix(n):
        if n < 0:
            return 0
        if k == 0:
            return n + 1
        elif k == 1:
            return n * (n + 1) // 2
        elif k == 2:
            return n * (n + 1) * (2 * n + 1) // 6
        elif k == 3:
            t = n * (n + 1) // 2
            return t * t
        elif k == 4:
            return n * (n + 1) * (2 * n + 1) * (3 * n * n + 3 * n - 1) // 30

    return prefix(b) - prefix(a - 1)

def sum_P_linear(v1, v2, coeff, offset):
    """Compute sum_{v=v1}^{v2} P(coeff*v + offset) exactly.

    P(w) = (w^4 + 2w^3 + 11w^2 + 34w + 24) / 24.
    Substituting w = coeff*v + offset and expanding gives a degree-4
    polynomial in v whose partial sums are computable in O(1).
    """
    if v1 > v2:
        return 0
    c, d = coeff, offset

    def wk_coeff(k, j):
        if j > k:
            return 0
        return comb(k, j) * c ** j * d ** (k - j)

    total = 0
    for j in range(5):
        cj = (wk_coeff(4, j) + 2 * wk_coeff(3, j) +
              11 * wk_coeff(2, j) + 34 * wk_coeff(1, j) +
              24 * (1 if j == 0 else 0))
        total += cj * sum_pow(v1, v2, j)

    return total // 24

# ---------------------------------------------------------------------------
# Brute-force solver (for verification, small N only)
# ---------------------------------------------------------------------------

def L_bruteforce(a, b):
    """Count lattice points (x,y) with x^2 <= y <= ax + b."""
    disc = a * a + 4 * b
    if disc < 0:
        return 0
    sq = int(math.isqrt(disc))
    if sq * sq != disc:
        return 0
    r1 = (a - sq) // 2
    r2 = (a + sq) // 2
    count = 0
    for x in range(r1, r2 + 1):
        gap = a * x + b - x * x
        if gap >= 0:
            count += gap + 1
    return count

def S_bruteforce(N):
    """Brute-force S(N): enumerate all (a,b) with |a|,|b| <= N."""
    total = 0
    for a in range(-N, N + 1):
        for b in range(-N, N + 1):
            total += L_bruteforce(a, b)
    return total

# ---------------------------------------------------------------------------
# O(N) solver
# ---------------------------------------------------------------------------

def S_linear(N):
    """O(N) algorithm using the u-parametrization."""
    total = 0
    sq = int(math.isqrt(N))

    # u = 0: s in [0, N]
    total += P(N)

    # u >= 1: a = s + 2u, b = u(s+u), need |a| <= N and |b| <= N
    for u in range(1, sq + 1):
        s_hi = min((N - u * u) // u, N - 2 * u)
        if s_hi >= 0:
            total += P(s_hi)

    # v >= 1 (u = -v): a = s - 2v, b = -v(s-v)
    for v in range(1, sq + 1):
        s_hi = v + N // v
        s_lo = max(0, 2 * v - N)
        if s_hi >= s_lo:
            total += P(s_hi) - P(s_lo - 1)

    for v in range(sq + 1, N + 1):
        q = N // v
        s_hi = v + q
        s_lo = max(v - q, 2 * v - N, 0)
        if s_lo <= s_hi:
            total += P(s_hi) - P(s_lo - 1)

    return total

# ---------------------------------------------------------------------------
# O(sqrt(N)) solver
# ---------------------------------------------------------------------------

def S_fast(N):
    """O(sqrt(N)) algorithm using floor(N/v) grouping."""
    total = 0
    sq = int(math.isqrt(N))

    # Part 1: u = 0 => s in [0, N]
    total += P(N)

    # Part 2: u = 1 .. sqrt(N)
    for u in range(1, sq + 1):
        s_hi = min((N - u * u) // u, N - 2 * u)
        if s_hi >= 0:
            total += P(s_hi)

    # Part 3: v = 1 .. sqrt(N) (u = -v)
    for v in range(1, sq + 1):
        s_hi = v + N // v
        s_lo = max(0, 2 * v - N)
        if s_hi >= s_lo:
            total += P(s_hi) - P(s_lo - 1)

    # Part 4: v = sqrt(N)+1 .. N, grouped by q = floor(N/v)
    v = sq + 1
    while v <= N:
        q = N // v
        if q == 0:
            break
        v_end = min(N // q, N)
        v_cut = N - q

        v_start = v
        v_end_A = min(v_cut, v_end)
        v_start_B = max(v_start, v_cut + 1)
        v_end_B = v_end

        # Sub-range A: s_lo = v - q, s_hi = v + q
        if v_start <= v_end_A:
            total += (sum_P_linear(v_start, v_end_A, 1, q) -
                      sum_P_linear(v_start, v_end_A, 1, -q - 1))

        # Sub-range B: s_lo = 2v - N, s_hi = v + q
        if v_start_B <= v_end_B:
            total += (sum_P_linear(v_start_B, v_end_B, 1, q) -
                      sum_P_linear(v_start_B, v_end_B, 2, -N - 1))

        v = v_end + 1

    return total

# ---------------------------------------------------------------------------
# Verification and main computation
# ---------------------------------------------------------------------------

def verify():
    """Verify all implementations against known values."""
    print("=" * 60)
    print("Verification")
    print("=" * 60)

    # Check L examples
    assert L_bruteforce(1, 2) == 8, f"L(1,2) = {L_bruteforce(1,2)}, expected 8"
    assert L_bruteforce(2, -1) == 1, f"L(2,-1) = {L_bruteforce(2,-1)}, expected 1"
    print("L(1,2) = 8, L(2,-1) = 1  [OK]")

    # Check formula L(s) = (s+1)(s^2-s+6)/6
    assert L(3) == 8   # L(1,2): disc=9, s=3
    assert L(0) == 1   # L(2,-1): disc=0, s=0
    print("L formula matches brute force  [OK]")

    # Check S examples
    for N, expected in [(5, 344), (100, 26709528)]:
        bf = S_bruteforce(N) if N <= 100 else None
        lin = S_linear(N)
        fast = S_fast(N)
        if bf is not None:
            assert bf == expected, f"S_bruteforce({N}) = {bf}, expected {expected}"
        assert lin == expected, f"S_linear({N}) = {lin}, expected {expected}"
        assert fast == expected, f"S_fast({N}) = {fast}, expected {expected}"
        print(f"S({N}) = {expected}  [OK]")

    # Cross-check larger values
    for N in [1000, 5000]:
        lin = S_linear(N)
        fast = S_fast(N)
        assert lin == fast, f"Mismatch at N={N}: linear={lin}, fast={fast}"
        print(f"S({N}) = {fast}  [OK, linear and fast agree]")

    print()

def main():
    verify()

    print("=" * 60)
    print("Computing S(10^12) mod 10^8")
    print("=" * 60)

    N = 10 ** 12
    MOD = 10 ** 8

    t0 = time.time()
    result = S_fast(N) % MOD
    elapsed = time.time() - t0

    print(f"S(10^12) mod 10^8 = {result}")
    print(f"Time: {elapsed:.2f}s")

# ---------------------------------------------------------------------------