Problem 462: Permutation of 3-smooth Numbers

Project Euler Problem 462 Solution

EulerSolve provides an optimized solution for Project Euler Problem 462, Permutation of 3-smooth Numbers, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary A 3-smooth number is a positive integer of the form $2^a3^b$ with $a,b\ge 0$. For a given $n$, define $$\mathcal{S}(n)=\{2^a3^b \le n : a,b \in \mathbb{Z}_{\ge 0}\}.$$ We must count the permutations of $\mathcal{S}(n)$ with the rule that whenever $x\mid y$, the number $x$ must appear earlier than $y$. If $F(n)$ denotes that count, the checkpoints used by the implementations are $$F(6)=5,\qquad F(8)=9,\qquad F(20)=450,\qquad F(1000)=8.8521816557\times 10^{21}.$$ The key simplification is that divisibility among numbers $2^a3^b$ is just coordinatewise order on the exponent pair $(a,b)$. That turns the problem into counting linear extensions of a finite Ferrers-shaped poset. Mathematical Approach Represent each element of $\mathcal{S}(n)$ by its exponent pair. Then $$2^a3^b \mid 2^c3^d \iff a\le c,\ b\le d.$$ So the admissible permutations are exactly the linear extensions of the set of lattice points $(a,b)$ satisfying $2^a3^b\le n$. Step 1: Convert the exponent set into a partition Fix a value of $b$. The admissible values of $a$ are those with $2^a\le n/3^b$. Therefore the row length is $$\lambda_{b+1}=1+\left\lfloor \log_2\left(\frac{n}{3^b}\right)\right\rfloor,$$ for every $b$ such that $3^b\le n$. As $b$ increases, $n/3^b$ decreases, so the row lengths are nonincreasing....

Detailed mathematical approach

Problem Summary

A 3-smooth number is a positive integer of the form $2^a3^b$ with $a,b\ge 0$. For a given $n$, define

$$\mathcal{S}(n)=\{2^a3^b \le n : a,b \in \mathbb{Z}_{\ge 0}\}.$$

We must count the permutations of $\mathcal{S}(n)$ with the rule that whenever $x\mid y$, the number $x$ must appear earlier than $y$. If $F(n)$ denotes that count, the checkpoints used by the implementations are

$$F(6)=5,\qquad F(8)=9,\qquad F(20)=450,\qquad F(1000)=8.8521816557\times 10^{21}.$$

The key simplification is that divisibility among numbers $2^a3^b$ is just coordinatewise order on the exponent pair $(a,b)$. That turns the problem into counting linear extensions of a finite Ferrers-shaped poset.

Mathematical Approach

Represent each element of $\mathcal{S}(n)$ by its exponent pair. Then

$$2^a3^b \mid 2^c3^d \iff a\le c,\ b\le d.$$

So the admissible permutations are exactly the linear extensions of the set of lattice points $(a,b)$ satisfying $2^a3^b\le n$.

Step 1: Convert the exponent set into a partition

Fix a value of $b$. The admissible values of $a$ are those with $2^a\le n/3^b$. Therefore the row length is

$$\lambda_{b+1}=1+\left\lfloor \log_2\left(\frac{n}{3^b}\right)\right\rfloor,$$

for every $b$ such that $3^b\le n$. As $b$ increases, $n/3^b$ decreases, so the row lengths are nonincreasing. Hence they form a partition

$$\lambda=(\lambda_1,\lambda_2,\dots,\lambda_R),\qquad R=1+\left\lfloor \log_3 n\right\rfloor.$$

Geometrically, the cell in row $b+1$ and column $a+1$ corresponds to the number $2^a3^b$.

Step 2: Interpret valid permutations as standard tableaux

Assign to every cell the position of its number inside the permutation. The divisibility rule says that if one cell is weakly north-west of another, its label must be smaller. Therefore a valid permutation is equivalent to a filling of the Ferrers diagram with $1,2,\dots,N$ that increases from left to right and from top to bottom, where

$$N=|\mathcal{S}(n)|=\sum_{r=1}^{R}\lambda_r.$$

That is exactly a standard Young tableau of shape $\lambda$.

Step 3: Apply the hook-length formula

The number of standard Young tableaux of shape $\lambda$ is

$$F(n)=\frac{N!}{\prod_{(r,c)\in\lambda} h_{r,c}}.$$

If $\lambda'_c$ is the height of column $c$, then the hook length of cell $(r,c)$ is

$$h_{r,c}=(\lambda_r-c)+(\lambda'_c-r)+1=\lambda_r+\lambda'_c-r-c+1.$$

In words, the hook counts the cell itself, every cell to its right in the same row, and every cell below it in the same column.

Step 4: Why this formula matches the divisibility poset

The exponent pairs form a left-justified diagram, and moving one step right multiplies by $2$ while moving one step down multiplies by $3$. Both moves increase the number, so both moves must also increase the permutation label. The partial order is therefore the usual row-and-column order of a Ferrers diagram. Hook-length theory applies directly, which removes any need to enumerate permutations.

Step 5: Worked example for $n=8$

The 3-smooth numbers up to $8$ are

$$1,\ 2,\ 3,\ 4,\ 6,\ 8.$$

The corresponding Ferrers shape is

$$\lambda=(4,2),$$

because the row for $b=0$ contains $1,2,4,8$ and the row for $b=1$ contains $3,6$. The hook lengths are

$$\begin{array}{cccc} 5 & 4 & 2 & 1\\ 2 & 1 \end{array}$$

$$F(8)=\frac{6!}{5\cdot 4\cdot 2\cdot 1\cdot 2\cdot 1}=\frac{720}{80}=9.$$

This matches the checkpoint and shows the full method on a small instance.

How the Code Works

The C++, Python, and Java implementations use the same pipeline. First they enumerate the rows by stepping through powers of $3$. For each row they determine the largest admissible power of $2$ with integer operations, so the row length is computed exactly without floating-point roundoff.

Next they sum all row lengths to obtain $N$, then compute the height of every column by counting how many rows reach that column. With these two arrays available, the implementation visits every cell once, evaluates its hook length, and multiplies all hooks into one exact denominator.

After that the implementation forms $N!$ as an arbitrary-precision integer and divides by the hook product to obtain $F(n)$. The C++ implementation also includes an exact-divisibility guard before the final division, which is a useful sanity check even though the theorem guarantees that the quotient is an integer. Finally the exact integer is rendered in scientific notation with 10 digits after the decimal point, using decimal-string rounding instead of floating-point arithmetic.

Complexity Analysis

Let $R=1+\lfloor \log_3 n\rfloor$ and $W=1+\lfloor \log_2 n\rfloor$. Constructing the row lengths costs $O(R)$. Computing column heights by scanning the row lengths costs $O(RW)$. Multiplying all hook lengths touches each cell once, which is $O(N)$ with

$$N=\sum_{r=1}^{R}\lambda_r=\Theta(\log^2 n).$$

Therefore the combinatorial part of the algorithm is $O(RW)=O(\log^2 n)$. The extra array storage is $O(R+W)$, while the dominant arithmetic cost comes from multiplying and dividing arbitrary-precision integers whose size grows with the final answer.

Footnotes and References

Problem page: https://projecteuler.net/problem=462
3-smooth numbers / regular numbers: Wikipedia — Regular number
Ferrers and Young diagrams: Wikipedia — Young diagram
Hook-length formula: Wikipedia — Hook-length formula
Linear extension of a poset: Wikipedia — Linear extension

Problem 462 source code

C++

#include <cstdint>
#include <iomanip>
#include <iostream>
#include <sstream>
#include <string>
#include <vector>

#include <boost/multiprecision/cpp_int.hpp>

namespace {

using u64 = std::uint64_t;
using boost::multiprecision::cpp_int;

struct Options {
    u64 n = 1'000'000'000'000'000'000ULL;
    bool run_checkpoints = true;
};

bool parse_u64_after_prefix(const std::string& arg, const std::string& prefix, u64& out) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }
    try {
        out = static_cast<u64>(std::stoull(tail));
    } catch (...) {
        return false;
    }
    return true;
}

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (parse_u64_after_prefix(arg, "--n=", options.n)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    if (options.n == 0ULL) {
        std::cerr << "--n must be positive.\n";
        return false;
    }
    return true;
}

std::vector<int> build_rows(const u64 n) {
    std::vector<int> rows;
    u64 p3 = 1ULL;
    while (p3 <= n) {
        const u64 bound = n / p3;
        const int max_a = 63 - __builtin_clzll(bound);
        rows.push_back(max_a + 1);
        if (p3 > n / 3ULL) {
            break;
        }
        p3 *= 3ULL;
    }
    return rows;
}

cpp_int count_exact(const u64 n) {
    const std::vector<int> rows = build_rows(n);
    const int max_width = rows.empty() ? 0 : rows.front();
    int cell_count = 0;
    for (const int w : rows) {
        cell_count += w;
    }

    std::vector<int> col_heights(static_cast<std::size_t>(max_width), 0);
    for (int c = 0; c < max_width; ++c) {
        int h = 0;
        for (const int w : rows) {
            if (w > c) {
                ++h;
            }
        }
        col_heights[static_cast<std::size_t>(c)] = h;
    }

    cpp_int numerator = 1;
    for (int x = 2; x <= cell_count; ++x) {
        numerator *= x;
    }

    cpp_int denominator = 1;
    for (int r = 0; r < static_cast<int>(rows.size()); ++r) {
        const int w = rows[static_cast<std::size_t>(r)];
        for (int c = 0; c < w; ++c) {
            const int right = w - c - 1;
            const int below = col_heights[static_cast<std::size_t>(c)] - r - 1;
            const int hook = right + below + 1;
            denominator *= hook;
        }
    }

    if ((numerator % denominator) != 0) {
        std::cerr << "Internal error: non-integral hook-length ratio.\n";
        return 0;
    }
    return numerator / denominator;
}

std::string format_scientific_10(const cpp_int& value) {
    std::string digits = value.convert_to<std::string>();
    long long exponent = static_cast<long long>(digits.size()) - 1LL;

    if (digits.size() < 12U) {
        digits.append(12U - digits.size(), '0');
    }

    u64 leading = 0ULL;
    for (int i = 0; i < 11; ++i) {
        leading = leading * 10ULL + static_cast<u64>(digits[static_cast<std::size_t>(i)] - '0');
    }
    const int round_digit = digits[11] - '0';
    if (round_digit >= 5) {
        ++leading;
    }
    if (leading == 100'000'000'000ULL) {
        leading = 10'000'000'000ULL;
        ++exponent;
    }

    const u64 integer_part = leading / 10'000'000'000ULL;
    const u64 fractional_part = leading % 10'000'000'000ULL;

    std::ostringstream out;
    out << integer_part << '.'
        << std::setw(10) << std::setfill('0') << fractional_part
        << 'e' << exponent;
    return out.str();
}

std::string solve_formatted(const u64 n) {
    return format_scientific_10(count_exact(n));
}

bool run_checkpoints() {
    if (count_exact(6ULL) != cpp_int(5)) {
        std::cerr << "Checkpoint failed: F(6)\n";
        return false;
    }
    if (count_exact(8ULL) != cpp_int(9)) {
        std::cerr << "Checkpoint failed: F(8)\n";
        return false;
    }
    if (count_exact(20ULL) != cpp_int(450)) {
        std::cerr << "Checkpoint failed: F(20)\n";
        return false;
    }
    if (solve_formatted(1000ULL) != "8.8521816557e21") {
        std::cerr << "Checkpoint failed: F(1000)\n";
        return false;
    }
    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }

    if (options.run_checkpoints && !run_checkpoints()) {
        return 1;
    }

    std::cout << solve_formatted(options.n) << '\n';
    return 0;
}

Python

import math

def count_exact(n):
    rows = []
    p3 = 1
    while p3 <= n:
        bound = n // p3
        max_a = bound.bit_length() - 1
        rows.append(max_a + 1)
        if p3 > n // 3:
            break
        p3 *= 3
    
    if not rows:
        return 0
    max_width = rows[0]
    cell_count = sum(rows)
    
    col_heights = [0] * max_width
    for c in range(max_width):
        h = sum(1 for w in rows if w > c)
        col_heights[c] = h
        
    numerator = math.factorial(cell_count)
    denominator = 1
    for r in range(len(rows)):
        w = rows[r]
        for c in range(w):
            right = w - c - 1
            below = col_heights[c] - r - 1
            hook = right + below + 1
            denominator *= hook
            
    return numerator // denominator

def format_scientific_10(value):
    digits = str(value)
    exponent = len(digits) - 1
    if len(digits) < 12:
        digits += '0' * (12 - len(digits))
        
    leading = int(digits[:11])
    round_digit = int(digits[11])
    if round_digit >= 5:
        leading += 1
        
    if leading == 100000000000:
        leading = 10000000000
        exponent += 1
        
    integer_part = leading // 10000000000
    fractional_part = leading % 10000000000
    return f"{integer_part}.{fractional_part:010d}e{exponent}"

def solve():
    n = 1000000000000000000
    return format_scientific_10(count_exact(n))

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

Java

import java.math.BigInteger;
import java.util.ArrayList;
import java.util.List;

public class Euler462 {
    private static BigInteger countExact(long n) {
        List<Integer> rows = new ArrayList<>();
        long p3 = 1;
        while (p3 <= n) {
            long bound = n / p3;
            int maxA = 63 - Long.numberOfLeadingZeros(bound);
            rows.add(maxA + 1);
            if (p3 > n / 3)
                break;
            p3 *= 3;
        }

        if (rows.isEmpty())
            return BigInteger.ZERO;

        int maxWidth = rows.get(0);
        int cellCount = 0;
        for (int w : rows)
            cellCount += w;

        int[] colHeights = new int[maxWidth];
        for (int c = 0; c < maxWidth; ++c) {
            int h = 0;
            for (int w : rows) {
                if (w > c)
                    h++;
            }
            colHeights[c] = h;
        }

        BigInteger numerator = BigInteger.ONE;
        for (int x = 2; x <= cellCount; ++x) {
            numerator = numerator.multiply(BigInteger.valueOf(x));
        }

        BigInteger denominator = BigInteger.ONE;
        for (int r = 0; r < rows.size(); ++r) {
            int w = rows.get(r);
            for (int c = 0; c < w; ++c) {
                int right = w - c - 1;
                int below = colHeights[c] - r - 1;
                int hook = right + below + 1;
                denominator = denominator.multiply(BigInteger.valueOf(hook));
            }
        }

        return numerator.divide(denominator);
    }

    private static String formatScientific10(BigInteger value) {
        String digits = value.toString();
        long exponent = digits.length() - 1;

        if (digits.length() < 12) {
            StringBuilder sb = new StringBuilder(digits);
            while (sb.length() < 12)
                sb.append('0');
            digits = sb.toString();
        }

        long leading = 0;
        for (int i = 0; i < 11; ++i) {
            leading = leading * 10 + (digits.charAt(i) - '0');
        }
        int roundDigit = digits.charAt(11) - '0';
        if (roundDigit >= 5) {
            leading++;
        }
        if (leading == 100000000000L) {
            leading = 10000000000L;
            exponent++;
        }

        long integerPart = leading / 10000000000L;
        long fractionalPart = leading % 10000000000L;

        return String.format("%d.%010de%d", integerPart, fractionalPart, exponent);
    }

    public static String solve() {
        return formatScientific10(countExact(1000000000000000000L));
    }

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