Problem 647: Linear Transformations of Polygonal Numbers

Project Euler Problem 647 Solution

EulerSolve provides an optimized solution for Project Euler Problem 647, Linear Transformations of Polygonal Numbers, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The arithmetic reduction used by the implementation rewrites the polygonal-number question as a sum over odd shape parameters $k\ge 3$. For each odd $k$, define $$\delta=k-2,\qquad \gamma=(k-4)^2.$$ Every admissible configuration is then indexed by an integer $t\ge 1$ and contributes two nonnegative quantities $$X_{k,t}=(2\delta t+1)^2,\qquad Y_{k,t}=\frac{\gamma}{2}\left(\delta t^2+t\right).$$ A pair $(k,t)$ is valid only when both $X_{k,t}\le N$ and $Y_{k,t}\le N$. The total quantity to compute is therefore $$S(N)=\sum_{\substack{k\ge 3\\ k\equiv 1 \!\!\!\pmod 2}}\ \sum_{\substack{t\ge 1\\ X_{k,t}\le N,\ Y_{k,t}\le N}}\left(X_{k,t}+Y_{k,t}\right).$$ The challenge is that $N$ is enormous, so the useful task is to collapse the inner summation to a closed form and to stop the outer loop as soon as no larger odd $k$ can contribute. Mathematical Approach For fixed odd $k$, let $C_k(N)$ be the total contribution of that single value of $k$. The derivation below matches exactly the arithmetic carried out by the implementations. Step 1: Fix One Odd $k$ and Isolate Its Contribution Once $k$ is fixed, the only free parameter is $t\ge 1$. Because $\delta=k-2>0$ and $\gamma=(k-4)^2>0$, the sequences $X_{k,t}$ and $Y_{k,t}$ both grow with $t$. Also, $k$ is odd, so $\delta$ is odd....

Detailed mathematical approach

Problem Summary

The arithmetic reduction used by the implementation rewrites the polygonal-number question as a sum over odd shape parameters $k\ge 3$. For each odd $k$, define

$$\delta=k-2,\qquad \gamma=(k-4)^2.$$

Every admissible configuration is then indexed by an integer $t\ge 1$ and contributes two nonnegative quantities

$$X_{k,t}=(2\delta t+1)^2,\qquad Y_{k,t}=\frac{\gamma}{2}\left(\delta t^2+t\right).$$

A pair $(k,t)$ is valid only when both $X_{k,t}\le N$ and $Y_{k,t}\le N$. The total quantity to compute is therefore

$$S(N)=\sum_{\substack{k\ge 3\\ k\equiv 1 \!\!\!\pmod 2}}\ \sum_{\substack{t\ge 1\\ X_{k,t}\le N,\ Y_{k,t}\le N}}\left(X_{k,t}+Y_{k,t}\right).$$

The challenge is that $N$ is enormous, so the useful task is to collapse the inner summation to a closed form and to stop the outer loop as soon as no larger odd $k$ can contribute.

Mathematical Approach

For fixed odd $k$, let $C_k(N)$ be the total contribution of that single value of $k$. The derivation below matches exactly the arithmetic carried out by the implementations.

Step 1: Fix One Odd $k$ and Isolate Its Contribution

Once $k$ is fixed, the only free parameter is $t\ge 1$. Because $\delta=k-2>0$ and $\gamma=(k-4)^2>0$, the sequences $X_{k,t}$ and $Y_{k,t}$ both grow with $t$.

Also, $k$ is odd, so $\delta$ is odd. Hence

$$\delta t^2+t=t(\delta t+1)$$

is always even: if $t$ is odd then $\delta t$ is odd and $\delta t+1$ is even, while if $t$ is even then the factor $t$ itself is even. Therefore $Y_{k,t}$ is always an integer.

Since both expressions increase with $t$, the feasible values form an initial interval

$$1\le t\le T_k(N),$$

and the contribution for this $k$ is

$$C_k(N)=\sum_{t=1}^{T_k(N)}\left(X_{k,t}+Y_{k,t}\right).$$

Step 2: Derive the First Bound for $t$

Let

$$R=\left\lfloor\sqrt{N}\right\rfloor.$$

The inequality $X_{k,t}\le N$ becomes

$$\left(2\delta t+1\right)^2\le N.$$

All quantities are nonnegative, so this is equivalent to

$$2\delta t+1\le R,$$

which yields the exact upper bound

$$T_A(k,N)=\left\lfloor\frac{R-1}{2\delta}\right\rfloor.$$

Step 3: Derive the Second Bound for $t$

Now use the second feasibility condition $Y_{k,t}\le N$:

$$\frac{\gamma}{2}\left(\delta t^2+t\right)\le N.$$

Since $\gamma>0$, divide by $\gamma$ and define

$$M=\left\lfloor\frac{2N}{\gamma}\right\rfloor.$$

Then $t$ must satisfy

$$\delta t^2+t\le M.$$

This is a quadratic inequality. Solving $\delta t^2+t-M=0$ gives the positive root

$$\frac{-1+\sqrt{1+4\delta M}}{2\delta}.$$

Therefore the exact integer bound is

$$T_B(k,N)=\left\lfloor\frac{\left\lfloor\sqrt{1+4\delta M}\right\rfloor-1}{2\delta}\right\rfloor.$$

The actual usable range is the intersection of the two bounds:

$$T_k(N)=\min\!\bigl(T_A(k,N),T_B(k,N)\bigr).$$

Step 4: Replace the Inner Loop by Closed Forms

Let $T=T_k(N)$. Introduce the standard sums

$$S_1(T)=\sum_{t=1}^{T} t=\frac{T(T+1)}{2},\qquad S_2(T)=\sum_{t=1}^{T} t^2=\frac{T(T+1)(2T+1)}{6}.$$

Expanding $X_{k,t}$ gives

$$X_{k,t}=4\delta^2 t^2+4\delta t+1,$$

hence

$$\sum_{t=1}^{T} X_{k,t}=4\delta^2S_2(T)+4\delta S_1(T)+T.$$

For the second term,

$$\sum_{t=1}^{T} Y_{k,t}=\frac{\gamma}{2}\sum_{t=1}^{T}\left(\delta t^2+t\right)=\frac{\gamma}{2}\left(\delta S_2(T)+S_1(T)\right).$$

So the whole contribution of one odd $k$ is

$$\boxed{C_k(N)=T+4\delta S_1(T)+4\delta^2S_2(T)+\frac{\gamma}{2}\left(\delta S_2(T)+S_1(T)\right).}$$

This is the key optimization: after computing $T_k(N)$, the contribution of that $k$ is obtained in constant time.

Step 5: Prove the Outer Loop Can Stop Early

The smallest admissible parameter is always $t=1$. Therefore a necessary condition for any contribution from a given odd $k$ is

$$X_{k,1}=(2k-3)^2\le N,\qquad Y_{k,1}=\frac{(k-4)^2(k-1)}{2}\le N.$$

If either inequality already fails at $t=1$, then no larger $t$ can work, because both $X_{k,t}$ and $Y_{k,t}$ increase with $t$.

Moreover, both $X_{k,1}$ and $Y_{k,1}$ increase as odd $k$ increases. So once one odd $k$ fails this test, every larger odd $k$ fails as well. That is why the outer loop can terminate immediately instead of checking all remaining values.

Step 6: Worked Example with $N=100$

This small case is large enough to show every ingredient of the method but still easy to check by hand.

For $k=3$, we have $\delta=1$ and $\gamma=1$. Since $\lfloor\sqrt{100}\rfloor=10$,

$$T_A(3,100)=\left\lfloor\frac{10-1}{2}\right\rfloor=4.$$

Also

$$M=\left\lfloor\frac{2\cdot 100}{1}\right\rfloor=200,\qquad 1+4\delta M=801,\qquad \left\lfloor\sqrt{801}\right\rfloor=28,$$

$$T_B(3,100)=\left\lfloor\frac{28-1}{2}\right\rfloor=13.$$

Thus $T_3(100)=4$. With $S_1(4)=10$ and $S_2(4)=30$,

$$C_3(100)=4+4\cdot 10+4\cdot 30+\frac{1}{2}(30+10)=164+20=184.$$

For $k=5$, we have $\delta=3$ and $\gamma=1$. Then

$$T_A(5,100)=\left\lfloor\frac{10-1}{6}\right\rfloor=1,$$

and

$$1+4\delta M=1+4\cdot 3\cdot 200=2401,\qquad \left\lfloor\sqrt{2401}\right\rfloor=49,$$

$$T_B(5,100)=\left\lfloor\frac{49-1}{6}\right\rfloor=8.$$

Hence $T_5(100)=1$, and the contribution is simply

$$C_5(100)=(2\cdot 3\cdot 1+1)^2+\frac{1}{2}(3\cdot 1^2+1)=49+2=51.$$

For $k=7$, the very first square term already fails:

$$X_{7,1}=(2\cdot 7-3)^2=11^2=121>100.$$

So no larger odd $k$ can contribute, and the final total is

$$S(100)=184+51=235.$$

How the Code Works

The C++, Python, and Java implementations all follow the same plan. They iterate only over odd $k\ge 3$, compute the two exact bounds $T_A$ and $T_B$ using integer square roots, and keep the minimum as the last admissible $t$.

If that minimum is $0$, the current odd $k$ contributes nothing. Otherwise the implementation evaluates $S_1(T)$ and $S_2(T)$, substitutes them into the closed-form expressions above, and adds the result to the running total.

No implementation performs a full inner loop over $t$. The only loop that remains is the outer loop over odd $k$, and it stops as soon as the $t=1$ test proves that the current $k$ and all larger odd $k$ are impossible.

Complexity Analysis

Each admissible odd $k$ is processed in $O(1)$ arithmetic time after a constant number of integer-square-root evaluations, and the memory usage is $O(1)$.

Let $K_{\max}$ be the largest odd $k$ that can pass the $t=1$ test. Since

$$Y_{k,1}=\frac{(k-4)^2(k-1)}{2}=\Theta(k^3),$$

we have $K_{\max}=O(N^{1/3})$. The bound coming from $X_{k,1}=(2k-3)^2$ is only $O(\sqrt{N})$, so the cubic term is the one that controls the asymptotic stopping point. Therefore the full running time is

$$O(N^{1/3}),$$

with constant extra memory.

Footnotes and References

Project Euler problem page: https://projecteuler.net/problem=647
Polygonal numbers: Wikipedia — Polygonal number
Quadratic equations and discriminants: Wikipedia — Quadratic equation
Arithmetic progressions and triangular sums: Wikipedia — Arithmetic progression
Sum of squares: Wikipedia — Square pyramidal number

Problem 647 source code

C++

#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <string>

namespace {

using u64 = std::uint64_t;
using i128 = __int128_t;
using u128 = __uint128_t;

u64 isqrt_u64(const u64 n) {
    u64 r = static_cast<u64>(std::sqrt(static_cast<long double>(n)));
    while ((r + 1) * (r + 1) <= n) ++r;
    while (r * r > n) --r;
    return r;
}

u64 floor_sqrt_u128(const u128 n) {
    u64 r = static_cast<u64>(std::sqrt(static_cast<long double>(n)));
    while ((u128)(r + 1) * (r + 1) <= n) ++r;
    while ((u128)r * r > n) --r;
    return r;
}

i128 Fk(const int k, const u64 N) {
    if ((k & 1) == 0 || k < 3) return 0;
    const u64 d = static_cast<u64>(k - 2);
    const u64 c = static_cast<u64>(k - 4) * static_cast<u64>(k - 4);

    const u64 sqrtN = isqrt_u64(N);
    if (sqrtN <= 1) return 0;

    const u64 tA = (sqrtN - 1) / (2 * d);

    const u64 M = (2 * N) / c;
    const u128 disc = 1 + (u128)4 * d * (u128)M;
    const u64 s = floor_sqrt_u128(disc);
    const u64 tB = (s > 0) ? (static_cast<u64>((s - 1) / (2 * d))) : 0;

    const u64 T = std::min(tA, tB);
    if (T == 0) return 0;

    const i128 TT = static_cast<i128>(T);
    const i128 S1 = TT * (TT + 1) / 2;
    const i128 S2 = TT * (TT + 1) * (2 * TT + 1) / 6;

    const i128 dd = static_cast<i128>(d);
    const i128 cc = static_cast<i128>(c);

    const i128 sumA = TT + 4 * dd * S1 + 4 * dd * dd * S2;
    const i128 sum_dt2_t = dd * S2 + S1;
    const i128 sumB = cc * (sum_dt2_t / 2);
    return sumA + sumB;
}

i128 sum_all(const u64 N) {
    i128 ans = 0;
    for (int k = 3;; k += 2) {
        const u64 d = static_cast<u64>(k - 2);
        const u64 c = static_cast<u64>(k - 4) * static_cast<u64>(k - 4);
        const u128 A1 = (u128)(2ULL * (u64)k - 3) * (2ULL * (u64)k - 3);
        const u128 B1 = (u128)c * (u64)(k - 1) / 2;
        if (A1 > N || B1 > N) break;
        ans += Fk(k, N);
    }
    return ans;
}

std::string to_string_i128(i128 x) {
    if (x == 0) return "0";
    bool neg = false;
    if (x < 0) {
        neg = true;
        x = -x;
    }
    std::string s;
    while (x > 0) {
        const int digit = static_cast<int>(x % 10);
        s.push_back(static_cast<char>('0' + digit));
        x /= 10;
    }
    if (neg) s.push_back('-');
    std::reverse(s.begin(), s.end());
    return s;
}

}  // namespace

int main() {
    assert(Fk(3, 100) == 184);
    assert(sum_all(1000) == 14993);

    std::cout << to_string_i128(sum_all(1'000'000'000'000ULL)) << "\n";
    return 0;
}

Python

import math

def isqrt_floor(n):
    if n < 0: return 0
    return math.isqrt(n)

def floor_sqrt(n):
    if n < 0: return 0
    r = math.isqrt(n)
    while (r + 1)**2 <= n: r += 1
    while r**2 > n: r -= 1
    return r

def Fk(k, N):
    if (k & 1) == 0 or k < 3: return 0
    d = k - 2
    c = (k - 4) * (k - 4)
    
    sqrtN = isqrt_floor(N)
    if sqrtN <= 1: return 0
    
    tA = (sqrtN - 1) // (2 * d)
    
    M = (2 * N) // c
    disc = 1 + 4 * d * M
    s = floor_sqrt(disc)
    tB = (s - 1) // (2 * d) if s > 0 else 0
    
    T = min(tA, tB)
    if T == 0: return 0
    
    TT = T
    S1 = TT * (TT + 1) // 2
    S2 = TT * (TT + 1) * (2 * TT + 1) // 6
    
    sumA = TT + 4 * d * S1 + 4 * d * d * S2
    sum_dt2_t = d * S2 + S1
    sumB = c * (sum_dt2_t // 2)
    return sumA + sumB

def sum_all(N):
    ans = 0
    k = 3
    while True:
        d = k - 2
        c = (k - 4) * (k - 4)
        A1 = (2 * k - 3) * (2 * k - 3)
        B1 = c * (k - 1) // 2
        if A1 > N or B1 > N: break
        ans += Fk(k, N)
        k += 2
    return ans

def solve():
    N = 1000000000000
    ans = sum_all(N)
    return str(ans)

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

Java

import java.math.BigInteger;

public class Euler647 {

    static long isqrtU64(long n) {
        if (n < 0)
            return 0;
        long r = (long) Math.sqrt(n);
        while ((r + 1) * (r + 1) <= n)
            r++;
        while (r * r > n)
            r--;
        return r;
    }

    static long floorSqrtU128(BigInteger n) {
        if (n.signum() < 0)
            return 0;
        // BigInteger sqrt
        BigInteger rBig = n.sqrt();
        long r = rBig.longValue();
        return r;
    }

    static BigInteger Fk(int k, long N) {
        if ((k & 1) == 0 || k < 3)
            return BigInteger.ZERO;
        long d = k - 2;
        long c = (long) (k - 4) * (k - 4);

        long sqrtN = isqrtU64(N);
        if (sqrtN <= 1)
            return BigInteger.ZERO;

        long tA = (sqrtN - 1) / (2 * d);

        long M = (2 * N) / c;
        BigInteger bigD = BigInteger.valueOf(d);
        BigInteger bigM = BigInteger.valueOf(M);
        BigInteger disc = BigInteger.ONE.add(BigInteger.valueOf(4).multiply(bigD).multiply(bigM));

        long s = floorSqrtU128(disc);
        long tB = (s > 0) ? (s - 1) / (2 * d) : 0;

        long T = Math.min(tA, tB);
        if (T == 0)
            return BigInteger.ZERO;

        BigInteger TT = BigInteger.valueOf(T);
        BigInteger S1 = TT.multiply(TT.add(BigInteger.ONE)).divide(BigInteger.valueOf(2));
        BigInteger S2 = TT.multiply(TT.add(BigInteger.ONE))
                .multiply(TT.multiply(BigInteger.valueOf(2)).add(BigInteger.ONE)).divide(BigInteger.valueOf(6));

        BigInteger bigC = BigInteger.valueOf(c);

        BigInteger dS1 = bigD.multiply(S1);
        BigInteger d2S2 = bigD.multiply(bigD).multiply(S2);

        BigInteger sumA = TT.add(dS1.multiply(BigInteger.valueOf(4))).add(d2S2.multiply(BigInteger.valueOf(4)));
        BigInteger sum_dt2_t = bigD.multiply(S2).add(S1);
        BigInteger sumB = bigC.multiply(sum_dt2_t.divide(BigInteger.valueOf(2)));

        return sumA.add(sumB);
    }

    static BigInteger sumAll(long N) {
        BigInteger ans = BigInteger.ZERO;
        for (int k = 3;; k += 2) {
            long d = k - 2;
            long c = (long) (k - 4) * (k - 4);
            long A1 = (2L * k - 3) * (2L * k - 3);
            long B1 = c * (long) (k - 1) / 2;

            // In Java, we should check overflow. A1 and B1 fit in long easily since N is up
            // to 10^12
            // If they overflow long, they are definitely > N
            boolean brk = false;
            if (2L * k - 3 > 3000000 || A1 > N)
                brk = true;
            if (B1 > N)
                brk = true;
            if (brk)
                break;

            ans = ans.add(Fk(k, N));
        }
        return ans;
    }

    public static String solve() {
        long N = 1000000000000L;
        BigInteger ans = sumAll(N);
        return ans.toString();
    }

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