BZOJ 4174 - tty的求助

#数学 #莫比乌斯反演 #数论
Published on 2016-08-18

题目地址

描述

给定整数 n,m(n,m≤500000)n, m(n, m\le 500000)，给定实数 x(x≤100000)x(x\le 100000)，请你求：

∑i=1n∑j=1m∑0≤k<m⌊ik+xj⌋\sum_{i = 1}^n\sum_{j = 1}^m\sum_{0\le k< m}\left\lfloor\frac {ik + x} j\right\rfloor

分析

首先一看这个式子显然不是什么好求的东西，我们试图对 ∑0≤k<m⌊ik+xj⌋\sum_{0\le k< m}\left\lfloor\frac {ik + x} j\right\rfloor 求一个封闭形式，看看有什么发现：

中间的省略号略去了推导过程，关于详尽的推导过程，参见具体数学 P75。最关键的一步推导是这样的。对于 knmodmkn \bmod m 这个式子，让 k∈[0,m)k \in [0, m)，最终得到的数按照某种方式排列之后一定是 0,d,2d,...,m−d0, d, 2d, ..., m - d，后面跟着这个式子的 d−1d - 1 次重复，其中 d=gcd(n,m)d = \gcd(n, m)。所以第一项和第三项都可以使用等差数列求和来搞定。

于是这个封闭形式中，我们发现可以用莫比乌斯反演来求：

由于只有一组，所以可以不用继续推导了，复杂度为 O(n+n⋅n)=O(n)O(n + \sqrt n\cdot\sqrt n) = O(n)。
注意乘逆元！

代码

//  Created by Sengxian on 8/16/16.
//  Copyright (c) 2016年 Sengxian. All rights reserved.
//  BZOJ 4174 莫比乌斯反演 + 推公式
#include <bits/stdc++.h>
using namespace std;

typedef long long ll;
const int maxn = 500000 + 3, modu = 998244353;
int primes[maxn], mu[maxn], num;
bool isPrime[maxn];
ll s[maxn], ss[maxn];

void sieve() {
    fill(isPrime, isPrime + maxn, 1);
    mu[1] = 1, num = 0;
    for (int i = 2; i < maxn; ++i) {
        if (isPrime[i]) primes[num++] = i, mu[i] = -1;
        static int d;
        for (int j = 0; j < num && (d = i * primes[j]) < maxn; ++j) {
            isPrime[d] = false;
            if (i % primes[j] == 0) {
                mu[d] = 0;
                break;
            } else mu[d] = -mu[i];
        }
    }
    for (int i = 1; i < maxn; ++i) s[i] = s[i - 1] + mu[i];
}

inline ll S(ll n) {
    return ((n + 1) * n >> 1) % modu;
}

inline ll G(ll n, ll m) {
    ll ans = 0;
    for (int i = 1, j = 1; i <= n; i = j + 1) {
        j = min(n /    (n / i), m / (m / i));
        (ans += (ll)(s[j] - s[i - 1]) * (n / i) % modu * (m / i) % modu) %= modu;
    }
    return ans;
}

inline int F(ll n, ll m, double x) {
    if (n > m) swap(n, m);
    ss[0] = 0;
    for (int i = 1; i <= n; ++i) ss[i] = 2 * i * (ll)floor(x / i) + i, (ss[i] += ss[i - 1]) %= modu;
    ll ans = (S(n - 1) * S(m - 1) % modu - n * m) % modu, tmp = 0;
    for (int i = 1, j = 1; i <= n; i = j + 1) {
        j = min(n / (n / i), m / (m / i));
        tmp = G(n / i, m / i);
        (ans += (ss[j] - ss[i - 1]) * tmp % modu) %= modu;
    }
    ans *= 499122177LL;
    return (ans % modu + modu) % modu;
}

int main() {
    sieve();
    int n, m; double x;
    scanf("%d%d%lf", &n, &m, &x);
    printf("%d\n", F(n, m, x));    
    return 0;
}
莫比乌斯反演简要笔记