Problem 581: $47$-smooth Triangular Numbers

Project Euler Problem 581 Solution

EulerSolve provides an optimized solution for Project Euler Problem 581, $47$-smooth Triangular Numbers, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary A positive integer is called $47$-smooth when every prime divisor belongs to $$P=\{2,3,5,7,11,13,17,19,23,29,31,37,41,43,47\}.$$ For the triangular numbers $$T(n)=\frac{n(n+1)}{2},$$ we must compute the sum of all indices $n$ for which $T(n)$ is $47$-smooth. The implementations do not factor every triangular number directly. Instead they rewrite each valid case in terms of one smooth core $x$ and one neighboring value $2x-1$ or $2x+1$. Mathematical Approach The search is built around the fixed prime set $P$ above and around the coprimality of consecutive integers. Step 1: Split the triangular number by parity If $n$ is even, write $n=2x$. Then $$T(n)=\frac{(2x)(2x+1)}{2}=x(2x+1).$$ If $n$ is odd, write $n=2x-1$. Then $$T(n)=\frac{(2x-1)(2x)}{2}=x(2x-1).$$ In both cases the two factors are coprime, because $$\gcd(x,2x+1)=1,\qquad \gcd(x,2x-1)=1.$$ Step 2: Turn smooth triangular numbers into smooth neighbors A product of two coprime integers is $47$-smooth if and only if each factor is $47$-smooth. Therefore $$T(2x)\text{ is }47\text{-smooth} \iff x\text{ and }2x+1\text{ are }47\text{-smooth},$$ $$T(2x-1)\text{ is }47\text{-smooth} \iff x\text{ and }2x-1\text{ are }47\text{-smooth}.$$ This converts the original problem into two tests around the same smooth core $x$....

Detailed mathematical approach

Problem Summary

A positive integer is called $47$-smooth when every prime divisor belongs to

$$P=\{2,3,5,7,11,13,17,19,23,29,31,37,41,43,47\}.$$

For the triangular numbers

$$T(n)=\frac{n(n+1)}{2},$$

we must compute the sum of all indices $n$ for which $T(n)$ is $47$-smooth. The implementations do not factor every triangular number directly. Instead they rewrite each valid case in terms of one smooth core $x$ and one neighboring value $2x-1$ or $2x+1$.

Mathematical Approach

The search is built around the fixed prime set $P$ above and around the coprimality of consecutive integers.

Step 1: Split the triangular number by parity

If $n$ is even, write $n=2x$. Then

$$T(n)=\frac{(2x)(2x+1)}{2}=x(2x+1).$$

If $n$ is odd, write $n=2x-1$. Then

$$T(n)=\frac{(2x-1)(2x)}{2}=x(2x-1).$$

In both cases the two factors are coprime, because

$$\gcd(x,2x+1)=1,\qquad \gcd(x,2x-1)=1.$$

Step 2: Turn smooth triangular numbers into smooth neighbors

A product of two coprime integers is $47$-smooth if and only if each factor is $47$-smooth. Therefore

$$T(2x)\text{ is }47\text{-smooth} \iff x\text{ and }2x+1\text{ are }47\text{-smooth},$$

$$T(2x-1)\text{ is }47\text{-smooth} \iff x\text{ and }2x-1\text{ are }47\text{-smooth}.$$

This converts the original problem into two tests around the same smooth core $x$. Every valid index $n$ determines $x$ uniquely from its parity, so there is no double counting between the even and odd branches.

Step 3: Enumerate every possible smooth core $x$

Every $47$-smooth integer has a unique factorization

$$x=\prod_{p\in P} p^{e_p},\qquad e_p\ge 0.$$

The implementations recurse over the primes in $P$, choose each exponent $e_p$ in turn, and multiply the running product by successive powers of that prime while the product stays below the fixed search cutoff

$$B=10^{13}.$$

Because of unique factorization, this depth-first enumeration visits each $47$-smooth $x\le B$ exactly once.

Step 4: Test the two neighbors $2x-1$ and $2x+1$

For each enumerated $x$, the implementations form

$$y_-=2x-1,\qquad y_+=2x+1.$$

To check whether one of these numbers is $47$-smooth, the algorithm repeatedly divides it by every prime in $P$. If the remaining residue is $1$, then all prime factors were allowed and the number is $47$-smooth. If the residue is larger than $1$, then some prime factor exceeds $47$, so the candidate must be rejected.

Whenever $y_+$ passes, the even index $n=2x$ is added to the total. Whenever $y_-$ passes, the odd index $n=2x-1$ is added.

Worked Example

Take $x=3$. This $x$ is $47$-smooth because its only prime factor is $3$.

Then

$$2x-1=5,\qquad 2x+1=7,$$

and both numbers are also $47$-smooth. So the algorithm accepts both

$$n=5\qquad\text{and}\qquad n=6.$$

Indeed,

$$T(5)=\frac{5\cdot 6}{2}=15=3\cdot 5,$$

$$T(6)=\frac{6\cdot 7}{2}=21=3\cdot 7,$$

and every prime factor in both triangular numbers is at most $47$.

Step 5: Why the enumeration matches the required sum

Every accepted candidate produces a valid index, because the parity split reconstructs

$$T(2x)=x(2x+1),\qquad T(2x-1)=x(2x-1),$$

with both factors already verified to be $47$-smooth. Conversely, every valid index $n$ belongs to exactly one of the two parity cases, so it yields one smooth core $x$ and one successful neighbor test. The implementations therefore compute the desired sum by exhaustively enumerating the smooth cores in their chosen search range and adding precisely the matching indices.

How the Code Works

The C++, Python, and Java implementations all follow the same structure. First they store the fifteen allowed primes up to $47$. Next they perform a recursive search over exponent choices, so the current product always stays $47$-smooth by construction.

When the recursion has processed all fifteen primes, the current product is one completed candidate $x$. At that point the implementation evaluates $2x-1$ and $2x+1$, checks each one by repeated division with the same prime list, and adds the corresponding index or indices to a running total.

The recursion only branches while another multiplication still fits under the fixed cutoff $10^{13}$. Since the arithmetic never leaves this range, ordinary 64-bit integer types are enough for the search and for the accumulated sum.

Complexity Analysis

Let $S(B)$ denote the number of $47$-smooth integers $x\le B$, where $B=10^{13}$ is the cutoff used by the implementations. The recursive generator visits each such $x$ once, and the recursion depth is the fixed constant $15$.

For every generated $x$, the code performs two smoothness checks, each over the same fixed prime set $P$. Repeated divisions by prime powers add only a small extra factor in practice, so the total running time is essentially proportional to the number of enumerated smooth candidates. With $47$ fixed, the memory usage is $O(1)$ apart from the recursion stack, which also has constant depth.

Footnotes and References

Problem page: https://projecteuler.net/problem=581
Smooth number: Wikipedia — Smooth number
Triangular number: Wikipedia — Triangular number
Coprime integers: Wikipedia — Coprime integers
Fundamental theorem of arithmetic: Wikipedia — Fundamental theorem of arithmetic

Problem 581 source code

C++

#include <cstdint>
#include <iostream>
#include <vector>

// Project Euler 581: 47-smooth triangular numbers
//
// T(n) = n(n+1)/2 is 47-smooth  <=>  its prime factors are all <= 47.
// Since gcd(n, n+1)=1, after dividing by 2 the two coprime factors must each be 47-smooth:
// - If n is even:  T(n) = (n/2) * (n+1)  =>  n/2 and n+1 are 47-smooth.
// - If n is odd :  T(n) = n * ((n+1)/2) =>  n and (n+1)/2 are 47-smooth.
//
// Let x be the 47-smooth factor:
// - n = 2x      requires 2x+1 to be 47-smooth.
// - n = 2x-1    requires 2x-1 to be 47-smooth.
//
// So we enumerate all 47-smooth x up to a safe bound and test 2x±1.

using u64 = std::uint64_t;

static constexpr u64 LIMIT_X = 10'000'000'000'000ULL;  // 1e13

static const int PRIMES[] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47};

static inline bool is_47_smooth(u64 x) {
    for (int p : PRIMES) {
        while (x % (u64)p == 0) x /= (u64)p;
    }
    return x == 1;
}

static u64 sum_n = 0;

static void dfs(int idx, u64 cur) {
    if (idx == (int)(sizeof(PRIMES) / sizeof(PRIMES[0]))) {
        const u64 x = cur;
        const u64 y_plus = 2 * x + 1;
        const u64 y_minus = 2 * x - 1;

        if (is_47_smooth(y_plus)) sum_n += 2 * x;
        if (is_47_smooth(y_minus)) sum_n += 2 * x - 1;
        return;
    }
    const u64 p = (u64)PRIMES[idx];
    u64 v = cur;
    while (v <= LIMIT_X) {
        dfs(idx + 1, v);
        if (v > LIMIT_X / p) break;
        v *= p;
    }
}

int main() {
    dfs(0, 1);
    std::cout << sum_n << "\n";
    return 0;
}

Python

def solve():
    PRIMES = [2,3,5,7,11,13,17,19,23,29,31,37,41,43,47]
    LIMIT = 10000000000000
    total = 0

    def is_smooth(x):
        for p in PRIMES:
            while x % p == 0: x //= p
        return x == 1

    def dfs(idx, cur):
        nonlocal total
        if idx == len(PRIMES):
            x = cur
            if is_smooth(2*x+1): total += 2*x
            if is_smooth(2*x-1): total += 2*x-1
            return
        p = PRIMES[idx]; v = cur
        while v <= LIMIT:
            dfs(idx+1, v)
            if v > LIMIT // p: break
            v *= p

    dfs(0, 1)
    return str(total)

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

Java

public class Euler581 {
    static final long LIMIT_X = 10_000_000_000_000L;
    static final int[] PRIMES = { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47 };

    static boolean is_47_smooth(long x) {
        for (int p : PRIMES) {
            while (x % p == 0)
                x /= p;
        }
        return x == 1;
    }

    static long sum_n = 0;

    static void dfs(int idx, long cur) {
        if (idx == PRIMES.length) {
            long y_plus = 2 * cur + 1;
            long y_minus = 2 * cur - 1;
            if (is_47_smooth(y_plus))
                sum_n += 2 * cur;
            if (is_47_smooth(y_minus))
                sum_n += 2 * cur - 1;
            return;
        }
        long p = PRIMES[idx];
        long v = cur;
        while (v <= LIMIT_X) {
            dfs(idx + 1, v);
            if (v > LIMIT_X / p)
                break;
            v *= p;
        }
    }

    public static String solve() {
        sum_n = 0;
        dfs(0, 1);
        return Long.toString(sum_n);
    }

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