Project Euler 464

发表于 2026-01-21 分类于 Project Euler 阅读次数：本文字数： 2.5k 阅读时长 ≈ 2 分钟

Project Euler 464

题目

Möbius function and intervals

The Möbius function, denoted \(mu(n)\), is defined as:

\(\mu(n) = (-1)^{\omega(n)}\) if \(n\) is squarefree (where \(\omega(n)\) is the number of distinct prime factors of \(n\))
\(\mu(n) = 0\) if \(n\) is not squarefree.

Let \(P(a,b)\) be the number of integers \(n\) in the interval \([a,b]\) such that \(\mu(n) = 1\).

Let \(N(a,b)\) be the number of integers \(n\) in the interval \([a,b]\) such that \(\mu(n) = -1\).

For example, \(P(2,10) = 2\) and \(N(2,10) = 4\).

Let \(C(n)\) be the number of integer pairs \((a,b)\) such that:

\(1\le a\le b\le n\),
\(99\cdot N(a,b)\le100\cdot P(a,b)\), and
\(99\cdot P(a,b)\le100\cdot N(a,b)\).

For example, \(C(10)=13, C(500)=16676\) and \(C(10000)=20155319\).

Find \(C(20000000)\).

解决方案

令 \(N=20000000\)。定义两个前缀函数：

\(p(m)=\#\{1\le k\le m:\mu(k)=1\}\)，并约定 \(p(0)=0\)
\(n(m)=\#\{1\le k\le m:\mu(k)=-1\}\)，并约定 \(n(0)=0\)

则对任意区间 \([a,b]\)（\(1\le a\le b\le N\)）有：\(P(a,b)=p(b)-p(a-1), N(a,b)=n(b)-n(a-1).\)

把题目中的两条不等式分别展开，移项后得到：

\(100p(b)-99n(b)\ge100p(a-1)-99n(a-1).\)
\(100n(b)-99p(b)\ge100n(a-1)-99p(a-1).\)

因此定义两个前缀势能：\(X(t)=100p(t)-99n(t), Y(t)=100n(t)-99p(t).\)

则区间 \([a,b]\) 同时满足两条条件，当且仅当\(X(b)\ge X(a-1)\land Y(b)\ge Y(a-1).\)

把 \((a,b)\) 映射成前缀端点 \((i,j)=(a-1,b)\)，则 \(0\le i<j\le N\)，并且需要\(X(i)\le X(j), Y(i)\le Y(j).\)

到这里看起来是二维偏序计数问题。

对任意非负整数 \(P_0,N_0\)（这里就是某个区间内的 \(P_0=P(a,b),N_0=N(a,b)\)），考虑两条条件：

\(99N_0\le 100P_0\)
\(99P_0\le 100N_0\)

若二者都不成立，则必须同时有严格不等式：\(99N_0>100P_0\land 99P_0>100N_0.\)

当 \(PN>0\) 时，两式相乘得到\(9801PN>10000PN,\)矛盾。因此，当 \(PN=0\) 时：

若 \(P=0,N>0\)，则第二条 \(99P\le 100N\) 必成立；
若 \(N=0,P>0\)，则第一条 \(99N\le 100P\) 必成立；
若 \(P=N=0\)，两条都成立。

因此对每个区间 \([a,b]\)，两条条件中至少有一条成立。设

\(I_1(a,b)\) 为第一条是否成立的指示变量（成立为 \(1\)，否则为 \(0\)）
\(I_2(a,b)\) 为第二条是否成立的指示变量

则对任意区间都有 \(I_1+I_2\ge 1\)，并且：

当且仅当两条都成立时，\(I_1+I_2-1=1\)
若恰有一条成立，则 \(I_1+I_2-1=0\)

于是

\[ C(N)=\sum_{1\le a\le b\le N}\bigl(I_1(a,b)+I_2(a,b)-1\bigr) = C_1(N)+C_2(N)-\frac{N(N+1)}2, \]

其中

\(C_1(N)=\#\{(a,b):99N(a,b)\le 100P(a,b)\}\)
\(C_2(N)=\#\{(a,b):99P(a,b)\le 100N(a,b)\}\)

也就是说，只要分别统计单条不等式成立的区间数，再减去总区间数即可。

接下来针对单条不等式进行统计，变成一维前缀非逆序对计数。以第一条为例，等价于\(X(b)-X(a-1)=100P(a,b)-99N(a,b)\ge 0\Longleftrightarrow X(a-1)\le X(b).\)

所以 \(C_1(N)\) 等于前缀序列 \(X(0),X(1),\dots,X(N)\) 中满足\(0\le i<j\le N, X(i)\le X(j)\)的配对数量。

这就是经典的统计非逆序对（等价于总对数减逆序对数）。做法是按时间顺序扫描 \(j\)，维护一个数据结构，支持：

插入已出现的 \(X(i)\)
查询已出现的值中 \(\le X(j)\) 的个数

用树状数组在离散化后的坐标上即可做到 \(O(\log M)\) 每次。

第二条同理，只是把序列换成 \(Y(t)\)。

我们无需显式存 \(p(t),n(t)\)以获得序列\(X(t),Y(t)\)。只需按 \(\mu(t)\) 更新即可：

若 \(\mu(t)=1\)：\(p\) 增 \(1\)，所以 \(X\) 增 \(100\)，\(Y\) 减 \(99\)
若 \(\mu(t)=-1\)：\(n\) 增 \(1\)，所以 \(X\) 减 \(99\)，\(Y\) 增 \(100\)
若 \(\mu(t)=0\)：都不变

因此扫描 \(t=1\ldots N\) 就能递推得到整个步行轨迹，并顺便得到其最小/最大值用于坐标平移。

代码

#include <bits/stdc++.h>
using namespace std;

// 0-index Fenwick (BIT) for counts
struct Fenwick
{
    int n;
    vector<uint32_t> bit; // counts fit in 32-bit (<= N+1)
    Fenwick(int n_ = 0) : n(n_), bit(n_, 0) {}
    inline void add(int i, uint32_t v = 1)
    {
        while (i < n)
        {
            bit[i] += v;
            i |= i + 1;
        }
    }
    inline uint64_t sumPrefix(int i) const
    { // sum [0..i]
        uint64_t s = 0;
        while (i >= 0)
        {
            s += bit[i];
            i = (i & (i + 1)) - 1;
        }
        return s;
    }
};

static const int N = 20000000;

// Linear sieve for Möbius mu[0..N], mu in {-1,0,1}
static vector<int8_t> mobius_sieve(int n)
{
    vector<int8_t> mu(n + 1, 0);
    vector<int> lp(n + 1, 0);
    vector<int> primes;
    primes.reserve(n / 10);

    mu[0] = 0;
    if (n >= 1)
        mu[1] = 1;

    for (int i = 2; i <= n; ++i)
    {
        if (lp[i] == 0)
        {
            lp[i] = i;
            primes.push_back(i);
            mu[i] = -1;
        }
        int li = lp[i];
        int8_t mi = mu[i];
        for (int p : primes)
        {
            long long ip = 1LL * i * p;
            if (ip > n)
                break;
            lp[(int)ip] = p;
            if (p == li)
            {
                mu[(int)ip] = 0;
                break;
            }
            else
            {
                mu[(int)ip] = (int8_t)(-mi);
            }
        }
    }
    return mu;
}

// Walk to get min/max of prefix with step sizes for mu==1 and mu==-1
static pair<int, int> walk_range(const vector<int8_t> &mu, int n, int d_pos, int d_neg)
{
    int x = 0, mn = 0, mx = 0;
    for (int t = 1; t <= n; ++t)
    {
        int8_t m = mu[t];
        if (m == 1)
            x += d_pos;
        else if (m == -1)
            x += d_neg;
        mn = min(mn, x);
        mx = max(mx, x);
    }
    return {mn, mx};
}

// Count pairs (i,j), 0<=i<j<=n, X(i) <= X(j)
static uint64_t count_prefix_pairs(const vector<int8_t> &mu, int n, int d_pos, int d_neg)
{
    auto [mn, mx] = walk_range(mu, n, d_pos, d_neg);
    int size = mx - mn + 1;

    Fenwick fw(size);

    int x = 0;
    uint64_t ans = 0;

    // prefixes t=0..n
    for (int t = 0; t <= n; ++t)
    {
        int idx = x - mn;
        ans += fw.sumPrefix(idx); // count previous i < t with X(i) <= X(t)
        fw.add(idx, 1);

        if (t == n)
            break;
        int8_t m = mu[t + 1];
        if (m == 1)
            x += d_pos;
        else if (m == -1)
            x += d_neg;
    }
    return ans;
}

int main()
{
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    auto mu = mobius_sieve(N);

    // C1: 99*N(a,b) <= 100*P(a,b)  <=>  X(b) >= X(a-1)
    // X step: mu==1 => +100, mu==-1 => -99
    uint64_t C1 = count_prefix_pairs(mu, N, +100, -99);

    // C2: 99*P(a,b) <= 100*N(a,b)  <=>  Y(b) >= Y(a-1)
    // Y step: mu==1 => -99, mu==-1 => +100
    uint64_t C2 = count_prefix_pairs(mu, N, -99, +100);

    uint64_t total = 1ULL * N * (N + 1) / 2;
    uint64_t ans = C1 + C2 - total;

    cout << ans << "\n";
    return 0;
}