题目背景
数学题，无背景。
题目描述
求
$sumlimits_{i=1}^{n} sumlimits_{j=1}^{m} (n mod i) imes (m mod j), i eq j ;mod;19940417$ 的值
输入输出格式
输入格式：
两个整数n m
输出格式：
答案 mod 19940417
输入输出样例
说明
30%: n,m <= 1000
60%: n,m <= 10^6
100% n,m <= 10^9
Solution:
本题实在是太贼有意思了。。。
开始没有发现条件$i eq j$，结果一直以为取模错误，搞了半天。
首先，我们先忽略$i eq j$的条件，直接求$sumlimits_{i=1}^{n} sumlimits_{j=1}^{m} (n mod i) imes (m mod j)$。
对原式化简：原式$=sumlimits_{i=1}^{n}{(n-i imeslfloor{n/i} floor)sumlimits_{j=1}^{m}{(m-j imeslfloor{m/j} floor)}}=(n^2-sumlimits_{i=1}^{n}{(i imeslfloor{n/i} floor))(m^2-sumlimits_{j=1}^{m}{(j imeslfloor{m/j} floor)})}$，然后对这个式子两边各自一遍数列分块套上等差数列求和，求出$ans$值并取模。
由于多算了$i==j$的情况，所以我们还要从$ans$中减去$i==j$的情况。
对于$i==j$的情况累加的值$tot$，容易得出$tot=sumlimits_{i=1}^{min(n,m)}{(n-i imeslfloor{n/i} floor) imes(m-j imeslfloor{m/i} floor)}=sumlimits_{i=1}^{min(n,m)}{n imes m+i^2 imeslfloor{n/i imes m/i} floor-n imes i imes lfloor{m/i} floor-m imes i imeslfloor{n/i} floor}$
然后对于$n imes m$直接累加，对于$n imes i imes lfloor{m/i} floor;,;m imes i imeslfloor{n/i} floor$还是数论分块套上等差数列求和。
难的是求$i^2 imeslfloor{n/i imes m/i} floor$，此时应用一个数学公式：$1^2+2^2+…+x^2=frac{x imes(x+1) imes(2 imes x+1)}{6}$，设$c=min(n/(n/i),m/(m/i))$，那么$i^2+{(i+1)}^2+…+c^2=frac{c imes(c+1) imes(2 imes c+1)}{6}-frac{(i-1) imes((i-1)+1) imes(2 imes (i-1)+1)}{6}$，这样原式就能直接公式求了，但是因为存在取模的条件，所以此时骚操作是处理出$6$的因子并约掉，或者直接预先算出$6$关于模数的逆元$inv$。
最后输出$ans$就好了。
代码：
#include<bits/stdc++.h> #define il inline #define ll long long #define For(i,a,b) for(int (i)=(a);(i)<=(b);(i)++) #define Bor(i,a,b) for(int (i)=(b);(i)>=(a);(i)--) #define Max(a,b) ((a)>(b)?(a):(b)) #define Min(a,b) ((a)>(b)?(b):(a)) using namespace std; const int mod = 19940417 , inv = 3323403; ll n,m; il ll solve(ll x){ ll ans=(x%mod*x%mod)%mod,p,c; for(ll i=1;i<=x;i=p+1){ p=x/(x/i); ans=(ans-(p+i)*(p-i+1)/2%mod*(x/i)%mod+mod)%mod; } return ans; } il ll get(ll x){return x*(x+1)%mod*(x<<1|1)%mod*inv%mod;} int main(){ ios::sync_with_stdio(0); cin>>n>>m; ll p,sum1,sum2,sum3,ans=solve(n)*solve(m)%mod; if(n>m)swap(n,m); for(ll i=1;i<=n;i=p+1){ p=Min(n/(n/i),m/(m/i)); sum1=(m*n%mod*(p-i+1))%mod; sum2=(n/i)*(m/i)%mod*(get(p)-get(i-1)+mod)%mod; sum3=(p-i+1)*(p+i)/2%mod*(n/i*m%mod+m/i*n%mod); ans=(ans-(sum1+sum2-sum3)%mod+mod)%mod; } cout<<ans%mod; return 0; }