1. 程式人生 > >【刷題記錄】BZOJ2154 crash的數字表格 莫比烏斯反演

【刷題記錄】BZOJ2154 crash的數字表格 莫比烏斯反演

治好了我多年的公式恐懼症

https://www.luogu.org/problemnew/show/P1829

(放洛谷題號的原因是BZOJ上很多內容是限制級的(霧))

題目大意:已知M,N,求ans =\sum_{i=1}^n\sum_{j=1}^m{\frac {ij} {(i,j)} }

(我採用數論教材上的方法表示最大公因數,即(i,j)表示i,j的最大公約數)

那麼,我們開始吧。

不妨設n<m

d = (i,j),有ans = \sum _{d=1}^{n} \sum_{\begin{matrix} (i,j)=d \\ 1\leqslant i\leqslant n, \\ 1\leqslant j \leqslant m, \end{matrix}}\frac {ij} d

ans = \sum _{d=1}^{n} d\sum_{(i,j)=1}^{1\leq i\leq [\frac nd],1\leq j\leq [\frac md]}{ij}

接下來是莫比烏斯函式的性質:\sum _{x|d} \mu (x) = 1當且僅當d = 1

ans = \sum _{d=1}^{n} d\sum_{i=1}^{[\frac nd]}\sum_{j=1}^{[\frac md]}\sum_{x|(i,j)}\mu (x){ij}

轉化為列舉x

ans = \sum_{d=1}^{ n}d\sum_{x=1}^{[\frac nd]}\mu(x)x^2\sum_{j=1}^{[\frac n {dx}]}\sum_{i=1}^{[\frac m {dx}]}ij

列舉t=dx

ans = \sum_{t=1}^{ n}t\sum_{i=1}^{[\frac nt]}\sum_{j=1}^{[\frac mt]}ij\sum_{x|t}\mu(x)x=\sum_{t=1}^{ n}t(\sum_{i=1}^{[\frac nt]}i)(\sum_{j=1}^{[\frac mt]}j)\sum_{x|t}\mu(x)x

中間兩個和式為等差數列求和,可以O(1)得到結果

 

經過完全歸納法證明,f(i)=\sum_{x|i}\mu(x)x為積性函式

證明如下:

當i為質數時,顯然f(i)=1-i

記p為任意質數

(i,p)=1

f(ip)=\sum_{x|ip}x\mu(x)=f(i)+\sum _{xp|ip}(xp)\mu(xp)=f(i)-pf(i)=f(p)*f(i)

證完

接下來考慮線性篩的實現

考慮當(i,p)=p

f(ip)=\sum_{x|ip}x\mu(x)=f(i)+\sum _{xp|ip}(xp)\mu(xp)=f(i)-0=f(i)

所以線性篩完成

依靠線性篩預處理後,可以O(n)時間預處理,再用O(n)時間得出結果

然而此題還可以再次優化

考慮到[\frac nt][\frac mt]可以數論分塊

時間複雜度可以優化至O(n+sqrt(n))

這樣可以解決多組資料的情況

程式碼

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <algorithm>
#include <cstring>
using namespace std;
int n,m;
#define ll long long
#define mod 20101009
#define inv2 10050505
#define N 10000000
int prime[N+9],mu[N+9];
bool notprime[N+8];
ll sum(int x)
{
    return (((((ll)x*(x+1))%mod)*inv2)%mod);
}
int main()
{
    scanf("%d%d",&n,&m);
    if(n>m)swap(n,m);
    notprime[1]=mu[1]=1;int tot = 0;
    for(int i = 2; i <= n; i ++)
    {
        if(!notprime[i])
        {
            prime[++tot]=i;
            mu[i]=(1-i+mod)%mod;
        }
        for(int j = 1; j <= tot &&prime[j]*i<=n; j ++)
        {
            notprime[prime[j]*i]=1;
            if(i%prime[j]==0)
            {
                mu[i*prime[j]]=mu[i];
                break;
            }
            else 
            {
                mu[i*prime[j]]=(((ll)mu[i]*(1-prime[j])%mod)+mod)%mod;
            }
        }
        mu[i]=((((ll)i*mu[i])%mod+mu[i-1])%mod+mod)%mod;
    }
    int ans = 0;
    for(int i = 1,last; i <= n; i =last +1)
    {
        last = min(n/(n/i),m/(m/i));
        ans = (((((sum(n/i)*sum(m/i))%mod)*(mu[last]-mu[i-1]))%mod+mod)%mod+ans)%mod;
    }
    printf("%d",ans);
}