Luogu P2260 模积和

题目描述

\sum^{n}_{i=1}\sum^{m}_{j=1}(n \bmod i) \times(m \bmod j) \quad i\neq j

\bmod 19940417 的值

Solution

考虑容斥, 先计算出原式 无 i \neq j 限制 的答案,然后再减去 原式 限制为 i=j 的答案,就是原式 附加i \neq j 的答案。

我们令 n\leq m

即求
\sum^{n}_{i=1}\sum^{m}_{j=1}(n \bmod i) \times(m \bmod j)-\sum^{n}_{i=1}(n \bmod i) \times(m \bmod i)
我们知道
n \bmod i = n- \lfloor \frac {n} {i} \rfloor*i
于是代入得到
\sum^{n}_{i=1}\sum^{m}_{j=1}(n- \lfloor \frac {n} {i} \rfloor*i) \times(m- \lfloor \frac {m} {j} \rfloor*j)
括号拆开,得到
\sum^{n}_{i=1}\sum^{m}_{j=1}n*m-(m*\lfloor \frac {n} {i} \rfloor*i+n*\lfloor \frac {m} {j} \rfloor*j)+\lfloor\frac {n} {i} \rfloor*i*\lfloor \frac {m} {j} \rfloor*j
利用整除分块优化复杂度


N=\sum^n_{i=1}\lfloor \frac {n} {i} \rfloor*i \quad M=\sum^m_{j=1}\lfloor \frac {m} {j} \rfloor*j
原式化简为

n^2*m^2-(m^2*N+n*n*M)+N*M

同样化简 \sum^{n}_{i=1}(n \bmod i) \times (m \bmod i),可得到
\sum^n_{i=1}n*m-(n*\lfloor \frac{m}{i} \rfloor+m*\lfloor \frac{n}{i} \rfloor)*i+i^2*\lfloor \frac{m}{i} \rfloor*\lfloor \frac{n}{i} \rfloor
\sum^n_{i=1} i^2 =\frac{n*(n+1)*(2*n+1)}{6}

此处除以 6 需要计算逆元

Code:

#include<bits/stdc++.h>
#define LL long long
LL mod = 19940417;
LL n,m;
LL inv[10];
LL calc_s(LL x){
    return (((x*(x+1)%mod)*(2*x+1)%mod)*inv[6])%mod;
}
LL calc_s(LL l,LL r){
    return ((calc_s(r)-calc_s(l-1))%mod+mod)%mod;
}
LL calc(LL n,LL x){
    LL ans=0;
    for(LL l=1,r;l<=n;l=r+1){
        r=x/(x/l);
        r=std::min(r,n);
        ans+=(((x/l)*(r-l+1)%mod*(l+r)%mod)*inv[2])%mod;
        ans%=mod;
    }
    return ans;
}
int main(){
    scanf("%lld%lld",&n,&m);
    if(n>m) std::swap(n,m);
    LL ans=(((((n*m)%mod)*n)%mod)*m)%mod,mans=0,nans=0,m_nans=0;
    inv[1]=1;
    for(int i=2;i<=6;++i){
        inv[i]=((1ll*(mod-mod/i))*inv[mod%i])%mod; 
    }
    nans=calc(n,n);
    mans=calc(m,m);
    ans=(ans-(((m*m%mod)*nans)%mod)+mod)%mod;
    ans=(ans-((n*n%mod)*mans%mod)%mod+mod)%mod;
    ans=(ans+mans*nans%mod)%mod;
    ans=(ans-((n*m%mod)*n%mod)%mod+mod)%mod;
    ans+=m*nans%mod,ans%=mod;
    for(LL l=1,r;l<=n;l=r+1){
        r=m/(m/l);
        r=std::min(r,n);
        m_nans+=((((m/l)*(r-l+1)%mod)*(l+r)%mod)*inv[2])%mod;
        m_nans%=mod;
    }
    ans+=m_nans*n%mod;ans%=mod;
    for(LL l=1,r;l<=n;l=r+1){
        r=std::min(n/(n/l),m/(m/l));
        ans=(ans-((calc_s(l,r)*(n/l)%mod)*(m/l)%mod)+mod)%mod;
    }
    printf("%lld",ans);
    return 0;
}
分类: 数论

0 条评论

发表评论

Avatar placeholder

您的电子邮箱地址不会被公开。 必填项已用*标注