P2260 [清华集训2012]模积和（数论分块）

#include <iostream>
#include <cstdio>
#include <algorithm>
#include <vector>
#include<set>
#include <map>
#include <string>
using namespace std;

#define ll long long
#define pb push_back

const ll MOD = 19940417;
const int N = 1e6;


ll inv2, inv6;

// ll gcd( ll a, ll b)
// {
//     return b == 0 ? a : gcd(b, a % b);
// }

ll ksm(ll a, ll b, ll mod)
{
    ll res = 1;
    while(b){
        if(b & 1) res = (res * a) % mod;
        a = a * a % mod;
        b >>= 1;
    }
    return res;
}

inline ll fx2(ll n)
{
    return n * (n + 1) % MOD * inv2 % MOD;
}

ll fun1(ll n)
{
    ll l, r, sum1, sum2;
    sum1 = sum2 = 0;
    sum1 = n * n % MOD;
    for(l = 1, r = 0; l <= n; l = r + 1){
        r = n / (n / l);
        sum2 = (sum2 + (fx2(r) - fx2(l - 1) + MOD) * (n / l) % MOD) % MOD;
    }
    return (sum1 - sum2 + MOD) % MOD;
}

ll fun2(ll n, ll m)
{
    ll l, r, sum1, sum2;
    sum1 = sum2 = 0;
    //cout << n << " " << m << endl;
    for(l = 1, r = 0; l <= n; l = r + 1){
        r = n / (n / l);
        sum1 = (sum1 + (fx2(r) - fx2(l - 1) + MOD) * (n / l)) % MOD;
    }
    sum1 = m % MOD * sum1 % MOD;

    for(l = 1, r = 0; l <= n; l = r + 1){
        r = min(n, m / (m / l));
        sum2 = (sum2 + (fx2(r) - fx2(l - 1) + MOD) * (m / l)) % MOD;
    }
    sum2 = n % MOD * sum2 % MOD;

    return (sum1 + sum2) % MOD;
}

inline ll fx6(ll n)
{
    return n * (n + 1) % MOD * (2 * n + 1) % MOD * inv6 % MOD;
}

ll fun3(ll n, ll m)
{   
    ll l, r, sum;
    sum = 0;
    for(l = 1, r = 0; l <= n; l = r + 1){
        //取右边小的使得值正确，如果取到了右边大的，会导致a或者b变小
        r = min(n / (n / l), m / (m / l));
        // /r = n / (n / l);
        ll a = n / l;
        ll b = m / l;
        sum = (sum + (fx6(r) - fx6(l - 1) + MOD) * a % MOD * b % MOD) % MOD;     
    }
    return sum % MOD;
}



ll exgcd(ll a, ll b, ll &x, ll &y)
{

    if(b == 0){//推理1，终止条件
        x = 1;
        y = 0;
        return a;
    }
    ll r = exgcd(b, a%b, x, y);
    //先得到更底层的x2,y2,再根据计算好的x2,y2计算x1，y1。
    //推理2，递推关系
    ll t = y;
    y = x - (a/b) * y;
    x = t;
    return r;
}

void inv()
{

    //ax同余1(mod prime)  求最小正整数解
    // ax + by = n ---> ax' + by' = gcd(a, b)
    // k = gcd(a, b)  k*(ax' + by') = n
    //以下 a, b, x, y对应上面 a, b, x', y'
    ll x, y;
    //以下求出mod不为质数的逆元
    // a / b % mod c
    // a / b * b^-1 ≡ x * b^-1(mod c)
    // b * b^-1 ≡ 1 (mod c)
    // b * b^-1 + y * c  = 1
    // a * x + b * y = 1
    exgcd(2, MOD, x, y);
    inv2 = (x + MOD) % MOD;
    exgcd(6, MOD, x, y);
    inv6 = (x + MOD) % MOD;

    // cout << "inv2 = " << (inv2 + MOD) % MOD << endl;
    // cout << "inv6 = " << (inv6 + MOD) % MOD << endl;
}

void solve()
{   

    //cout << gcd(2, MOD) << " " << gcd(6, MOD) << endl;
    //cout << inv2 << " " << inv6 << endl;
    inv();
   // cout << inv2 << " " << inv6 << endl;
    ll n, m;
    cin >> n >> m;
    if(n > m) swap(n, m);
    ll ans1 = (fun1(n) * fun1(m)) % MOD;
    ll ans2 = (n * n % MOD * m % MOD - fun2(n, m) + fun3(n, m) + MOD) % MOD;
    ll ans = (ans1 - ans2 + MOD) % MOD;

    cout << ans << endl;
}

int main()
{   
    ios::sync_with_stdio(false);
    cin.tie(0);
    cout.tie(0);
    solve();
    //int ok = 0;
    //cout << "ok" << endl;
    return 0;
}