BZOJ2154/BZOJ2693/Luogu1829 Crash的数字表格/JZPFAR 莫比乌斯反演

BZOJ2154/BZOJ2693/Luogu1829 Crash的数字表格/JZPFAR 莫比乌斯反演

传送门——Luogu

传送门——BZOJ2154

BZOJ2693是权限题

其中JZPFAR是多组询问，Crash的数字表格是单组询问

先推式子(默认(N leq M)，所有分数下取整)

(begin{align*} sumlimits_{i=1}^N sumlimits_{j=1}^M lcm(i,j) & = sumlimits_{i=1}^N sumlimits_{j=1}^M frac{ij}{gcd(i,j)} \ & = sumlimits_{d=1}^N dsumlimits_{i=1}^frac{N}{d} sumlimits_{j=1}^frac{M}{d} ij[gcd(i,j) == 1] \ & = sumlimits_{d=1}^N dsumlimits_{i=1}^frac{N}{d} sumlimits_{j=1}^frac{M}{d} ij sumlimits_{p mid gcd(i,j)} mu(p) \ & = sumlimits_{d=1}^N d sumlimits_{p=1}^frac{N}{d} p^2 mu(p) sumlimits_{i=1}^frac{N}{dp} sumlimits_{j=1}^frac{M}{dp} ij \ & = sumlimits_{T=1}^N (sumlimits_{i=1}^frac{N}{T}sumlimits_{j=1}^frac{M}{T} ij) sumlimits_{p | T} p^2 times frac{T}{p} times mu(p) end{align*})

推到这里开始做

首先(sumlimits_{i=1}^frac{N}{T} sumlimits_{j=1}^frac{M}{T} ij = frac{frac{N}{T}(frac{N}{T} + 1) times frac{M}{T}(frac{M}{T} + 1)}{4})，可以数论分块，那么要算(sumlimits_{p | T} p^2 times frac{T}{p} times mu(p) = T sumlimits_{p|T} p mu(p))的前缀和

首先(sumlimits_{p|T} p mu(p) = (id · mu) * I)是一个积性函数，而(N leq 10^7)还没有大到用杜教筛的程度，考虑线性筛。（PS：亲测直接写枚举倍数在Luogu的神机上是能跑过的）

设(f(i) = sumlimits_{p | i} p mu (p))，首先(f(1) = 1 , f(p)(p in Prime) = (1-p) )。现在假设已计算出了(f(i))，要计算(f(i times j))，其中(j in Prime)。

若(i = j^k times x(k,x geq 1 && x notmid j))，那么相比于(i)，(i times j)新产生的因数都可写成(j^{k+1} times y(y mid x))的形式，而(mu (j^{k+1} times y) = 0)，所以(f(i times j) = f(i) )

否则，对于(i)的一个因数(x(mu(x) neq 0))，在(i times j)中对应两个因数(x)和(x times j)，它们会在(f(i times j))中产生(x times j times (-mu(x)) + x times mu(x) = x times mu(x) times (1 - j))的贡献。

所以(f(i times j) = sumlimits_{p | i}p mu(p) times (1 - j) = (1-j)f(i) = f(i)f(j))

根据上面的式子线性筛即可。总复杂度(O(n + tsqrt{n}))

#include<bits/stdc++.h>
//This code is written by Itst
using namespace std;inline int read(){int a = 0;char c = getchar();bool f = 0;while(!isdigit(c)){if(c == '-')f = 1;c = getchar();}while(isdigit(c)){a = (a << 3) + (a << 1) + (c ^ '0');c = getchar();}return f ? -a : a;
}const int MOD = 20101009 , MAXN = 1e7 + 7;
int prime[MAXN] , xs[MAXN];
bool nprime[MAXN];
int cnt , N , M;void init(){xs[1] = 1;for(int i = 2 ; i <= N ; ++i){if(!nprime[i]){prime[++cnt] = i;xs[i] = MOD - i + 1;}for(int j = 1 ; j <= cnt && prime[j] * i <= N ; ++j){nprime[prime[j] * i] = 1;if(i % prime[j] == 0){xs[i * prime[j]] = xs[i];break;}xs[i * prime[j]] = 1ll * xs[i] * xs[prime[j]] % MOD;}}for(int i = 1 ; i <= N ; ++i)xs[i] = (1ll * xs[i] * i + xs[i - 1]) % MOD;
}int main(){N = read();M = read();if(N > M)swap(N , M);init();int sum = 0;for(int i = 1 , pi ; i <= N ; i = pi + 1){pi = min(N / (N / i) , M / (M / i));sum = (sum + (xs[pi] - xs[i - 1] + MOD) % MOD * (1ll * (N / i) * (N / i + 1) / 2 % MOD) % MOD * (1ll * (M / i) * (M / i + 1) / 2 % MOD)) % MOD;}cout << sum;return 0;
}

转载于:.html