Problem 478: Mixtures

Problem 478 source code

C++

#include <algorithm>
#include <array>
#include <atomic>
#include <cstdint>
#include <cstdlib>
#include <iostream>
#include <thread>
#include <vector>

using namespace std;

namespace {

constexpr int64_t MOD = 214358881LL;     // 11^8
constexpr int64_t PHI = 194871710LL;     // 10 * 11^7
constexpr int64_t PHI2 = 389743420LL;    // 2 * PHI

struct Pow2Mod {
    int64_t mod;
    array<int64_t, 32> table{};
    explicit Pow2Mod(int64_t mod_) : mod(mod_) {
        table[0] = 2 % mod;
        for (size_t i = 1; i < table.size(); ++i) {
            table[i] = (table[i - 1] * table[i - 1]) % mod;
        }
    }

    int64_t pow(int64_t exp) const {
        int64_t res = 1;
        size_t bit = 0;
        while (exp > 0) {
            if (exp & 1LL) res = (res * table[bit]) % mod;
            exp >>= 1;
            ++bit;
        }
        return res;
    }
};

void sieve_mu_phi(int n, vector<int>& mu, vector<int>& phi) {
    mu.assign(n + 1, 0);
    phi.assign(n + 1, 0);
    vector<int> primes;
    vector<uint8_t> is_comp(n + 1, 0);

    mu[1] = 1;
    phi[1] = 1;

    for (int i = 2; i <= n; ++i) {
        if (!is_comp[i]) {
            primes.push_back(i);
            mu[i] = -1;
            phi[i] = i - 1;
        }
        for (int p : primes) {
            int64_t v = 1LL * i * p;
            if (v > n) break;
            is_comp[v] = 1;
            if (i % p == 0) {
                mu[v] = 0;
                phi[v] = phi[i] * p;
                break;
            }
            mu[v] = -mu[i];
            phi[v] = phi[i] * (p - 1);
        }
    }
}

int64_t count_M_mod2phi(int n, const vector<int>& mu, const vector<int>& n_div) {
    int64_t total = 0;
    int64_t mertens = 0;
    for (int d = 1; d <= n; ++d) {
        int mu_d = mu[d];
        if (mu_d == 0) continue;
        int64_t k = n_div[d];
        __int128 cube = (k + 1);
        cube = cube * cube * cube;
        int64_t term = static_cast<int64_t>(cube % PHI2);
        total += static_cast<int64_t>(mu_d) * term;
        mertens += mu_d;
        if (total >= PHI2 || total <= -PHI2) total %= PHI2;
        if (mertens >= PHI2 || mertens <= -PHI2) mertens %= PHI2;
    }
    total = (total - mertens) % PHI2;
    if (total < 0) total += PHI2;
    return total;
}

int64_t compute_E(int n, unsigned threads) {
    vector<int> mu;
    vector<int> phi;
    sieve_mu_phi(n, mu, phi);

    vector<int> n_div(n + 1, 0);
    for (int d = 1; d <= n; ++d) n_div[d] = n / d;

    int64_t sum_phi_mod = 0;
    for (int i = 1; i <= n; ++i) {
        sum_phi_mod += phi[i];
        if (sum_phi_mod >= MOD) sum_phi_mod -= MOD;
    }
    int64_t D_mod = (6LL * sum_phi_mod) % MOD;

    int64_t N_mod2phi = count_M_mod2phi(n, mu, n_div);
    int64_t N_mod_phi = N_mod2phi % PHI;
    int64_t Nprime_mod2phi = (N_mod2phi + PHI2 - 1) % PHI2;
    if (Nprime_mod2phi & 1LL) {
        cerr << "Validation failed: N' not even." << '\n';
        return -1;
    }
    int64_t Nhalf_mod_phi = (Nprime_mod2phi / 2) % PHI;

    Pow2Mod pow2(MOD);
    int64_t pow2_N = pow2.pow(N_mod_phi);
    int64_t pow2_half = pow2.pow(Nhalf_mod_phi);

    if (threads == 0) threads = 1;
    if (n < 100000) threads = 1;
    threads = min<unsigned>(threads, static_cast<unsigned>(n));

    vector<int64_t> partial(threads, 0);
    atomic<int> nextL{1};
    const int chunk = 64;

    auto worker = [&](unsigned idx) {
        int64_t local = 0;
        while (true) {
            int start = nextL.fetch_add(chunk, memory_order_relaxed);
            if (start > n) break;
            int end = min(n, start + chunk - 1);

            for (int L = start; L <= end; ++L) {
                int m = n / L;
                int64_t total = 1;
                for (int d = 1; d <= m; ++d) {
                    int mu_d = mu[d];
                    if (mu_d == 0) continue;
                    int md = m / d;
                    int64_t term = 1LL * md * n_div[d] - 1LL * L * md * (md + 1) / 2;
                    total += static_cast<int64_t>(mu_d) * term;
                }
                int64_t g_mod = total % PHI;
                if (g_mod < 0) g_mod += PHI;
                int64_t exp = Nhalf_mod_phi - g_mod;
                if (exp < 0) exp += PHI;
                int64_t pow = pow2.pow(exp);
                int64_t mL = (6LL * phi[L]) % MOD;
                local += (mL * pow) % MOD;
                if (local >= MOD) local %= MOD;
            }
        }
        partial[idx] = local % MOD;
    };

    vector<thread> pool;
    pool.reserve(threads);
    for (unsigned t = 0; t < threads; ++t) pool.emplace_back(worker, t);
    for (auto& th : pool) th.join();

    int64_t S = 0;
    for (int64_t v : partial) {
        S += v;
        if (S >= MOD) S %= MOD;
    }

    int64_t F = (1 + (pow2_half * D_mod) % MOD - S) % MOD;
    if (F < 0) F += MOD;
    int64_t E = (pow2_N - F) % MOD;
    if (E < 0) E += MOD;
    return E;
}

bool validate() {
    struct Test {
        int n;
        int64_t expected;
    };
    const Test tests[] = {
        {1, 103},
        {2, 520447},
        {10, 82608406},
        {500, 13801403},
    };

    for (const auto& test : tests) {
        int64_t got = compute_E(test.n, 1);
        if (got != test.expected) {
            cerr << "Validation failed for n=" << test.n
                 << ": got " << got << ", expected " << test.expected << '\n';
            return false;
        }
    }
    return true;
}

} // namespace

int main(int argc, char** argv) {
    if (!validate()) return 1;

    int n = 10000000;
    if (argc > 1) {
        n = max(1, atoi(argv[1]));
    }
    unsigned threads = thread::hardware_concurrency();
    int64_t result = compute_E(n, threads);
    if (result < 0) return 1;

    cout << result << '\n';
    return 0;
}

Python

def solve():
    MOD = 214358881  # 11^8
    PHI = 194871710  # 10 * 11^7
    PHI2 = 2 * PHI
    n = 10000000

    # Sieve mu and phi
    mu = [0]*(n+1); phi = [0]*(n+1); mu[1] = 1; phi[1] = 1
    is_comp = bytearray(n+1); primes = []
    for i in range(2, n+1):
        if not is_comp[i]: primes.append(i); mu[i] = -1; phi[i] = i-1
        for p in primes:
            v = i*p
            if v > n: break
            is_comp[v] = 1
            if i % p == 0: mu[v] = 0; phi[v] = phi[i]*p; break
            mu[v] = -mu[i]; phi[v] = phi[i]*(p-1)

    n_div = [0]*(n+1)
    for d in range(1, n+1): n_div[d] = n // d

    def pow2(exp):
        r, b = 1, 2 % MOD
        while exp > 0:
            if exp & 1: r = r * b % MOD
            b = b * b % MOD; exp >>= 1
        return r

    # Count M mod 2*PHI
    total = 0; mertens = 0
    for d in range(1, n+1):
        if mu[d] == 0: continue
        k = n_div[d]
        cube = (k+1)**3 % PHI2
        total += mu[d] * cube; mertens += mu[d]
        if abs(total) >= PHI2: total %= PHI2
        if abs(mertens) >= PHI2: mertens %= PHI2
    N_mod2phi = (total - mertens) % PHI2
    N_mod_phi = N_mod2phi % PHI
    Np = (N_mod2phi - 1) % PHI2
    Nhalf_mod_phi = (Np // 2) % PHI

    pow2_N = pow2(N_mod_phi)
    pow2_half = pow2(Nhalf_mod_phi)

    sum_phi = 0
    for i in range(1, n+1):
        sum_phi += phi[i]
        if sum_phi >= MOD: sum_phi -= MOD
    D_mod = 6 * sum_phi % MOD

    S = 0
    for L in range(1, n+1):
        m = n // L
        total_g = 1
        for d in range(1, m+1):
            if mu[d] == 0: continue
            md = m // d
            term = md * n_div[d] - L * md * (md+1) // 2
            total_g += mu[d] * term
        g_mod = total_g % PHI
        exp = (Nhalf_mod_phi - g_mod) % PHI
        pw = pow2(exp)
        mL = 6 * phi[L] % MOD
        S = (S + mL * pw) % MOD

    F = (1 + pow2_half * D_mod - S) % MOD
    E = (pow2_N - F) % MOD
    return str(E)

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

Java

import java.util.ArrayList;
import java.util.List;
import java.util.concurrent.Callable;
import java.util.concurrent.ExecutionException;
import java.util.concurrent.ExecutorService;
import java.util.concurrent.Executors;
import java.util.concurrent.Future;
import java.util.concurrent.atomic.AtomicInteger;

public class Euler478 {

    private static final long MOD = 214358881L;
    private static final long PHI = 194871710L;
    private static final long PHI2 = 389743420L;

    static class Pow2Mod {
        long mod;
        long[] table = new long[32];

        Pow2Mod(long mod) {
            this.mod = mod;
            table[0] = 2 % mod;
            for (int i = 1; i < 32; i++) {
                table[i] = (table[i - 1] * table[i - 1]) % mod;
            }
        }

        long pow(long exp) {
            long res = 1;
            int bit = 0;
            while (exp > 0) {
                if ((exp & 1) == 1)
                    res = (res * table[bit]) % mod;
                exp >>= 1;
                bit++;
            }
            return res;
        }
    }

    private static void sieveMuPhi(int n, int[] mu, int[] phi) {
        byte[] isComp = new byte[n + 1];
        int[] primes = new int[n + 1]; // Size overhead is fine
        int primeCount = 0;

        mu[1] = 1;
        phi[1] = 1;

        for (int i = 2; i <= n; i++) {
            if (isComp[i] == 0) {
                primes[primeCount++] = i;
                mu[i] = -1;
                phi[i] = i - 1;
            }
            for (int j = 0; j < primeCount; j++) {
                int p = primes[j];
                long v = (long) i * p;
                if (v > n)
                    break;
                isComp[(int) v] = 1;
                if (i % p == 0) {
                    mu[(int) v] = 0;
                    phi[(int) v] = phi[i] * p;
                    break;
                }
                mu[(int) v] = -mu[i];
                phi[(int) v] = phi[i] * (p - 1);
            }
        }
    }

    private static long countMMod2Phi(int n, int[] mu, int[] nDiv) {
        long total = 0;
        long mertens = 0;
        for (int d = 1; d <= n; d++) {
            int muD = mu[d];
            if (muD == 0)
                continue;
            long k = nDiv[d];
            // BigInteger equivalent for cube
            long rem = (k + 1) % PHI2;
            long term = (rem * rem) % PHI2;
            term = (term * rem) % PHI2;
            // Better to do this correctly without overflow:
            // k+1 <= 10^7, so cube can be 10^21, which overflows long.
            // Using 128-bit or multiple steps with modulo:
            // But PHI2 is ~3.8 * 10^8, so (k+1)^3 mod PHI2 doesn't immediately fit into
            // long multiplication stages implicitly.
            long q1 = ((k + 1) * (k + 1)) % PHI2;
            long q2 = (q1 * ((k + 1) % PHI2)) % PHI2;

            total += (long) muD * q2;
            mertens += muD;
            if (total >= PHI2 || total <= -PHI2)
                total %= PHI2;
            if (mertens >= PHI2 || mertens <= -PHI2)
                mertens %= PHI2;
        }
        total = (total - mertens) % PHI2;
        if (total < 0)
            total += PHI2;
        return total;
    }

    private static long computeE(int n) throws InterruptedException, ExecutionException {
        int[] mu = new int[n + 1];
        int[] phi = new int[n + 1];
        sieveMuPhi(n, mu, phi);

        int[] nDiv = new int[n + 1];
        for (int d = 1; d <= n; d++)
            nDiv[d] = n / d;

        long sumPhiMod = 0;
        for (int i = 1; i <= n; i++) {
            sumPhiMod += phi[i];
            if (sumPhiMod >= MOD)
                sumPhiMod -= MOD;
        }
        long dMod = (6L * sumPhiMod) % MOD;

        long nMod2Phi = countMMod2Phi(n, mu, nDiv);
        long nModPhi = nMod2Phi % PHI;
        long nPrimeMod2Phi = (nMod2Phi + PHI2 - 1) % PHI2;
        if ((nPrimeMod2Phi & 1) != 0) {
            System.err.println("Validation failed: N' not even.");
            return -1;
        }
        long nHalfModPhi = (nPrimeMod2Phi / 2) % PHI;

        Pow2Mod pow2 = new Pow2Mod(MOD);
        long pow2N = pow2.pow(nModPhi);
        long pow2Half = pow2.pow(nHalfModPhi);

        int threads = Math.min(Runtime.getRuntime().availableProcessors(), n);
        if (n < 100000)
            threads = 1;

        AtomicInteger nextL = new AtomicInteger(1);
        int chunk = 64;

        ExecutorService executor = Executors.newFixedThreadPool(threads);
        List<Future<Long>> futures = new ArrayList<>();

        for (int t = 0; t < threads; t++) {
            futures.add(executor.submit(() -> {
                long local = 0;
                while (true) {
                    int start = nextL.getAndAdd(chunk);
                    if (start > n)
                        break;
                    int end = Math.min(n, start + chunk - 1);

                    for (int L = start; L <= end; L++) {
                        int m = n / L;
                        long total = 1;
                        for (int d = 1; d <= m; d++) {
                            int muD = mu[d];
                            if (muD == 0)
                                continue;
                            int md = m / d;
                            long term = 1L * md * nDiv[d] - 1L * L * md * (md + 1) / 2;
                            total += (long) muD * term;
                        }
                        long gMod = total % PHI;
                        if (gMod < 0)
                            gMod += PHI;
                        long exp = nHalfModPhi - gMod;
                        if (exp < 0)
                            exp += PHI;
                        long p2 = pow2.pow(exp);
                        long mL = (6L * phi[L]) % MOD;
                        local += (mL * p2) % MOD;
                        if (local >= MOD)
                            local %= MOD;
                    }
                }
                return local % MOD;
            }));
        }

        long s = 0;
        for (Future<Long> f : futures) {
            s += f.get();
            if (s >= MOD)
                s %= MOD;
        }
        executor.shutdown();

        long f = (1 + (pow2Half * dMod) % MOD - s) % MOD;
        if (f < 0)
            f += MOD;
        long e = (pow2N - f) % MOD;
        if (e < 0)
            e += MOD;
        return e;
    }

    public static void main(String[] args) throws Exception {
        int n = 10_000_000;
        System.out.println(computeE(n));
    }
}