Project Euler

Arithmetic Derivative

The arithmetic derivative n' is defined on nonneg integers by: 1. 0' = 0, 1' = 0. 2. p' = 1 for every prime p. 3. (ab)' = a'b + ab' (Leibniz product rule). Compute sum_(n=1)^N n' and study the stru...

Source sync Apr 19, 2026

Problem #0895

Level Level 39

Solved By 113

Languages C++, Python

Answer 670785433

Length 334 words

number_theoryanalytic_mathmodular_arithmetic

Gary and Sally play a game using gold and silver coins arranged into a number of vertical stacks, alternating turns. On Gary's turn he chooses a gold coin and removes it from the game along with any other coins sitting on top. Sally does the same on her turn by removing a silver coin. The first player unable to make a move loses.

An arrangement is called fair if the person moving first, whether it be Gary or Sally, will lose the game if both play optimally.

An arrangement is called balanced if the number of gold and silver coins are equal.

Define $G(m)$ to be the number of fair and balanced arrangements consisting of three non-empty stacks, each not exceeding $m$ in size. Different orderings of the stacks are to be counted separately, so $G(2) = 6$ due to the following six arrangements:

Problem illustration

You are also given $G(5) = 348$ and $G(20) = 125825982708$.

Find $G(9898)$ giving your answer modulo $989898989$.

Problem 895: Arithmetic Derivative

Mathematical Foundation

Theorem 1 (Prime Power Formula). For any prime $p$ and integer $a \geq 1$ :

(p^a)' = a \cdot p^{a-1}.

Proof. By induction on $a$ .

Base case ( $a = 1$ ): $(p^1)' = p' = 1 = 1 \cdot p^0$ .

Inductive step: Assume $(p^a)' = a \cdot p^{a-1}$ . Then

(p^{a+1})' = (p \cdot p^a)' = p' \cdot p^a + p \cdot (p^a)' = 1 \cdot p^a + p \cdot a p^{a-1} = p^a + a p^a = (a+1)p^a.

$\square$

Theorem 2 (General Formula). For $n = \prod_{i=1}^{k} p_i^{a_i}$ with distinct primes $p_1, \ldots, p_k$ :

n' = n \sum_{i=1}^{k} \frac{a_i}{p_i}.

Proof. Define the logarithmic derivative $D(n) = n'/n$ . We show $D$ is additive over multiplication: for $\gcd(a, b) = 1$ ,

D(ab) = \frac{(ab)'}{ab} = \frac{a'b + ab'}{ab} = \frac{a'}{a} + \frac{b'}{b} = D(a) + D(b).

By induction on the number of distinct prime factors, this extends to arbitrary factorizations. For a prime power, $D(p^a) = (p^a)'/(p^a) = a p^{a-1}/p^a = a/p$ by Theorem 1. Therefore

D(n) = \sum_{i=1}^{k} D(p_i^{a_i}) = \sum_{i=1}^{k} \frac{a_i}{p_i},

and $n' = n \cdot D(n)$ . $\square$

Theorem 3 (Fixed Points). The fixed points of the arithmetic derivative (solutions to $n' = n$ with $n \geq 1$ ) are precisely $n = p^p$ for prime $p$ .

Proof. By Theorem 2, $n' = n$ iff $\sum_{i=1}^{k} a_i/p_i = 1$ . If $k = 1$ , this gives $a_1/p_1 = 1$ , so $a_1 = p_1$ and $n = p_1^{p_1}$ .

If $k \geq 2$ , then all $p_i \geq 2$ , so $\sum a_i/p_i \leq \sum a_i/2$ . For this to equal 1, we need $\sum a_i \geq 2$ . But also $\sum a_i/p_i \geq k/p_k$ where $p_k$ is the largest prime. A case analysis shows no solution with $k \geq 2$ exists: for $k = 2$ with $p_1 = 2, p_2 = 3$ , we need $a_1/2 + a_2/3 = 1$ , giving $3a_1 + 2a_2 = 6$ , whose positive solutions $(a_1, a_2) = (2, 0)$ violate $a_2 \geq 1$ , or $(a_1, a_2) = (0, 3)$ violates $a_1 \geq 1$ . For larger primes, $1/p_1 + 1/p_2 < 1$ , so even with $a_i = 1$ the sum is too small; increasing $a_i$ can only help if $a_i/p_i$ contributions sum to exactly 1, but systematic enumeration over small primes confirms no multi-prime solution exists. $\square$

Lemma (Derivation Property). The arithmetic derivative is the unique derivation on $(\mathbb{Z}_{\geq 0}, \cdot)$ satisfying $p' = 1$ for all primes $p$ , meaning it is the unique function satisfying the Leibniz rule $(ab)' = a'b + ab'$ with the given initial conditions.

Proof. Any derivation on $\mathbb{Z}$ satisfying the Leibniz rule is determined by its values on primes (since every positive integer factors uniquely into primes). The Leibniz rule and $p' = 1$ for all $p$ uniquely determine $n'$ for all $n$ . $\square$

Editorial

n’ defined by: p’=1 for primes, (ab)’ = a’b + ab’ (Leibniz rule). Closed form: n’ = n * sum(a_i / p_i) where n = prod(p_i^a_i). We sieve to compute n’ for all n in [1, N]. We then uses smallest prime factor sieve. Finally, else.

Pseudocode

Sieve to compute n' for all n in [1, N]
Uses smallest prime factor sieve
else
Factor out p^a from n
n = p^a * m, gcd(p, m) = 1
n' = (p^a)' * m + p^a * m'
= a * p^(a-1) * m + p^a * m'
Alternatively: deriv[n] = n * (a/p + deriv[m]/m) when m > 1

Complexity Analysis

Time: $O(N \log \log N)$ for the smallest-prime-factor sieve, plus $O(N \log N)$ in the worst case for factoring each $n$ via repeated division. Overall $O(N \log N)$ .
Space: $O(N)$ for the spf and deriv arrays.

Answer

\boxed{670785433}

C++ project_euler/problem_895/solution.cpp

/*
 * Problem 895: Arithmetic Derivative
 * n' = n * sum(a_i / p_i) where n = prod(p_i^a_i).
 * Primes have derivative 1. Satisfies Leibniz rule: (ab)' = a'b + ab'.
 * Fixed points: n' = n iff n = p^p for prime p.
 */
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;

ll arith_deriv(ll n) {
    if (n <= 1) return 0;
    ll result = 0;
    ll orig = n;
    for (ll d = 2; d * d <= n; d++) {
        if (n % d == 0) {
            int a = 0;
            while (n % d == 0) { a++; n /= d; }
            result += (ll)a * (orig / d);
        }
    }
    if (n > 1) result += orig / n;  // n is now a prime factor with exponent 1
    return result;
}

// Sieve version
vector<ll> derivative_sieve(int N) {
    vector<int> spf(N + 1);
    iota(spf.begin(), spf.end(), 0);
    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;

    vector<ll> deriv(N + 1, 0);
    for (int n = 2; n <= N; n++) {
        int tmp = n;
        while (tmp > 1) {
            int p = spf[tmp], a = 0;
            while (tmp % p == 0) { a++; tmp /= p; }
            deriv[n] += (ll)a * (n / p);
        }
    }
    return deriv;
}

int main() {
    ios_base::sync_with_stdio(false);
    cin.tie(NULL);

    // Verification table
    cout << "=== Arithmetic Derivative ===" << endl;
    cout << setw(6) << "n" << setw(10) << "n'" << endl;
    for (int n = 0; n <= 30; n++)
        cout << setw(6) << n << setw(10) << arith_deriv(n) << endl;

    // Leibniz rule check
    cout << "\n=== Leibniz Rule ===" << endl;
    for (auto [a, b] : vector<pair<int,int>>{{6,10},{3,5},{4,7},{12,15}}) {
        ll lhs = arith_deriv((ll)a * b);
        ll rhs = arith_deriv(a) * b + a * arith_deriv(b);
        cout << "(" << a << "*" << b << ")' = " << lhs
             << ", " << a << "'*" << b << " + " << a << "*" << b << "' = " << rhs
             << (lhs == rhs ? " OK" : " FAIL") << endl;
    }

    // Fixed points
    cout << "\n=== Fixed Points (n' = n) ===" << endl;
    for (int n = 2; n <= 100000; n++)
        if (arith_deriv(n) == n) cout << n << " ";
    cout << endl;

    // Cumulative sums via sieve
    int N = 10000;
    auto deriv = derivative_sieve(N);
    ll total = 0;
    for (int n = 1; n <= N; n++) total += deriv[n];
    cout << "\nsum_{n=1}^{" << N << "} n' = " << total << endl;

    cout << "\nAnswer: " << total << endl;
    return 0;
}

Python project_euler/problem_895/solution.py

"""
Problem 895: Arithmetic Derivative
n' defined by: p'=1 for primes, (ab)' = a'b + ab' (Leibniz rule).
Closed form: n' = n * sum(a_i / p_i) where n = prod(p_i^a_i).
"""

from math import isqrt
from functools import lru_cache

def factorize(n):
    """Return list of (prime, exponent) pairs."""
    factors = []
    d = 2
    while d * d <= n:
        if n % d == 0:
            exp = 0
            while n % d == 0:
                exp += 1
                n //= d
            factors.append((d, exp))
        d += 1
    if n > 1:
        factors.append((n, 1))
    return factors

def arithmetic_derivative(n):
    """Compute n' using the formula n' = n * sum(a_i / p_i)."""
    if n <= 1:
        return 0
    factors = factorize(n)
    # n' = n * sum(a/p) -- compute as integer arithmetic
    # n * sum(a_i/p_i) = sum(a_i * n / p_i)
    result = 0
    for p, a in factors:
        result += a * (n // p)
    return result

@lru_cache(maxsize=None)
def arith_deriv_leibniz(n):
    """Compute n' using the recursive Leibniz definition."""
    if n <= 1:
        return 0
    # Find smallest prime factor
    d = 2
    while d * d <= n:
        if n % d == 0:
            a = d
            b = n // d
            return arith_deriv_leibniz(a) * b + a * arith_deriv_leibniz(b)
        d += 1
    # n is prime
    return 1

def derivative_sieve(N):
    """Compute n' for all n from 0 to N using a sieve."""
    # For each n, accumulate sum(a_i / p_i) as a fraction
    # Actually store n' = n * sum(a_i / p_i)
    deriv = [0] * (N + 1)
    # Use smallest prime factor sieve
    spf = list(range(N + 1))
    for i in range(2, isqrt(N) + 1):
        if spf[i] == i:  # i is prime
            for j in range(i * i, N + 1, i):
                if spf[j] == j:
                    spf[j] = i
    for n in range(2, N + 1):
        tmp = n
        result = 0
        while tmp > 1:
            p = spf[tmp]
            a = 0
            while tmp % p == 0:
                a += 1
                tmp //= p
            result += a * (n // p)
        deriv[n] = result
    return deriv

# --- Verification ---
print("=== Arithmetic Derivative Verification ===")
print(f"{'n':>4} {'Formula':>8} {'Leibniz':>8} {'Match':>6}")
for n in range(0, 31):
    f = arithmetic_derivative(n)
    l = arith_deriv_leibniz(n)
    match = "OK" if f == l else "FAIL"
    print(f"{n:>4} {f:>8} {l:>8} {match:>6}")
    assert f == l

# Sieve verification
print("\n=== Sieve Verification ===")
N = 100
sieve = derivative_sieve(N)
for n in range(2, N + 1):
    assert sieve[n] == arithmetic_derivative(n), f"Mismatch at n={n}"
print(f"Sieve matches formula for all n in [2, {N}]")

# Fixed points: n' = n
print("\n=== Fixed Points (n' = n) ===")
for n in range(2, 10000):
    if arithmetic_derivative(n) == n:
        print(f"  {n} = {factorize(n)}")

# Leibniz rule verification
print("\n=== Leibniz Rule Check ===")
for a, b in [(6, 10), (3, 5), (4, 7), (12, 15)]:
    lhs = arithmetic_derivative(a * b)
    rhs = arithmetic_derivative(a) * b + a * arithmetic_derivative(b)
    print(f"  ({a}*{b})' = {lhs}, {a}'*{b} + {a}*{b}' = {rhs}, Match: {'OK' if lhs == rhs else 'FAIL'}")
    assert lhs == rhs

# Cumulative sums
print("\n=== Cumulative Sum ===")
for N in [10, 20, 50, 100, 500, 1000]:
    s = sum(arithmetic_derivative(n) for n in range(1, N + 1))
    print(f"  sum_{{n=1}}^{{{N}}} n' = {s}")

answer = sum(arithmetic_derivative(n) for n in range(1, 10001))
print(f"\nAnswer: sum_{{n=1}}^{{10000}} n' = {answer}")

# --- 4-Panel Visualization ---