【莫比乌斯反演】BZOJ2154Crash的数字表格&BZOJ2693jzptab

首先是BZOJ2154

题目大意

给出n,m求

∑i=1n∑j=1mlcm(i,j)

满足（n,m<=107）

分析

因为比起lcm，我们更常用到gcd，因此首先将题目转换为求gcd的形式

ans=∑i=1n∑j=1mlcm(i,j)=∑i=1n∑j=1mi∗jgcd(i,j)

接下来，一扯到gcd，本能驱使我们这么设（说白了就是做多了的经验）

设f(n,m,k)=∑i=1n∑j=1m[(i,j)=k]i∗j

接下来才开始推（以上基本没什么营养）

ans=∑d=1min(n,m)f(n,m,d)d=∑d=1min(n,m)d2f(nd,md,1)d=∑d=1min(n,m)d∗f(nd,md,1)

推到这里，不难想到求出f函数，再分块求和，求和方法类似bzoj2301problem b
接下来的问题就是求出f

设F(n,m,d)=∑i=1n∑j=1m[d|(i,j)]i∗j=∑i=1nd∑j=1mdd2∗i∗j=d2(n/d)∗(n/d+1)2(m/d)∗(m/d+1)2=∑d|kf(n,m,k)

所以

F(n,m,d)=∑d|kf(n,m,k)

这就是典型的可以做莫比乌斯反演的形式，那么。。反演呗

f(n,m,k)=∑k|dμ(dk)F(n,m,d)

因为我们要求的k=1，所以

f(n,m,1)=∑d=1min(n,m)μ(d)F(n,m,d)=∑d=1min(n,m)μ(d)d2(n/d)∗(n/d+1)2(m/d)∗(m/d+1)2

因此再用老方法，分块，n−√内求出f
因此总复杂度为O(n−√∗n−√)=O(n)

#include<cstdio>#include<cstring>#include<algorithm>#include<cmath>#include<vector>#define SF scanf#define PF printf#define MAXN 10000010int MAXPRIME;#define MOD 20101009using namespace std;vector<long long> primes;int isprime[MAXN],mu[MAXN];long long sum[MAXN],n,m;void prepare(){    mu[1]=1;    for(long long i=2;i<=MAXPRIME;i++){        if(isprime[i]==0){            mu[i]=-1;            primes.push_back(i);        }        for(int j=0;j<primes.size()&&i*primes[j]<=MAXPRIME;j++){            isprime[i*primes[j]]=1;            if(i%primes[j]==0){                mu[i*primes[j]]=0;                break;            }            mu[i*primes[j]]=-mu[i];        }    }    for(long long i=1;i<=MAXPRIME;i++)        sum[i]=(1ll*mu[i]*i*i+sum[i-1]+MOD)%MOD;}long long que(long long x,long long y){    long long res=0,last;    for(long long i=1;i<=min(x,y);i=last+1){        last=min(x/(x/i),y/(y/i));        res=(res+((sum[last]-sum[i-1])%MOD*((x/i)*((x/i)+1)/2%MOD*((y/i)*((y/i)+1)/2%MOD)%MOD))%MOD+MOD)%MOD;        res%=MOD;    }    return res;}int main(){    //freopen("data.in","r",stdin);    //freopen("data.out","w",stdout);    SF("%lld%lld",&n,&m);    MAXPRIME=min(n,m);    prepare();    long long ans=0,last;    for(long long i=1;i<=min(n,m);i=last+1){        last=min(n/(n/i),m/(m/i));        ans=(ans+(i+last)*(last-i+1)/2%MOD*que(n/i,m/i))%MOD;        //PF("%I64d\n",ans);    }    PF("%lld",ans);}

以下为BZOJ2693jzptap的内容

题目大意

与上题完全相同，但改为多组数据，最多10000组

分析

10000组数据，明显上题O(n)的复杂度是不可能的了
因此考虑优化，为了方便起见，将答案的转移试整理一下。

ans=∑d=1min(n,m)d∗∑k=1min(nd,md)μ(k)∗k2(n/dk)∗(n/dk+1)2(m/dk)∗(m/dk+1)2

BZOJ2820的经验告诉我们，当有多个枚举量，且这多个枚举量多以乘积的形式出现在公式中，一种比较常用的优化方式，就是枚举这几个枚举量的乘积，即:

ans=∑t=1min(n,m)(n/t)∗(n/t+1)2(m/t)∗(m/t+1)2∗∑d=1min(n/t,m/t)td∗μ(d)∗d2

设S(t)=∑d|ttd∗μ(d)∗d2

因此，只要预处理出S函数的的前缀和，那么就可以在O(tn−√)内求解，可以接受。
接下来就要用到一直以来从未用过的积性函数的特征。
因为积性函数的因子和倍数都是积性函数，所以这个函数必定为积性函数。那么就可以用线性筛的方式快速求出函数。
线性筛的特征是：任何一个数都被它最小的质因子筛到，因此我们可以利用这一点求值（类似于筛莫比乌斯函数）设一个数i,其最小质因子为primej
若i  mod  primej≠0即i与primej互质，那么S(i∗primej)=S(i)∗S(primej)
若i  mod  primej=0即i唯一分解后primej的次数大于1,
这种情况积性函数不能直接求值，于是观察表达式，
不难发现，在这种情况下，μ(i∗primej)的值为0，
所以相当于这一部分对答案无贡献，所以S(i)的值与S(i∗primej)的值必定成倍数关系，
于是可以表示为为S(i∗primej)=S(i)∗primej
至此，本题就可以解决了。
说说一些题外话，相比上题，本题我的代码不多出15行，更重要的是还是1A，但上题我却因为取模的问题WA了近十次…

#include<cstdio>#include<cstring>#include<algorithm>#include<cmath>#include<vector>#define SF scanf#define PF printf#define MAXN 10000010#define MAXPRIME  10000000#define MOD 100000009using namespace std;int mu[MAXN],isprime[MAXN];vector<int> primes;long long sum[MAXN];void prepare(){    mu[1]=1;    sum[1]=1;    for(long long i=2;i<=MAXPRIME;i++){        if(isprime[i]==0){            primes.push_back(i);            mu[i]=-1;            sum[i]=-i*i+i;        }        sum[i]%=MOD;        //PF("(%I64d)\n",sum[i]);        for(int j=0;j<primes.size()&&i*primes[j]<=MAXPRIME;j++){            isprime[primes[j]*i]=1;            if(i%primes[j]==0){                sum[i*primes[j]]=(primes[j]*sum[i])%MOD;                break;            }            mu[i*primes[j]]=-mu[i];            sum[i*primes[j]]=(sum[i]*sum[primes[j]])%MOD;        }    }    for(int i=1;i<=MAXPRIME;i++){        sum[i]+=sum[i-1];        sum[i]%=MOD;       }}int t;long long n,m;long long s(long long x){    return ((x*(x+1))/2)%MOD;}int main(){    prepare();    SF("%d",&t);    while(t--){        SF("%I64d%I64d",&n,&m);        long long last=0,res=0;        for(long long i=1;i<=min(n,m);i=last+1){            last=min(n/(n/i),m/(m/i));            res+=((s(n/i)*s(m/i))%MOD*(sum[last]-sum[i-1]+MOD+MOD))%MOD;            //PF("(%lld %lld)",last,res);            res%=MOD;        }        PF("%lld\n",res);    }}