Problem 464: Möbius Function and Intervals

Project Euler Problem 464 Solution

EulerSolve provides an optimized solution for Project Euler Problem 464, Möbius Function and Intervals, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Let $C(n)$ be the number of intervals $I=[l,r]\subseteq [1,n]$ for which the values $\mu(k)=1$ and $\mu(k)=-1$ stay in almost perfect balance. If we write $$P_I=\#\{k\in I:\mu(k)=1\},\qquad N_I=\#\{k\in I:\mu(k)=-1\},$$ then the counted intervals are exactly those satisfying $$99P_I\le 100N_I,\qquad 99N_I\le 100P_I.$$ Values with $\mu(k)=0$ are neutral: they lengthen the interval but do not change either count. The implementations verify the derivation with $$C(10)=13,\qquad C(500)=16676,\qquad C(10000)=20155319,$$ and then evaluate the same method at $n=20000000$. Mathematical Approach The key reduction is to replace interval counting by a pair-counting problem on weighted prefix sums. Step 1: Generate the Möbius Values Efficiently We need $\mu(1),\mu(2),\dots,\mu(n)$. The Möbius function takes the values $$\mu(m)=\begin{cases} 1, & m=1,\\ (-1)^t, & m\text{ is a product of }t\text{ distinct primes},\\ 0, & m\text{ is divisible by }p^2\text{ for some prime }p. \end{cases}$$ The implementations build these values with a linear sieve, so the preprocessing is essentially one forward pass through the integers....

Detailed mathematical approach

Problem Summary

Let $C(n)$ be the number of intervals $I=[l,r]\subseteq [1,n]$ for which the values $\mu(k)=1$ and $\mu(k)=-1$ stay in almost perfect balance. If we write

$$P_I=\#\{k\in I:\mu(k)=1\},\qquad N_I=\#\{k\in I:\mu(k)=-1\},$$

then the counted intervals are exactly those satisfying

$$99P_I\le 100N_I,\qquad 99N_I\le 100P_I.$$

Values with $\mu(k)=0$ are neutral: they lengthen the interval but do not change either count. The implementations verify the derivation with

$$C(10)=13,\qquad C(500)=16676,\qquad C(10000)=20155319,$$

and then evaluate the same method at $n=20000000$.

Mathematical Approach

The key reduction is to replace interval counting by a pair-counting problem on weighted prefix sums.

Step 1: Generate the Möbius Values Efficiently

We need $\mu(1),\mu(2),\dots,\mu(n)$. The Möbius function takes the values

$$\mu(m)=\begin{cases} 1, & m=1,\\ (-1)^t, & m\text{ is a product of }t\text{ distinct primes},\\ 0, & m\text{ is divisible by }p^2\text{ for some prime }p. \end{cases}$$

The implementations build these values with a linear sieve, so the preprocessing is essentially one forward pass through the integers.

Step 2: Express the Interval Condition with Two Counts

For an interval $I=[l,r]$, define

$$P_I=\#\{k\in I:\mu(k)=1\},\qquad N_I=\#\{k\in I:\mu(k)=-1\}.$$

The interval is counted precisely when neither sign dominates the other by more than the factor $100/99$:

$$99P_I\le 100N_I,\qquad 99N_I\le 100P_I.$$

This formulation makes clear that entries with $\mu(k)=0$ do not affect the balance at all. In particular, an interval containing only zeros is valid because then $P_I=N_I=0$.

Step 3: Pass to the Complementary Bad Families

Instead of counting good intervals directly, count the two ways an interval can fail:

$$\mathcal{A}=\{I:99P_I-100N_I\gt 0\},\qquad \mathcal{B}=\{I:99N_I-100P_I\gt 0\}.$$

These two families are disjoint. Indeed, if both inequalities held, then adding them would give

$$-(P_I+N_I)\gt 0,$$

which is impossible. Therefore

$$C(n)=\frac{n(n+1)}{2}-\#\mathcal{A}-\#\mathcal{B}.$$

The total number of intervals appears because there are $n(n+1)/2$ choices of $[l,r]\subseteq [1,n]$.

Step 4: Turn Each Bad Family into a Weighted Prefix-Sum Problem

Define two weight systems:

$$a(k)=\begin{cases} 99, & \mu(k)=1,\\ -100, & \mu(k)=-1,\\ 0, & \mu(k)=0, \end{cases}\qquad b(k)=\begin{cases} -100, & \mu(k)=1,\\ 99, & \mu(k)=-1,\\ 0, & \mu(k)=0. \end{cases}$$

Now define prefix sums

$$A_j=\sum_{k=1}^{j}a(k),\qquad B_j=\sum_{k=1}^{j}b(k),\qquad A_0=B_0=0.$$

For an interval $I=[l,r]$,

$$A_r-A_{l-1}=99P_I-100N_I,\qquad B_r-B_{l-1}=99N_I-100P_I.$$

$$I\in\mathcal{A}\iff A_r\gt A_{l-1},\qquad I\in\mathcal{B}\iff B_r\gt B_{l-1}.$$

This converts interval counting into the problem of counting ordered pairs $(i,j)$ with $0\le i\lt j\le n$ and a strict prefix inequality.

Step 5: Count Increasing Prefix Pairs

For the first bad family we must count

$$\#\mathcal{A}=\#\{(i,j):0\le i\lt j\le n,\ A_j\gt A_i\},$$

and similarly

$$\#\mathcal{B}=\#\{(i,j):0\le i\lt j\le n,\ B_j\gt B_i\}.$$

Scanning the prefixes from left to right, we only need to know how many earlier prefix values are strictly smaller than the current one. A Fenwick tree supports exactly this operation after the prefix values are shifted into a dense range or compressed into sorted coordinates.

Hence the full problem becomes a pair of standard order-statistics passes over two prefix arrays.

Worked Example: $n=10$

The Möbius values are

$$\bigl(\mu(1),\dots,\mu(10)\bigr)=(1,-1,-1,0,-1,1,-1,0,0,1).$$

Some representative intervals are easy to classify:

$$[1,1]:\ P_I=1,\ N_I=0\ \Rightarrow\ 99P_I-100N_I=99\gt 0,$$

so $[1,1]\in\mathcal{A}$.

$$[2,3]:\ P_I=0,\ N_I=2\ \Rightarrow\ 99N_I-100P_I=198\gt 0,$$

so $[2,3]\in\mathcal{B}$.

$$[1,2]:\ P_I=1,\ N_I=1,$$

hence both balance inequalities hold and $[1,2]$ is counted. Also $[8,9]$ contains only zeros, so $P_I=N_I=0$ and it is counted as well.

For the first weight system, the prefix sums are

$$A=(0,99,-1,-101,-101,-201,-102,-202,-202,-202,-103),$$

and for the second they are

$$B=(0,-100,-1,98,98,197,97,196,196,196,96).$$

The number of increasing pairs in the first sequence is $6$, and in the second sequence it is $36$. Since the total number of intervals is

$$\frac{10\cdot 11}{2}=55,$$

we obtain

$$C(10)=55-6-36=13,$$

which matches the checkpoint used by the implementations.

How the Code Works

The C++, Python, and Java implementations all follow the same pipeline. First they compute the Möbius values up to $n$ with a linear sieve. Then they run the interval-counting routine twice: once with weights $(99,-100)$ and once with weights $(-100,99)$. In each pass they build the weighted prefix sums, map those prefix values to Fenwick-tree indices, and accumulate how many earlier prefixes are strictly smaller than the current prefix. After obtaining the two bad counts, they subtract them from $n(n+1)/2$.

The C++ and Java implementations also include a coordinate-compression fallback when the raw prefix-value span would make a dense Fenwick tree unnecessarily large. The mathematical quantity being counted is the same in every language.

Complexity Analysis

Computing all Möbius values up to $n$ by linear sieve costs $O(n)$ time and $O(n)$ memory. Each weighted pass constructs one prefix array in $O(n)$ time and performs Fenwick queries and updates in $O(n\log n)$ time in the worst case; when the prefix span is already dense, the logarithm is taken over that span instead of over the number of distinct compressed values. Therefore the overall method runs in $O(n\log n)$ time with $O(n)$ memory, a dramatic improvement over naive $O(n^2)$ interval enumeration.

Footnotes and References

Problem page: https://projecteuler.net/problem=464
Möbius function: Wikipedia - Möbius function
Fenwick tree: Wikipedia - Fenwick tree
Prefix sum: Wikipedia - Prefix sum

Problem 464 source code

C++

#include <algorithm>
#include <cstdint>
#include <iostream>
#include <string>
#include <vector>

namespace {

using i64 = std::int64_t;

struct Options {
    int n = 20'000'000;
    bool run_checkpoints = true;
};

class Fenwick {
public:
    explicit Fenwick(const int n) : bit_(static_cast<std::size_t>(n + 1), 0) {}

    void add(int idx, const int delta) {
        const int n = static_cast<int>(bit_.size()) - 1;
        while (idx <= n) {
            bit_[static_cast<std::size_t>(idx)] += delta;
            idx += idx & -idx;
        }
    }

    int sum(int idx) const {
        int result = 0;
        while (idx > 0) {
            result += bit_[static_cast<std::size_t>(idx)];
            idx -= idx & -idx;
        }
        return result;
    }

private:
    std::vector<int> bit_;
};

bool parse_int_after_prefix(const std::string& arg, const std::string& prefix, int& out) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }
    try {
        out = std::stoi(tail);
    } catch (...) {
        return false;
    }
    return true;
}

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (parse_int_after_prefix(arg, "--n=", options.n)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    if (options.n <= 0) {
        std::cerr << "--n must be positive.\n";
        return false;
    }
    return true;
}

std::vector<std::int8_t> mobius_values(const int n) {
    std::vector<int> lp(static_cast<std::size_t>(n + 1), 0);
    std::vector<int> primes;
    primes.reserve(static_cast<std::size_t>(n / 10));
    std::vector<std::int8_t> mu(static_cast<std::size_t>(n + 1), 0);
    mu[1] = 1;

    for (int i = 2; i <= n; ++i) {
        if (lp[static_cast<std::size_t>(i)] == 0) {
            lp[static_cast<std::size_t>(i)] = i;
            primes.push_back(i);
            mu[static_cast<std::size_t>(i)] = -1;
        }
        for (const int p : primes) {
            const int64_t x = static_cast<int64_t>(p) * static_cast<int64_t>(i);
            if (x > n || p > lp[static_cast<std::size_t>(i)]) {
                break;
            }
            lp[static_cast<std::size_t>(x)] = p;
            if (p == lp[static_cast<std::size_t>(i)]) {
                mu[static_cast<std::size_t>(x)] = 0;
                break;
            }
            mu[static_cast<std::size_t>(x)] =
                static_cast<std::int8_t>(-mu[static_cast<std::size_t>(i)]);
        }
    }
    return mu;
}

i64 count_positive_weighted_segments(
    const std::vector<std::int8_t>& mu,
    const int weight_pos,
    const int weight_neg) {
    const int n = static_cast<int>(mu.size()) - 1;
    std::vector<int> prefix(static_cast<std::size_t>(n + 1), 0);

    int cur = 0;
    int min_prefix = 0;
    int max_prefix = 0;
    for (int i = 1; i <= n; ++i) {
        const int m = static_cast<int>(mu[static_cast<std::size_t>(i)]);
        if (m == 1) cur += weight_pos;
        if (m == -1) cur += weight_neg;
        prefix[static_cast<std::size_t>(i)] = cur;
        if (cur < min_prefix) min_prefix = cur;
        if (cur > max_prefix) max_prefix = cur;
    }

    const i64 span = static_cast<i64>(max_prefix) - static_cast<i64>(min_prefix) + 1LL;
    constexpr i64 DENSE_LIMIT = 50'000'000LL;

    i64 positives = 0;
    if (span <= DENSE_LIMIT) {
        Fenwick bit(static_cast<int>(span) + 2);
        for (int i = 0; i <= n; ++i) {
            const int idx = prefix[static_cast<std::size_t>(i)] - min_prefix + 1;
            positives += static_cast<i64>(bit.sum(idx - 1));
            bit.add(idx, 1);
        }
        return positives;
    }

    std::vector<int> coords = prefix;
    std::sort(coords.begin(), coords.end());
    coords.erase(std::unique(coords.begin(), coords.end()), coords.end());

    Fenwick bit(static_cast<int>(coords.size()) + 2);
    for (int i = 0; i <= n; ++i) {
        const int v = prefix[static_cast<std::size_t>(i)];
        const int idx = static_cast<int>(
            std::lower_bound(coords.begin(), coords.end(), v) - coords.begin()) + 1;
        positives += static_cast<i64>(bit.sum(idx - 1));
        bit.add(idx, 1);
    }
    return positives;
}

i64 solve(const int n) {
    const std::vector<std::int8_t> mu = mobius_values(n);
    const i64 total_pairs = static_cast<i64>(n) * static_cast<i64>(n + 1) / 2LL;
    const i64 bad_a = count_positive_weighted_segments(mu, 99, -100);
    const i64 bad_b = count_positive_weighted_segments(mu, -100, 99);
    return total_pairs - bad_a - bad_b;
}

bool run_checkpoints() {
    if (solve(10) != 13LL) {
        std::cerr << "Checkpoint failed: C(10)\n";
        return false;
    }
    if (solve(500) != 16'676LL) {
        std::cerr << "Checkpoint failed: C(500)\n";
        return false;
    }
    if (solve(10'000) != 20'155'319LL) {
        std::cerr << "Checkpoint failed: C(10000)\n";
        return false;
    }
    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }

    if (options.run_checkpoints && !run_checkpoints()) {
        return 1;
    }

    std::cout << solve(options.n) << '\n';
    return 0;
}

Python

def solve():
    n = 20000000
    # Mobius sieve
    lp = [0] * (n + 1)
    primes = []
    mu = [0] * (n + 1)
    mu[1] = 1
    for i in range(2, n + 1):
        if lp[i] == 0:
            lp[i] = i; primes.append(i); mu[i] = -1
        for p in primes:
            v = i * p
            if v > n or p > lp[i]: break
            lp[v] = p
            if p == lp[i]: mu[v] = 0; break
            mu[v] = -mu[i]

    def count_pos_seg(wp, wn):
        prefix = [0] * (n + 1)
        cur = 0; mn = 0; mx = 0
        for i in range(1, n + 1):
            m = mu[i]
            if m == 1: cur += wp
            if m == -1: cur += wn
            prefix[i] = cur
            if cur < mn: mn = cur
            if cur > mx: mx = cur

        span = mx - mn + 1
        # Fenwick tree
        bit = [0] * (span + 2)
        def bit_add(idx, delta):
            while idx <= span + 1:
                bit[idx] += delta
                idx += idx & (-idx)
        def bit_sum(idx):
            s = 0
            while idx > 0:
                s += bit[idx]
                idx -= idx & (-idx)
            return s

        positives = 0
        for i in range(n + 1):
            idx = prefix[i] - mn + 1
            positives += bit_sum(idx - 1)
            bit_add(idx, 1)
        return positives

    total_pairs = n * (n + 1) // 2
    bad_a = count_pos_seg(99, -100)
    bad_b = count_pos_seg(-100, 99)
    return str(total_pairs - bad_a - bad_b)

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

Java

import java.util.ArrayList;
import java.util.List;
import java.util.Arrays;

public class Euler464 {

    static class Fenwick {
        int[] bit;

        Fenwick(int n) {
            bit = new int[n + 1];
        }

        void add(int idx, int delta) {
            while (idx < bit.length) {
                bit[idx] += delta;
                idx += idx & -idx;
            }
        }

        int sum(int idx) {
            int result = 0;
            while (idx > 0) {
                result += bit[idx];
                idx -= idx & -idx;
            }
            return result;
        }
    }

    static byte[] mobiusValues(int n) {
        int[] lp = new int[n + 1];
        List<Integer> primes = new ArrayList<>(n / 10);
        byte[] mu = new byte[n + 1];
        mu[1] = 1;

        for (int i = 2; i <= n; i++) {
            if (lp[i] == 0) {
                lp[i] = i;
                primes.add(i);
                mu[i] = -1;
            }
            for (int p : primes) {
                long x = (long) p * i;
                if (x > n || p > lp[i])
                    break;
                lp[(int) x] = p;
                if (p == lp[i]) {
                    mu[(int) x] = 0;
                    break;
                }
                mu[(int) x] = (byte) -mu[i];
            }
        }
        return mu;
    }

    static long countPositiveWeightedSegments(byte[] mu, int weightPos, int weightNeg) {
        int n = mu.length - 1;
        int[] prefix = new int[n + 1];

        int cur = 0;
        int minPrefix = 0;
        int maxPrefix = 0;

        for (int i = 1; i <= n; i++) {
            int m = mu[i];
            if (m == 1)
                cur += weightPos;
            else if (m == -1)
                cur += weightNeg;
            prefix[i] = cur;
            if (cur < minPrefix)
                minPrefix = cur;
            if (cur > maxPrefix)
                maxPrefix = cur;
        }

        long span = (long) maxPrefix - minPrefix + 1;
        long DENSE_LIMIT = 50000000L;
        long positives = 0;

        if (span <= DENSE_LIMIT) {
            Fenwick bit = new Fenwick((int) span + 2);
            for (int i = 0; i <= n; i++) {
                int idx = prefix[i] - minPrefix + 1;
                positives += bit.sum(idx - 1);
                bit.add(idx, 1);
            }
            return positives;
        }

        int[] coords = prefix.clone();
        Arrays.sort(coords);

        int uniqueCount = 0;
        if (coords.length > 0) {
            uniqueCount = 1;
            for (int i = 1; i < coords.length; i++) {
                if (coords[i] != coords[i - 1]) {
                    coords[uniqueCount++] = coords[i];
                }
            }
        }

        Fenwick bit = new Fenwick(uniqueCount + 2);
        for (int i = 0; i <= n; i++) {
            int v = prefix[i];
            int idx = Arrays.binarySearch(coords, 0, uniqueCount, v) + 1;
            positives += bit.sum(idx - 1);
            bit.add(idx, 1);
        }

        return positives;
    }

    public static String solve() {
        int n = 20000000;
        byte[] mu = mobiusValues(n);
        long totalPairs = (long) n * (n + 1) / 2L;

        long badA = countPositiveWeightedSegments(mu, 99, -100);
        long badB = countPositiveWeightedSegments(mu, -100, 99);

        return Long.toString(totalPairs - badA - badB);
    }

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