题解 P1829 【[国家集训队]Crash的数字表格 / JZPTAB】

· · 题解

看到没有写杜教筛的十分开心.

n,m\leq 1000?

我会求lcm!

n,m\leq 10000000?

首先一定可以套上这里的做法.

不过还是来推一遍(跳跃性可能略大，详细步骤请看上面的链接).以下除法均向下取整.不失一般性，设n\leq m，那么

\begin{array}{lcc}\sum\limits_{i=1}^n\sum\limits_{j=1}^m lcm(i,j)=\sum\limits_{i=1}^n\sum\limits_{j=1}^m\dfrac{ij}{\gcd(i,j)}\\=\sum\limits_{d=1}^n\frac{1}{d}\sum\limits_{i=1}^n\sum\limits_{j=1}^mij[\gcd(i,j)=d]\\=\sum\limits_{d=1}^nd\sum\limits_{i=1}^{\frac{n}{d}}\sum\limits_{j=1}^\frac{m}{d}ij[gcd(i,j)=1]\\=\sum\limits_{d=1}^nd\sum\limits_{k}\mu(k)\sum\limits_{i=1}^\frac{n}{d}\sum\limits_{j=1}^\frac{m}{d}ij[k|i][k|j]\\=\sum\limits_{d=1}^nd\sum\limits_{k=1}^\frac{n}{d}\mu(k)k^2\sum\limits_{i=1}^\frac{n}{kd}\sum\limits_{j=1}^\frac{m}{kd}ij\end{array}

设S(n)=\sum\limits_{i=1}^ni=\dfrac{n(n+1)}{2}，那么可以写成

\sum\limits_{d=1}^nd\sum\limits_{k=1}^\frac{n}{d}\mu(k)k^2S(\frac{n}{kd})S(\frac{m}{kd})

如果我们再设一个F(n,m)=\sum\limits_{k=1}^n\mu(k)k^2S(\frac{n}{k})S(\frac{m}{k})，那么显然F(n,m)可以预处理\mu后在O(\sqrt{n})时间内数论分块算出来

原式可以写成\sum\limits_{d=1}^ndF(\frac{n}{d},\frac{m}{d})也可以数论分块.

总的复杂度其实是O(\sum\limits_{i=1}^{\sqrt{n}}\sqrt{i}+\sum\limits_{i=1}^{\sqrt{n}}\sqrt{\frac{n}{i}})=O(n^\frac{3}{4})

但是要预处理\mu，所以是O(n)的.

n,m\leq 1000000,T$组询问，$T\leq 10000?

回到设F之前的那个位置，换一个角度.

令T=kd，枚举T，那么一定有d|T

原式=\sum\limits_{T=1}^nS(\frac{n}{T})S(\frac{m}{T})\sum\limits_{d|T}d\mu(\frac{T}{d})\left(\frac{T}{d}\right)^2=\sum\limits_{T=1}^nS(\frac{n}{T})S(\frac{m}{T})\sum\limits_{d|T}Td\mu(d)

那么我们只需要预处理\sum\limits_{d|T}Td\mu(d)这个东西即可做到\sqrt{n}回答一个询问.

怎么预处理呢?显然可以直接O(n\log n)，但是我们还是考虑线筛.

先把T拿出来，然后思考\mu的意义.我们只需要考虑d无平方因子的情况.如果n=\prod p_i，那么\sum\limits_{d|T}d\mu(d)=1-\sum p_i+\sum p_ip_j-\cdots.经过一番思考我们可以发现这个式子=\prod(1-p_i).(每个质因子有取和不取两种选择，并且如果取就会有一个-1乘上去).于是我们就可以快乐线筛了.

void make(int n)
{
    f[1]=p[1]=1;
    for(int i=2;i<=n;i++)
    {
        if(!p[i])prime[++cnt]=i,f[i]=i*(1LL-i);
        for(int j=1;j<=cnt&&i*prime[j]<=n;j++)
        {
            int x=i*prime[j];p[x]=1;
            if(i%prime[j])f[x]=f[i]*prime[j]*(1-prime[j]);
            else f[x]=f[i]*prime[j];
        }
    }
    for(int i=2;i<=n;i++)f[i]=(f[i]+f[i-1])%mod;//处理一个前缀和
}

现在是O(n+T\sqrt{n})

单组询问，n,m\leq 10^9?

我刚才看到某个神仙写完O(n+T\sqrt{n})忘记了需要线筛这个事实从而认为自己做到了O(\sqrt{n})

显然杜教筛了.设f(n)=n\sum\limits_{d|n}d\mu(d)，我们需要杜教筛出它的前缀和.

换一换形式我们得到f(n)=\sum\limits_{d|n}d^2\mu(d)\frac{T}{d}，即f(n)=(id^2\cdot \mu)\ast id.(定义f\cdot g(n)=g(n)g(n))

我们来试试找一个g使得f\ast g和g的前缀和都可以快速计算.因为(id^k\cdot\mu)\ast id^k=\epsilon，所以我们取g=id^2，于是f\ast g=id

杜教筛就好啦.

卡了卡常数，现在是换新评测姬之后跑的最快的了(

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#define mod 20101009
#define INF 1000000000000000000LL
#define inv2 10050505
#define inv6 16750841
using namespace std;
const int N=1e7+10;
int p[N],prime[N],cnt,mu[N],N2;
long long f[N],n,m;
struct Hash_map
{
    struct pair{long long u,v;}val[N];
    int fst[N],nxt[N],node_cnt;
    long long find(long long n)
    {
        int key=n&8388607;//秘技·手写哈希表·模数2的幂
        for(int i=fst[key];i;i=nxt[i])if(val[i].u==n)return val[i].v;
        return 1e18;
    }
    void ins(long long u,long long v)
    {
        int key=u&8388607;
        val[++node_cnt]=(pair){u,v},nxt[node_cnt]=fst[key],fst[key]=node_cnt;
    }
}mpf;
void make(int n)
{
    f[1]=p[1]=1;
    for(int i=2;i<=n;i++)
    {
        if(!p[i])prime[++cnt]=i,f[i]=i*(1LL-i);
        for(int j=1;j<=cnt&&i*prime[j]<=n;j++)
        {
            int x=i*prime[j];p[x]=1;
            if(i%prime[j])f[x]=f[i]*f[prime[j]];
            else {f[x]=f[i]*prime[j],mu[x]=0;break;}
        }
    }
    for(int i=2;i<=n;i++)f[i]=(f[i]%mod+f[i-1]+mod)%mod;
}
long long S(long long n){return n*(n+1)%mod*inv2%mod;}
long long S2(long long n){return n*(n+1)%mod*(2*n+1)%mod*inv6%mod;}
long long Sf(long long n)
{
    if(n<=N2)return f[n];
    long long x=mpf.find(n);
    if(x!=INF)return x;
    long long ans=S(n);long long lstans=1;//这里记录上一次的答案来避免重复计算.大约能快那么一点点???
    for(long long i=2,lt;i<=n;i=lt+1)
    {
        long long fn=n/i;lt=n/fn;long long tmp=S2(lt);
        ans-=(tmp-lstans)*Sf(fn)%mod;
        lstans=tmp;
    }
    ans=(ans%mod+mod)%mod;mpf.ins(n,ans);return ans;
}
int main()
{
    scanf("%lld%lld",&n,&m);
    if(n>m)swap(n,m);N2=pow(n,0.7)+1;make(N2);
    long long ans=0;long long lstans=0;
    for(int i=1,lt;i<=n;i=lt+1)
    {
        int fn=n/i,fm=m/i;lt=min(n/fn,m/fm);
        long long tmp=Sf(lt);//同上
        ans=(ans+(tmp-lstans)*S(fn)%mod*S(fm)%mod)%mod;
        lstans=tmp;
    }
    cout<<(ans+mod)%mod<<endl;
}

真是快乐的卡常数题