【BZOJ-2956】(清华集训2012)模积和

题目链接

题目描述

求$\sum\sum((n\;mod\;i)*(m\;mod\;j))其中1<=i<=n,1<=j<=m,i≠j$。

例如对于 n=3,m=4 :
答案为(3 mod 1)(4 mod 2)+(3 mod 1) (4 mod 3)+(3 mod 1) * (4 mod 4) + (3 mod 2) * (4 mod 1) + (3 mod 2) * (4 mod 3) + (3 mod 2) * (4 mod 4) + (3 mod 3) * (4 mod 1) + (3 mod 3) * (4 mod 2) + (3 mod 3) * (4 mod 4) = 1

题解

可以先去做做余数求和 ，我的题解。

一开始看漏了$i\ne j$，以为这和余数求和一样…
但是既然又和取模的和有关，与余数求和的思想是差不多的。
我们有$a\;mod\;b=a−floor(a/b)∗b$
于是我们再来推一波式子:

$$\sum_{i=1}^{n}\sum_{j=1,j\ne i}^{m}((n\;mod\;i)*(m\;mod\;j))$$
$=\sum\sum((n\;mod\;i)*(m\;mod\;j))-\sum^{min(n,m)}_{i=1}((n\;mod\;i)*(m\;mod\;i))$

前面那段就是余数求和的方法，不多讲，我们来看后面一段。
不妨设$n\le m$
那么:
$\sum^{n}_{i=1}((n\;mod\;i)*(m\;mod\;i))$
$=\sum((n-n/i*i)*(m-m/i*i))$
令$n/i=x,m/i=y$
原式化为:
$\sum(mn-mxi-nyi+xy*i^2)$
$\sum(mn+xy*i^2-(mx+ny)*i)$

到这里式子已经推的差不多了。
我们发现可以对$x,y$进行除法分块，且前面的mn可以预先处理掉。
于是我们着眼于求出一个块内的$(xy*i^2-(mx+ny)*i)$
假定当前块内 $xy=S_1,mx+ny=S_2$
那么要求的就是:
$S_1*i^2-S_2*i$
前面这一坨要用到平方和公式即:$\sum^n_{i=1}(i^2)=n*(n+1)*(2n+1)/6$
注意这里要求个逆元。
后面的直接等差数列求和即可，那么本题就做完了。

代码如下:

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<algorithm>
#include<cstring>
#include<set>
using namespace std;
typedef long long ll;
const int N=1e9+10;
const ll mod=19940417;
const ll inv=3323403;//6的逆元
inline ll Mod_sum(ll k,ll n)
{
    ll ans=n*k;
    register ll l=1,r=0;
    for(l=1;l<=n;l=r+1)
    {
        if(k/l!=0) r=min(k/(k/l),n);
        else r=n;
        ans-=((k/l)*(r-l+1)*(l+r)/2)%mod;
        if(ans<0) ans+=mod;
    }
    return ans%mod;
}
inline ll sum(ll l,ll r)
{
    return ((r+l)*(r-l+1)/2)%mod;
}
inline ll sum2(ll n)
{
    if(n==0) return 0;
    if(n==1) return 1;
    return ((n*(n+1)%mod)*(2*n%mod+1)%mod)*inv%mod;
}
int main()
{
    register ll ans=0;ll n,m;
    scanf("%lld %lld",&n,&m);
    ans=Mod_sum(n,n);
    ans*=Mod_sum(m,m);
    if(n>m) swap(n,m);
    ans-=(n*(m*n%mod))%mod;
    ans=(ans+mod)%mod;
    ll l;ll r=0;
    for(l=1;l<=n;l=r+1){
        r=min(m/(m/l),n/(n/l));//取两个都相同的块
        register ll s1=((m/l)*(n/l)%mod)*((sum2(r)-sum2(l-1)+mod)%mod)%mod;//平方和公式
        register ll s2=(((n*(m/l)%mod+m*(n/l)%mod)%mod)*sum(l,r))%mod;
        ans=(ans-(s1-s2+mod)+mod)%mod;
    }
    while(ans<0) ans+=mod;
    printf("%lld\n",ans);
    return 0;
}