模積和 洛谷2260 bzoj2956 整除分塊、逆元

題意：

求

∑i=1n∑j=1m(n mod i)∗(m mod j),i≠j ∑ i = 1 n ∑ j = 1 m ( n   m o d   i ) ∗ ( m   m o d   j ) , i ≠ j

mod 19940417 m o d   19940417 的值

首先我們直接求是 O(n∗m) O ( n ∗ m ) 的，那麼我們考慮推一下這個式子。我們發現如果沒有 i≠j i ≠ j 的限制，那麼我們隻需要分别求出 ∑ni=1n mod i ∑ i = 1 n n   m o d   i 和 ∑mj=1m mod j ∑ j = 1 m m   m o d   j ，因為 n mod i=n−⌊ni⌋∗i n   m o d   i = n − ⌊ n i ⌋ ∗ i ，是以可以整數分塊來算每一個的結果，可以用然後做乘積即可。但是有了這個限制我們應該怎麼辦呢？我們考慮用包含 i=j i = j 的總和減去 i=j i = j 的情況。 ∑i=1n∑j=1m(n mod i)∗(m mod j),i≠j ∑ i = 1 n ∑ j = 1 m ( n   m o d   i ) ∗ ( m   m o d   j ) , i ≠ j

∑i=1n∑j=1m(n mod i)∗(m mod j)−∑i=1min(n,m)(n mod i)∗(m mod i) ∑ i = 1 n ∑ j = 1 m ( n   m o d   i ) ∗ ( m   m o d   j ) − ∑ i = 1 m i n ( n , m ) ( n   m o d   i ) ∗ ( m   m o d   i )

∑i=1n(n mod i)∗∑j=1m(m mod j)−∑i=1min(n,m)(n−⌊ni⌋∗i)∗(m−⌊mi⌋∗i) ∑ i = 1 n ( n   m o d   i ) ∗ ∑ j = 1 m ( m   m o d   j ) − ∑ i = 1 m i n ( n , m ) ( n − ⌊ n i ⌋ ∗ i ) ∗ ( m − ⌊ m i ⌋ ∗ i )

∑i=1n(n−⌊ni⌋∗i)∗∑j=1m(m−⌊mi⌋∗i)−∑i=1min(n,m)n∗m−m∗⌊ni⌋∗i−n∗⌊mi⌋∗i+⌊ni⌋∗⌊mi⌋∗i2 ∑ i = 1 n ( n − ⌊ n i ⌋ ∗ i ) ∗ ∑ j = 1 m ( m − ⌊ m i ⌋ ∗ i ) − ∑ i = 1 m i n ( n , m ) n ∗ m − m ∗ ⌊ n i ⌋ ∗ i − n ∗ ⌊ m i ⌋ ∗ i + ⌊ n i ⌋ ∗ ⌊ m i ⌋ ∗ i 2

好了，推導完畢，我們發現這個式子的前半部分已經講過可以數論分塊來用 O(n−−√+m−−√) O ( n + m ) 的時間求了，後半部分同理也可以。

這裡需要補充一個公式：

∑i=1ni2=16n(n+1)(2n+1) ∑ i = 1 n i 2 = 1 6 n ( n + 1 ) ( 2 n + 1 )

這裡我就不給出證明了，有興趣的話可以自行搜尋。

最後吐槽一句，這種取模是真的煩，式子本來就長，再加上不停地取模操作，還不能忘了 (x+mod)%mod ( x + m o d ) % m o d ，不然還會出現負數，真的很讨厭。

代碼：

#include <bits/stdc++.h>
using namespace std;

long long n,m,ans,mod=,ans1,ans2,ni;
void exgcd(long long a,long long b,long long &x,long long &y)
{
    if(!b)
    {
        x=;
        y=;
    }
    else
    {
        exgcd(b,a%b,y,x);
        y-=a/b*x;
    }
}
long long calc(long long x)
{
    return (x*(x+)%mod*(*x+)%mod*ni%mod)%mod;
}
int main()
{
    scanf("%lld%lld",&n,&m);
    long long r;
    long long x,y;
    exgcd(,mod,x,y);
    x=(x%mod+mod)%mod;
    if(!x)
    x+=mod;
    ni=x;
    ans1=(n*n)%mod;
    ans2=(m*m)%mod;
    for(long long l=;l<=n;l=r+)
    {
        if(n/l!=)
        r=min(n/(n/l),n);
        else
        r=n;
        ans1=(ans1-(((r-l+)*(l+r)/%mod*(n/l))%mod+mod)%mod)%mod;
        ans1=(ans1+mod)%mod;
    }
    for(long long l=;l<=m;l=r+)
    {
        if(m/l!=)
        r=min(m/(m/l),m);
        else
        r=m;
        ans2=(ans2-(((r-l+)*(l+r)/%mod*(m/l))%mod+mod)%mod)%mod;
        ans2=(ans2+mod)%mod;
    }
    ans=(ans1*ans2)%mod;
    for(int l=;l<=min(n,m);l=r+)
    {
        r=min(n/(n/l),m/(m/l));
        ans=(ans-((r-l+)*((m*n)%mod)%mod+(calc(r)-calc(l-))%mod*((n/l)*(m/l)%mod)%mod-((r-l+)*(l+r)/)%mod*m%mod*(n/l)%mod-((r-l+)*(l+r)/)%mod*n%mod*(m/l)%mod))%mod;
        ans=(ans+mod)%mod;
    }
    printf("%lld\n",ans);
    return ;
}