整数分块毒瘤题（非素数模）--模积和（洛谷P2260&BZOJ2956）详解

题目描述

求

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