BZOJ2956：模积和

传送门

• 题意
  求

  $$\sum\sum((n\%i)*(m\%j))(1\le i\le n,1\le j\le m,i!=j)$$

• 题解

如果没有$i!=j$这一条件这道题就简单了。

$$\begin{align}\sum(n\%i)&=\sum(n-\lfloor \frac{n}{i}\rfloor*i)\\&=n^2-\sum(\lfloor \frac{n}{i}\rfloor*i)\end{align}$$

对于$\lfloor \frac{n}{i}\rfloor$相同的区间一起处理即可。

证明：对于任意$n$，$\lfloor \frac{n}{i}\rfloor$的数量级为$\sqrt{n}$。
因为对于$i<\sqrt{n}$，i的取值最多$\sqrt{n}$个，此时区间也不超过$\sqrt{n}$个，对于$i\ge \sqrt{n}$,$\lfloor \frac{n}{i}\rfloor$的取值最多$\sqrt{n}$个(因为不超过$\sqrt{n}$)。所以$\lfloor \frac{n}{i}\rfloor$的数量级为$\sqrt{n}$。

加上$i!=j$这一条件后，将原式变形，得：

$$\begin{align}\sum_{i=1}^{n}\sum_{j=1\lor i\not =j}^{m}((n\%i)*(m\%j))&=\sum_{i=1}^{n}(n\%i)*\sum_{j=1\lor i\not =j}^{m}(m\%j)\\&=\sum_{i=1}^{n}(n\%i)*\sum_{j=1}^{m}(m\%j)-\sum_{i=1}^{min(n,m)}(n\%i)*(m\%i)\end{align}$$

对于前面一部分很好处理，对于后面一部分，再进一步变形得

$$\begin{align}\sum_{i=1}^{min(n,m)}(n\%i)*(m\%i)&=\sum_{i=1}^{min(n,m)}(n-\lfloor \frac{n}{i}\rfloor*i )*(m-\lfloor \frac{m}{i}\rfloor*i )\\&=\sum_{i=1}^{min(n,m)}(n*m+\lfloor \frac{n}{i}\rfloor\lfloor \frac{m}{i}\rfloor*i^2-m*\lfloor \frac{n}{i}\rfloor*i-n*\lfloor \frac{m}{i}\rfloor*i)\end{align}$$

对于$deg(i)=1$的部分处理方法同上，对于$deg(i)=2$的部分不能再用前面使用的等差数列求值法，这时容易推出：

$$n^2=\frac{n^3-(n-1)^3+3n-1}{3}$$
那么前缀和相抵消可以快速求出 $$1^2+2^2+3^2+...+n^2=\frac{n*(n+1)*(2*n+1)}{6}$$
要求得一段连续的平方和，只需用前缀和相减即可。

那么问题是，在模意义下，除法不能直接运算，这里就需要用到数论中逆元的知识：

$$a/b=a*b_1(mod \,m)(b*b1\%m=1,b|a)$$

要注意的是，对于一个数采取逆元运算前不能对其进行任何操作（必须保证这个数能被b整除）。

分析到这里代码也就很好写：

• Code:

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<string>
#include<algorithm>
#include<cmath>
using namespace std;
typedef long long ll;
const int Mod=19940417;
const int inv=3323403;

ll n,m;

inline ll calc(ll N)
{
    ll ans=0,pos=0;
    ans=N*N;
    for(ll bg=1;bg<=N;bg=pos+1)
    {
        pos=min(N/(N/bg),N);
        ans-=((pos-bg+1)*(pos+bg)*(N/bg))/2;
    }
    return ans%Mod;
}

inline ll sum2(ll x,ll y)
{
    return ((y-x+1)*(x+y)/2)%Mod;
}
inline ll sum1(ll x)
{
    return (x*(x+1)%Mod*(2*x+1)%Mod*inv)%Mod;
}

inline ll calc2(ll N)
{
    ll ans=0,pos=0;
    ans=n*m%Mod*N%Mod;
    for(ll bg=1;bg<=N;bg=pos+1)
    {
        pos=min(n/(n/bg),m/(m/bg));
        ll t1=n/bg,t2=m/bg;
        ll s1=t1*t2%Mod*((sum1(pos)-sum1(bg-1)+Mod)%Mod)%Mod;
        ll s2=(t1*m%Mod+t2*n%Mod)%Mod*sum2(bg,pos)%Mod;
        ans=(ans+(s1-s2+Mod)%Mod)%Mod;
    }
    return ans%Mod;
}

int main()
{
    cin>>n>>m;
    cout<<(calc(n)*calc(m)-calc2(min(n,m))+Mod)%Mod;
}