Problem 821: 123-Separable

Project Euler Problem 821 Solution

EulerSolve provides an optimized solution for Project Euler Problem 821, 123-Separable, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The problem asks for the largest total number of integers in $1,\dots,N$ that can be arranged into pairwise non-overlapping local $123$-blocks. Every positive integer has a unique representation $$n=m\,2^a3^b,\qquad \gcd(m,6)=1,$$ so the search splits into independent components indexed by the part $m$ that is coprime to $6$. Inside one component, only $2^a3^b$-smooth multipliers matter, and the final answer is the sum of the optimal local contributions. Mathematical Approach Write $A(N)$ for the required optimum. The fast solution has two layers: first solve one component exactly, then sum those component values over all possible $m$. Step 1: Split by the part coprime to 6 Fix $m$ with $\gcd(m,6)=1$ and let $$L=\left\lfloor\frac{N}{m}\right\rfloor.$$ Every integer in $1,\dots,N$ whose $6$-free part is $m$ can be written uniquely as $m s$ with $$\mathcal{S}(L)=\{2^a3^b\le L: a,b\ge 0\}.$$ Different values of $m$ never mix, so the global optimum is a sum of independent component optima: $$A(N)=\sum_{\substack{1\le m\le N\\ \gcd(m,6)=1}} G\!\left(\left\lfloor\frac{N}{m}\right\rfloor\right),$$ where $G(L)$ denotes the best contribution of one component with smooth bound $L$....

Detailed mathematical approach

Problem Summary

The problem asks for the largest total number of integers in $1,\dots,N$ that can be arranged into pairwise non-overlapping local $123$-blocks. Every positive integer has a unique representation

$$n=m\,2^a3^b,\qquad \gcd(m,6)=1,$$

so the search splits into independent components indexed by the part $m$ that is coprime to $6$. Inside one component, only $2^a3^b$-smooth multipliers matter, and the final answer is the sum of the optimal local contributions.

Mathematical Approach

Write $A(N)$ for the required optimum. The fast solution has two layers: first solve one component exactly, then sum those component values over all possible $m$.

Step 1: Split by the part coprime to 6

Fix $m$ with $\gcd(m,6)=1$ and let

$$L=\left\lfloor\frac{N}{m}\right\rfloor.$$

Every integer in $1,\dots,N$ whose $6$-free part is $m$ can be written uniquely as $m s$ with

$$\mathcal{S}(L)=\{2^a3^b\le L: a,b\ge 0\}.$$

Different values of $m$ never mix, so the global optimum is a sum of independent component optima:

$$A(N)=\sum_{\substack{1\le m\le N\\ \gcd(m,6)=1}} G\!\left(\left\lfloor\frac{N}{m}\right\rfloor\right),$$

where $G(L)$ denotes the best contribution of one component with smooth bound $L$.

Step 2: Turn one component into a weighted graph

For a smooth base $s\in \mathcal{S}(L)$, the local block is

$$B_L(s)=\{s,2s,3s\}\cap [1,L].$$

After multiplying by $m$, this becomes $\{ms,2ms,3ms\}\cap [1,N]$. Its size is

$$|B_L(s)|=1+\mathbf{1}_{2s\le L}+\mathbf{1}_{3s\le L}.$$

So inside one component we want a collection of smooth bases whose blocks do not overlap, maximizing the total weight contributed by those block sizes.

Step 3: Why only six local conflicts matter

Two local blocks $B_L(s)$ and $B_L(t)$ overlap exactly when

$$i s=j t\quad \text{for some } i,j\in\{1,2,3\}.$$

For distinct $s$ and $t$, this is equivalent to

$$\frac{t}{s}\in \left\{2,3,\frac{3}{2},\frac{1}{2},\frac{1}{3},\frac{2}{3}\right\}.$$

If we write $s=2^a3^b$, these six ratios become the six neighbor moves

$$ (a,b)\leftrightarrow(a\pm 1,b),\qquad (a,b)\leftrightarrow(a,b\pm 1),\qquad (a,b)\leftrightarrow(a+1,b-1),\qquad (a,b)\leftrightarrow(a-1,b+1). $$

Thus $G(L)$ is a maximum-weight independent-set value on the triangular lattice formed by the exponent pairs $(a,b)$ with $2^a3^b\le L$.

Step 4: Closed form for the component value

The key structural fact used by the implementations is that this graph optimization has a simple exact value:

$$G(L)=|\mathcal{S}(L)|-|\mathcal{E}\cap [1,L]|,$$

where the exceptional smooth numbers are

$$\mathcal{E}=\{6,24,54\}\cup \{384\cdot 8^t:t\ge 0\}\cup \{243\cdot 27^t:t\ge 0\}.$$

So every $2^a3^b$-smooth number normally contributes one unit to the optimum, and only the explicit exceptional values subtract one unit. All three implementations are based on this identity.

Step 5: Reassemble the global answer

Substituting the component formula into the outer sum gives

$$A(N)=\sum_{\substack{1\le m\le N\\ \gcd(m,6)=1}}\left(|\mathcal{S}(\lfloor N/m\rfloor)|-|\mathcal{E}\cap [1,\lfloor N/m\rfloor]|\right).$$

This already removes the hard combinatorial search. The remaining work is purely arithmetic: count smooth numbers, count exceptional smooth numbers, and aggregate them over the quotient values $\lfloor N/m\rfloor$.

Step 6: Group equal quotients into intervals

List the smooth numbers in increasing order. Between two consecutive smooth breakpoints, neither $|\mathcal{S}(L)|$ nor $|\mathcal{E}\cap [1,L]|$ changes, because every exceptional number is itself smooth. Hence $G(L)$ is constant on each interval $[L,R]$, where $L$ is a smooth number and $R$ is the next smooth number minus $1$, or $R=N$ for the final interval.

For such an interval, the condition

$$L\le \left\lfloor\frac{N}{m}\right\rfloor\le R$$

is equivalent to

$$m_{\mathrm{low}}=\left\lfloor\frac{N}{R+1}\right\rfloor+1,\qquad m_{\mathrm{high}}=\left\lfloor\frac{N}{L}\right\rfloor.$$

The count of integers up to $x$ that are coprime to $6$ is

$$C(x)=x-\left\lfloor\frac{x}{2}\right\rfloor-\left\lfloor\frac{x}{3}\right\rfloor+\left\lfloor\frac{x}{6}\right\rfloor.$$

Therefore one quotient interval contributes

$$G(L)\left(C(m_{\mathrm{high}})-C(m_{\mathrm{low}}-1)\right).$$

Worked Example: $N=20$

The integers up to $20$ that are $2^a3^b$-smooth are

$$1,2,3,4,6,8,9,12,16,18,$$

so $|\mathcal{S}(20)|=10$. At this scale the only exceptional value is $6$, so

$$G(20)=10-1=9.$$

The integers $m\le 20$ with $\gcd(m,6)=1$ are

$$1,5,7,11,13,17,19.$$

The corresponding quotient values are

$$\left\lfloor\frac{20}{1}\right\rfloor=20,\quad \left\lfloor\frac{20}{5}\right\rfloor=4,\quad \left\lfloor\frac{20}{7}\right\rfloor=2,\quad \left\lfloor\frac{20}{11}\right\rfloor=\left\lfloor\frac{20}{13}\right\rfloor=\left\lfloor\frac{20}{17}\right\rfloor=\left\lfloor\frac{20}{19}\right\rfloor=1.$$

Since $G(4)=4$, $G(2)=2$, and $G(1)=1$, we get

$$A(20)=9+4+2+1+1+1+1=19,$$

which matches the checkpoint used by the implementations.

How the Code Works

The C++, Python, and Java implementations never build the conflict graph explicitly. Instead they first generate all $2^a3^b$-smooth numbers up to $N$, sort them, and remove duplicates. They also generate the exceptional subset up to the same limit.

Next, they scan the smooth list once in increasing order. At each smooth breakpoint they update the current component value $G(L)$: ordinary smooth numbers increase it by one, while exceptional smooth numbers leave it unchanged because they represent the one-unit loss encoded by the closed form.

Finally, the implementation walks through consecutive breakpoint intervals $[L,R]$, converts each interval into the matching range of $m$-values, counts how many of those $m$ are coprime to $6$ with the inclusion-exclusion formula for $C(x)$, and adds the resulting block contribution. Python gets arbitrary-precision arithmetic automatically, while the compiled implementations use a wider accumulator for the final sum.

Complexity Analysis

Let $S(N)=|\mathcal{S}(N)|$. Since $2^a3^b\le N$ implies $0\le a\le \lfloor \log_2 N\rfloor$ and $0\le b\le \lfloor \log_3 N\rfloor$, we have $S(N)=O((\log N)^2)$. Generating and sorting the smooth list costs $O(S(N)\log S(N))$, the exceptional list has only $O(\log N)$ terms, and the final sweep over the breakpoint intervals is $O(S(N))$. Memory usage is $O(S(N))=O((\log N)^2)$.

Footnotes and References

Problem page: Project Euler 821
Smooth numbers: Wikipedia — Smooth number
Independent sets in graphs: Wikipedia — Independent set
Inclusion-exclusion principle: Wikipedia — Inclusion-exclusion principle
Floor function: Wikipedia — Floor and ceiling functions

Problem 821 source code

C++

#include <algorithm>
#include <array>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <cstdlib>
#include <future>
#include <iomanip>
#include <iostream>
#include <limits>
#include <map>
#include <numeric>
#include <queue>
#include <set>
#include <string>
#include <thread>
#include <tuple>
#include <unordered_map>
#include <unordered_set>
#include <utility>
#include <vector>
using namespace std;

using int64 = long long;
using i128 = __int128_t;

static inline int64 count_coprime6(int64 x) {
    if (x <= 0) return 0;
    return x - x / 2 - x / 3 + x / 6;
}

static vector<int64> generate_smooth(int64 N) {
    vector<int64> v;
    for (int64 p2 = 1; p2 <= N; ) {
        for (int64 cur = p2; cur <= N; ) {
            v.push_back(cur);
            if (cur > N / 3) break;
            cur *= 3;
        }
        if (p2 > N / 2) break;
        p2 *= 2;
    }
    sort(v.begin(), v.end());
    v.erase(unique(v.begin(), v.end()), v.end());
    return v;
}

static vector<int64> generate_bad(int64 N) {
    vector<int64> b;
    auto push_if = [&](int64 x) { if (x <= N) b.push_back(x); };

    push_if(6);
    push_if(24);
    push_if(54);

    {
        int64 x = 384;
        while (x <= N) {
            b.push_back(x);
            if (x > N / 8) break;
            x *= 8;
        }
    }
    {
        int64 x = 243;
        while (x <= N) {
            b.push_back(x);
            if (x > N / 27) break;
            x *= 27;
        }
    }

    sort(b.begin(), b.end());
    b.erase(unique(b.begin(), b.end()), b.end());
    return b;
}

static int64 solve_fast(int64 N) {
    vector<int64> smooth = generate_smooth(N);
    if (smooth.empty() || smooth.front() != 1) {
        throw runtime_error("smooth list must start with 1");
    }

    vector<int64> bad = generate_bad(N);

    for (int64 x : bad) {
        if (!binary_search(smooth.begin(), smooth.end(), x)) {
            cerr << "Bad number not found in smooth list: " << x << "\n";
            throw runtime_error("bad number missing from smooth list");
        }
    }

    vector<int64> fAt(smooth.size(), 0);
    size_t ib = 0;
    int64 f = 0;
    for (size_t i = 0; i < smooth.size(); i++) {
        if (ib < bad.size() && smooth[i] == bad[ib]) {
            ++ib;
        } else {
            ++f;
        }
        fAt[i] = f;
        if (fAt[i] > (int64)(i + 1)) {
            throw runtime_error("fAt exceeded i+1 (logic error)");
        }
    }
    if (ib != bad.size()) {
        throw runtime_error("Not all bad numbers were consumed in scan");
    }

    int numThreads = 1;
#ifdef USE_THREADS
    numThreads = max(1u, thread::hardware_concurrency());
#endif

    vector<i128> partial(numThreads, 0);

    auto worker = [&](int tid) {
        size_t n = smooth.size();
        size_t lo = (n * (size_t)tid) / (size_t)numThreads;
        size_t hi = (n * (size_t)(tid + 1)) / (size_t)numThreads;

        i128 sum = 0;
        for (size_t i = lo; i < hi; i++) {
            int64 L = smooth[i];
            int64 R = (i + 1 < n ? smooth[i + 1] - 1 : N);
            if (R > N) R = N;

            int64 mLow  = N / (R + 1) + 1;
            int64 mHigh = N / L;
            if (mLow > mHigh) continue;

            int64 cnt = count_coprime6(mHigh) - count_coprime6(mLow - 1);
            sum += (i128)fAt[i] * (i128)cnt;
        }
        partial[tid] = sum;
    };

    if (numThreads == 1) {
        worker(0);
    } else {
        vector<thread> th;
        th.reserve(numThreads);
        for (int t = 0; t < numThreads; t++) th.emplace_back(worker, t);
        for (auto &tt : th) tt.join();
    }

    i128 ans = 0;
    for (int t = 0; t < numThreads; t++) ans += partial[t];

    return (int64)ans;
}

#ifdef VALIDATE
static int64 brute_component_f_small(int64 L) {
    vector<pair<int,int>> pts;
    vector<int64> vals;
    for (int a = 0, p2 = 1; (int64)p2 <= L; a++) {
        for (int b = 0, cur = p2; (int64)cur <= L; b++) {
            pts.push_back({a,b});
            vals.push_back(cur);
            if (cur > L/3) break;
            cur *= 3;
        }
        if (p2 > L/2) break;
        p2 *= 2;
    }
    int n = (int)vals.size();
    if (n > 64) {
        throw runtime_error("Validation brute_component_f_small exceeded 64 nodes");
    }

    vector<int64> w(n);
    for (int i = 0; i < n; i++) {
        int64 v = vals[i];
        w[i] = 1 + (2*v <= L) + (3*v <= L);
    }

    unordered_map<long long, long long> idx;
    idx.reserve(n*2);
    for (int i = 0; i < n; i++) {
        long long key = ((long long)pts[i].first << 32) ^ (unsigned)pts[i].second;
        idx[key] = i;
    }

    vector<uint64_t> adj(n, 0);
    auto get = [&](int a,int b)->int {
        long long key = ((long long)a << 32) ^ (unsigned)b;
        auto it = idx.find(key);
        return (it == idx.end() ? -1 : (int)it->second);
    };

    for (int i = 0; i < n; i++) {
        int a = pts[i].first, b = pts[i].second;
        const int da[6] = {1,-1,0,0,1,-1};
        const int db[6] = {0,0,1,-1,-1,1};
        for (int k = 0; k < 6; k++) {
            int j = get(a + da[k], b + db[k]);
            if (j >= 0) adj[i] |= (1ULL << j);
        }
    }

    unordered_map<uint64_t, int64> memo;
    memo.reserve(1<<16);

    function<int64(uint64_t)> dp = [&](uint64_t mask)->int64 {
        if (!mask) return 0;
        auto it = memo.find(mask);
        if (it != memo.end()) return it->second;

        int i = __builtin_ctzll(mask);
        uint64_t m0 = mask & ~(1ULL << i);
        int64 best = dp(m0);
        int64 take = w[i] + dp(m0 & ~adj[i]);
        if (take > best) best = take;
        memo[mask] = best;
        return best;
    };

    uint64_t full = (n == 64 ? ~0ULL : ((1ULL << n) - 1ULL));
    return dp(full);
}

static int64 brute_F_small(int64 N) {
    unordered_map<int64, int64> memo_f;
    memo_f.reserve(N);
    auto f = [&](int64 L)->int64 {
        auto it = memo_f.find(L);
        if (it != memo_f.end()) return it->second;
        int64 val = brute_component_f_small(L);
        memo_f[L] = val;
        return val;
    };

    int64 ans = 0;
    for (int64 m = 1; m <= N; m++) {
        if (std::gcd<int64>(m, 6) != 1) continue;
        ans += f(N / m);
    }
    return ans;
}

static void run_validation() {
    {
        int64 N = 6;
        int64 fast = solve_fast(N);
        int64 brute = brute_F_small(N);
        assert(fast == brute);
        assert(fast == 5);
    }
    {
        int64 N = 20;
        int64 fast = solve_fast(N);
        int64 brute = brute_F_small(N);
        assert(fast == brute);
        assert(fast == 19);
    }

    vector<int64> tests = {1,2,3,4,5,7,8,9,10,50,100,200,400,1000,2000,5000};
    for (auto N : tests) {
        int64 fast = solve_fast(N);
        int64 brute = brute_F_small(N);
        if (fast != brute) {
            cerr << "Validation failed at N=" << N
                 << " fast=" << fast << " brute=" << brute << "\n";
            assert(false);
        }
    }
}
#endif

int main() {
#ifdef VALIDATE
    run_validation();
#endif

    const int64 N = 10000000000000000LL;
    cout << solve_fast(N) << "\n";
    return 0;
}

Python

def count_coprime6(x):
    if x <= 0:
        return 0
    return x - x // 2 - x // 3 + x // 6

def generate_smooth(N):
    v = []
    p2 = 1
    while p2 <= N:
        cur = p2
        while cur <= N:
            v.append(cur)
            if cur > N // 3:
                break
            cur *= 3
        if p2 > N // 2:
            break
        p2 *= 2
    return sorted(list(set(v)))

def generate_bad(N):
    b = []
    def push_if(x):
        if x <= N:
            b.append(x)
            
    push_if(6)
    push_if(24)
    push_if(54)
    
    x = 384
    while x <= N:
        b.append(x)
        if x > N // 8:
            break
        x *= 8
        
    x = 243
    while x <= N:
        b.append(x)
        if x > N // 27:
            break
        x *= 27
        
    return sorted(list(set(b)))

def solve_fast(N):
    smooth = generate_smooth(N)
    bad = generate_bad(N)
    
    fAt = [0] * len(smooth)
    ib = 0
    f = 0
    for i in range(len(smooth)):
        if ib < len(bad) and smooth[i] == bad[ib]:
            ib += 1
        else:
            f += 1
        fAt[i] = f
        
    ans = 0
    n = len(smooth)
    for i in range(n):
        L = smooth[i]
        R = smooth[i + 1] - 1 if i + 1 < n else N
        if R > N:
            R = N
            
        mLow = N // (R + 1) + 1
        mHigh = N // L
        if mLow > mHigh:
            continue
            
        cnt = count_coprime6(mHigh) - count_coprime6(mLow - 1)
        ans += fAt[i] * cnt
        
    return ans

def solve():
    ans = solve_fast(10000000000000000)
    return str(ans)

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

Java

import java.util.ArrayList;
import java.util.Collections;

public class Euler821 {

    static long countCoprime6(long x) {
        if (x <= 0)
            return 0;
        return x - x / 2 - x / 3 + x / 6;
    }

    static ArrayList<Long> generateSmooth(long N) {
        ArrayList<Long> v = new ArrayList<>();
        for (long p2 = 1; p2 <= N;) {
            for (long cur = p2; cur <= N;) {
                v.add(cur);
                if (cur > N / 3)
                    break;
                cur *= 3;
            }
            if (p2 > N / 2)
                break;
            p2 *= 2;
        }
        Collections.sort(v);
        ArrayList<Long> res = new ArrayList<>();
        for (long x : v) {
            if (res.isEmpty() || res.get(res.size() - 1) != x) {
                res.add(x);
            }
        }
        return res;
    }

    static ArrayList<Long> generateBad(long N) {
        ArrayList<Long> b = new ArrayList<>();
        if (6 <= N)
            b.add(6L);
        if (24 <= N)
            b.add(24L);
        if (54 <= N)
            b.add(54L);

        long x = 384;
        while (x <= N) {
            b.add(x);
            if (x > N / 8)
                break;
            x *= 8;
        }

        x = 243;
        while (x <= N) {
            b.add(x);
            if (x > N / 27)
                break;
            x *= 27;
        }

        Collections.sort(b);
        ArrayList<Long> res = new ArrayList<>();
        for (long val : b) {
            if (res.isEmpty() || res.get(res.size() - 1) != val) {
                res.add(val);
            }
        }
        return res;
    }

    public static String solve() {
        long N = 10000000000000000L;
        ArrayList<Long> smooth = generateSmooth(N);
        ArrayList<Long> bad = generateBad(N);

        long[] fAt = new long[smooth.size()];
        int ib = 0;
        long f = 0;
        for (int i = 0; i < smooth.size(); i++) {
            if (ib < bad.size() && smooth.get(i).equals(bad.get(ib))) {
                ib++;
            } else {
                f++;
            }
            fAt[i] = f;
        }

        // use BigInteger or java.math if it overflows, but maximum is sum of fAt * cnt
        // The max of fAt is smooth.size() which is around ~1200
        // max of cnt is N/L. So max sum happens at large N/L * fAt.
        // N = 10^16. Max product is roughly 1200 * 10^16 = 1.2 * 10^19.
        // Needs java.math.BigInteger, wait. The answer in C++ used `__int128_t`!
        // So we must use BigInteger for the `ans`.
        java.math.BigInteger ans = java.math.BigInteger.ZERO;

        int n = smooth.size();
        for (int i = 0; i < n; i++) {
            long L = smooth.get(i);
            long R = (i + 1 < n) ? smooth.get(i + 1) - 1 : N;
            if (R > N)
                R = N;

            long mLow = N / (R + 1) + 1;
            long mHigh = N / L;
            if (mLow > mHigh)
                continue;

            long cnt = countCoprime6(mHigh) - countCoprime6(mLow - 1);

            java.math.BigInteger term = java.math.BigInteger.valueOf(fAt[i])
                    .multiply(java.math.BigInteger.valueOf(cnt));
            ans = ans.add(term);
        }

        return ans.toString();
    }

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