題目描述
求
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;
}