[bzoj2154][积性函数]Crash的数字表格

2154: Crash的数字表格

Time Limit: 20 Sec Memory Limit: 259 MB
Description

今天的数学课上，Crash小朋友学习了最小公倍数(Least Common Multiple)。对于两个正整数a和b，LCM(a, b)表示能同时被a和b整除的最小正整数。例如，LCM(6, 8) = 24。回到家后，Crash还在想着课上学的东西，为了研究最小公倍数，他画了一张N*M的表格。每个格子里写了一个数字，其中第i行第j列的那个格子里写着数为LCM(i, j)。一个4*5的表格如下： 1 2 3 4 5 2 2 6 4 10 3 6 3 12 15 4 4 12 4 20 看着这个表格，Crash想到了很多可以思考的问题。不过他最想解决的问题却是一个十分简单的问题：这个表格中所有数的和是多少。当N和M很大时，Crash就束手无策了，因此他找到了聪明的你用程序帮他解决这个问题。由于最终结果可能会很大，Crash只想知道表格里所有数的和mod 20101009的值。

Input

输入的第一行包含两个正整数，分别表示N和M。

Output

输出一个正整数，表示表格中所有数的和mod 20101009的值。

Sample Input

4 5

Sample Output

122

【数据规模和约定】

100%的数据满足N, M ≤ 10^7。

sol1.0

由于特殊原因，回头又做了一次这道题。因此题解有2份，第一份是当初很naive的时候膜行止出来的。第二份是后来做的。
这里我设n < m,不然就swap
题目要求
$\sum_{i=1}^n\sum_{j=1}^mlcm(i,j)$

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

发现这题的贡献和i，j都有关，并且和gcd有关系。那么我们发现难点主要在gcd，所以我们把gcd拆开做。
那么我们设
f(n,m,d)=$\sum_{i=1}^n\sum_{j=1}^mi*j*[gcd(i,j)=d]$

这么设的优点在哪里呢？
我们可以方便的表示出ans的式子
ans=$\sum_{d=1}^n\frac {f(n,m,d)}{d}$

=$\sum_{d=1}^n \frac {f(\frac {n}{d},\frac {m}{d},1)*d*d}{d}$

=$\sum_{d=1}^n f(\frac {n}{d},\frac {m}{d},1)*d$

所以我们现在的问题是怎么求f
根据经验，我们知道如果f中存在[xxx=d]的话，那么F中应当有一个[d|xxx]的状况，然后就能倍数和乱搞。
构造F(n,m,d)=$\sum_{i=1}^n\sum_{j=1}^m i*j*[d|gcd(i,j)]$使得我们能够反演

现在我们对F进行操作
F(n,m,d)=$\sum_{i=1}^{n/d}\sum_{j=1}^{m/d} i*j*d^2$

=$d^2*\sum_{i=1}^{n/d}\sum_{j=1}^{m/d} i*j$
通过这个操作，我们去掉了d|gcd(i,j)的条件

对于$\sum_{i=1}^{n/d}\sum_{j=1}^{m/d} i*j$

=$\sum_{i=1}^{n/d} \frac {(1+ \frac {m}{d})* \frac {m}{d}}{2}$

=$\frac {(1+ \frac {n}{d})* \frac {n}{d}}{2} \frac {(1+ \frac {m}{d})* \frac {m}{d}}{2}$

那么显然这个数据范围大F是不能出现O(n)求的情况，因为ans外面已经套了个$\sum$了，所以必须分段求和优化
设sum(n/d,m/d)=这个式子

所以F(n,m,d)的式子就出来了

F(n,m,d)=$d^2*sum(n/d,m/d)$这就能分段求和了

那么可以真正上反演了
F(n,m,d)=$\sum_{d|n} f(n)$
f(n,m,d)=$\sum_{d|n} \mu(n/d)F(d)$

f(n,m,1)=$\sum_{d=1}^{n} \mu(d)F(d)$

=f(n,m,1)=$\sum_{d=1}^{n} \mu(d)*d^2*sum(n/d,m/d)$

这里ans还是不能直接求，发现ans式子中f有除号。那这里分段求和一下就好了。

#include<iostream>
#include<cstring>
#include<string>
#include<cstdlib>
#include<cmath>
#include<algorithm>
#include<cstdio>
#define ll long long
using namespace std;

const int N=11000000;
const int pyz=20101009;

int n,m;
int prime[N],mu[N],sqri[N];
bool is[N];
inline int read()
{
    char c;
    bool flag=0;
    while((c=getchar())>'9'||c<'0') if(c=='-')flag=1;
    int res=c-'0';
    while((c=getchar())>='0'&&c<='9') res=(res<<3)+(res<<1)+c-'0';
    return flag?-res:res;
}
inline int sum(ll x,ll y)
{
    return (((x+1)*x/2)%pyz)*(y*(y+1)/2%pyz)%pyz;
}
inline ll F(int x,int y)
{
    int j=0;
    ll ans=0;
    for(int i=1;i<=x;i=j+1)
    {
        j=min(x/(x/i),y/(y/i));
        ans=(ans+(ll)(sqri[j]-sqri[i-1]+pyz)%pyz*sum((ll)x/i,(ll)y/i)%pyz)%pyz;
    }
    return ans;
}
int main()
{
//  freopen("2154.txt","r",stdin);
    n=read();m=read();
    if(n>m) swap(n,m);
    mu[1]=1;
    for(int i=2;i<=n;++i)
    {
        if(!is[i]) prime[++prime[0]]=i,mu[i]=-1;
        for(int j=1;j<=prime[0]&&i*prime[j]<=n;++j)
        {
            is[i*prime[j]]=1;
            if(!(i%prime[j]))
            {
                mu[i*prime[j]]=0;;
                break;
            }
            mu[i*prime[j]]=-mu[i];
        }
    }
    for(int i=1;i<=n;++i) sqri[i]=(sqri[i-1]+(ll)i*i%pyz*mu[i])%pyz; 
    int j=0;
    ll ans=0;
    for(int d=1;d<=n;d=j+1)
    {
        j=min(n/(n/d),m/(m/d));
        ans=(ans+(ll)(d+j)*(j-d+1)/2%pyz*F(n/d,m/d)%pyz)%pyz;
    }
    printf("%lld",ans);
}

sol2.0

在学习杜教筛，积性函数求和的时候，又看到了这道题，可以说是很怀念了，当初可把我做个半死

$\sum_i^n\sum_j^m \frac {ij}{(i,j)}$
$\sum_{d=1}^nd\sum_i^{n/d}\sum_j^{m/d}ij,[(i,j)==1]$
记
f(x,y)=$\sum_i^x\sum_j^yij,[(i,j)==1]$
$\sum_i^x\sum_j^yij\sum_{d|(i,j)}\mu(d)$
$\sum_d^x\mu(d)\sum_i^{x/d}\sum_j^{y/d}ij$
$\sum_d^x\mu(d)(1+y)y/2(1+x)x/2$
带回去就行了。分段求和即可。