Problem 603: Substring Sums of Prime Concatenations
View on Project EulerProject Euler Problem 603 Solution
EulerSolve provides an optimized solution for Project Euler Problem 603, Substring Sums of Prime Concatenations, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary For a decimal string \(w\), let \(S(w)\) denote the sum of all contiguous decimal substrings of \(w\), counted with multiplicity and interpreted as ordinary integers. Let \(P(t)\) be the decimal concatenation of the first \(t\) primes, and let \(C(t,k)\) be \(k\) consecutive copies of \(P(t)\). The target is $$S\big(C(10^6,10^{12})\big)\bmod(10^9+7).$$ The direct string is far too large to build explicitly: even one copy of \(P(10^6)\) already contains millions of digits, and the problem asks for \(10^{12}\) copies. The solution therefore compresses the effect of appending digits into a fixed linear transformation and then raises one block transformation to a huge power. Mathematical Approach Write any decimal string as \(d_1d_2\cdots d_L\), where each \(d_j\in\{0,\dots,9\}\). We first derive an exact recurrence for substring sums, then package one digit-append step into a \(4\times4\) matrix, and finally exploit the fact that the large input is just one repeated prime-concatenation block. Step 1: Sum All Substrings Ending at One Position Let \(F_j\) be the sum of all substrings that end exactly at position \(j\). Then the total answer for the whole string is $$S(d_1d_2\cdots d_L)=\sum_{j=1}^{L} F_j.$$ Why is there a simple recurrence? Every substring ending at position \(j-1\) can be extended by digit \(d_j\), which multiplies its previous value by \(10\)....
Detailed mathematical approach
Problem Summary
For a decimal string \(w\), let \(S(w)\) denote the sum of all contiguous decimal substrings of \(w\), counted with multiplicity and interpreted as ordinary integers. Let \(P(t)\) be the decimal concatenation of the first \(t\) primes, and let \(C(t,k)\) be \(k\) consecutive copies of \(P(t)\). The target is
$$S\big(C(10^6,10^{12})\big)\bmod(10^9+7).$$
The direct string is far too large to build explicitly: even one copy of \(P(10^6)\) already contains millions of digits, and the problem asks for \(10^{12}\) copies. The solution therefore compresses the effect of appending digits into a fixed linear transformation and then raises one block transformation to a huge power.
Mathematical Approach
Write any decimal string as \(d_1d_2\cdots d_L\), where each \(d_j\in\{0,\dots,9\}\). We first derive an exact recurrence for substring sums, then package one digit-append step into a \(4\times4\) matrix, and finally exploit the fact that the large input is just one repeated prime-concatenation block.
Step 1: Sum All Substrings Ending at One Position
Let \(F_j\) be the sum of all substrings that end exactly at position \(j\). Then the total answer for the whole string is
$$S(d_1d_2\cdots d_L)=\sum_{j=1}^{L} F_j.$$
Why is there a simple recurrence? Every substring ending at position \(j-1\) can be extended by digit \(d_j\), which multiplies its previous value by \(10\). In addition, the new digit \(d_j\) appears as the last digit in each of the \(j\) substrings ending at \(j\): the one-digit substring \(d_j\), the two-digit substring ending at \(j\), and so on up to the full prefix \(d_1\cdots d_j\). Therefore
$$F_j=10F_{j-1}+j\,d_j,\qquad F_0=0.$$
This recurrence is exact and already reduces the problem from quadratic substring enumeration to a single left-to-right scan.
Step 2: Turn the Recurrence into a Linear State
To compose many digit steps efficiently, track four quantities after processing \(\ell\) digits:
$$v_\ell=\begin{bmatrix}F_\ell\\ S_\ell\\ \ell\\ 1\end{bmatrix},\qquad S_\ell=\sum_{j=1}^{\ell} F_j.$$
If the next digit is \(d\), then the new position is \(\ell+1\), so the recurrence becomes
$$F_{\ell+1}=10F_\ell+(\ell+1)d,$$
$$S_{\ell+1}=S_\ell+F_{\ell+1},$$
$$\ell'=\ell+1.$$
That update is linear in the state vector:
$$v_{\ell+1}=T(d)\,v_\ell,$$
with
$$T(d)=\begin{bmatrix} 10 & 0 & d & d\\ 10 & 1 & d & d\\ 0 & 0 & 1 & 1\\ 0 & 0 & 0 & 1 \end{bmatrix}.$$
So appending one digit is no longer a special-case arithmetic step; it is just multiplication by a fixed-size matrix.
Step 3: A Whole String Becomes One Matrix
For a block \(w=d_1d_2\cdots d_L\), define its transform \(M(w)\) by feeding the digits from left to right. Starting from
$$v_0=\begin{bmatrix}0\\0\\0\\1\end{bmatrix},$$
we get
$$v_L=M(w)\,v_0.$$
Because matrix multiplication composes successive digit steps, concatenation of strings corresponds to multiplication of their block transforms. If \(u\) is processed first and then \(v\), then
$$M(uv)=M(v)\,M(u).$$
This is the key structural fact: instead of thinking about billions or trillions of individual digits, we can think about a much smaller number of block transformations.
Step 4: Repeating One Prime Block Means Matrix Powers
The problem input is not an arbitrary string; it is the same block repeated many times. Let
$$B=M\big(P(10^6)\big).$$
Then one copy of \(P(10^6)\) acts as \(B\), two copies act as \(B^2\), and in general
$$M\big(C(10^6,k)\big)=B^k.$$
Therefore the required state is
$$v_{\text{final}}=B^{10^{12}}v_0,$$
and the answer to the problem is the second component of \(v_{\text{final}}\), namely the accumulated substring sum.
This replaces an astronomically long concatenation with binary exponentiation on a \(4\times4\) matrix.
Step 5: Build the Prime Block Without Building the Giant String
To obtain \(B\), we only need the decimal digits of the first \(10^6\) primes in order. A sieve gives those primes up to a safe upper bound for the millionth prime, and each prime is then decomposed into decimal digits from most significant to least significant. Every digit updates the current \(4\times4\) block transform once.
No copy of \(P(10^6)\) is stored as one long string, and certainly no copy of \(C(10^6,10^{12})\) is stored. The entire method keeps only a tiny state matrix, a few fixed-size vectors, and the sieve structure needed to enumerate primes.
Worked Example
Take the ordinary string \(2024\). The recurrence gives
$$F_1=2,$$
$$F_2=10\cdot 2+2\cdot 0=20,$$
$$F_3=10\cdot 20+3\cdot 2=206,$$
$$F_4=10\cdot 206+4\cdot 4=2076.$$
Hence
$$S(2024)=2+20+206+2076=2304.$$
This matches the direct substring sum \(2+0+2+4+20+2+24+202+24+2024\). The same recurrence is what the full solution uses; the only difference is that the actual problem feeds in the digits of prime concatenations and then compresses repeated blocks with matrix powers.
How the Code Works
The C++, Python, and Java implementations all follow the same structure. First, they define modular arithmetic under \(10^9+7\) and implement \(4\times4\) matrix multiplication together with matrix-vector multiplication for the state described above.
Next, they generate the first million primes with a sieve up to a safe analytic upper bound. The digits of each prime are processed from most significant to least significant, and each digit left-composes the current block transform by the corresponding digit matrix. After the millionth prime has been processed, the code has the single matrix \(B\) for one block \(P(10^6)\).
Finally, they apply exponentiation by squaring to compute the action of \(B^{10^{12}}\) on the initial state vector. Because the matrix dimension is fixed at \(4\), this stage is tiny compared with prime generation and digit streaming. The printed output is the second state component after the exponentiation loop finishes.
Complexity Analysis
Let \(n=10^6\), let \(U\) be a safe upper bound for the \(n\)-th prime, and let \(D\) be the number of decimal digits in \(P(n)\). The sieve costs \(O(U\log\log U)\) time and \(O(U)\) memory, or \(O(U/2)\) logical storage in an odd-only representation. Feeding the digits of all primes into the block matrix costs \(O(D)\) time. Exponentiating a \(4\times4\) matrix to the power \(10^{12}\) costs \(O(\log 10^{12})\) matrix multiplications and only \(O(1)\) extra memory beyond the fixed matrices and vectors.
So the overall complexity is
$$O(U\log\log U + D + \log 10^{12})$$
time and \(O(U)\) memory. In practice, the dominant work is prime generation plus the single pass over the digits of the first million primes.
Footnotes and References
- Problem page: Project Euler 603
- Exponentiation by squaring: Wikipedia - Exponentiation by squaring
- Sieve of Eratosthenes: Wikipedia - Sieve of Eratosthenes
- Prime number theorem: Wikipedia - Prime number theorem
- Matrix multiplication: Wikipedia - Matrix multiplication
Problem 603 source code
C++
#include <array>
#include <cstdint>
#include <cmath>
#include <iostream>
#include <string>
#include <vector>
// Project Euler 603: Substring Sums
// Use the recurrence for substrings ending at j:
// F_j = 10*F_{j-1} + d_j * j, and S = sum_j F_j.
// Model the update on (F,S,j,1) with a 4x4 matrix; build the matrix for P(10^6) and
// exponentiate it to apply 10^12 repetitions modulo 1e9+7.
using u64 = std::uint64_t;
static constexpr int MOD = 1000000007;
static inline int addmod(int a, int b) {
int s = a + b;
if (s >= MOD) s -= MOD;
return s;
}
static inline int submod(int a, int b) {
int s = a - b;
if (s < 0) s += MOD;
return s;
}
static inline int mulmod(long long a, long long b) {
return (int)((a * b) % MOD);
}
struct Mat {
int a[4][4];
};
static Mat mat_identity() {
Mat m{};
for (int i = 0; i < 4; ++i) m.a[i][i] = 1;
return m;
}
static Mat mat_mul(const Mat& x, const Mat& y) {
Mat z{};
for (int i = 0; i < 4; ++i) {
for (int k = 0; k < 4; ++k) {
const int xik = x.a[i][k];
if (!xik) continue;
for (int j = 0; j < 4; ++j) {
z.a[i][j] = addmod(z.a[i][j], mulmod(xik, y.a[k][j]));
}
}
}
return z;
}
static std::array<int, 4> mat_vec_mul(const Mat& m, const std::array<int, 4>& v) {
std::array<int, 4> r{0, 0, 0, 0};
for (int i = 0; i < 4; ++i) {
long long acc = 0;
for (int k = 0; k < 4; ++k) acc += (long long)m.a[i][k] * v[k];
r[i] = (int)(acc % MOD);
}
return r;
}
static void feed_digit(Mat& m, int d) {
// Left-multiply the current matrix by the per-digit transform.
// Update columns using the state update rules.
for (int col = 0; col < 4; ++col) {
const int f = m.a[0][col];
const int s = m.a[1][col];
const int t = m.a[2][col];
const int one = m.a[3][col];
const int dt = mulmod(d, t);
const int done = mulmod(d, one);
const int tenf = mulmod(10, f);
const int f2 = addmod(addmod(tenf, dt), done);
const int s2 = addmod(s, f2);
const int t2 = addmod(t, one);
m.a[0][col] = f2;
m.a[1][col] = s2;
m.a[2][col] = t2;
m.a[3][col] = one;
}
}
static void feed_number(Mat& m, int x) {
// Feed decimal digits of x in most-significant-first order, no allocations.
int buf[16];
int len = 0;
while (x > 0) {
buf[len++] = x % 10;
x /= 10;
}
for (int i = len - 1; i >= 0; --i) feed_digit(m, buf[i]);
}
static Mat matrix_from_string(const std::string& s) {
Mat m = mat_identity();
for (char ch : s) feed_digit(m, ch - '0');
return m;
}
static int substring_sum_mod(const std::string& s) {
long long f = 0;
long long ans = 0;
for (int i = 0; i < (int)s.size(); ++i) {
const int d = s[i] - '0';
const long long j = i + 1; // 1-indexed
f = (10 * f + (long long)d * j) % MOD;
ans += f;
if (ans >= MOD) ans %= MOD;
}
return (int)(ans % MOD);
}
static int upper_bound_nth_prime(int n) {
// For n >= 6: p_n < n (log n + log log n) (Rosser-Schoenfeld).
if (n < 6) return 15;
const double nd = (double)n;
const double bound = nd * (std::log(nd) + std::log(std::log(nd)));
return (int)std::ceil(bound + 100.0); // generous margin for rounding
}
static Mat matrix_for_P_primes(int n) {
// Build the per-block matrix for P(n): concatenation of the first n primes.
int limit = upper_bound_nth_prime(n);
for (;;) {
std::vector<std::uint8_t> comp((limit >> 1) + 1, 0); // only odds; idx = x>>1
const int r = (int)std::sqrt((double)limit);
for (int p = 3; p <= r; p += 2) {
if (comp[p >> 1]) continue;
const int step = p << 1;
for (int x = p * p; x <= limit; x += step) comp[x >> 1] = 1;
}
Mat m = mat_identity();
int count = 0;
// prime 2
++count;
feed_number(m, 2);
if (count == n) return m;
for (int x = 3; x <= limit; x += 2) {
if (comp[x >> 1]) continue;
++count;
feed_number(m, x);
if (count == n) return m;
}
// Not enough primes; increase the limit and retry (shouldn't happen with the bound).
limit = (int)(limit * 1.2) + 1000;
}
}
int main() {
// Validations from the statement / definition.
if (substring_sum_mod("2024") != 2304) {
std::cerr << "Validation failed: S(2024)\n";
return 1;
}
// Cross-check matrix method on a small case: P(7) repeated 3 times.
const std::string p7 = "2357111317";
const std::string c73 = p7 + p7 + p7;
const int naive = substring_sum_mod(c73);
Mat m7 = matrix_from_string(p7);
std::array<int, 4> v{0, 0, 0, 1};
u64 k = 3;
while (k) {
if (k & 1) v = mat_vec_mul(m7, v);
m7 = mat_mul(m7, m7);
k >>= 1;
}
if (v[1] != naive) {
std::cerr << "Validation failed: matrix repetition\n";
return 1;
}
const int N = 1000000;
const u64 K = 1000000000000ULL;
Mat base = matrix_for_P_primes(N);
std::array<int, 4> state{0, 0, 0, 1};
u64 e = K;
while (e) {
if (e & 1) state = mat_vec_mul(base, state);
base = mat_mul(base, base);
e >>= 1;
}
std::cout << state[1] << "\n";
return 0;
}
Python
import math
def solve():
MOD = 1000000007
def addmod(a, b):
s = a + b
return s - MOD if s >= MOD else s
def mulmod(a, b): return a * b % MOD
# 4x4 matrix operations
def mat_id():
m = [[0]*4 for _ in range(4)]
for i in range(4): m[i][i] = 1
return m
def mat_mul(x, y):
z = [[0]*4 for _ in range(4)]
for i in range(4):
for k in range(4):
if x[i][k] == 0: continue
for j in range(4):
z[i][j] = addmod(z[i][j], mulmod(x[i][k], y[k][j]))
return z
def mat_vec(m, v):
r = [0]*4
for i in range(4):
acc = 0
for k in range(4): acc += m[i][k] * v[k]
r[i] = acc % MOD
return r
def feed_digit(m, d):
for col in range(4):
f, s, t, one = m[0][col], m[1][col], m[2][col], m[3][col]
f2 = addmod(addmod(mulmod(10, f), mulmod(d, t)), mulmod(d, one))
s2 = addmod(s, f2); t2 = addmod(t, one)
m[0][col] = f2; m[1][col] = s2; m[2][col] = t2
def feed_number(m, x):
digits = []
while x > 0: digits.append(x % 10); x //= 10
for d in reversed(digits): feed_digit(m, d)
# Sieve primes up to bound for first N primes
N = 1000000; K = 1000000000000
bound = int(N * (math.log(N) + math.log(math.log(N)))) + 200
sieve = bytearray(b'\x01') * ((bound >> 1) + 1)
r = int(bound**0.5)
for p in range(3, r+1, 2):
if sieve[p >> 1]:
for x in range(p*p, bound+1, 2*p): sieve[x >> 1] = 0
# Build per-block matrix
base = mat_id(); count = 0
feed_number(base, 2); count += 1
for x in range(3, bound+1, 2):
if count >= N: break
if sieve[x >> 1]:
feed_number(base, x); count += 1
# Matrix exponentiation: base^K applied to [0,0,0,1]
state = [0, 0, 0, 1]; e = K
while e:
if e & 1: state = mat_vec(base, state)
base = mat_mul(base, base); e >>= 1
return str(state[1])
if __name__ == '__main__':
print(solve())
Java
public class Euler603 {
static final int MOD = 1000000007;
static int addmod(int a, int b) {
int s = a + b;
if (s >= MOD)
s -= MOD;
return s;
}
static int mulmod(long a, long b) {
return (int) ((a * b) % MOD);
}
static long[][] matIdentity() {
long[][] m = new long[4][4];
for (int i = 0; i < 4; i++)
m[i][i] = 1;
return m;
}
static long[][] matMul(long[][] x, long[][] y) {
long[][] z = new long[4][4];
for (int i = 0; i < 4; i++) {
for (int k = 0; k < 4; k++) {
if (x[i][k] == 0)
continue;
for (int j = 0; j < 4; j++) {
z[i][j] = (z[i][j] + x[i][k] * y[k][j]) % MOD;
}
}
}
return z;
}
static long[] matVecMul(long[][] m, long[] v) {
long[] r = new long[4];
for (int i = 0; i < 4; i++) {
long acc = 0;
for (int k = 0; k < 4; k++) {
acc += m[i][k] * v[k];
acc %= MOD;
}
r[i] = acc;
}
return r;
}
static void feedDigit(long[][] m, int d) {
for (int col = 0; col < 4; col++) {
int f = (int) m[0][col];
int s = (int) m[1][col];
int t = (int) m[2][col];
int one = (int) m[3][col];
int dt = mulmod(d, t);
int done = mulmod(d, one);
int tenf = mulmod(10, f);
int f2 = addmod(addmod(tenf, dt), done);
int s2 = addmod(s, f2);
int t2 = addmod(t, one);
m[0][col] = f2;
m[1][col] = s2;
m[2][col] = t2;
m[3][col] = one;
}
}
static void feedNumber(long[][] m, int x) {
if (x == 0) {
feedDigit(m, 0);
return;
}
int[] buf = new int[16];
int len = 0;
while (x > 0) {
buf[len++] = x % 10;
x /= 10;
}
for (int i = len - 1; i >= 0; i--) {
feedDigit(m, buf[i]);
}
}
static int upperBoundNthPrime(int n) {
if (n < 6)
return 15;
double nd = (double) n;
double bound = nd * (Math.log(nd) + Math.log(Math.log(nd)));
return (int) Math.ceil(bound + 100.0);
}
static long[][] matrixForPPrimes(int n) {
int limit = upperBoundNthPrime(n);
byte[] comp = new byte[(limit >> 1) + 1];
int r = (int) Math.sqrt(limit);
for (int p = 3; p <= r; p += 2) {
if (comp[p >> 1] == 1)
continue;
int step = p << 1;
for (int x = p * p; x <= limit; x += step) {
comp[x >> 1] = 1;
}
}
long[][] m = matIdentity();
int count = 1;
feedNumber(m, 2);
if (count == n)
return m;
for (int x = 3; x <= limit; x += 2) {
if (comp[x >> 1] == 0) {
count++;
feedNumber(m, x);
if (count == n)
return m;
}
}
return m;
}
public static String solve() {
int N = 1000000;
long K = 1000000000000L;
long[][] base = matrixForPPrimes(N);
long[] state = { 0, 0, 0, 1 };
long e = K;
while (e > 0) {
if ((e & 1) == 1)
state = matVecMul(base, state);
base = matMul(base, base);
e >>= 1;
}
return Long.toString(state[1]);
}
public static void main(String[] args) {
System.out.println(solve());
}
}