Problem 350: Constraining the Least Greatest and the Greatest Least

View on Project Euler

Project Euler Problem 350 Solution

EulerSolve provides an optimized solution for Project Euler Problem 350, Constraining the Least Greatest and the Greatest Least, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary We must count positive integer \(N\)-tuples \((x_1,\dots,x_N)\) such that $$\gcd(x_1,\dots,x_N)\ge G,\qquad \operatorname{lcm}(x_1,\dots,x_N)\le L,$$ and return the result modulo \(101^4\). A direct search is impossible for the actual parameters, so the solution reorganizes each tuple by separating its common gcd from its normalized lcm structure. Mathematical Approach Let $$F(G,L,N)=\#\left\{(x_1,\dots,x_N)\in \mathbb{Z}_{>0}^N:\gcd(x_1,\dots,x_N)\ge G,\ \operatorname{lcm}(x_1,\dots,x_N)\le L\right\}.$$ The key idea is to classify tuples by their common gcd and by the lcm that remains after dividing that gcd out. Step 1: Normalize by the Common GCD For any admissible tuple define $$g=\gcd(x_1,\dots,x_N),\qquad x_i=g\,y_i.$$ Then the normalized tuple satisfies $$\gcd(y_1,\dots,y_N)=1.$$ Now define $$m=\operatorname{lcm}(y_1,\dots,y_N).$$ Because lcm scales linearly when every coordinate is multiplied by the same factor, we obtain $$\operatorname{lcm}(x_1,\dots,x_N)=g\,m.$$ So once the normalized tuple is fixed, the only remaining freedom is the choice of \(g\)....

Detailed mathematical approach

Problem Summary

We must count positive integer \(N\)-tuples \((x_1,\dots,x_N)\) such that

$$\gcd(x_1,\dots,x_N)\ge G,\qquad \operatorname{lcm}(x_1,\dots,x_N)\le L,$$

and return the result modulo \(101^4\). A direct search is impossible for the actual parameters, so the solution reorganizes each tuple by separating its common gcd from its normalized lcm structure.

Mathematical Approach

Let

$$F(G,L,N)=\#\left\{(x_1,\dots,x_N)\in \mathbb{Z}_{>0}^N:\gcd(x_1,\dots,x_N)\ge G,\ \operatorname{lcm}(x_1,\dots,x_N)\le L\right\}.$$

The key idea is to classify tuples by their common gcd and by the lcm that remains after dividing that gcd out.

Step 1: Normalize by the Common GCD

For any admissible tuple define

$$g=\gcd(x_1,\dots,x_N),\qquad x_i=g\,y_i.$$

Then the normalized tuple satisfies

$$\gcd(y_1,\dots,y_N)=1.$$

Now define

$$m=\operatorname{lcm}(y_1,\dots,y_N).$$

Because lcm scales linearly when every coordinate is multiplied by the same factor, we obtain

$$\operatorname{lcm}(x_1,\dots,x_N)=g\,m.$$

So once the normalized tuple is fixed, the only remaining freedom is the choice of \(g\). The constraints become

$$G\le g\le \left\lfloor\frac{L}{m}\right\rfloor.$$

Hence the number of admissible gcd values for this normalized tuple is

$$w(m)=\left\lfloor\frac{L}{m}\right\rfloor-G+1,$$

and this is positive only when

$$m\le K=\left\lfloor\frac{L}{G}\right\rfloor.$$

Step 2: Define the Normalized Counting Function

Let

$$H_N(m)=\#\left\{(y_1,\dots,y_N)\in \mathbb{Z}_{>0}^N:\gcd(y_1,\dots,y_N)=1,\ \operatorname{lcm}(y_1,\dots,y_N)=m\right\}.$$

This counts exactly the normalized tuples whose lcm is \(m\). Summing over all possible normalized lcm values gives

$$F(G,L,N)=\sum_{m=1}^{\lfloor L/G\rfloor} H_N(m)\left(\left\lfloor\frac{L}{m}\right\rfloor-G+1\right).$$

Therefore the hard part of the problem is reduced to computing \(H_N(m)\) efficiently for every \(m\le K\).

Step 3: Count One Prime Power

Write the prime factorization of \(m\) as

$$m=\prod_{p} p^{e_p}.$$

Fix one prime power \(p^e\parallel m\). For each coordinate \(y_i\), let \(\alpha_i=v_p(y_i)\). Since \(y_i\mid m\), each exponent lies in

$$\alpha_i\in\{0,1,\dots,e\}.$$

The normalized conditions translate into simple conditions on these exponents:

$$\gcd(y_1,\dots,y_N)=1 \iff \min(\alpha_1,\dots,\alpha_N)=0,$$

$$\operatorname{lcm}(y_1,\dots,y_N)=m \iff \max(\alpha_1,\dots,\alpha_N)=e.$$

So for each prime \(p^e\parallel m\), we only need to count exponent vectors in \(\{0,\dots,e\}^N\) whose minimum is \(0\) and whose maximum is \(e\).

Step 4: Inclusion-Exclusion for the Local Factor

The total number of exponent vectors is \((e+1)^N\).

Vectors with no zero exponent use only \(\{1,\dots,e\}\), so there are \(e^N\) of them.

Vectors with no exponent equal to \(e\) also number \(e^N\).

Vectors having neither \(0\) nor \(e\) use only \(\{1,\dots,e-1\}\), so there are \((e-1)^N\).

Therefore inclusion-exclusion gives the local count

$$C_e=(e+1)^N-2e^N+(e-1)^N.$$

This is the number of admissible \(p\)-adic exponent patterns for one prime power. For example, when \(e=1\), we get \(C_1=2^N-2\): all binary exponent vectors except the all-zero and all-one vectors.

Step 5: Reconstruct \(H_N(m)\)

Exponent choices for distinct primes are independent. Once an admissible exponent vector is chosen for each prime \(p^{e_p}\parallel m\), multiplying the prime-power contributions reconstructs one unique normalized tuple \((y_1,\dots,y_N)\).

Hence \(H_N(m)\) is multiplicative and factors as

$$H_N(m)=\prod_{p^{e_p}\parallel m} C_{e_p}=\prod_{p^{e_p}\parallel m}\left((e_p+1)^N-2e_p^N+(e_p-1)^N\right).$$

Step 6: Final Closed Form

Substituting the normalized count into the outer sum yields

$$\boxed{F(G,L,N)=\sum_{m=1}^{\lfloor L/G\rfloor}\left(\prod_{p^{e_p}\parallel m}\left((e_p+1)^N-2e_p^N+(e_p-1)^N\right)\right)\left(\left\lfloor\frac{L}{m}\right\rfloor-G+1\right)\pmod{101^4}.}$$

This is exactly the formula implemented by the program.

Worked Example: \((G,L,N)=(10,100,2)\)

Here

$$K=\left\lfloor\frac{100}{10}\right\rfloor=10.$$

For \(N=2\), the local factor simplifies to

$$C_e=(e+1)^2-2e^2+(e-1)^2=2.$$

So every prime dividing \(m\) contributes a factor \(2\), and therefore

$$H_2(m)=2^{\omega(m)},$$

where \(\omega(m)\) is the number of distinct prime divisors of \(m\).

Now evaluate the first few terms:

$$\begin{aligned} m=1&:&&H_2(1)=1,\quad w(1)=100-10+1=91,\quad c=91,\\ m=2&:&&H_2(2)=2,\quad w(2)=50-10+1=41,\quad c=82,\\ m=3&:&&H_2(3)=2,\quad w(3)=33-10+1=24,\quad c=48,\\ m=4&:&&H_2(4)=2,\quad w(4)=25-10+1=16,\quad c=32,\\ m=5&:&&H_2(5)=2,\quad w(5)=20-10+1=11,\quad c=22,\\ m=6&:&&H_2(6)=4,\quad w(6)=16-10+1=7,\quad c=28. \end{aligned}$$

The remaining terms are \(10,6,4,4\) for \(m=7,8,9,10\). Therefore

$$91+82+48+32+22+28+10+6+4+4=327,$$

which matches the checkpoint used by the implementation.

How the Code Works

The program first sets \(K=\lfloor L/G\rfloor\) and builds an SPF table (smallest prime factor) up to \(K\). It then precomputes

$$\texttt{local}[e]=C_e=(e+1)^N-2e^N+(e-1)^N \pmod{101^4}$$

using fast modular exponentiation. For each \(m\le K\), the SPF table gives the factorization of \(m\), so the code multiplies the relevant \(\texttt{local}[e]\) values to obtain

$$\texttt{h}[m]=H_N(m)\pmod{101^4}.$$

Finally it accumulates

$$\texttt{h}[m]\left(\left\lfloor\frac{L}{m}\right\rfloor-G+1\right)$$

for all \(1\le m\le K\), always reducing modulo \(101^4\).

Complexity Analysis

Let \(K=\lfloor L/G\rfloor\). Building the SPF sieve costs \(O(K\log\log K)\) time and \(O(K)\) memory. Factoring every \(m\le K\) via the SPF table touches each prime-power exponent once; the total divisor-stripping work has average order \(O(K\log\log K)\). Thus the overall method is near-linear in \(K\) and uses \(O(K)\) memory.

Footnotes and References

  1. Problem page: https://projecteuler.net/problem=350
  2. Least common multiple: Wikipedia — Least common multiple
  3. Greatest common divisor: Wikipedia — Greatest common divisor
  4. Inclusion-exclusion principle: Wikipedia — Inclusion-exclusion principle
  5. Multiplicative function: Wikipedia — Multiplicative function

Problem 350 source code

C++

#include <algorithm>
#include <cstdint>
#include <iostream>
#include <string>
#include <vector>

namespace {

using i64 = std::int64_t;
using u64 = std::uint64_t;

constexpr i64 MOD = 104060401;  // 101^4

i64 mod_add(const i64 a, const i64 b) {
    i64 x = a + b;
    if (x >= MOD) {
        x -= MOD;
    }
    return x;
}

i64 mod_sub(const i64 a, const i64 b) {
    i64 x = a - b;
    if (x < 0) {
        x += MOD;
    }
    return x;
}

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

i64 mod_pow(i64 base, u64 exp) {
    base %= MOD;
    i64 result = 1;
    while (exp > 0) {
        if ((exp & 1ULL) != 0ULL) {
            result = mod_mul(result, base);
        }
        base = mod_mul(base, base);
        exp >>= 1ULL;
    }
    return result;
}

std::vector<int> build_spf(const int n) {
    std::vector<int> spf(static_cast<std::size_t>(n + 1));
    for (int i = 0; i <= n; ++i) {
        spf[static_cast<std::size_t>(i)] = i;
    }

    for (int i = 2; static_cast<i64>(i) * i <= n; ++i) {
        if (spf[static_cast<std::size_t>(i)] != i) {
            continue;
        }
        for (int j = i * i; j <= n; j += i) {
            if (spf[static_cast<std::size_t>(j)] == j) {
                spf[static_cast<std::size_t>(j)] = i;
            }
        }
    }

    return spf;
}

i64 solve(const u64 G, const u64 L, const u64 N) {
    const int K = static_cast<int>(L / G);
    const std::vector<int> spf = build_spf(K);

    int max_exp = 0;
    {
        int x = K;
        while (x > 0) {
            x /= 2;
            ++max_exp;
        }
    }

    std::vector<i64> pows(static_cast<std::size_t>(max_exp + 2), 0);
    for (int b = 0; b <= max_exp + 1; ++b) {
        if (b == 0) {
            pows[static_cast<std::size_t>(b)] = 0;  // N > 0 here
        } else {
            pows[static_cast<std::size_t>(b)] = mod_pow(b, N);
        }
    }

    std::vector<i64> local(static_cast<std::size_t>(max_exp + 1), 0);
    for (int a = 1; a <= max_exp; ++a) {
        // (a+1)^N - 2*a^N + (a-1)^N
        i64 value = mod_sub(pows[static_cast<std::size_t>(a + 1)], mod_mul(2, pows[static_cast<std::size_t>(a)]));
        value = mod_add(value, pows[static_cast<std::size_t>(a - 1)]);
        local[static_cast<std::size_t>(a)] = value;
    }

    std::vector<i64> h(static_cast<std::size_t>(K + 1), 0);
    h[1] = 1;

    for (int n = 2; n <= K; ++n) {
        int x = n;
        i64 value = 1;

        while (x > 1) {
            const int p = spf[static_cast<std::size_t>(x)];
            int e = 0;
            while (x % p == 0) {
                x /= p;
                ++e;
            }
            value = mod_mul(value, local[static_cast<std::size_t>(e)]);
        }

        h[static_cast<std::size_t>(n)] = value;
    }

    i64 answer = 0;
    for (int m = 1; m <= K; ++m) {
        const u64 weight = L / static_cast<u64>(m) - G + 1ULL;
        const i64 w = static_cast<i64>(weight % static_cast<u64>(MOD));
        answer = mod_add(answer, mod_mul(h[static_cast<std::size_t>(m)], w));
    }

    return answer;
}

bool run_checkpoints() {
    if (solve(10ULL, 100ULL, 1ULL) != 91) {
        std::cerr << "Checkpoint failed: f(10,100,1)\n";
        return false;
    }
    if (solve(10ULL, 100ULL, 2ULL) != 327) {
        std::cerr << "Checkpoint failed: f(10,100,2)\n";
        return false;
    }
    if (solve(10ULL, 100ULL, 3ULL) != 1135) {
        std::cerr << "Checkpoint failed: f(10,100,3)\n";
        return false;
    }
    if (solve(10ULL, 100ULL, 1000ULL) != 3286053) {
        std::cerr << "Checkpoint failed: f(10,100,1000) mod 101^4\n";
        return false;
    }
    return true;
}

}  // namespace

int main(int argc, char** argv) {
    bool skip_checkpoints = false;
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            skip_checkpoints = true;
        } else {
            std::cerr << "Unknown argument: " << arg << '\n';
            return 1;
        }
    }

    if (!skip_checkpoints && !run_checkpoints()) {
        return 2;
    }

    const i64 answer = solve(1000000ULL, 1000000000000ULL, 1000000000000000000ULL);
    std::cout << answer << '\n';
    return 0;
}

Python

MOD = 104060401

def mod_pow(base, exp):
    return pow(base, exp, MOD)

def build_spf(n):
    spf = list(range(n + 1))
    for i in range(2, int(n ** 0.5) + 1):
        if spf[i] == i:
            for j in range(i * i, n + 1, i):
                if spf[j] == j:
                    spf[j] = i
    return spf

def solve_350(G, L, N):
    K = L // G
    spf = build_spf(K)
    
    max_exp = 0
    x = K
    while x > 0:
        x //= 2
        max_exp += 1
        
    pows = [0] * (max_exp + 2)
    for b in range(max_exp + 2):
        if b == 0:
            pows[b] = 0
        else:
            pows[b] = mod_pow(b, N)
            
    local = [0] * (max_exp + 1)
    for a in range(1, max_exp + 1):
        val = (pows[a + 1] - 2 * pows[a] + pows[a - 1]) % MOD
        local[a] = (val + MOD) % MOD
        
    h = [0] * (K + 1)
    h[1] = 1
    
    for n in range(2, K + 1):
        x = n
        value = 1
        while x > 1:
            p = spf[x]
            e = 0
            while x % p == 0:
                x //= p
                e += 1
            value = (value * local[e]) % MOD
        h[n] = value
        
    answer = 0
    for m in range(1, K + 1):
        weight = L // m - G + 1
        w = weight % MOD
        answer = (answer + h[m] * w) % MOD
        
    return answer

def solve():
    G = 1000000
    L = 1000000000000
    N = 1000000000000000000
    ans = solve_350(G, L, N)
    return str(ans)

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

Java

import java.util.*;

public class Euler350 {

    static final long MOD = 104060401L;

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

    static int[] buildSpf(int n) {
        int[] spf = new int[n + 1];
        for (int i = 0; i <= n; i++)
            spf[i] = i;

        for (int i = 2; i * i <= n; i++) {
            if (spf[i] == i) {
                for (int j = i * i; j <= n; j += i) {
                    if (spf[j] == j) {
                        spf[j] = i;
                    }
                }
            }
        }
        return spf;
    }

    static long solve(long G, long L, long N) {
        int K = (int) (L / G);
        int[] spf = buildSpf(K);

        int maxExp = 0;
        int x = K;
        while (x > 0) {
            x /= 2;
            maxExp++;
        }

        long[] pows = new long[maxExp + 2];
        for (int b = 0; b <= maxExp + 1; b++) {
            if (b == 0) {
                pows[b] = 0;
            } else {
                pows[b] = modPow(b, N);
            }
        }

        long[] local = new long[maxExp + 1];
        for (int a = 1; a <= maxExp; a++) {
            long value = (pows[a + 1] - 2 * pows[a] % MOD + MOD) % MOD;
            value = (value + pows[a - 1]) % MOD;
            local[a] = value;
        }

        long[] h = new long[K + 1];
        h[1] = 1;

        for (int n = 2; n <= K; n++) {
            int cx = n;
            long value = 1;

            while (cx > 1) {
                int p = spf[cx];
                int e = 0;
                while (cx % p == 0) {
                    cx /= p;
                    e++;
                }
                value = (value * local[e]) % MOD;
            }
            h[n] = value;
        }

        long answer = 0;
        for (int m = 1; m <= K; m++) {
            long weight = L / m - G + 1;
            long w = weight % MOD;
            answer = (answer + h[m] * w) % MOD;
        }

        return answer;
    }

    public static String solve() {
        long ans = solve(1000000L, 1000000000000L, 1000000000000000000L);
        return String.valueOf(ans);
    }

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