P2260 [清华集训2012]模积和

题目背景

数学题，无背景。

题目描述

求

$sumlimits_{i=1}^{n} sumlimits_{j=1}^{m} (n mod i) imes (m mod j), i eq j ;mod;19940417$ 的值

输入输出格式

输入格式：

 

两个整数n m

 

输出格式：

 

答案 mod 19940417

 

输入输出样例

输入样例#1： 

3 4

输出样例#1： 

1

输入样例#2： 

123456 654321

输出样例#2： 

116430

说明

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;
}