bzoj 4176: Lucas的数论 （反演）

题目描述

传送门

题目大意：
设f(ij)表示i*j的约束个数，求

∑i=1n∑j=1nf(ij) mod 1e9+7

n<=1e9

题解

要化简上面的式子，就必须科学的表示出f(ij)
如果你做过SDOI的约数个数和，那么就应该知道

f(nm)=∑i|n∑j|m[gcd(i,j)=1]

这个怎么证？其实可以感受一下。
当j=1的时候，我们加入了n的所有约数
当i=1的时候，我们加入了m的所有约数
但是n,m的约数可能是有交集的，但是交集中的元素都可以用一个之前未出现的约数替换，并且在后面计算中这个约数不会再加入。比如说n有一个质因子为pi指数为x,m也有pi这个质因子指数为y，那么m中的pi^{1..y}其实都可以换成pi^{x+1…y}，其他的也是同样的道理。
剩下的数当且仅当gcd(i,j)=1，可以形成新的约数。总之如果i,j不互质，那么一定可以用一对互质的i’,j’替换掉。

那么利用上面的式子，我们进行化简

∑a=1n∑b=1n∑i|a∑j|b[gcd(i,j)=1]

∑i=1n∑j=1n⌊ni⌋⌊nj⌋[gcd(i,j)=1]

∑i=1n⌊ni⌋∑j=1n⌊nj⌋∑d|(i,j)μ(d)

∑d=1nμ(d)∑i=1⌊ni⌋⌊ni∗d⌋∑j=1⌊nj⌋⌊nj∗d⌋

∑d=1nμ(d)(∑i=1⌊ni⌋⌊ni∗d⌋)2

h(x)=∑xi=1⌊xi⌋ 可以在O(x√)的时间内求解，如果我们可以预处理μ 的前缀和，那么就可以在O(n34)的时间内求解。
n的范围是1e9，所以肯定不能O(n)的预处理，而且我们也做不到处理处所有数的前缀和，所以我们只考虑对我们有贡献的数。
这里的话应该可以想到μ要用杜教筛。
设M(n)=∑i=1nμ(i)
1=∑i=1n[i=1]=∑i=1n∑d|iμ(d)=∑i=1n∑d=1⌊ni⌋μ(d)=∑i=1nM(⌊ni⌋)
M(n)=1−∑i=2nM(⌊ni⌋) 对于这个式子我们只要预处理出n34的前缀和，然后计算n√个M就可以了，总的时间复杂度应该是O(n34)，与计算上面的式子是一样的。
这是要算一个M的时候，那么我们可不可以预处理出一些关键的值啊，其实是可以的。
考虑哪些值对我们的计算有用处，n/(n/i)如果这个值在n34内我们可以直接调用我们预处理的前缀和，那么如果大于的话，我们找到所有的n/t>n34，然后从大到小枚举t，计算M(n/t),这样每次计算的时候需要用到的值都已经预处理好了，直接O(n/t−−−√)的计算即可。

代码

#include<iostream>#include<cstdio>#include<cstring>#include<cmath>#include<algorithm>#define N 10000000#define LL long long #define p 1000000007using namespace std;int n,n1,pd[N],prime[N];LL mu[N],sum[N];void init(){    mu[1]=1;    for (int i=2;i<=n1;i++) {        if (!pd[i]) {            prime[++prime[0]]=i;            mu[i]=-1;        }        for (int j=1;j<=prime[0];j++) {            if (i*prime[j]>n1) break;            pd[i*prime[j]]=1;            if (i%prime[j]==0) {                mu[i*prime[j]]=0;                break;            }            mu[i*prime[j]]=-mu[i];        }    }    for (int i=1;i<=n1;i++) mu[i]=mu[i-1]+mu[i];}LL calc(int n){    LL ans=0;    for (int i=1,j;i<=n;i=j+1) {        j=n/(n/i);        ans=(ans+(LL)(j-i+1)*(n/i))%p;    }    return ans*ans%p;}LL get_sum(int x){    if (x<=n1) return mu[x];    return sum[n/x];}void summu(){    int t;    for (t=1;n/t>n1;t++);    for (int k=t;k;k--) {        int m=n/k;        sum[k]=1;        for (int i=2,j;i<=m;i=j+1) {            j=m/(m/i);            sum[k]-=(LL)(j-i+1)*get_sum(m/i)%p;            sum[k]%=p;        }    }}int main(){    freopen("a.in","r",stdin);    freopen("my.out","w",stdout);    scanf("%d",&n); n1=ceil(pow(n,0.75));    init(); summu();    LL ans=0;    for (int i=1,j;i<=n;i=j+1) {        j=n/(n/i);        ans=(ans+(get_sum(j)-get_sum(i-1))*calc(n/i)%p)%p;    }    printf("%lld\n",(ans%p+p)%p);}