洛谷P2260 [清华集训2012]模积和

@Pelom  October 12, 2021

题意


$$\sum_{i=1}^n\sum_{j=1}^m (n\%i) \cdot (m\%j),i \neq j$$
$\bmod \ 19940417$的值

数据范围:$n,m \le 10^9$

题解

算法:数论分块
先不考虑条件$i \neq j$,化简原式

$$ \begin{aligned} \sum_{i=1}^n\sum_{j=1}^m (n\%i) \cdot (m\%j) &= \sum_{i=1}^n (n\%i) \cdot \sum_{j=1}^m (m\%j) \\ &= \sum_{i=1}^n (n- i\lfloor \frac{n}{i} \rfloor) \cdot \sum_{j=1}^m (m- j\lfloor \frac{m}{j} \rfloor) \\ &= (n^2-\sum_{i=1}^n i\lfloor \frac{n}{i} \rfloor) \cdot (m^2-\sum_{j=1}^m j\lfloor \frac{m}{j} \rfloor) \end{aligned} $$

对两项分别使用数论分块计算即可

inline int calc(int n){
    int res=n%mod*n%mod;
    for(int l=1,r;l<=n;l=r+1){
        r=n/(n/l);
        (res-=(r-l+1)*(l+r)/2%mod*(n/l)%mod)%=mod;
    }
    return res;
}
int main(){
    ans=calc(n)*calc(m)%mod;
    return 0;
}

由容斥原理,现在需要减去满足$i=j$的项

$$ \begin{aligned} \sum_{i=1}^{Min(n,m)} (n\%i) \cdot (m\%i) &= \sum_{i=1}^{Min(n,m)} (n-\lfloor \frac{n}{i} \rfloor) \cdot (m-\lfloor \frac{m}{i} \rfloor) \\ &= \sum_{i=1}^{Min(n,m)} (nm-in\lfloor \frac{m}{i} \rfloor-im\lfloor \frac{n}{i} \rfloor+i^2 \lfloor \frac{n}{i} \rfloor \lfloor \frac{m}{i} \rfloor) \\ &= nm \cdot Min(n,m)-\sum_{i=1}^{Min(n,m)}(in\lfloor \frac{m}{i} \rfloor+im\lfloor \frac{n}{i} \rfloor-i^2 \lfloor \frac{n}{i} \rfloor \lfloor \frac{m}{i} \rfloor) \end{aligned} $$

对于求和中的前两项可以进行同样的数论分块,但对于第三项,在数论分块时会遇到累计的$i^2$,于是需要求$\sum\limits_{i=l}^r i^2$,此时需要用到公式
$$\sum_{i=1}^n i^2=\frac{n(n+1)(2n+1)}{6}$$
对于$\sum\limits_{i=l}^r i^2$,运用前缀和思想即可
而对于此公式的求值,由于取模的限制,需要先求出$6$关于模数$19940417$的逆元。此处模数为合数,不可使用费马小定理,但$(6,19940417)=1$,可以使用欧拉定理或者使用扩展欧几里得算法求得

注意:
1.数据较大,需要使用$long \ long$
2.随时取模,以免溢出
3.由于代码中有需要$/2$的地方,请将其改为乘以$2$关于模数的逆元,或将(r-l+1)*(l+r)/2写在一起以保证整除($/6$同理,但由于分子已过大,必须先取模,因此必须使用逆元)

代码:

#include<iostream>
#include<cstdio>
using namespace std;
#define int long long
const int mod=19940417;
const int inv=3323403;
int n,m,p;
int s1,s2,s3,ans;
template <typename T>
inline T Min(T a,T b){
    return a<b?a:b;
}
inline int calc(int n){
    int res=n%mod*n%mod;
    for(int l=1,r;l<=n;l=r+1){
        r=n/(n/l);
        (res-=(r-l+1)*(l+r)/2%mod*(n/l)%mod)%=mod;
    }
    return res;
}
inline int S2(int n){
    return n*(n+1)%mod*(2*n+1)%mod*inv%mod;
}
signed main(){
    scanf("%lld%lld",&n,&m);
    ans=calc(n)*calc(m)%mod;
    p=Min(n,m);
    (ans-=p*n%mod*m%mod)%=mod;
    for(int l=1,r;l<=p;l=r+1){
        r=Min(n/(n/l),m/(m/l));
        s1=(r-l+1)*(l+r)/2%mod*((n/l)*m%mod+(m/l)*n%mod)%mod;
        s2=(n/l)*(m/l)%mod*(S2(r)-S2(l-1))%mod;
        (ans+=s1-s2)%=mod;
    }
    (ans+=mod)%=mod;
    printf("%lld",ans);
    return 0;
}

添加新评论