Problem 593: Fleeting Medians
View on Project EulerProject Euler Problem 593 Solution
EulerSolve provides an optimized solution for Project Euler Problem 593, Fleeting Medians, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary The sequence is defined by $$S(k)=p_k^k \bmod 10007,\qquad S_2(k)=S(k)+S\left(\left\lfloor\frac{k}{10000}\right\rfloor+1\right),$$ where \(p_k\) is the \(k\)-th prime. For each contiguous block \(S_2(i),S_2(i+1),\dots,S_2(i+K-1)\), let \(M(i,i+K-1)\) be its median, using the average of the two middle values when \(K\) is even. The target quantity is $$F(n,K)=\sum_{i=1}^{n-K+1} M(i,i+K-1).$$ A direct recomputation of every window median would be far too slow. The implementations succeed because the sequence can be generated online, every value lies in a very small fixed range, and a frequency structure can therefore answer median queries without re-sorting each window. Mathematical Approach The solution combines two observations: generating each term of \(S_2\) is cheap once the primes are streamed in order, and medians of a multiset can be recovered from prefix frequencies. Step 1: Generate \(S(k)\) efficiently from the prime stream Since \(10007\) is prime, every prime \(p_k\neq 10007\) is invertible modulo \(10007\). Therefore Fermat's little theorem allows the exponent to be reduced modulo $$\varphi(10007)=10006.$$ So, except for the single case \(p_k=10007\), we may compute $$p_k^k \bmod 10007 = p_k^{\,k \bmod 10006} \bmod 10007.$$ This turns each term into one fast modular exponentiation once the \(k\)-th prime is known....
Detailed mathematical approach
Problem Summary
The sequence is defined by
$$S(k)=p_k^k \bmod 10007,\qquad S_2(k)=S(k)+S\left(\left\lfloor\frac{k}{10000}\right\rfloor+1\right),$$
where \(p_k\) is the \(k\)-th prime. For each contiguous block \(S_2(i),S_2(i+1),\dots,S_2(i+K-1)\), let \(M(i,i+K-1)\) be its median, using the average of the two middle values when \(K\) is even. The target quantity is
$$F(n,K)=\sum_{i=1}^{n-K+1} M(i,i+K-1).$$
A direct recomputation of every window median would be far too slow. The implementations succeed because the sequence can be generated online, every value lies in a very small fixed range, and a frequency structure can therefore answer median queries without re-sorting each window.
Mathematical Approach
The solution combines two observations: generating each term of \(S_2\) is cheap once the primes are streamed in order, and medians of a multiset can be recovered from prefix frequencies.
Step 1: Generate \(S(k)\) efficiently from the prime stream
Since \(10007\) is prime, every prime \(p_k\neq 10007\) is invertible modulo \(10007\). Therefore Fermat's little theorem allows the exponent to be reduced modulo
$$\varphi(10007)=10006.$$
So, except for the single case \(p_k=10007\), we may compute
$$p_k^k \bmod 10007 = p_k^{\,k \bmod 10006} \bmod 10007.$$
This turns each term into one fast modular exponentiation once the \(k\)-th prime is known. The C++, Python, and Java implementations all estimate a safe upper bound for the \(n\)-th prime and use an odd-only sieve to enumerate primes in increasing order.
Step 2: Exploit the very small value range of \(S_2(k)\)
By construction,
$$0\le S(k)\le 10006.$$
Hence
$$0\le S_2(k)=S(k)+S\left(\left\lfloor\frac{k}{10000}\right\rfloor+1\right)\le 20012.$$
This bound is the key simplification. Every window value lies in the fixed domain
$$\{0,1,2,\dots,20012\},$$
whose size is only \(V=20013\). For the target instance \(n=10^7\), the shifted index \(\left\lfloor k/10000\right\rfloor+1\) never exceeds \(1001\), so the implementation only needs to remember those early sequence values to form the second summand.
Step 3: Rewrite the median as an order-statistics problem
For one window of length \(K\), let \(f_v\) be the number of occurrences of value \(v\), and define the prefix counts
$$P(v)=\sum_{x=0}^{v} f_x.$$
The \(t\)-th smallest value in the window is the unique \(v\) satisfying
$$P(v-1) \lt t \le P(v).$$
If the sorted window values are \(a_1\le a_2\le \cdots \le a_K\), then setting
$$t_1=\left\lfloor\frac{K+1}{2}\right\rfloor,\qquad t_2=\left\lfloor\frac{K+2}{2}\right\rfloor$$
gives the compact formula
$$M=\frac{a_{t_1}+a_{t_2}}{2}.$$
So every window median reduces to finding one or two order statistics inside a frequency table.
Step 4: Maintain those frequencies under a sliding window
Moving the window by one position changes only two elements: one outgoing value disappears and one incoming value appears. In frequency form this is
$$f_{v_{\text{out}}}\leftarrow f_{v_{\text{out}}}-1,\qquad f_{v_{\text{in}}}\leftarrow f_{v_{\text{in}}}+1.$$
A Fenwick tree stores the array \((f_0,f_1,\dots,f_{20012})\), supports both updates in \(O(\log V)\), and can also search for the smallest value whose prefix sum reaches a target rank \(t\). That makes median extraction \(O(\log V)\) per window as well.
Step 5: Keep the total integral by doubling every median
When \(K\) is even, a median may end in \(0.5\). Instead of using floating-point arithmetic, the implementation accumulates the doubled contribution
$$2M=a_{t_1}+a_{t_2}.$$
Equivalently, it sums
$$2F(n,K)=\sum_{i=1}^{n-K+1} \bigl(a_{t_1}^{(i)}+a_{t_2}^{(i)}\bigr),$$
where \(a_{t_1}^{(i)}\) and \(a_{t_2}^{(i)}\) are the two middle order statistics of the \(i\)-th window. Only at the very end is the result divided by \(2\), which guarantees an exact decimal ending in .0 or .5.
Worked Example: the first window of length \(10\)
For \(1\le k\le 10\), we have \(\left\lfloor k/10000\right\rfloor+1=1\), and \(S(1)=2\). Therefore the first ten values are
$$\{S_2(1),\dots,S_2(10)\}=\{4,11,127,2403,941,3437,1640,2867,6554,9242\}.$$
After sorting, the window becomes
$$\{4,11,127,941,1640,2403,2867,3437,6554,9242\}.$$
The middle positions are \(5\) and \(6\), so
$$M(1,10)=\frac{1640+2403}{2}=\frac{4043}{2}=2021.5.$$
The doubled contribution is \(4043\), exactly the value used by the implementation when checking this sample window.
How the Code Works
The C++, Python, and Java implementations follow the same pipeline. They first estimate a bound for the \(n\)-th prime and sieve only odd numbers, which is enough to stream primes in order without storing a full list of all candidates. As each prime arrives, the implementation computes the corresponding residue \(S(k)\), caches the early residues needed by the shifted term, and immediately forms the next \(S_2(k)\).
The current window is maintained by combining a circular buffer with a Fenwick tree over the value range \(0\) through \(20012\). Filling the first \(K\) values creates the initial multiset. Each later step removes the outgoing value, inserts the new one, and asks the Fenwick tree for the middle rank or the two middle ranks. The running sum is kept doubled, so no rounding is needed. The C++ implementation also verifies the sample medians and sample sums from the problem statement before evaluating the full instance; the Python and Java versions implement the same algorithmic idea with the same mathematics.
Complexity Analysis
Let \(B\) be the sieve bound used for the \(n\)-th prime, and let \(V=20013\). Building the odd-only sieve costs \(O(B\log\log B)\) time and uses sieve storage proportional to \(B\). Streaming the sequence contributes one modular exponentiation per prime, and because the modulus is fixed and the exponent is reduced modulo \(10006\), this behaves as constant-cost arithmetic per term. The sliding-window part performs two Fenwick updates and one or two rank searches per window, each in \(O(\log V)\) time. Since \(V\) is fixed, the median-maintenance cost is \(O(n\log V)\) with very small constants, and the additional working memory is \(O(V+K)\) besides the sieve.
Footnotes and References
- Problem page: https://projecteuler.net/problem=593
- Median: Wikipedia — Median
- Fenwick tree: Wikipedia — Fenwick tree
- Modular exponentiation: Wikipedia — Modular exponentiation
- Sieve of Eratosthenes: Wikipedia — Sieve of Eratosthenes
Problem 593 source code
C++
#include <algorithm>
#include <array>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <numeric>
#include <string>
#include <vector>
// Project Euler 593: Fleeting Medians
//
// S(k) = p_k^k mod 10007, where p_k is the k-th prime.
// S2(k) = S(k) + S(floor(k/10000) + 1).
//
// Values are small: S(k) in [0,10006], so S2(k) in [0,20012].
// For sliding-window medians we maintain a Fenwick tree of frequencies and query order
// statistics in O(log V) time.
//
// Sum of medians can be half-integer; we accumulate twice the answer to stay integral.
using u64 = std::uint64_t;
using u128 = unsigned __int128;
static void print_u128(u128 x) {
if (x == 0) {
std::cout << '0';
return;
}
char buf[64];
int n = 0;
while (x > 0) {
const u64 digit = (u64)(x % 10);
buf[n++] = (char)('0' + digit);
x /= 10;
}
while (n--) std::cout << buf[n];
}
static u64 nth_prime_upper_bound(u64 n) {
if (n < 6) return 15;
const long double nn = (long double)n;
const long double bound = nn * (logl(nn) + logl(logl(nn)));
return (u64)(bound + 16.0L);
}
struct Fenwick {
int n = 0;
std::vector<int> bit;
explicit Fenwick(int n_) : n(n_), bit(n_ + 1, 0) {}
void add(int idx0, int delta) {
// idx0 is 0-based value; Fenwick is 1-based.
int i = idx0 + 1;
for (; i <= n; i += i & -i) bit[i] += delta;
}
int kth(int k) const {
// Smallest value v (0-based) such that prefix sum >= k, assuming 1 <= k <= total.
int idx = 0;
int step = 1;
while (step << 1 <= n) step <<= 1;
for (; step; step >>= 1) {
const int nxt = idx + step;
if (nxt <= n && bit[nxt] < k) {
idx = nxt;
k -= bit[nxt];
}
}
return idx; // idx is 0-based because we return (idx+1)-1
}
};
static int powmod_int(int a, int e, int mod) {
int r = 1;
int x = a % mod;
while (e > 0) {
if (e & 1) r = (int)((1LL * r * x) % mod);
x = (int)((1LL * x * x) % mod);
e >>= 1;
}
return r;
}
struct OddCompositeBits {
// Composite flags for odd numbers up to limit.
// Index is x>>1 (so 3 -> 1, 5 -> 2, ...). Bit=1 means composite.
u64 limit;
std::vector<u64> bits;
explicit OddCompositeBits(u64 limit_) : limit(limit_) {
const u64 sz_bits = (limit >> 1) + 1;
bits.assign((sz_bits + 63) >> 6, 0ULL);
}
inline bool get(u64 half) const {
return (bits[half >> 6] >> (half & 63ULL)) & 1ULL;
}
inline void set(u64 half) {
bits[half >> 6] |= 1ULL << (half & 63ULL);
}
};
static void sieve_odd(OddCompositeBits& comp) {
const u64 limit = comp.limit;
const u64 r = (u64)std::sqrt((long double)limit);
for (u64 p = 3; p <= r; p += 2) {
if (comp.get(p >> 1)) continue;
const u64 step = p << 1;
for (u64 x = p * p; x <= limit; x += step) comp.set(x >> 1);
}
}
static std::vector<int> generate_S2(u64 n) {
constexpr int MOD = 10007;
constexpr int PHI = MOD - 1;
const u64 limit = nth_prime_upper_bound(n);
OddCompositeBits comp(limit);
sieve_odd(comp);
std::vector<int> Sfirst(1002, 0);
std::vector<int> S2(n + 1, 0);
u64 k = 0;
auto handle_prime = [&](u64 p) {
++k;
const int base = (int)(p % MOD);
int s = 0;
if (base != 0) {
const int e = (int)(k % PHI);
s = powmod_int(base, e, MOD);
}
if (k <= 1001) Sfirst[(size_t)k] = s;
const int idx = (int)(k / 10000 + 1);
const int s2 = s + Sfirst[(size_t)idx];
S2[(size_t)k] = s2;
};
handle_prime(2);
for (u64 x = 3; x <= limit && k < n; x += 2) {
if (!comp.get(x >> 1)) handle_prime(x);
}
assert(k == n);
return S2;
}
static u128 median2_from_sorted(const std::vector<int>& v) {
const size_t len = v.size();
if (len & 1) return (u128)2 * (u128)v[len / 2];
return (u128)v[len / 2 - 1] + (u128)v[len / 2];
}
static u128 median2_range(const std::vector<int>& S2, u64 l, u64 r) {
std::vector<int> tmp;
tmp.reserve((size_t)(r - l + 1));
for (u64 i = l; i <= r; ++i) tmp.push_back(S2[(size_t)i]);
std::sort(tmp.begin(), tmp.end());
return median2_from_sorted(tmp);
}
static u128 F2(u64 n, int K) {
constexpr int MOD = 10007;
constexpr int PHI = MOD - 1;
constexpr int VMAX = 20013; // values 0..20012
const u64 limit = nth_prime_upper_bound(n);
OddCompositeBits comp(limit);
sieve_odd(comp);
std::vector<int> Sfirst(1002, 0);
Fenwick fw(VMAX);
std::vector<std::uint16_t> ring((size_t)K);
int ptr = 0;
u128 sum2 = 0;
auto add_median = [&]() {
if ((K & 1) == 1) {
const int v = fw.kth((K + 1) / 2);
sum2 += (u128)2 * (u128)v;
} else {
const int a = fw.kth(K / 2);
const int b = fw.kth(K / 2 + 1);
sum2 += (u128)a + (u128)b;
}
};
u64 k = 0;
auto handle_prime = [&](u64 p) {
++k;
const int base = (int)(p % MOD);
int s = 0;
if (base != 0) {
const int e = (int)(k % PHI);
s = powmod_int(base, e, MOD);
}
if (k <= 1001) Sfirst[(size_t)k] = s;
const int idx = (int)(k / 10000 + 1);
const int v = s + Sfirst[(size_t)idx];
if ((int)k <= K) {
ring[(size_t)ptr++] = (std::uint16_t)v;
fw.add(v, +1);
if ((int)k == K) {
ptr = 0;
add_median();
}
} else {
const int out = (int)ring[(size_t)ptr];
fw.add(out, -1);
ring[(size_t)ptr] = (std::uint16_t)v;
fw.add(v, +1);
++ptr;
if (ptr == K) ptr = 0;
add_median();
}
};
handle_prime(2);
for (u64 x = 3; x <= limit && k < n; x += 2) {
if (!comp.get(x >> 1)) handle_prime(x);
}
assert(k == n);
return sum2;
}
int main() {
// Median validations from the statement.
{
const auto S2 = generate_S2(1000);
if (median2_range(S2, 1, 10) != 4043) {
std::cerr << "Validation failed: M(1,10)\n";
return 1;
}
if (median2_range(S2, 100, 1000) != 9430) {
std::cerr << "Validation failed: M(100,1000)\n";
return 1;
}
}
// F validations from the statement.
if (F2(100, 10) != 927257) {
std::cerr << "Validation failed: F(100,10)\n";
return 1;
}
if (F2(100000, 10000) != 1350696415ULL) {
std::cerr << "Validation failed: F(1e5,1e4)\n";
return 1;
}
const u128 ans2 = F2(10000000ULL, 100000);
print_u128(ans2 / 2);
std::cout << ((ans2 & 1) ? ".5\n" : ".0\n");
return 0;
}
Python
import math
def solve():
n = 10000000; K = 100000
MOD_S = 10007; PHI = MOD_S - 1; VMAX = 20013
def powmod(a, e, m):
r = 1; a %= m
while e > 0:
if e & 1: r = r * a % m
a = a * a % m; e >>= 1
return r
# Sieve primes up to nth_prime upper bound
bound = int(n * (math.log(n) + math.log(math.log(n)))) + 200
sieve = bytearray(b'\x01') * ((bound >> 1) + 1)
for p in range(3, int(bound**0.5)+1, 2):
if sieve[p >> 1]:
for x in range(p*p, bound+1, 2*p): sieve[x >> 1] = 0
# Fenwick tree for order statistics
bit = [0] * (VMAX + 1)
def fw_add(idx, d):
idx += 1
while idx <= VMAX: bit[idx] += d; idx += idx & (-idx)
def fw_kth(k):
idx = 0; step = 1
while step * 2 <= VMAX: step *= 2
while step:
nxt = idx + step
if nxt <= VMAX and bit[nxt] < k:
idx = nxt; k -= bit[nxt]
step >>= 1
return idx
Sfirst = [0] * 1002
ring = [0] * K; ptr = 0; k = 0
sum2 = 0
def add_median():
nonlocal sum2
if K & 1:
sum2 += 2 * fw_kth((K+1)//2)
else:
sum2 += fw_kth(K//2) + fw_kth(K//2+1)
def handle_prime(p):
nonlocal k, ptr
k += 1
base = p % MOD_S
s = powmod(base, k % PHI, MOD_S) if base else 0
if k <= 1001: Sfirst[k] = s
idx = k // 10000 + 1; v = s + Sfirst[idx]
if k <= K:
ring[ptr] = v; ptr += 1; fw_add(v, 1)
if k == K: ptr = 0; add_median()
else:
out = ring[ptr]; fw_add(out, -1)
ring[ptr] = v; fw_add(v, 1)
ptr += 1
if ptr == K: ptr = 0
add_median()
handle_prime(2)
for x in range(3, bound+1, 2):
if k >= n: break
if sieve[x >> 1]: handle_prime(x)
whole = sum2 // 2; frac = sum2 & 1
return f"{whole}.{'5' if frac else '0'}"
if __name__ == '__main__':
print(solve())
Java
public class Euler593 {
static long nthPrimeUpperBound(long n) {
if (n < 6)
return 15;
double nn = (double) n;
double bound = nn * (Math.log(nn) + Math.log(Math.log(nn)));
return (long) (bound + 16.0);
}
static class OddCompositeBits {
long limit;
long[] bits;
OddCompositeBits(long limit) {
this.limit = limit;
long szBits = (limit >> 1) + 1;
bits = new long[(int) ((szBits + 63) >> 6)];
}
boolean get(long half) {
return ((bits[(int) (half >> 6)] >> (half & 63L)) & 1L) == 1L;
}
void set(long half) {
bits[(int) (half >> 6)] |= (1L << (half & 63L));
}
}
static void sieveOdd(OddCompositeBits comp) {
long limit = comp.limit;
long r = (long) Math.sqrt((double) limit);
for (long p = 3; p <= r; p += 2) {
if (comp.get(p >> 1))
continue;
long step = p << 1;
for (long x = p * p; x <= limit; x += step) {
comp.set(x >> 1);
}
}
}
static int powmodInt(int a, int e, int mod) {
int r = 1;
int x = a % mod;
while (e > 0) {
if ((e & 1) != 0)
r = (int) ((1L * r * x) % mod);
x = (int) ((1L * x * x) % mod);
e >>= 1;
}
return r;
}
static class Fenwick {
int n;
int[] bit;
Fenwick(int n) {
this.n = n;
bit = new int[n + 1];
}
void add(int idx0, int delta) {
int i = idx0 + 1;
for (; i <= n; i += i & -i)
bit[i] += delta;
}
int kth(int k) {
int idx = 0;
int step = 1;
while (step << 1 <= n)
step <<= 1;
for (; step > 0; step >>= 1) {
int nxt = idx + step;
if (nxt <= n && bit[nxt] < k) {
idx = nxt;
k -= bit[nxt];
}
}
return idx;
}
}
public static String solve() {
int MOD = 10007;
int PHI = MOD - 1;
int VMAX = 20013;
long n = 10000000L;
int K = 100000;
long limit = nthPrimeUpperBound(n);
OddCompositeBits comp = new OddCompositeBits(limit);
sieveOdd(comp);
int[] Sfirst = new int[1002];
Fenwick fw = new Fenwick(VMAX);
int[] ring = new int[K];
int ptr = 0;
long sum2Long = 0;
long k_final = 0;
int kHalf1 = (K % 2 == 1) ? (K + 1) / 2 : K / 2;
int kHalf2 = (K % 2 == 1) ? kHalf1 : K / 2 + 1;
// Process 2
k_final++;
int baseVal = 2 % MOD;
int s = 0;
if (baseVal != 0) {
int e = (int) (k_final % PHI);
s = powmodInt(baseVal, e, MOD);
}
Sfirst[1] = s;
int v2 = s + Sfirst[1];
ring[ptr++] = v2;
fw.add(v2, 1);
for (long x = 3; x <= limit && k_final < n; x += 2) {
if (!comp.get(x >> 1)) {
k_final++;
baseVal = (int) (x % MOD);
s = 0;
if (baseVal != 0) {
int e = (int) (k_final % PHI);
if (e != 0) {
s = powmodInt(baseVal, e, MOD);
} else {
s = 1;
}
}
if (k_final <= 1001)
Sfirst[(int) k_final] = s;
int idx = (int) (k_final / 10000 + 1);
int v = s + Sfirst[idx];
if (k_final <= K) {
ring[ptr++] = v;
fw.add(v, 1);
if (k_final == K) {
ptr = 0;
if (K % 2 == 1) {
sum2Long += 2L * fw.kth(kHalf1);
} else {
sum2Long += (long) fw.kth(kHalf1) + fw.kth(kHalf2);
}
}
} else {
int out = ring[ptr];
fw.add(out, -1);
ring[ptr] = v;
fw.add(v, 1);
ptr++;
if (ptr == K)
ptr = 0;
if (K % 2 == 1) {
sum2Long += 2L * fw.kth(kHalf1);
} else {
sum2Long += (long) fw.kth(kHalf1) + fw.kth(kHalf2);
}
}
}
}
String out = Long.toString(sum2Long / 2);
out += (sum2Long % 2 == 1) ? ".5" : ".0";
return out;
}
public static void main(String[] args) {
System.out.println(solve());
}
}