题目描述
求
mod 19940417 的值
输入格式
两个整数n m
输出格式
答案 mod 19940417
输入输出样例
输入
3 4
输出
1
输入
123456 654321
输出 #2复制
116430
说明/提示
30%: n,m <= 1000
60%: n,m <= 10^6
100% n,m <= 10^9
。。。。公式很快就可以推出来。。。就是被傻逼BUG卡了好久QAQ。
我们可以用不考虑i==j的情况扣去i==j的情况,于是等式就变成了:
那么前面一坨的东西很好求,就是一个比较裸的整数分块,我们可以直接将拆开:那么这一坨就宣告结束了,如果不会的话请看戳这里
重点是后面的一坨东西。我们可以将它化简一下:(假设k=min(n,m))
然后我智障了。。。中间两项的前面多加了一个k,找了半天BUG。。。
那么前面的三项挺好求的,就不多说了,主要是最后的一项。
对于最后一项我们同样可以分块,我们知道是可以分块的,那么同样的
也可以分块,后者是前者的整数倍,那么后者的分块绝对也是前者的整数倍,也就是说虽然我们无法一次性求得后者的块的左端点和右端点,但我们可以通过前者的左右端点来一块块加上:
ll ans=0;
for (ll l=1,r=0; l<=k; l=r+1){
r=min(n/(n/l),m/(m/l));
ll p1=n/l,p2=m/l;
ans=(ans+((p1*p2%mod)*sum2(l,r)%mod)*inv6%mod)%mod;
}
这里还有一个傻逼BUG。。。乘了两次(r-l+1)。。。sum2里面已经求过了,我又再求了一次。。。
这里说明一下,由于p模数不是素数。。。。所以我们不能用费马定理求逆元。。sum2是i^2的求和公式具体的如下:
。。。至于怎么推的,哈哈哈哈,我也不知道。(说真的,傻逼BUG真的找得挺自卑的QAQ)
以下是AC代码:
#include <bits/stdc++.h>
using namespace std;
#define ll long long
const ll mod=19940417;
ll inv2,inv6;
ll exgcd(ll a,ll b,ll &x,ll &y)
{
if (b==0) {
x=1;y=0;
return a;
}
ll ret=exgcd(b,a%b,y,x);
y-=a/b*x;
return ret;
}
ll getinv(ll a,ll mod)
{
ll x,y;
ll d=exgcd(a,mod,x,y);
return d==1?(x%mod+mod)%mod:-1;
}
ll solve(ll n)
{
ll ans=0;
for (ll l=1,r=0; l<=n; l=r+1){
r=n/(n/l);
ans=(ans+((((r+l)%mod)*(r-l+1)%mod)*(n/l)%mod)*inv2%mod)%mod;
}
return ans;
}
ll work(ll k,ll n)
{
ll ans=0;
for (ll l=1,r=0; l<=k; l=r+1){
r=min(k,n/(n/l));//注意
ans=(ans+((((r+l)%mod)*(r-l+1)%mod)*(n/l)%mod)*inv2%mod)%mod;
}
return ans;
}
ll sum2(ll l,ll r)
{
l--;
ll ans1=0,ans2=0;
ans1=(l*(l+1)%mod)*(2*l+1)%mod;
ans2=(r*(r+1)%mod)*(2*r+1)%mod;
return (ans2-ans1+mod)%mod;
}
ll deal(ll k,ll n,ll m)
{
ll ans=0;
for (ll l=1,r=0; l<=k; l=r+1){
r=min(n/(n/l),m/(m/l));//注意
ll p1=n/l,p2=m/l;
ans=(ans+((p1*p2%mod)*sum2(l,r)%mod)*inv6%mod)%mod;
}
return ans;
}
int main()
{
ll n,m;
cin>>n>>m;
inv2=getinv(2ll,mod);
inv6=getinv(6ll,mod);
ll ans=0,ans1=0,ans2=0,ans3=0;
ans=((n*n%mod-solve(n)+mod)%mod)*((m*m%mod-solve(m)+mod)%mod)%mod;
// cout<<ans<<endl;
ll k=min(n,m);
ans1=(k*n%mod)*m%mod;
//cout<<(ans-ans1+mod)%mod<<endl;
//cout<<endl;
ans2=(m*work(k,n)%mod+n*work(k,m)%mod)%mod;
ans3=deal(k,n,m);
//cout<<ans2<<endl<<ans3<<endl;
cout<<(((ans-ans1+mod)%mod+ans2)%mod-ans3+mod)%mod<<endl;
return 0;
}