[BZOJ 2693]jzptab 莫比乌斯反演

这道题什么玩意啊!鬼畜啊!内网OD标程4900ms而他时限开了5s!!

卡了一个多小时常,什么鬼东西啊!!

下面开始正常推导:

$$ans=\sum_{i=1}^n \sum_{j=1}^m lcm(i,j)$$

$$=\sum_{i=1}^n \sum_{j=1}^m \frac{i*j}{gcd(i,j)}$$

枚举$d=gcd$得

$$ans=\sum_{d=1}^n d  \sum_{i=1}^{\lfloor \frac{n}{d} \rfloor} \sum_{j=1}^{\lfloor \frac{m}{d} \rfloor} [gcd(i,j)==1]\frac{i*d*j*d}{d}$$

即算出除以$d$后互质的对数,两个数都乘$d$就得到了两个数的乘积,在除以$d$就是他们的最小公倍数。

然后把$d$提到前面来,后边的东西可以反演得到

$$ans=\sum_{d=1}^n d \sum_{t=1}^{\lfloor \frac{n}{d} \rfloor} \mu(t)  \sum_{i=1}^{\lfloor \frac{n}{td} \rfloor} \sum_{j=1}^{\lfloor \frac{m}{td} \rfloor} i*t*j*t$$

$$=\sum_{d=1}^n d \sum_{t=1}^{\lfloor \frac{n}{d} \rfloor} \mu(t) t^2  \sum_{i=1}^{\lfloor \frac{n}{td} \rfloor}i \sum_{j=1}^{\lfloor \frac{m}{td} \rfloor}j$$

设$$sum(x,y)=\sum_{i=1}^x\sum_{j=1}^yi*j=\frac{x*(x+1)}{2}*\frac{y*(y+1)}{2}$$

于是后边的$ \sum_{i=1}^{\lfloor \frac{n}{td} \rfloor}i \sum_{j=1}^{\lfloor \frac{m}{td} \rfloor}j$是可以$O(1)$求出来的

于是

$$ans=\sum_{d=1}^n d \sum_{t=1}^{\lfloor \frac{n}{d} \rfloor} \mu(t) t^2  sum({\lfloor \frac{n}{td} \rfloor},{\lfloor \frac{m}{td} \rfloor})$$

这样就可以$O(n)$求出来啦!可以过掉 $BZOJ 2154 Crash$的数字表格 ,但是这道题不行,需要再推导

令$D=kd$,稍加变形可得

$$ans=\sum_{D=1}^{min(n,m)}sum(\lfloor \frac{n}{D}\rfloor,\lfloor\frac{m}{D}\rfloor)\sum_{i|D}\mu(i)*i^2*\frac Di$$

然后就变成喜闻乐见的形式啦!

$$g(x)=\sum_{i|D}\mu(i)*i^2*\frac Di$$

我们要线性筛出所有$g(i)$

还是那个套路,分类

当$i$为质数时

$g(p)=\mu(1)*p+\mu(p)*p*p=p-p^2$

当$i$与$p$互质时,把约数分成含有$p$的和不含$p$的

$$g(i*p)=\sum_{d|i}\mu(d)*d*i*p+\sum_{d*p|i*p}\mu(d*p)*d*p*i*p$$

$$=p*g(i)-p^2*g(i)$$

$$=g(i)*g(p)$$

可见它是一个积性函数

当$i$与$p$不互质时

$i=i_1*p^y$

那么我们可以把$i*p$的约数分为含有$p^{y+1}$的和其他次方的

由于$y+1>=2$,所以含有$p^{y+1}$的约数的$\mu$值为$0$,所以对答案带来的贡献为$0$

所以

$$g(i*p)=\sum_{d|i}\mu(d)*d*i*p$$

$$=g(i)*p$$

于是就可以$\sqrt{n}$搞出来啦!

#include <cstdio>
#include <cstring>
using namespace std;
#define N 10000010
#define LL long long
#define mod 100000009
#define pos(i,a,b) for(LL i=(a);i<=(b);i++)
char B[1<<15],*S=B,*T=B;
#define getc (S==T&&(T=(S=B)+fread(B,1,1<<15,stdin),S==T)?0:*S++)
inline int read()
{
    int x=0;register char c=getc;
    while(c<'0'||c>'9')c=getc;
    while(c>='0'&&c<='9')x=10*x+(c^48),c=getc;
    return x;
}//快读 
int prime[N/10],notprime[N];
LL g[N];
inline int min(int a,int b){return a<b?a:b;}
inline LL gsum(LL a,LL b){
	return ( (a*(a+1)>>1)%mod )*( (b*(b+1)>>1)%mod )%mod;
}
LL ans;
int t,n,m;
int main()
{
	g[1]=1;
    pos(i,2,N-10){
        if(!notprime[i])
            prime[++prime[0]]=i,g[i]=(i-(LL)i*i)%mod;
        for(int j=1;j<=prime[0]&&(i*prime[j])<=N-10;++j){
            notprime[i*prime[j]]=1;
            if(i%prime[j]==0){g[i*prime[j]]=g[i]*prime[j]%mod;break;}
            g[i*prime[j]]=g[i]*g[prime[j]]%mod;
        }
    }
    pos(i,1,N-10) g[i]=(g[i-1]+g[i])%mod;
    t=read();
    while(t--){
        n=read(),m=read();
        if(n>m)n^=m,m^=n,n^=m;
        int last=0;ans=0;
        for(int i=1;i<=n;i=last+1)
            last=min(n/(n/i),m/(m/i)),ans=ans+gsum(n/i,m/i)*(g[last]-g[i-1])%mod;
        printf("%lld\n",(ans%mod+mod)%mod);
    }
}

  

posted @ 2017-12-10 18:02 Hallmeow 阅读(...) 评论(...) 编辑 收藏