Problem 963: Removing Trits

Project Euler Problem 963 Solution

EulerSolve provides an optimized solution for Project Euler Problem 963, Removing Trits, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For each integer $1\le n\le N$, the implementations do not work with $n$ directly. Instead they attach to $n$ a canonical ternary signature $V(n)$ that preserves the quantity relevant to the problem. After that compression, the task becomes a counting problem: how many unordered pairs $\{a,b\}$ and $\{c,d\}$ with entries in $[1,N]$ have the same combined signature? Brute force over all quadruples is impossible. Even brute force over all $\binom{N+1}{2}$ unordered pairs would be far too large for $N=100000$. The key observation in the code is that many integers share the same signature, the signature is built recursively from the ternary expansion, and pair equality can therefore be counted class-by-class instead of number-by-number. Mathematical Approach Write the signature as $$V(n)=\bigl(f(n),\varepsilon(n),\mathbf{u}(n)\bigr),$$ where $f(n)$ is a nonpositive dyadic rational, $\varepsilon(n)\in\{0,1\}$ is a parity bit, and $\mathbf{u}(n)$ is a finite vector of nonnegative integers. The whole solution is built around the fact that this signature can be computed recursively from the ternary digits of $n$, and that pair signatures combine componentwise. The Ternary Signature Recurrence Let $n=3q+r$ with $r\in\{0,1,2\}$....

Detailed mathematical approach

Problem Summary

For each integer $1\le n\le N$, the implementations do not work with $n$ directly. Instead they attach to $n$ a canonical ternary signature $V(n)$ that preserves the quantity relevant to the problem. After that compression, the task becomes a counting problem: how many unordered pairs $\{a,b\}$ and $\{c,d\}$ with entries in $[1,N]$ have the same combined signature?

Brute force over all quadruples is impossible. Even brute force over all $\binom{N+1}{2}$ unordered pairs would be far too large for $N=100000$. The key observation in the code is that many integers share the same signature, the signature is built recursively from the ternary expansion, and pair equality can therefore be counted class-by-class instead of number-by-number.

Mathematical Approach

Write the signature as

$$V(n)=\bigl(f(n),\varepsilon(n),\mathbf{u}(n)\bigr),$$

where $f(n)$ is a nonpositive dyadic rational, $\varepsilon(n)\in\{0,1\}$ is a parity bit, and $\mathbf{u}(n)$ is a finite vector of nonnegative integers. The whole solution is built around the fact that this signature can be computed recursively from the ternary digits of $n$, and that pair signatures combine componentwise.

The Ternary Signature Recurrence

Let $n=3q+r$ with $r\in\{0,1,2\}$. Start from

$$V(0)=\bigl(0,0,\mathbf{0}\bigr).$$

Then the recurrence used by all three implementations is

$$V(3q+1)=\bigl(f(q)-1,\varepsilon(q),\mathbf{u}(q)\bigr),$$

$$V(3q+2)=\bigl(f(q),\varepsilon(q)\oplus 1,\mathbf{u}(q)\bigr),$$

and

$$V(3q)=\begin{cases} \bigl(f(q)+1,\varepsilon(q),\mathbf{u}(q)\bigr), & f(q)\le -2,\\ \bigl(f(q)/2,\varepsilon(q),\mathbf{u}(q)\bigr), & -2<f(q)<0,\\ Z(3q), & f(q)=0. \end{cases}$$

The code stores $f(n)$ after multiplying by a fixed power of two, so that every subtraction, increment, and halving stays exact in integer arithmetic. This works because the first component is always built from the operations $-1$, $+1$, and division by 2, so only dyadic rational numbers can appear.

When the Special Branch Appears

The branch $Z(3q)$ is not an arbitrary exception. It is reached exactly when the ternary expansion of $n$ contains only the digits $0$ and $2$.

The proof is a short induction. From a state with $f=0$, appending a ternary digit $2$ keeps $f=0$, while appending a digit $0$ also keeps the first component at 0 via the special branch. Appending a digit $1$ makes the first component negative immediately. Once $f<0$, none of the three transitions can return it to 0: the $r=1$ branch decreases it, the $r=2$ branch leaves it unchanged, and the $r=0$ branch either adds 1 to a value at most $-2$ or halves a value in $(-2,0)$, both of which remain negative. Therefore

Equivalently, $f(n)=0$ if and only if every ternary digit of $n$ belongs to $\{0,2\}$.

This is why the solution splits the problem into two regimes: ordinary numbers with $f(n)<0$, and the special $0/2$-digit numbers where the scalar part vanishes and a more detailed combinatorial profile is needed.

The Zero-Count Profile on $0/2$-Digit Numbers

Suppose $n=\sum_{i\ge 0} t_i 3^i$ and every ternary digit satisfies $t_i\in\{0,2\}$. Then the implementations encode the remaining information as a parity bit together with a vector that records how zeros are distributed to the right of each digit $2$.

The parity bit is simply

$$\varepsilon(n)\equiv \#\{i:t_i=2\}\pmod 2.$$

For $k\ge 1$, define

$$u_k(n)=\#\Bigl\{i:t_i=2,\ \#\{j<i:t_j=0\}=k\Bigr\}.$$

So $u_k(n)$ counts how many digits $2$ have exactly $k$ zeros somewhere to their right. Digits $2$ with no zero to their right affect only the parity bit and do not contribute to the stored vector. The implementations reconstruct this profile by scanning the ternary expansion from least significant trit to most significant trit, keeping a running count of zeros already seen.

A few concrete examples make the rule transparent. For $18=200_3$, the parity bit is 1 and the only nonzero vector entry is $u_2=1$. For $20=202_3$, the parity bit is 0 and the only nonzero entry is $u_1=1$, because one digit $2$ has one zero to its right and the other has none. For $24=220_3$, the parity bit is 0 and the only nonzero entry is $u_1=2$, because both digits $2$ have exactly one zero to their right.

Combining Two Signatures

Once a number has been replaced by its signature, a pair $\{a,b\}$ is represented by the componentwise combination

$$V(a,b)=\bigl(f(a)+f(b),\ \varepsilon(a)\oplus\varepsilon(b),\ \mathbf{u}(a)+\mathbf{u}(b)\bigr).$$

If a signature class $K$ occurs $F(K)$ times among $1,2,\dots,N$, then there is no reason to enumerate all members of that class separately. The number of unordered pairs contributed by one class on the diagonal is

$$\binom{F(K)+1}{2},$$

while two distinct classes $K_1$ and $K_2$ contribute

$$F(K_1)F(K_2).$$

Accumulating these weights by the resulting pair signature produces a distribution $C(T)$:

$$C(T)=\#\bigl\{\{a,b\}:1\le a\le b\le N,\ V(a,b)=T\bigr\}.$$

The final answer is then

$$A(N)=\sum_T C(T)^2.$$

In other words, the program counts ordered pairs of unordered pairs that land in the same combined signature class.

Worked Example: The Validation Case $N=5$

For the first five integers, the signatures are

$$V(1)=(-1,0,\mathbf{0}),\quad V(2)=(0,1,\mathbf{0}),\quad V(3)=(-1/2,0,\mathbf{0}),$$

$$V(4)=(-2,0,\mathbf{0}),\quad V(5)=(-1,1,\mathbf{0}).$$

There are 15 unordered pairs $\{a,b\}$ with $1\le a\le b\le 5$, but only 12 distinct combined signatures. The collisions are

$$V(1,1)=V(5,5),\qquad V(1,5)=V(2,4),\qquad V(2,5)=V(3,3).$$

So three signature classes occur twice and nine classes occur once. Therefore

$$A(5)=3\cdot 2^2+9\cdot 1^2=21,$$

which is exactly the value used by the implementations as their built-in sanity check.

How the Code Works

Building the Signature Table

The C++, Python, and Java implementations compute $V(n)$ for every $0\le n\le N$ in increasing order. Because each transition only depends on the quotient $\lfloor n/3\rfloor$, this is a simple dynamic-programming pass over the integers. The dyadic scalar is stored as a scaled integer, the parity bit is stored explicitly, and the zero-count profile is kept in a short fixed-length vector that is comfortably large enough for $N=100000$.

The only non-constant-time update is the special $0/2$-digit branch for multiples of 3 with $f=0$. In that case the implementation rescans the ternary digits of the current number to rebuild the zero-count profile. This is deliberate: once the scalar part has collapsed to zero, appending a trailing zero can shift every recorded zero count, so rebuilding from the digits is simpler and safer than trying to update the compressed profile locally.

Compressing Equal Signatures and Counting Pair Classes

After the table of $V(n)$ values is available, the implementations build a frequency map $F(K)$ for the distinct signature classes. This compression is substantial: for the target input, the 100000 integers collapse to only a few thousand distinct signatures, so the expensive part of the search moves from numbers to classes.

The next pass loops over unordered pairs of signature classes. For each diagonal class it adds $\binom{F(K)+1}{2}$, and for each off-diagonal class pair it adds $F(K_1)F(K_2)$, to the bucket indexed by the combined signature. Once that map has been built, the answer is the sum of the squares of all bucket sizes.

Because the final count is much larger than 64-bit signed range, the implementations use an extended integer type for the last accumulation stage: native arbitrary precision in Python, a big integer type in Java, and a 128-bit unsigned integer in C++.

Complexity Analysis

Let $U$ be the number of distinct signatures among $1,\dots,N$. Building the main table is essentially linear in $N$, with occasional extra $O(\log_3 N)$ rescans on the special $0/2$-digit branch. The compression phase is $O(N)$. The dominant counting phase is $O(U^2)$, because every unordered pair of signature classes is processed once, and each combination uses only fixed-size arithmetic on the three signature components.

Memory usage is $O(N)$ for the table of signatures together with $O(U)$ and $O(M)$ for the frequency map and the pair-signature map, where $M$ is the number of distinct combined signatures that actually occur. For the target input this is practical precisely because $U$ is far smaller than $N$, so the program never needs to examine all integer pairs directly.

Footnotes and References

Problem page: https://projecteuler.net/problem=963
Ternary numeral system: Wikipedia - Ternary numeral system
Dyadic rational: Wikipedia - Dyadic rational
Cantor set and ternary digits $0$ and $2$: Wikipedia - Cantor set
Discrete convolution as a counting pattern on class frequencies: Wikipedia - Convolution

Problem 963 source code

C++

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

namespace {

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

constexpr int kMaxK = 24;
constexpr i64 kScale = 1LL << 20;

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

struct ValueKey {
    i64 f_scaled = 0;
    std::uint8_t star = 0;
    std::array<std::int16_t, kMaxK + 1> up{};

    bool operator==(const ValueKey& other) const {
        return f_scaled == other.f_scaled && star == other.star && up == other.up;
    }
};

struct ValueHash {
    std::size_t operator()(const ValueKey& v) const noexcept {
        std::uint64_t h = 1469598103934665603ULL;
        const auto mix = [&](std::uint64_t x) {
            h ^= x;
            h *= 1099511628211ULL;
        };
        mix(static_cast<std::uint64_t>(v.f_scaled));
        mix(static_cast<std::uint64_t>(v.star));
        for (std::int16_t c : v.up) {
            mix(static_cast<std::uint16_t>(c));
        }
        return static_cast<std::size_t>(h);
    }
};

ValueKey add_values(const ValueKey& a, const ValueKey& b) {
    ValueKey out;
    out.f_scaled = a.f_scaled + b.f_scaled;
    out.star = static_cast<std::uint8_t>(a.star ^ b.star);
    for (int i = 0; i <= kMaxK; ++i) {
        out.up[static_cast<std::size_t>(i)] =
            static_cast<std::int16_t>(a.up[static_cast<std::size_t>(i)] + b.up[static_cast<std::size_t>(i)]);
    }
    return out;
}

ValueKey value_for_zero_fraction_case(int n3) {
    ValueKey out;
    int zeros_to_right = 0;
    int x = n3;
    while (x > 0) {
        const int d = x % 3;
        if (d == 0) {
            ++zeros_to_right;
        } else {
            assert(d == 2);
            out.star ^= 1U;
            if (zeros_to_right <= kMaxK) {
                ++out.up[static_cast<std::size_t>(zeros_to_right)];
            }
        }
        x /= 3;
    }
    return out;
}

std::vector<ValueKey> build_values(int limit) {
    std::vector<ValueKey> values(static_cast<std::size_t>(limit + 1));
    values[0] = ValueKey{};

    for (int n = 1; n <= limit; ++n) {
        const int q = n / 3;
        const int r = n % 3;
        ValueKey cur = values[static_cast<std::size_t>(q)];

        if (r == 1) {
            cur.f_scaled -= kScale;
        } else if (r == 2) {
            cur.star ^= 1U;
        } else {
            if (cur.f_scaled <= -2 * kScale) {
                cur.f_scaled += kScale;
            } else if (cur.f_scaled < 0) {
                cur.f_scaled /= 2;
            } else {
                assert(cur.f_scaled == 0);
                cur = value_for_zero_fraction_case(n);
            }
        }

        values[static_cast<std::size_t>(n)] = cur;
    }
    return values;
}

u128 solve_fast(int limit) {
    const std::vector<ValueKey> values = build_values(limit);

    std::unordered_map<ValueKey, u64, ValueHash> freq;
    freq.reserve(static_cast<std::size_t>(limit));
    for (int n = 1; n <= limit; ++n) {
        ++freq[values[static_cast<std::size_t>(n)]];
    }

    std::vector<ValueKey> uniq;
    uniq.reserve(freq.size());
    for (const auto& [key, _] : freq) {
        uniq.push_back(key);
    }

    std::vector<u64> f;
    f.reserve(uniq.size());
    for (const auto& key : uniq) {
        f.push_back(freq[key]);
    }

    std::unordered_map<ValueKey, u64, ValueHash> pair_count;
    pair_count.reserve(uniq.size() * 8ULL + 1ULL);

    for (std::size_t i = 0; i < uniq.size(); ++i) {
        const u64 fi = f[i];
        {
            const ValueKey s = add_values(uniq[i], uniq[i]);
            pair_count[s] += fi * (fi + 1ULL) / 2ULL;
        }
        for (std::size_t j = i + 1; j < uniq.size(); ++j) {
            const ValueKey s = add_values(uniq[i], uniq[j]);
            pair_count[s] += fi * f[j];
        }
    }

    u128 ans = 0U;
    for (const auto& [_, c] : pair_count) {
        ans += static_cast<u128>(c) * static_cast<u128>(c);
    }
    return ans;
}

}  // namespace

int main() {
    assert(solve_fast(5) == 21U);

    const u128 answer = solve_fast(100'000);
    std::cout << to_string_u128(answer) << '\n';
    return 0;
}

Python

def solve():
    MAX_K = 24; SCALE = 1 << 20
    limit = 100000

    def val_zero_frac(n3):
        star = 0; up = [0]*(MAX_K+1)
        zeros_right = 0; x = n3
        while x > 0:
            d = x % 3
            if d == 0: zeros_right += 1
            else:
                star ^= 1
                if zeros_right <= MAX_K: up[zeros_right] += 1
            x //= 3
        return (0, star, tuple(up))

    def build_values(lim):
        vals = [None]*(lim+1)
        vals[0] = (0, 0, tuple([0]*(MAX_K+1)))
        for n in range(1, lim+1):
            q, r = divmod(n, 3)
            fs, st, up = vals[q]; up = list(up)
            if r == 1: fs -= SCALE
            elif r == 2: st ^= 1
            else:
                if fs <= -2*SCALE: fs += SCALE
                elif fs < 0: fs //= 2
                else: vals[n] = val_zero_frac(n); continue
            vals[n] = (fs, st, tuple(up))
        return vals

    vals = build_values(limit)
    freq = {}
    for n in range(1, limit+1):
        v = vals[n]; freq[v] = freq.get(v, 0) + 1

    def add_vals(a, b):
        return (a[0]+b[0], a[1]^b[1], tuple(x+y for x,y in zip(a[2], b[2])))

    uniq = list(freq.keys())
    f = [freq[k] for k in uniq]
    pair_count = {}

    for i in range(len(uniq)):
        s = add_vals(uniq[i], uniq[i])
        pair_count[s] = pair_count.get(s, 0) + f[i]*(f[i]+1)//2
        for j in range(i+1, len(uniq)):
            s = add_vals(uniq[i], uniq[j])
            pair_count[s] = pair_count.get(s, 0) + f[i]*f[j]

    ans = sum(c*c for c in pair_count.values())
    return str(ans)

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

Java

import java.util.ArrayList;
import java.util.Arrays;
import java.util.HashMap;
import java.util.List;
import java.util.Map;

public class Euler963 {

    static final int K_MAX_K = 24;
    static final long K_SCALE = 1L << 20;

    static class ValueKey {
        long fScaled;
        int star;
        short[] up;

        ValueKey() {
            this.fScaled = 0;
            this.star = 0;
            this.up = new short[K_MAX_K + 1];
        }

        ValueKey(long fScaled, int star, short[] up) {
            this.fScaled = fScaled;
            this.star = star;
            this.up = up;
        }

        @Override
        public boolean equals(Object o) {
            if (this == o)
                return true;
            if (o == null || getClass() != o.getClass())
                return false;
            ValueKey valueKey = (ValueKey) o;
            return fScaled == valueKey.fScaled &&
                    star == valueKey.star &&
                    Arrays.equals(up, valueKey.up);
        }

        @Override
        public int hashCode() {
            long h = 1469598103934665603L;
            h ^= fScaled;
            h *= 1099511628211L;
            h ^= star;
            h *= 1099511628211L;
            for (short c : up) {
                h ^= c;
                h *= 1099511628211L;
            }
            return (int) (h ^ (h >>> 32));
        }

        public ValueKey copy() {
            return new ValueKey(this.fScaled, this.star, Arrays.copyOf(this.up, this.up.length));
        }
    }

    static ValueKey addValues(ValueKey a, ValueKey b) {
        short[] newUp = new short[K_MAX_K + 1];
        for (int i = 0; i <= K_MAX_K; ++i) {
            newUp[i] = (short) (a.up[i] + b.up[i]);
        }
        return new ValueKey(a.fScaled + b.fScaled, a.star ^ b.star, newUp);
    }

    static ValueKey valueForZeroFractionCase(int n3) {
        ValueKey out = new ValueKey();
        int zerosToRight = 0;
        int x = n3;
        while (x > 0) {
            int d = x % 3;
            if (d == 0) {
                ++zerosToRight;
            } else {
                out.star ^= 1;
                if (zerosToRight <= K_MAX_K) {
                    out.up[zerosToRight]++;
                }
            }
            x /= 3;
        }
        return out;
    }

    static ValueKey[] buildValues(int limit) {
        ValueKey[] values = new ValueKey[limit + 1];
        values[0] = new ValueKey();

        for (int n = 1; n <= limit; ++n) {
            int q = n / 3;
            int r = n % 3;
            ValueKey cur = values[q].copy();

            if (r == 1) {
                cur.fScaled -= K_SCALE;
            } else if (r == 2) {
                cur.star ^= 1;
            } else {
                if (cur.fScaled <= -2 * K_SCALE) {
                    cur.fScaled += K_SCALE;
                } else if (cur.fScaled < 0) {
                    cur.fScaled /= 2;
                } else {
                    cur = valueForZeroFractionCase(n);
                }
            }
            values[n] = cur;
        }
        return values;
    }

    static String solveFast(int limit) {
        ValueKey[] values = buildValues(limit);

        Map<ValueKey, Long> freq = new HashMap<>();
        for (int n = 1; n <= limit; ++n) {
            freq.put(values[n], freq.getOrDefault(values[n], 0L) + 1L);
        }

        List<ValueKey> uniq = new ArrayList<>(freq.keySet());
        long[] f = new long[uniq.size()];
        for (int i = 0; i < uniq.size(); i++) {
            f[i] = freq.get(uniq.get(i));
        }

        Map<ValueKey, Long> pairCount = new HashMap<>();

        for (int i = 0; i < uniq.size(); ++i) {
            long fi = f[i];
            ValueKey sSelf = addValues(uniq.get(i), uniq.get(i));
            pairCount.put(sSelf, pairCount.getOrDefault(sSelf, 0L) + fi * (fi + 1L) / 2L);

            for (int j = i + 1; j < uniq.size(); ++j) {
                ValueKey s = addValues(uniq.get(i), uniq.get(j));
                pairCount.put(s, pairCount.getOrDefault(s, 0L) + fi * f[j]);
            }
        }

        java.math.BigInteger ans = java.math.BigInteger.ZERO;
        for (long c : pairCount.values()) {
            java.math.BigInteger bc = java.math.BigInteger.valueOf(c);
            ans = ans.add(bc.multiply(bc));
        }
        return ans.toString();
    }

    public static String solve() {
        return solveFast(100000);
    }

    public static void main(String[] args) {
        if (!solveFast(5).equals("21")) {
            System.out.println("Validation failed");
            return;
        }
        System.out.println(solve());
    }
}