Problem 723: Pythagorean Quadrilaterals

View on Project Euler

Project Euler Problem 723 Solution

EulerSolve provides an optimized solution for Project Euler Problem 723, Pythagorean Quadrilaterals, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The implementations start from a reduced arithmetic form of the Pythagorean quadrilateral counting problem. For an integer $$n=\prod_{i=1}^{r} p_i^{e_i},$$ the number of admissible quadrilaterals associated with the exact parameter \(n\) depends only on the exponent vector \((e_1,\dots,e_r)\), not on the actual prime values. The quantity required by the problem is then the divisor-summatory version of that exact count, so the solution never enumerates quadrilaterals geometrically. Once the reduction has been made, everything becomes multiplicative. The whole task is therefore to evaluate a short list of prime-power factors and combine them with fixed coefficients. Mathematical Approach Let \(F(n)\) denote the exact-count function for one fixed arithmetic parameter \(n\), and let $$S(n)=\sum_{d\mid n} F(d)$$ be the divisor-summatory quantity that the implementations actually evaluate. Because \(F(n)\) depends only on the exponents in the prime factorization of \(n\), it is convenient to work directly with those exponents. Step 1: Keep Only the Prime-Exponent Profile If \(n=\prod p_i^{e_i}\), then every divisor \(d\mid n\) has the form $$d=\prod_{i=1}^{r} p_i^{a_i},\qquad 0\le a_i\le e_i.$$ The reduced counting formula ignores the actual primes \(p_i\) and uses only the exponents....

Detailed mathematical approach

Problem Summary

The implementations start from a reduced arithmetic form of the Pythagorean quadrilateral counting problem. For an integer

$$n=\prod_{i=1}^{r} p_i^{e_i},$$

the number of admissible quadrilaterals associated with the exact parameter \(n\) depends only on the exponent vector \((e_1,\dots,e_r)\), not on the actual prime values. The quantity required by the problem is then the divisor-summatory version of that exact count, so the solution never enumerates quadrilaterals geometrically.

Once the reduction has been made, everything becomes multiplicative. The whole task is therefore to evaluate a short list of prime-power factors and combine them with fixed coefficients.

Mathematical Approach

Let \(F(n)\) denote the exact-count function for one fixed arithmetic parameter \(n\), and let

$$S(n)=\sum_{d\mid n} F(d)$$

be the divisor-summatory quantity that the implementations actually evaluate. Because \(F(n)\) depends only on the exponents in the prime factorization of \(n\), it is convenient to work directly with those exponents.

Step 1: Keep Only the Prime-Exponent Profile

If \(n=\prod p_i^{e_i}\), then every divisor \(d\mid n\) has the form

$$d=\prod_{i=1}^{r} p_i^{a_i},\qquad 0\le a_i\le e_i.$$

The reduced counting formula ignores the actual primes \(p_i\) and uses only the exponents. That means two integers with the same exponent profile give the same value of \(F\) and the same value of \(S\).

For one prime-power exponent \(e\), define the five local factors

$$B_1(e)=e+1,$$

$$B_2(e)=(e+1)^2,$$

$$B_3(e)=\left\lfloor\frac{(e+1)^2+1}{2}\right\rfloor,$$

$$B_4(e)=(e+1)^3,$$

$$B_5(e)=\frac{(e+1)(2e^2+4e+3)}{3}.$$

These are exactly the prime-power building blocks used by the reduced formula.

Step 2: Build the Exact-Count Function \(F(n)\)

Since the local contributions factor over primes, the exact count is multiplicative and can be written as

$$F(n)=7\prod_{i=1}^{r} B_1(e_i)-14\prod_{i=1}^{r} B_2(e_i)-4\prod_{i=1}^{r} B_3(e_i)+8\prod_{i=1}^{r} B_4(e_i)+4\prod_{i=1}^{r} B_5(e_i).$$

The first factor \(B_1(e)=e+1\) is the familiar local contribution from the divisor-count function. The higher-degree factors encode the remaining symmetry classes that survive the reduction from geometry to arithmetic.

Nothing in this formula depends on the magnitudes of the primes themselves. Only the exponent vector matters.

Step 3: Lift from an Exact Parameter to All Divisors

The problem asks for the divisor-sum lift of \(F\), namely

$$S(n)=\sum_{d\mid n} F(d)=\sum_{0\le a_i\le e_i} F(a_1,\dots,a_r).$$

Because \(F\) is a fixed linear combination of products, the divisor sum separates prime by prime:

$$S(n)=7\prod_{i=1}^{r}\left(\sum_{a=0}^{e_i} B_1(a)\right)-14\prod_{i=1}^{r}\left(\sum_{a=0}^{e_i} B_2(a)\right)-4\prod_{i=1}^{r}\left(\sum_{a=0}^{e_i} B_3(a)\right)+8\prod_{i=1}^{r}\left(\sum_{a=0}^{e_i} B_4(a)\right)+4\prod_{i=1}^{r}\left(\sum_{a=0}^{e_i} B_5(a)\right).$$

So the entire divisor summation is absorbed into five closed-form cumulative factors.

Step 4: Evaluate the Cumulative Local Sums

Define

$$A_j(e)=\sum_{a=0}^{e} B_j(a).$$

Then the implementations use the closed forms

$$A_1(e)=\frac{(e+1)(e+2)}{2},$$

$$A_2(e)=\frac{(e+1)(e+2)(2e+3)}{6},$$

$$A_3(e)=\sum_{a=0}^{e}\left\lfloor\frac{(a+1)^2+1}{2}\right\rfloor=\frac{1}{2}\left(A_2(e)+\left\lfloor\frac{e}{2}\right\rfloor+1\right),$$

$$A_4(e)=\sum_{a=0}^{e}(a+1)^3=A_1(e)^2,$$

$$A_5(e)=\sum_{a=0}^{e}\frac{(a+1)(2a^2+4a+3)}{3}=\frac{(e+1)(e+2)(e^2+3e+3)}{6}.$$

The parity correction inside \(A_3(e)\) comes from the fact that odd squares contribute the extra \(+1\) when the integer division in \(B_3\) is written as a floor.

Step 5: Final Closed Form for the Required Quantity

Substituting the cumulative factors into the separated divisor sum gives

$$\boxed{S(n)=7\prod_{i=1}^{r} A_1(e_i)-14\prod_{i=1}^{r} A_2(e_i)-4\prod_{i=1}^{r} A_3(e_i)+8\prod_{i=1}^{r} A_4(e_i)+4\prod_{i=1}^{r} A_5(e_i).}$$

This is the whole reason the computation is so small: all geometric counting has already been compressed into five multiplicative prime-power terms and their divisor-sum lifts.

Worked Example: Exponent Vector \((2,1)\)

For a quick exact-count checkpoint, a single prime square has exponent vector \((2)\). Then

$$B_1(2)=3,\quad B_2(2)=9,\quad B_3(2)=5,\quad B_4(2)=27,\quad B_5(2)=19,$$

so

$$F=7\cdot 3-14\cdot 9-4\cdot 5+8\cdot 27+4\cdot 19=167.$$

Now lift to the divisor sum for exponent vector \((2,1)\). The cumulative factors are

$$A_1(2)=6,\quad A_2(2)=14,\quad A_3(2)=8,\quad A_4(2)=36,\quad A_5(2)=26,$$

$$A_1(1)=3,\quad A_2(1)=5,\quad A_3(1)=3,\quad A_4(1)=9,\quad A_5(1)=7.$$

Therefore

$$\begin{aligned} S&=7(6\cdot 3)-14(14\cdot 5)-4(8\cdot 3)+8(36\cdot 9)+4(26\cdot 7)\\ &=126-980-96+2592+728\\ &=2370. \end{aligned}$$

This matches the checkpoint embedded in the implementations. Another checkpoint is \(S(1,1,1)=5535\).

How the Code Works

The C++, Python, and Java implementations store the target exponent vector \((6,3,2,1,1,1,1,1)\) directly. No runtime factorization is needed, because the reduced arithmetic formula already depends only on these exponents.

Each implementation evaluates the five cumulative product terms \(\prod A_j(e_i)\), combines them with the coefficients \(7,-14,-4,8,4\), and prints the resulting exact integer. The structure is intentionally minimal: one generic product loop, five local formulas, and one final linear combination.

All arithmetic is exact. The C++ implementation uses unsigned 128-bit arithmetic, Python uses built-in arbitrary-precision integers, and Java uses arbitrary-precision integers as well. No floating-point calculations or geometric searches are involved.

Complexity Analysis

If the exponent vector has length \(r\), each of the five product terms requires one pass over that vector. The total running time is therefore \(O(r)\), and the extra memory usage is \(O(1)\) beyond the fixed storage for the exponents and a few accumulator variables. In the target instance, \(r=8\), so the computation is effectively constant time.

Footnotes and References

  1. Problem page: https://projecteuler.net/problem=723
  2. Divisor function: Wikipedia — Divisor function
  3. Multiplicative function: Wikipedia — Multiplicative function
  4. Dirichlet convolution: Wikipedia — Dirichlet convolution
  5. Faulhaber's formula: Wikipedia — Faulhaber's formula

Problem 723 source code

C++

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

using u64 = std::uint64_t;
using u128 = unsigned __int128;

static std::string to_string_u128(u128 v) {
    if (v == 0) return "0";
    std::string s;
    while (v > 0) {
        const int d = static_cast<int>(v % 10);
        s.push_back(static_cast<char>('0' + d));
        v /= 10;
    }
    std::reverse(s.begin(), s.end());
    return s;
}

template <typename F>
static u128 product_over_exps(const std::array<int, 8>& exps, F term) {
    u128 out = 1;
    for (const int e : exps) out *= static_cast<u128>(term(e));
    return out;
}

static u128 F_value(const std::array<int, 8>& exps) {
    const u128 t1 = 7 * product_over_exps(exps, [](int e) -> u64 {
        return static_cast<u64>(e + 1);
    });
    const u128 t2 = 14 * product_over_exps(exps, [](int e) -> u64 {
        const u64 x = static_cast<u64>(e + 1);
        return x * x;
    });
    const u128 t3 = 4 * product_over_exps(exps, [](int e) -> u64 {
        const u64 x = static_cast<u64>(e + 1);
        return (x * x + 1) / 2;
    });
    const u128 t4 = 8 * product_over_exps(exps, [](int e) -> u64 {
        const u64 x = static_cast<u64>(e + 1);
        return x * x * x;
    });
    const u128 t5 = 4 * product_over_exps(exps, [](int e) -> u64 {
        const u64 x = static_cast<u64>(e + 1);
        const u64 y = static_cast<u64>(2 * e * e + 4 * e + 3);
        return (x * y) / 3;
    });
    return t1 - t2 - t3 + t4 + t5;
}

static u128 S_value(const std::array<int, 8>& exps) {
    const u128 t1 = 7 * product_over_exps(exps, [](int e) -> u64 {
        return static_cast<u64>((e + 1) * (e + 2) / 2);
    });
    const u128 t2 = 14 * product_over_exps(exps, [](int e) -> u64 {
        return static_cast<u64>((e + 1) * (e + 2) * (2 * e + 3) / 6);
    });
    const u128 t3 = 4 * product_over_exps(exps, [](int e) -> u64 {
        const u64 a = static_cast<u64>((e + 1) * (e + 2) * (2 * e + 3) / 6);
        const u64 b = static_cast<u64>(e / 2 + 1);
        return (a + b) / 2;
    });
    const u128 t4 = 8 * product_over_exps(exps, [](int e) -> u64 {
        const u64 x = static_cast<u64>((e + 1) * (e + 2) / 2);
        return x * x;
    });
    const u128 t5 = 4 * product_over_exps(exps, [](int e) -> u64 {
        return static_cast<u64>((e + 1) * (e + 2) * (e * e + 3 * e + 3) / 6);
    });
    return t1 - t2 - t3 + t4 + t5;
}

int main() {
    {
        std::array<int, 8> exps{};
        assert(F_value(exps) == 1);
    }
    {
        std::array<int, 8> exps{};
        exps[0] = 1;
        assert(F_value(exps) == 38);
    }
    {
        std::array<int, 8> exps{};
        exps[0] = 2;
        assert(F_value(exps) == 167);
    }
    {
        std::array<int, 8> exps{};
        exps[0] = 2;
        exps[1] = 1;
        assert(S_value(exps) == 2370);
    }
    {
        std::array<int, 8> exps{};
        exps[0] = 1;
        exps[1] = 1;
        exps[2] = 1;
        assert(S_value(exps) == 5535);
    }

    const std::array<int, 8> target_exps = {6, 3, 2, 1, 1, 1, 1, 1};
    std::cout << to_string_u128(S_value(target_exps)) << '\n';
    return 0;
}

Python

def solve():
    exps = [6, 3, 2, 1, 1, 1, 1, 1]
    
    def prod(term_func):
        out = 1
        for e in exps:
            out *= term_func(e)
        return out
        
    t1 = 7 * prod(lambda e: (e + 1) * (e + 2) // 2)
    t2 = 14 * prod(lambda e: (e + 1) * (e + 2) * (2 * e + 3) // 6)
    
    def term3(e):
        a = (e + 1) * (e + 2) * (2 * e + 3) // 6
        b = e // 2 + 1
        return (a + b) // 2
    t3 = 4 * prod(term3)
    
    def term4(e):
        x = (e + 1) * (e + 2) // 2
        return x * x
    t4 = 8 * prod(term4)
    
    t5 = 4 * prod(lambda e: (e + 1) * (e + 2) * (e * e + 3 * e + 3) // 6)
    
    ans = t1 - t2 - t3 + t4 + t5
    return str(ans)

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

Java

import java.math.BigInteger;

public class Euler723 {
    public static String solve() {
        int[] exps = { 6, 3, 2, 1, 1, 1, 1, 1 };

        BigInteger t1Term = BigInteger.ONE;
        for (int e : exps) {
            long val = (long) (e + 1) * (e + 2) / 2;
            t1Term = t1Term.multiply(BigInteger.valueOf(val));
        }
        BigInteger t1 = BigInteger.valueOf(7).multiply(t1Term);

        BigInteger t2Term = BigInteger.ONE;
        for (int e : exps) {
            long val = (long) (e + 1) * (e + 2) * (2 * e + 3) / 6;
            t2Term = t2Term.multiply(BigInteger.valueOf(val));
        }
        BigInteger t2 = BigInteger.valueOf(14).multiply(t2Term);

        BigInteger t3Term = BigInteger.ONE;
        for (int e : exps) {
            long a = (long) (e + 1) * (e + 2) * (2 * e + 3) / 6;
            long b = e / 2 + 1;
            long val = (a + b) / 2;
            t3Term = t3Term.multiply(BigInteger.valueOf(val));
        }
        BigInteger t3 = BigInteger.valueOf(4).multiply(t3Term);

        BigInteger t4Term = BigInteger.ONE;
        for (int e : exps) {
            long x = (long) (e + 1) * (e + 2) / 2;
            long val = x * x;
            t4Term = t4Term.multiply(BigInteger.valueOf(val));
        }
        BigInteger t4 = BigInteger.valueOf(8).multiply(t4Term);

        BigInteger t5Term = BigInteger.ONE;
        for (int e : exps) {
            long val = (long) (e + 1) * (e + 2) * (e * e + 3 * e + 3) / 6;
            t5Term = t5Term.multiply(BigInteger.valueOf(val));
        }
        BigInteger t5 = BigInteger.valueOf(4).multiply(t5Term);

        BigInteger ans = t1.subtract(t2).subtract(t3).add(t4).add(t5);
        return ans.toString();
    }

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