Problem 944: Sum of Elevisors
View on Project EulerProject Euler Problem 944 Solution
EulerSolve provides an optimized solution for Project Euler Problem 944, Sum of Elevisors, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary For each subset \(A \subseteq \{1,2,\dots,n\}\), call an element \(x \in A\) an elevisor if \(x\) divides at least one other element of \(A\). The problem asks for the total sum of all elevisors over all subsets of \(\{1,\dots,n\}\) when \(n=10^{14}\), modulo \(M=1234567891\). Enumerating subsets is hopeless because there are \(2^n\) of them. The key move is to reverse the order of summation: instead of looping over subsets, fix a value \(x\) and count how many subsets make that specific \(x\) an elevisor. Mathematical Approach If \(E(A)\) denotes the set of elevisors of \(A\), then the desired quantity is \[ S(n)=\sum_{A\subseteq \{1,\dots,n\}} \sum_{x\in E(A)} x \pmod M. \] Counting subsets for one fixed value \(x\) Fix \(x\). Let \[ q=\left\lfloor \frac{n}{x} \right\rfloor. \] Then the multiples of \(x\) inside \(\{1,\dots,n\}\) are \(x,2x,\dots,qx\), so there are exactly \(q\) of them. For \(x\) to be an elevisor of a subset \(A\), two conditions are necessary and sufficient: \(x\) itself must be chosen, and at least one of the other \(q-1\) multiples of \(x\) must also be chosen. Every non-multiple of \(x\) may be chosen freely. Therefore the number of subsets in which \(x\) is an elevisor is \[ 2^{n-q}\left(2^{q-1}-1\right). \] This can be rewritten in the more convenient form \[ 2^{n-q}\left(2^{q-1}-1\right)=2^{n-1}-2^{n-q}....
Detailed mathematical approach
Problem Summary
For each subset \(A \subseteq \{1,2,\dots,n\}\), call an element \(x \in A\) an elevisor if \(x\) divides at least one other element of \(A\). The problem asks for the total sum of all elevisors over all subsets of \(\{1,\dots,n\}\) when \(n=10^{14}\), modulo \(M=1234567891\).
Enumerating subsets is hopeless because there are \(2^n\) of them. The key move is to reverse the order of summation: instead of looping over subsets, fix a value \(x\) and count how many subsets make that specific \(x\) an elevisor.
Mathematical Approach
If \(E(A)\) denotes the set of elevisors of \(A\), then the desired quantity is
\[ S(n)=\sum_{A\subseteq \{1,\dots,n\}} \sum_{x\in E(A)} x \pmod M. \]
Counting subsets for one fixed value \(x\)
Fix \(x\). Let
\[ q=\left\lfloor \frac{n}{x} \right\rfloor. \]
Then the multiples of \(x\) inside \(\{1,\dots,n\}\) are \(x,2x,\dots,qx\), so there are exactly \(q\) of them.
For \(x\) to be an elevisor of a subset \(A\), two conditions are necessary and sufficient: \(x\) itself must be chosen, and at least one of the other \(q-1\) multiples of \(x\) must also be chosen. Every non-multiple of \(x\) may be chosen freely.
Therefore the number of subsets in which \(x\) is an elevisor is
\[ 2^{n-q}\left(2^{q-1}-1\right). \]
This can be rewritten in the more convenient form
\[ 2^{n-q}\left(2^{q-1}-1\right)=2^{n-1}-2^{n-q}. \]
The interpretation is simple: \(2^{n-1}\) counts all subsets containing \(x\), and \(2^{n-q}\) removes the subsets where \(x\) is present but no other multiple of \(x\) is present.
Swapping the order of summation
Now sum the contribution of each possible elevisor separately. By linearity,
\[ S(n)=\sum_{x=1}^{n} x\left(2^{n-1}-2^{n-\left\lfloor n/x\right\rfloor}\right). \]
This is exactly the formula used by the C++, Python, and Java implementations. Splitting off the easy part gives
\[ S(n)=2^{n-1}\sum_{x=1}^{n} x-\sum_{x=1}^{n} x\,2^{n-\left\lfloor n/x\right\rfloor}. \]
Since \(\sum_{x=1}^{n}x=n(n+1)/2\), all of the real work is concentrated in the weighted sum
\[ W(n)=\sum_{x=1}^{n} x\,2^{n-\left\lfloor n/x\right\rfloor}. \]
Grouping by equal values of \(\left\lfloor n/x\right\rfloor\)
The quotient \(\left\lfloor n/x\right\rfloor\) stays constant on long intervals. If
\[ q=\left\lfloor \frac{n}{x} \right\rfloor, \]
then the corresponding values of \(x\) are exactly those in
\[ \left\lfloor \frac{n}{q+1} \right\rfloor+1 \le x \le \left\lfloor \frac{n}{q} \right\rfloor. \]
On such an interval the power of 2 is constant, so the whole block contributes
\[ 2^{n-q}\sum_{x=l}^{r}x = 2^{n-q}\frac{(l+r)(r-l+1)}{2}. \]
This turns many consecutive terms of \(W(n)\) into one arithmetic-progression calculation.
Worked example: \(n=10\)
Take \(x=2\). Then \(q=\lfloor 10/2\rfloor=5\), and the relevant multiples are \(2,4,6,8,10\). To make 2 an elevisor, we must include 2, include at least one of \(\{4,6,8,10\}\), and choose any subset of the five non-multiples \(\{1,3,5,7,9\}\). Hence 2 is an elevisor in
\[ 2^{5}(2^{4}-1)=480 \]
subsets, contributing \(2\cdot 480=960\) to \(S(10)\).
The quotient grouping is visible in the same example. For \(W(10)\), the quotient \(q=1\) occurs for \(x=6,7,8,9,10\), so that whole block contributes
\[ 2^{10-1}(6+7+8+9+10)=2^{9}\cdot 40. \]
This is the exact compression used on the large-\(x\) side of the final algorithm.
The two-phase split used by the implementations
The implementations choose a fixed cutoff \(B=3{,}000{,}000\). Values \(x\le B\) are handled directly, while values \(x>B\) are regrouped by the quotient \(q=\lfloor n/x\rfloor\).
Because \(x>B\) implies \(q\le \lfloor n/(B+1)\rfloor\), the remaining quotients are relatively small. The resulting decomposition is
\[ W(n) = \sum_{x=1}^{B} x\,2^{n-\left\lfloor n/x\right\rfloor} + \sum_{q=1}^{\left\lfloor n/(B+1)\right\rfloor} 2^{n-q} \sum_{x=\max\left(B+1,\left\lfloor \frac{n}{q+1}\right\rfloor+1\right)}^{\left\lfloor n/q\right\rfloor} x. \]
This is not merely a conceptual derivation; it is the exact formula evaluated by the code.
How the Code Works
Precomputing the fixed factors
The implementation first computes \(2^{n-1}\bmod M\), because it appears in the closed form of the easy summand. It also evaluates arithmetic-progression sums modulo \(M\) so that a whole interval \([l,r]\) can be handled without iterating through every integer in that interval.
Direct accumulation for small \(x\)
The first pass runs through \(x=1,2,\dots,\min(n,B)\). For each \(x\), it forms \(q=\lfloor n/x\rfloor\), computes \(2^{n-q}\bmod M\) by binary exponentiation, multiplies by \(x\), and adds the result to the running value of \(W(n)\).
Block accumulation for large \(x\)
The second pass runs through \(q=1,2,\dots,\lfloor n/(B+1)\rfloor\). For each quotient it reconstructs the interval of values \(x>B\) satisfying \(\lfloor n/x\rfloor=q\), evaluates the interval sum, and multiplies by the common factor \(2^{n-q}\).
At this stage the implementation maintains a useful invariant: at the beginning of the iteration for a given \(q\), the stored power is exactly \(2^{n-q}\bmod M\). Moving to the next quotient therefore only requires multiplication by \(2^{-1}\bmod M\). Since \(M\) is odd, \(2^{-1}\equiv (M+1)/2 \pmod M\).
Final subtraction
After \(W(n)\) has been computed, the programs evaluate \(2^{n-1}\cdot n(n+1)/2 \bmod M\) and subtract \(W(n)\) modulo \(M\). The C++ and Java implementations also perform small sanity checks, including \(S(1)=0\) and \(S(10)=4927\), before printing the value for \(n=10^{14}\).
Complexity Analysis
Let \(B=3{,}000{,}000\). The direct pass performs \(B\) iterations, and each one computes a modular power by repeated squaring, so this part costs \(O(B\log n)\).
The block pass performs \(\lfloor n/(B+1)\rfloor\) iterations, each with constant-time arithmetic, so it costs \(O(n/B)\). Therefore the implemented algorithm runs in
\[ O(B\log n+n/B) \]
time and uses \(O(1)\) extra memory. For the target value \(n=10^{14}\), this is practical while still computing the exact mathematical sum.
Footnotes and References
- Project Euler problem page: https://projecteuler.net/problem=944
- Divisibility and divisors: Wikipedia - Divisor
- Floor function: Wikipedia - Floor and ceiling functions
- Arithmetic progression: Wikipedia - Arithmetic progression
- Modular exponentiation: Wikipedia - Modular exponentiation
- Quotient grouping in divisor-style sums: cp-algorithms - Number of divisors / sum of divisors
Problem 944 source code
C++
#include <cassert>
#include <cstdint>
#include <iostream>
namespace {
using u64 = std::uint64_t;
using u128 = unsigned __int128;
constexpr u64 kMod = 1'234'567'891ULL;
constexpr u64 kSplit = 3'000'000ULL;
inline u64 add_mod(u64 a, u64 b) {
a += b;
if (a >= kMod) {
a -= kMod;
}
return a;
}
inline u64 sub_mod(u64 a, u64 b) {
return (a >= b) ? (a - b) : (a + kMod - b);
}
inline u64 mul_mod(u64 a, u64 b) {
return (a * b) % kMod;
}
u64 pow_mod(u64 base, u64 exp) {
u64 result = 1 % kMod;
base %= kMod;
while (exp > 0) {
if (exp & 1ULL) {
result = mul_mod(result, base);
}
base = mul_mod(base, base);
exp >>= 1ULL;
}
return result;
}
u64 sum_range_mod(u64 l, u64 r) {
if (l > r) {
return 0;
}
const u128 cnt = static_cast<u128>(r) - static_cast<u128>(l) + 1U;
const u128 s = static_cast<u128>(l) + static_cast<u128>(r);
const u128 total = (cnt * s) / 2U;
return static_cast<u64>(total % kMod);
}
u64 solve(u64 n) {
const u64 pow_n_minus_1 = pow_mod(2, n - 1);
const u64 split = (n < kSplit) ? n : kSplit;
u64 weighted_sum = 0;
for (u64 x = 1; x <= split; ++x) {
const u64 q = n / x;
const u64 p = pow_mod(2, n - q);
weighted_sum = add_mod(weighted_sum, mul_mod(x % kMod, p));
}
const u64 inv2 = (kMod + 1ULL) / 2ULL;
const u64 qmax = n / (split + 1ULL);
u64 p = pow_n_minus_1; // 2^(n-1)
for (u64 q = 1; q <= qmax; ++q) {
u64 l = n / (q + 1ULL) + 1ULL;
if (l <= split) {
l = split + 1ULL;
}
const u64 r = n / q;
if (l <= r) {
const u64 sum_x = sum_range_mod(l, r);
weighted_sum = add_mod(weighted_sum, mul_mod(sum_x, p));
}
p = mul_mod(p, inv2); // 2^(n-(q+1))
}
const u64 tri = sum_range_mod(1, n);
const u64 first = mul_mod(pow_n_minus_1, tri);
return sub_mod(first, weighted_sum);
}
u64 solve_naive(u64 n) {
u64 ans = 0;
const u64 pow_n_minus_1 = pow_mod(2, n - 1);
for (u64 x = 1; x <= n; ++x) {
const u64 q = n / x;
const u64 ways = sub_mod(pow_n_minus_1, pow_mod(2, n - q));
ans = add_mod(ans, mul_mod(x % kMod, ways));
}
return ans;
}
void run_validations() {
assert(solve(1) == 0ULL);
assert(solve(10) == 4'927ULL);
for (u64 n = 2; n <= 200; ++n) {
assert(solve(n) == solve_naive(n));
}
}
} // namespace
int main() {
run_validations();
constexpr u64 kN = 100'000'000'000'000ULL;
std::cout << solve(kN) << '\n';
return 0;
}
Python
def solve():
MOD = 1234567891
N = 100000000000000
SPLIT = 3000000
def pow_mod(base, exp):
r = 1; base %= MOD
while exp > 0:
if exp & 1: r = r * base % MOD
base = base * base % MOD; exp >>= 1
return r
def sum_range_mod(l, r):
if l > r: return 0
cnt = r - l + 1; s = l + r
return cnt * s // 2 % MOD
pow_n1 = pow_mod(2, N - 1)
split = min(N, SPLIT)
wsum = 0
for x in range(1, split + 1):
q = N // x
p = pow_mod(2, N - q)
wsum = (wsum + x % MOD * p) % MOD
inv2 = (MOD + 1) // 2
qmax = N // (split + 1)
p = pow_n1
for q in range(1, qmax + 1):
l = N // (q + 1) + 1
if l <= split: l = split + 1
r = N // q
if l <= r:
sx = sum_range_mod(l, r)
wsum = (wsum + sx * p) % MOD
p = p * inv2 % MOD
tri = sum_range_mod(1, N)
first = pow_n1 * tri % MOD
return str((first - wsum) % MOD)
if __name__ == '__main__':
print(solve())
Java
public class Euler944 {
static final long kMod = 1234567891L;
static final long kSplit = 3000000L;
static long addMod(long a, long b) {
a += b;
if (a >= kMod) {
a -= kMod;
}
return a;
}
static long subMod(long a, long b) {
return (a >= b) ? (a - b) : (a + kMod - b);
}
static long mulMod(long a, long b) {
return (a * b) % kMod;
}
static long powMod(long base, long exp) {
long result = 1 % kMod;
base %= kMod;
while (exp > 0) {
if ((exp & 1) != 0) {
result = mulMod(result, base);
}
base = mulMod(base, base);
exp >>= 1;
}
return result;
}
static long sumRangeMod(long l, long r) {
if (l > r) {
return 0;
}
long cnt = (r - l + 1) % kMod;
long s = (l + r) % kMod;
long total = (cnt * s) % kMod;
long inv2 = (kMod + 1L) / 2L;
return (total * inv2) % kMod;
}
public static String solve(long n) {
long powNMinus1 = powMod(2, n - 1);
long split = Math.min(n, kSplit);
long weightedSum = 0;
for (long x = 1; x <= split; ++x) {
long q = n / x;
long p = powMod(2, n - q);
weightedSum = addMod(weightedSum, mulMod(x % kMod, p));
}
long inv2 = (kMod + 1L) / 2L;
long qmax = n / (split + 1L);
long p = powNMinus1;
for (long q = 1; q <= qmax; ++q) {
long l = n / (q + 1L) + 1L;
if (l <= split) {
l = split + 1L;
}
long r = n / q;
if (l <= r) {
long sumX = sumRangeMod(l, r);
weightedSum = addMod(weightedSum, mulMod(sumX, p));
}
p = mulMod(p, inv2);
}
long tri = sumRangeMod(1, n);
long first = mulMod(powNMinus1, tri);
return Long.toString(subMod(first, weightedSum));
}
public static void main(String[] args) {
if (!solve(1).equals("0") || !solve(10).equals("4927")) {
System.out.println("Validation failed");
return;
}
System.out.println(solve(100000000000000L));
}
}