Problem 672: One More One

Project Euler Problem 672 Solution

EulerSolve provides an optimized solution for Project Euler Problem 672, One More One, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Define $g(n)$ as follows: starting from $n$, repeatedly divide by $7$ when possible; otherwise add $1$, and count only those $+1$ operations. The problem asks for the large cumulative quantity $$F(N)=\sum_{m=1}^{N} g(m),\qquad H(k)=F\left(\frac{7^k-1}{11}\right)\pmod{M},\qquad M=1{,}117{,}117{,}717.$$ The required value is $H(10^9)$. A direct simulation is impossible because the upper limit $\frac{7^k-1}{11}$ is astronomically large, so the solution rewrites the process in base $7$, turns it into a constant-size linear update, and then exploits the repeating digit pattern of $\frac{1}{11}$ in base $7$. Mathematical Approach Step 1: Rewrite the primitive rule If $n=7a+r$ with $0\le r\le 6$, then one application of the rule gives $$g(7a)=g(a),\qquad g(7a+r)=7-r+g(a+1)\quad (1\le r\le 6).$$ The second formula says that when the last base-$7$ digit is $r\neq 0$, we need exactly $7-r$ increments to reach the next multiple of $7$, after which the number becomes $a+1$....

Detailed mathematical approach

Problem Summary

Define $g(n)$ as follows: starting from $n$, repeatedly divide by $7$ when possible; otherwise add $1$, and count only those $+1$ operations. The problem asks for the large cumulative quantity

$$F(N)=\sum_{m=1}^{N} g(m),\qquad H(k)=F\left(\frac{7^k-1}{11}\right)\pmod{M},\qquad M=1{,}117{,}117{,}717.$$

The required value is $H(10^9)$. A direct simulation is impossible because the upper limit $\frac{7^k-1}{11}$ is astronomically large, so the solution rewrites the process in base $7$, turns it into a constant-size linear update, and then exploits the repeating digit pattern of $\frac{1}{11}$ in base $7$.

Mathematical Approach

Step 1: Rewrite the primitive rule

If $n=7a+r$ with $0\le r\le 6$, then one application of the rule gives

$$g(7a)=g(a),\qquad g(7a+r)=7-r+g(a+1)\quad (1\le r\le 6).$$

The second formula says that when the last base-$7$ digit is $r\neq 0$, we need exactly $7-r$ increments to reach the next multiple of $7$, after which the number becomes $a+1$.

Step 2: Shift by one and obtain a digit formula

Introduce

$$f(n)=g(n+1).$$

For every positive $n=7a+d$ with $0\le d\le 6$, the previous recurrence becomes

$$f(7a+d)=f(a)+6-d,\qquad f(0)=0.$$

Repeatedly stripping the last base-$7$ digit shows that if $n=(d_m d_{m-1}\cdots d_0)_7$ is positive, then

$$f(n)=\sum_{i=0}^{m}(6-d_i).$$

So $g(n+1)$ is just a weighted digit sum. For example, $124=(235)_7$, hence

$$g(125)=f(124)=(6-2)+(6-3)+(6-5)=4+3+1=8,$$

which matches the small checkpoint used by the implementations.

Step 3: Derive a recurrence for the cumulative sum

Now define the prefix sum

$$F(N)=\sum_{m=1}^{N} g(m)=\sum_{n=0}^{N-1} f(n).$$

For a positive prefix $n$ and a next base-$7$ digit $d$, split the range $0\le x<7n+d$ into full blocks $7a,7a+1,\dots,7a+6$ for $0\le a<n$ and one partial block for $a=n$. This yields

$$F(7n+d)=7F(n)+21n-6+d\,f(n)+P(d),$$

where

$$P(d)=\sum_{x=0}^{d-1}(6-x)=\frac{d(13-d)}{2}.$$

The term $-6$ is the only correction needed for the special value $f(0)=0$; otherwise the full-block contribution would be completely uniform.

Step 4: Encode the update as a $5\times 5$ matrix

For positive $n$, use the augmented state

$$\mathbf{v}(n)=\begin{pmatrix}1\\ n-1\\ F(n)\\ f(n)\\ 1\end{pmatrix}.$$

Then appending one digit $d\in\{0,\dots,6\}$ corresponds to $n\mapsto 7n+d$, and the state evolves linearly:

$$\mathbf{v}(7n+d)=M_d\,\mathbf{v}(n),$$

with

$$M_d=\begin{pmatrix} 1 & 0 & 0 & 0 & 0\\ 6 & 7 & 0 & 0 & d\\ 15 & 21 & 7 & d & P(d)\\ 0 & 0 & 0 & 1 & 6-d\\ 0 & 0 & 0 & 0 & 1 \end{pmatrix} \pmod{M}.$$

This is the exact matrix family used by the C++, Python, and Java implementations.

Step 5: Exploit the special upper limit $\frac{7^k-1}{11}$

Because $11\mid 7^{10}-1$, we have a purely periodic base-$7$ expansion for $\frac{1}{11}$, and in particular

$$\frac{7^{10}-1}{11}=(431162355)_7.$$

For $k=10t$,

$$\frac{7^k-1}{11}=\frac{7^{10}-1}{11}\left(1+7^{10}+7^{20}+\cdots+7^{10(t-1)}\right).$$

Therefore its base-$7$ digits are the block 4311623550 repeated $t-1$ times, followed by 431162355. Since every such number begins with the digit $4$, the implementations start directly from the state

$$\mathbf{v}(4)=\begin{pmatrix}1\\ 3\\ 12\\ 2\\ 1\end{pmatrix},$$

which already stores $4-1$, $F(4)=12$, and $f(4)=2$.

Worked Example

A small local check comes from appending the digit $3$ to the prefix $4$, giving $43_7=31$. Using $F(4)=12$, $f(4)=2$, and $P(3)=15$,

$$F(31)=7\cdot 12+21\cdot 4-6+3\cdot 2+15=183.$$

So the recurrence already reproduces $\sum_{m=1}^{31} g(m)=183$. For the first nontrivial target, $k=10$, the upper limit is $(431162355)_7$. Starting from $\mathbf{v}(4)$ and processing the remaining suffix 31162355 gives

$$H(10)=690409338 \pmod{M},$$

exactly the checkpoint embedded in the implementations.

How the Code Works

The implementations precompute the seven digit-dependent matrices $M_0,\dots,M_6$, always working modulo $M$. They also precompute the transformations of the few fixed digit strings that appear in the base-$7$ representation of $\frac{7^k-1}{11}$: one short suffix for $k=10$, one initial boundary block, one repeating 10-digit block, and one final block.

For $k=10$, the short suffix is applied once to the initial state $\mathbf{v}(4)$. For $k=10t\ge 20$, the implementation applies the initial boundary block, raises the repeating block transformation to the power $t-2$ by binary exponentiation, applies the final block, and then reads the third state component. That component is $F\!\left(\frac{7^k-1}{11}\right)$, so it is exactly $H(k)$.

Complexity Analysis

Each matrix has fixed dimension $5$, so one multiplication is constant-time in the asymptotic sense. The only nontrivial growth comes from raising the repeating block transformation to the power $t-2$, where $t=k/10$. Binary exponentiation therefore gives $O(\log t)=O(\log k)$ matrix multiplications and $O(1)$ extra memory.

Footnotes and References

Project Euler problem page: https://projecteuler.net/problem=672
Base-$7$ numerals and positional notation: Wikipedia — Radix
Geometric series: Wikipedia — Geometric series
Matrix exponentiation: Wikipedia — Exponentiation by squaring

Problem 672 source code

C++

#include <array>
#include <cassert>
#include <cstdint>
#include <iostream>
#include <string>

namespace {

using u64 = std::uint64_t;
using u128 = __uint128_t;

constexpr u64 MOD = 1'117'117'717ULL;
constexpr int DIM = 5;

struct Mat {
    std::array<std::array<u64, DIM>, DIM> a{};

    static Mat identity() {
        Mat m;
        for (int i = 0; i < DIM; ++i) {
            m.a[static_cast<std::size_t>(i)][static_cast<std::size_t>(i)] = 1ULL;
        }
        return m;
    }
};

Mat multiply(const Mat& x, const Mat& y) {
    Mat z;
    for (int i = 0; i < DIM; ++i) {
        for (int k = 0; k < DIM; ++k) {
            const u64 xik = x.a[static_cast<std::size_t>(i)][static_cast<std::size_t>(k)];
            if (xik == 0ULL) {
                continue;
            }
            for (int j = 0; j < DIM; ++j) {
                const u64 ykj = y.a[static_cast<std::size_t>(k)][static_cast<std::size_t>(j)];
                if (ykj == 0ULL) {
                    continue;
                }
                u64& cell = z.a[static_cast<std::size_t>(i)][static_cast<std::size_t>(j)];
                cell = static_cast<u64>((static_cast<u128>(cell) +
                                         static_cast<u128>(xik) * static_cast<u128>(ykj)) %
                                        MOD);
            }
        }
    }
    return z;
}

Mat power(Mat base, u64 exp) {
    Mat result = Mat::identity();
    while (exp > 0ULL) {
        if ((exp & 1ULL) != 0ULL) {
            result = multiply(base, result);
        }
        exp >>= 1ULL;
        if (exp > 0ULL) {
            base = multiply(base, base);
        }
    }
    return result;
}

std::array<u64, DIM> apply_matrix(const Mat& m, const std::array<u64, DIM>& v) {
    std::array<u64, DIM> out{};
    for (int i = 0; i < DIM; ++i) {
        u64 s = 0ULL;
        for (int j = 0; j < DIM; ++j) {
            s = static_cast<u64>((static_cast<u128>(s) +
                                  static_cast<u128>(m.a[static_cast<std::size_t>(i)]
                                                        [static_cast<std::size_t>(j)]) *
                                      static_cast<u128>(v[static_cast<std::size_t>(j)])) %
                                 MOD);
        }
        out[static_cast<std::size_t>(i)] = s;
    }
    return out;
}

u64 slow_g(u64 n) {
    u64 added = 0ULL;
    while (n > 1ULL) {
        if (n % 7ULL == 0ULL) {
            n /= 7ULL;
        } else {
            ++n;
            ++added;
        }
    }
    return added;
}

Mat digit_matrix(int d, const std::array<int, 7>& weight, const std::array<int, 7>& pref) {
    Mat m;

    m.a[0][0] = 1ULL;

    m.a[1][0] = 6ULL;
    m.a[1][1] = 7ULL;
    m.a[1][4] = static_cast<u64>(d);

    m.a[2][0] = 15ULL;
    m.a[2][1] = 21ULL;
    m.a[2][2] = 7ULL;
    m.a[2][3] = static_cast<u64>(d);
    m.a[2][4] = static_cast<u64>(pref[static_cast<std::size_t>(d)]);

    m.a[3][3] = 1ULL;
    m.a[3][4] = static_cast<u64>(weight[static_cast<std::size_t>(d)]);

    m.a[4][4] = 1ULL;

    return m;
}

Mat sequence_matrix(const std::string& seq, const std::array<int, 7>& weight,
                    const std::array<int, 7>& pref) {
    Mat all = Mat::identity();
    for (char ch : seq) {
        const int d = ch - '0';
        const Mat md = digit_matrix(d, weight, pref);
        all = multiply(md, all);
    }
    return all;
}

u64 solve_H(u64 k) {
    assert(k % 10ULL == 0ULL);
    assert(k >= 10ULL);

    const std::array<int, 7> weight{6, 5, 4, 3, 2, 1, 0};
    std::array<int, 7> pref{};
    int run = 0;
    for (int d = 0; d < 7; ++d) {
        pref[static_cast<std::size_t>(d)] = run;
        run += weight[static_cast<std::size_t>(d)];
    }

    const std::string blockA = "4311623550";
    const std::string blockB = "431162355";
    const std::string remSingle = "31162355";
    const std::string remFirst = "311623550";

    const Mat MA = sequence_matrix(blockA, weight, pref);
    const Mat MB = sequence_matrix(blockB, weight, pref);
    const Mat MSingle = sequence_matrix(remSingle, weight, pref);
    const Mat MFirst = sequence_matrix(remFirst, weight, pref);

    const u64 t = k / 10ULL;

    std::array<u64, DIM> v{};
    v[0] = 1ULL;
    v[1] = 3ULL;
    v[2] = 12ULL;
    v[3] = 2ULL;
    v[4] = 1ULL;

    if (t == 1ULL) {
        v = apply_matrix(MSingle, v);
    } else {
        v = apply_matrix(MFirst, v);
        if (t > 2ULL) {
            const Mat mid = power(MA, t - 2ULL);
            v = apply_matrix(mid, v);
        }
        v = apply_matrix(MB, v);
    }

    return v[2] % MOD;
}

}  // namespace

int main() {
    assert(slow_g(125ULL) == 8ULL);
    assert(slow_g(1'000ULL) == 9ULL);
    assert(slow_g(10'000ULL) == 21ULL);

    assert(solve_H(10ULL) == 690'409'338ULL);

    std::cout << solve_H(1'000'000'000ULL) << "\n";
    return 0;
}

Python

MOD = 1117117717
DIM = 5

def identity():
    return [[1 if i == j else 0 for j in range(DIM)] for i in range(DIM)]

def multiply(x, y):
    z = [[0] * DIM for _ in range(DIM)]
    for i in range(DIM):
        for k in range(DIM):
            xik = x[i][k]
            if not xik:
                continue
            for j in range(DIM):
                ykj = y[k][j]
                if ykj:
                    z[i][j] = (z[i][j] + xik * ykj) % MOD
    return z

def power(base, exp):
    result = identity()
    while exp > 0:
        if exp & 1:
            result = multiply(base, result)
        exp >>= 1
        if exp > 0:
            base = multiply(base, base)
    return result

def apply_matrix(m, v):
    out = [0] * DIM
    for i in range(DIM):
        s = 0
        for j in range(DIM):
            s = (s + m[i][j] * v[j]) % MOD
        out[i] = s
    return out

def digit_matrix(d, weight, pref):
    m = [[0] * DIM for _ in range(DIM)]
    m[0][0] = 1

    m[1][0] = 6
    m[1][1] = 7
    m[1][4] = d

    m[2][0] = 15
    m[2][1] = 21
    m[2][2] = 7
    m[2][3] = d
    m[2][4] = pref[d]

    m[3][3] = 1
    m[3][4] = weight[d]

    m[4][4] = 1

    return m

def sequence_matrix(seq, weight, pref):
    all_mat = identity()
    for ch in seq:
        d = int(ch)
        md = digit_matrix(d, weight, pref)
        all_mat = multiply(md, all_mat)
    return all_mat

def solve_H(k):
    weight = [6, 5, 4, 3, 2, 1, 0]
    pref = [0] * 7
    run = 0
    for d in range(7):
        pref[d] = run
        run += weight[d]

    blockA = "4311623550"
    blockB = "431162355"
    remSingle = "31162355"
    remFirst = "311623550"

    MA = sequence_matrix(blockA, weight, pref)
    MB = sequence_matrix(blockB, weight, pref)
    MSingle = sequence_matrix(remSingle, weight, pref)
    MFirst = sequence_matrix(remFirst, weight, pref)

    t = k // 10

    v = [1, 3, 12, 2, 1]

    if t == 1:
        v = apply_matrix(MSingle, v)
    else:
        v = apply_matrix(MFirst, v)
        if t > 2:
            mid = power(MA, t - 2)
            v = apply_matrix(mid, v)
        v = apply_matrix(MB, v)

    return v[2] % MOD

def solve():
    ans = solve_H(1000000000)
    return str(ans)

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

Java

public class Euler672 {

    static final long MOD = 1117117717L;
    static final int DIM = 5;

    static long[][] identity() {
        long[][] m = new long[DIM][DIM];
        for (int i = 0; i < DIM; ++i) {
            m[i][i] = 1L;
        }
        return m;
    }

    static long[][] multiply(long[][] x, long[][] y) {
        long[][] z = new long[DIM][DIM];
        for (int i = 0; i < DIM; ++i) {
            for (int k = 0; k < DIM; ++k) {
                long xik = x[i][k];
                if (xik == 0L)
                    continue;
                for (int j = 0; j < DIM; ++j) {
                    long ykj = y[k][j];
                    if (ykj == 0L)
                        continue;
                    z[i][j] = (z[i][j] + xik * ykj) % MOD;
                }
            }
        }
        return z;
    }

    static long[][] power(long[][] base, long exp) {
        long[][] result = identity();
        while (exp > 0L) {
            if ((exp & 1L) != 0L) {
                result = multiply(base, result);
            }
            exp >>= 1L;
            if (exp > 0L) {
                base = multiply(base, base);
            }
        }
        return result;
    }

    static long[] applyMatrix(long[][] m, long[] v) {
        long[] out = new long[DIM];
        for (int i = 0; i < DIM; ++i) {
            long s = 0L;
            for (int j = 0; j < DIM; ++j) {
                s = (s + m[i][j] * v[j]) % MOD;
            }
            out[i] = s;
        }
        return out;
    }

    static long[][] digitMatrix(int d, int[] weight, int[] pref) {
        long[][] m = new long[DIM][DIM];
        m[0][0] = 1L;
        m[1][0] = 6L;
        m[1][1] = 7L;
        m[1][4] = d;
        m[2][0] = 15L;
        m[2][1] = 21L;
        m[2][2] = 7L;
        m[2][3] = d;
        m[2][4] = pref[d];
        m[3][3] = 1L;
        m[3][4] = weight[d];
        m[4][4] = 1L;
        return m;
    }

    static long[][] sequenceMatrix(String seq, int[] weight, int[] pref) {
        long[][] all = identity();
        for (int i = 0; i < seq.length(); ++i) {
            int d = seq.charAt(i) - '0';
            long[][] md = digitMatrix(d, weight, pref);
            all = multiply(md, all);
        }
        return all;
    }

    static long solveH(long k) {
        int[] weight = { 6, 5, 4, 3, 2, 1, 0 };
        int[] pref = new int[7];
        int run = 0;
        for (int d = 0; d < 7; ++d) {
            pref[d] = run;
            run += weight[d];
        }

        String blockA = "4311623550";
        String blockB = "431162355";
        String remSingle = "31162355";
        String remFirst = "311623550";

        long[][] MA = sequenceMatrix(blockA, weight, pref);
        long[][] MB = sequenceMatrix(blockB, weight, pref);
        long[][] MSingle = sequenceMatrix(remSingle, weight, pref);
        long[][] MFirst = sequenceMatrix(remFirst, weight, pref);

        long t = k / 10L;

        long[] v = { 1L, 3L, 12L, 2L, 1L };

        if (t == 1L) {
            v = applyMatrix(MSingle, v);
        } else {
            v = applyMatrix(MFirst, v);
            if (t > 2L) {
                long[][] mid = power(MA, t - 2L);
                v = applyMatrix(mid, v);
            }
            v = applyMatrix(MB, v);
        }

        return v[2] % MOD;
    }

    public static String solve() {
        long ans = solveH(1000000000L);
        return Long.toString(ans);
    }

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