bzoj 2693: jzptab （反演）

题目描述

传送门

题目大意：
$\sum\limits_{i=1}^n\sum\limits_{j=1}^m lcm(i,j) \space mod \space 100000009$

题解

哎，虽然式子很恶心，但是还是要化的 ,设n<=m
$\sum\limits_{i=1}^n\sum\limits_{j=1}^m lcm(i,j)$
$\sum\limits_{i=1}^n\sum\limits_{j=1}^m {i*j \over gcd(i,j)}$
$\sum\limits_{k=1}^n \sum\limits_{i=1}^n\sum\limits_{j=1}^m {i*j \over k}[gcd(i,j)=k]$
$\sum\limits_{k=1}^n {1\over k} \sum\limits_{i=1}^{\lfloor {n\over k}\rfloor }\sum\limits_{j=1}^{\lfloor {m\over k}\rfloor }{i*k*j*k} [gcd(i,j)=1]$
$\sum\limits_{k=1}^n {k} \sum\limits_{i=1}^{\lfloor {n\over k}\rfloor }\sum\limits_{j=1}^{\lfloor {m\over k}\rfloor }{i*j} \sum_{d|(i,j)}\mu(d)$
$\sum\limits_{k=1}^n {k} \sum_{d=1}^n \mu(d) \sum\limits_{i=1}^{\lfloor {n\over k}\rfloor } [d|i]*i\sum\limits_{j=1}^{\lfloor {m\over k}\rfloor } [d|j]*j$
$\sum\limits_{k=1}^n {k} \sum_{d=1}^n \mu(d) \sum\limits_{i=1}^{\lfloor {n\over {k*d}}\rfloor }d*i\sum\limits_{j=1}^{\lfloor {m\over {k*d}}\rfloor }d*j$
设$sum(x,y)=\sum\limits_{i=1}^x i\sum\limits_{j=1}^y j$ 这个式子直接用等差数列求和公式计算即可。
那么上面的式子就是

$$\sum\limits_{k=1}^n {k} \sum_{d=1}^n \mu(d) *d^2sum({\lfloor {n\over {k*d}}\rfloor },{\lfloor {m\over {k*d}}\rfloor})$$
然后设 $T=k*d$ ,转变枚举顺序
$$\sum\limits_{T=1}^n sum({\lfloor {n\over T}\rfloor },{\lfloor {m\over T}\rfloor})*\sum\limits_{d|n}\mu({T\over d})*{T\over d}^2*d$$
那么问题的关键就转换成了如何快速的求解 $h(n)=\sum\limits_{d|n}\mu({T\over d})*{T\over d}^2*d$
这个函数是一个积性函数，可以线性筛啊。a,b互质的时候，直接相乘。
a,b不互质的时候，直接乘上质数。 $h(pi^{xi})=\mu(pi)*pi^2*pi^{xi-1}+\mu(1)*pi^xi$

代码

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
#define LL long long 
#define p 100000009
#define N 10000000
using namespace std;
int n,m,T,pd[N+3],prime[N+3],g[N+3];
LL f[N+3],inv;
void init()
{
    f[1]=1; 
    for (int i=2;i<=N;i++) {
        if (!pd[i]) {
            prime[++prime[0]]=i;
            f[i]=(LL)(i-(LL)i*i%p)%p; 
        }
        for (int j=1;prime[j]*i<=N;j++) {
            pd[i*prime[j]]=1;
            if (i%prime[j]==0) {
                f[i*prime[j]]=(LL)f[i]*prime[j]%p;
                break;
            }
            f[i*prime[j]]=(LL)f[i]*f[prime[j]]%p;
        }
    }
    for (int i=1;i<=N;i++) f[i]=(f[i-1]+f[i])%p;
}
LL calc(LL x,LL y)
{
    LL t1=((x+1)*x>>1)%p; LL t2=((y+1)*y>>1)%p;
    return t1*t2%p;
}
int main()
{
    freopen("a.in","r",stdin);
    freopen("my.out","w",stdout);
    init(); scanf("%d",&T);
    inv=g[4];
    while (T--) {
        scanf("%d%d",&n,&m);
        if (n>m) swap(n,m);
        int j=0; LL ans=0;
        for (int i=1;i<=n;i=j+1) {
            j=min(n/(n/i),m/(m/i));
            ans=(ans+calc(n/i,m/i)*(f[j]-f[i-1])%p)%p;
        }
        printf("%lld\n",(ans%p+p)%p);
    }
}