题意
给定 \(n,m,p\)，求
\[\sum_{i=1}^n \sum_{j=1}^m (i \times j \bmod p) \]分析
看到 \(1 \le n,m \le 10^{12}\)，知道 \(O(nm)\) 的暴力算法行不通。但 \(1 \le p \le 10^3\)，范围很小，是本题的突破口。
那怎么利用这个特点呢？这里需要用到一条同余的性质，即
\[(i \times j) \bmod p = (i \bmod p \times j \bmod p) \bmod p \]具体证明方法这里不详解。
由此我们可以发现，\(i\) 和 \(j\) 取余后都满足 \(0 \le i,j < p\)，而 \(p\) 的范围很小，因此考虑先预处理所有的 \(i \times j \bmod p(0 \le i,j < p)\)，然后统计不同余数出现的次数，最后运用乘法原理计数，同时取模即可。
注意输入的 \(p\) 和 \(10^9+7\) 两个模数用在不同的地方，预处理时用 \(p\)，计数时用 \(10^9+7\)。
Code
代码如下：
#include <iostream>
#include <cstdio>
#define int long long
using namespace std;
const int MOD=1000000007;
int cnt[1001][1001];
int a[1001],b[1001],x,y;
int ans;
signed main()
{
int n,m,p;
cin>>n>>m>>p;
x=n%p,y=m%p;
for(int i=0;i<p;i++)
for(int j=0;j<p;j++)
cnt[i][j]=i*j%p;
for(int i=0;i<p;i++)
{
a[i]=n/p;
b[i]=m/p;
if(i<=x) a[i]++;
if(i<=y) b[i]++;
}
for(int i=0;i<p;i++)
for(int j=0;j<p;j++)
ans=(ans+(cnt[i][j]*a[i]%MOD*b[j])%MOD)%MOD;
cout<<ans;
return 0;
}