洛谷P1829 - [国家集训队]Crash的数字表格

最新推荐文章于 2025-06-11 16:51:29 发布

转载最新推荐文章于 2025-06-11 16:51:29 发布 · 146 阅读

0 ·

CC 4.0 BY-SA版权

原文链接：http://www.cnblogs.com/VisJiao/p/LgP1829.html

本文详细解析了Crash的数字表格问题的算法解决方案，通过数学推导和代码实现，展示了如何在O(n)时间复杂度下计算特定数学表达式的模运算结果，涉及积性函数、莫比乌斯函数及数论等核心概念。

Portal

Description

给出\(n,m(n,m\leq10^7)\)，计算\[ \sum_{i=1}^n \sum_{j=1}^m lcm(i,j) \bmod 20101009\]
fdsa

Solution

推推推...
\[\begin{align} ans &= \sum_{i=1}^n \sum_{j=1}^m lcm(i,j) \\ &= \sum_{i=1}^n \sum_{j=1}^m \frac{ij}{gcd(i,j)} \\ &= \sum_{d=1}^n \sum_{i=1}^n \sum_{j=1}^m \frac{ij}{d}[gcd(i,j)=d] \\ &= \sum_{d=1}^n \sum_{i=1}^{⌊\frac{n}{d}⌋} \sum_{j=1}^{⌊\frac{m}{d}⌋} ijd[gcd(i,j)=1] \\ &= \sum_{d=1}^n \sum_{i=1}^{⌊\frac{n}{d}⌋} \sum_{j=1}^{⌊\frac{m}{d}⌋} ijd \sum_{d'|gcd(i,j)} \mu(d') \\ &= \sum_{d=1}^n \sum_{d'=1}^n d\mu(d') \sum_{d'|i}^{⌊\frac{n}{d}⌋} \sum_{d'|j}^{⌊\frac{m}{d}⌋} ij \\ &= \sum_{d=1}^n \sum_{d'=1}^n d\mu(d') \sum_{i=1}^{⌊\frac{n}{dd'}⌋} \sum_{j=1}^{⌊\frac{m}{dd'}⌋} ijd'^2 \\ &= \frac{1}{4} \sum_{d=1}^n \sum_{d'=1}^n dd'^2\mu(d') (1+⌊\frac{n}{dd'}⌋)⌊\frac{n}{dd'}⌋(1+⌊\frac{m}{dd'}⌋)⌊\frac{m}{dd'}⌋ \\ &= \frac{1}{4} \sum_{d=1}^n (1+⌊\frac{n}{d}⌋)⌊\frac{n}{d}⌋(1+⌊\frac{m}{d}⌋)⌊\frac{m}{d}⌋ \sum_{d'|d} dd'\mu(d') \end{align}\] \((8)\)到\((9)\)转换了枚举对象，\(d=dd',d'=d'\)。
前面的部分可以\(O(\sqrt n)\)搞定，那么我们只需要考虑求出\(g(d)=\sum_{d'|d} dd'\mu(d')\)的前缀和。由于积性函数的积是积性函数，积性函数的卷积也是积性函数，所以\(g(d)\)是积性函数，可以在\(O(n)\)的时间内预处理出来。

时间复杂度\(O(n)\)。

Code

//[国家集训队]Crash的数字表格
#include <algorithm>
#include <cstdio>
using std::min; using std::swap;
typedef long long lint;
const int N=1e7+10;
const int H=20101009;
int mu[N]; lint f[N];
int cntP,pr[N],p1[N]; bool notP[N];
void init(int n)
{
    mu[1]=1,f[1]=1;
    for(lint i=2;i<=n;i++)
    {
        if(!notP[i])
        {
            pr[++cntP]=i; p1[i]=i;
            mu[i]=-1,f[i]=(1-i+H)*i%H;
            for(lint j=i*i;j<=n;j*=i) p1[j]=j,mu[j]=0,f[j]=(1-i+H)*j%H;
        }
        for(int j=1;j<=cntP;j++)
        {
            lint x=pr[j]*i; if(x>n) break;
            notP[x]=true;
            if(i%pr[j]) p1[x]=pr[j],mu[x]=-mu[i],f[x]=f[pr[j]]*f[i]%H;
            else {p1[x]=p1[i]*pr[j],mu[x]=0,f[x]=f[p1[x]]*f[x/p1[x]]%H; break;}
        }
    }
    for(int i=1;i<=n;i++) f[i]+=f[i-1];
}
int main()
{
    int n,m; scanf("%d%d",&n,&m);
    if(n>m) swap(n,m); init(n);
    lint ans=0;
    for(int L=1,R;L<=n;L=R+1)
    {
        lint v1=n/L,v2=m/L; R=min(n/v1,m/v2);
        ans+=(1+v1)*v1%H*(1+v2)%H*v2%H*(f[R]-f[L-1]+H)%H;
        ans%=H;
    }
    printf("%lld\n",ans*15075757%H);    //inv(4)=15075757
    return 0;
}