bzoj2693: jzptab

链接

　　http://www.lydsy.com/JudgeOnline/problem.php?id=2693

题解

　　网上的题解都和popoqqq长得一样，我来发篇不太一样的。
　　令
　　

g(x)=x(x+1)2

　　

s(n,m,x)=x2g(⌊nx⌋)g(⌊mx⌋)

　　沿用上一道题中的结论
　　

ans=∑d=1n1d∑d|xnμ(xd)s(n,m,x)

　　来变一下形
　　

∑x=1ns(n,m,x)∑d|x1dμ(xd)　　=∑x=1ng(⌊nx⌋)g(⌊mx⌋)x2∑d|x1dμ(xd)

　　令inv(x)=1x易知其为积性函数，用∗表示卷积
　　

ans=∑x=1ng(⌊nx⌋)g(⌊mx⌋)x2(inv∗μ)(x)　　=∑x=1ng(⌊nx⌋)g(⌊mx⌋)x2(μ∗inv)(x)　　=∑x=1ng(⌊nx⌋)g(⌊mx⌋)x2∑d|xμ(d)dx　　=∑x=1ng(⌊nx⌋)g(⌊mx⌋)x∑d|xμ(d)d

　　令h(d)=μ(d)d，显然h(d)也是积性函数
　　

ans=∑x=1ng(⌊nx⌋)g(⌊mx⌋)x∑d|xh(d)　　=∑x=1ng(⌊nx⌋)g(⌊mx⌋)x(h∗1)(x)

　　令f(x)=(h∗1)(x)，显然积性函数的卷积还是积性函数，所以f(x)可以线性筛。当i mod primej=0时，可以通过对f(i)乘primej的完全积性函数，再加上i×primej的非完全积性函数来求。
　　最后每次询问直接分块O(N−−√)

代码

//线性筛+莫比乌斯反演#include <cstdio>#include <algorithm>#define maxn 10000001#define ll long long#define mod 100000009llusing namespace std;ll f[maxn], s2[maxn];int mu[maxn], prime[maxn/10], inv[maxn];bool mark[maxn];void shai(){    int i, j;    for(inv[1]=1,i=2;i<maxn;i++)inv[i]=(ll)(mod-mod/i)*(inv[mod%i])%mod;    mu[1]=f[1]=1;    for(i=2;i<maxn;i++)    {        f[i]%=mod;        if(!mark[i]){prime[++prime[0]]=i;mu[i]=-1;f[i]=inv[i]+mu[i];}        for(j=1;j<=prime[0] and i*prime[j]<maxn;j++)        {            mark[i*prime[j]]=1;            if(i%prime[j]==0)            {                mu[i*prime[j]]=0;                f[i*prime[j]]=f[i]*inv[prime[j]];                break;            }            mu[i*prime[j]]=-mu[i];            f[i*prime[j]]=f[i]*f[prime[j]];        }    }    for(i=1;i<maxn;i++)f[i]=(f[i]*i*i+f[i-1])%mod;}inline ll g(ll x){return x*(x+1)/2%mod;}inline ll s(ll n, ll m, ll x){return x*x%mod*g(n/x)%mod*g(m/x)%mod;}void solve(int N, int  M){    if(N>M)swap(N,M);    int last, i; ll ans=0;    for(i=1;i<=N;i=last+1)    {        last=min(N/(N/i),M/(M/i));        ans=(ans+g(N/i)%mod*g(M/i)%mod*(f[last]-f[i-1]))%mod;    }    printf("%lld\n",(ans+mod)%mod);}int main(){    int N, M, T;    shai();    for(scanf("%d",&T);T;T--)    {        scanf("%d%d",&N,&M);        solve(N,M);    }    return 0;}