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

 

 

 

  1 #include <iostream>
  2 #include <cstdio>
  3 #include <algorithm>
  4 #include <vector>
  5 #include<set>
  6 #include <map>
  7 #include <string>
  8 using namespace std;
  9 
 10 #define ll long long
 11 #define pb push_back
 12 
 13 const ll MOD = 19940417;
 14 const int N = 1e6;
 15 
 16 
 17 ll inv2, inv6;
 18 
 19 // ll gcd( ll a, ll b)
 20 // {
 21 //     return b == 0 ? a : gcd(b, a % b);
 22 // }
 23 
 24 ll ksm(ll a, ll b, ll mod)
 25 {
 26     ll res = 1;
 27     while(b){
 28         if(b & 1) res = (res * a) % mod;
 29         a = a * a % mod;
 30         b >>= 1;
 31     }
 32     return res;
 33 }
 34 
 35 inline ll fx2(ll n)
 36 {
 37     return n * (n + 1) % MOD * inv2 % MOD;
 38 }
 39 
 40 ll fun1(ll n)
 41 {
 42     ll l, r, sum1, sum2;
 43     sum1 = sum2 = 0;
 44     sum1 = n * n % MOD;
 45     for(l = 1, r = 0; l <= n; l = r + 1){
 46         r = n / (n / l);
 47         sum2 = (sum2 + (fx2(r) - fx2(l - 1) + MOD) * (n / l) % MOD) % MOD;
 48     }
 49     return (sum1 - sum2 + MOD) % MOD;
 50 }
 51 
 52 ll fun2(ll n, ll m)
 53 {
 54     ll l, r, sum1, sum2;
 55     sum1 = sum2 = 0;
 56     //cout << n << " " << m << endl;
 57     for(l = 1, r = 0; l <= n; l = r + 1){
 58         r = n / (n / l);
 59         sum1 = (sum1 + (fx2(r) - fx2(l - 1) + MOD) * (n / l)) % MOD;
 60     }
 61     sum1 = m % MOD * sum1 % MOD;
 62 
 63     for(l = 1, r = 0; l <= n; l = r + 1){
 64         r = min(n, m / (m / l));
 65         sum2 = (sum2 + (fx2(r) - fx2(l - 1) + MOD) * (m / l)) % MOD;
 66     }
 67     sum2 = n % MOD * sum2 % MOD;
 68 
 69     return (sum1 + sum2) % MOD;
 70 }
 71 
 72 inline ll fx6(ll n)
 73 {
 74     return n * (n + 1) % MOD * (2 * n + 1) % MOD * inv6 % MOD;
 75 }
 76 
 77 ll fun3(ll n, ll m)
 78 {   
 79     ll l, r, sum;
 80     sum = 0;
 81     for(l = 1, r = 0; l <= n; l = r + 1){
 82         //取右边小的使得值正确，如果取到了右边大的，会导致a或者b变小
 83         r = min(n / (n / l), m / (m / l));
 84         // /r = n / (n / l);
 85         ll a = n / l;
 86         ll b = m / l;
 87         sum = (sum + (fx6(r) - fx6(l - 1) + MOD) * a % MOD * b % MOD) % MOD;     
 88     }
 89     return sum % MOD;
 90 }
 91 
 92 
 93 
 94 ll exgcd(ll a, ll b, ll &x, ll &y)
 95 {
 96 
 97     if(b == 0){//推理1，终止条件
 98         x = 1;
 99         y = 0;
100         return a;
101     }
102     ll r = exgcd(b, a%b, x, y);
103     //先得到更底层的x2,y2,再根据计算好的x2,y2计算x1，y1。
104     //推理2，递推关系
105     ll t = y;
106     y = x - (a/b) * y;
107     x = t;
108     return r;
109 }
110 
111 void inv()
112 {
113 
114     //ax同余1(mod prime)  求最小正整数解
115     // ax + by = n ---> ax' + by' = gcd(a, b)
116     // k = gcd(a, b)  k*(ax' + by') = n
117     //以下 a, b, x, y对应上面 a, b, x', y'
118     ll x, y;
119     //以下求出mod不为质数的逆元
120     // a / b % mod c
121     // a / b * b^-1 ≡ x * b^-1(mod c)
122     // b * b^-1 ≡ 1 (mod c)
123     // b * b^-1 + y * c  = 1
124     // a * x + b * y = 1
125     exgcd(2, MOD, x, y);
126     inv2 = (x + MOD) % MOD;
127     exgcd(6, MOD, x, y);
128     inv6 = (x + MOD) % MOD;
129 
130     // cout << "inv2 = " << (inv2 + MOD) % MOD << endl;
131     // cout << "inv6 = " << (inv6 + MOD) % MOD << endl;
132 }
133 
134 void solve()
135 {   
136 
137     //cout << gcd(2, MOD) << " " << gcd(6, MOD) << endl;
138     //cout << inv2 << " " << inv6 << endl;
139     inv();
140    // cout << inv2 << " " << inv6 << endl;
141     ll n, m;
142     cin >> n >> m;
143     if(n > m) swap(n, m);
144     ll ans1 = (fun1(n) * fun1(m)) % MOD;
145     ll ans2 = (n * n % MOD * m % MOD - fun2(n, m) + fun3(n, m) + MOD) % MOD;
146     ll ans = (ans1 - ans2 + MOD) % MOD;
147 
148     cout << ans << endl;
149 }
150 
151 int main()
152 {   
153     ios::sync_with_stdio(false);
154     cin.tie(0);
155     cout.tie(0);
156     solve();
157     //int ok = 0;
158     //cout << "ok" << endl;
159     return 0;
160 }

 

posted @ 2020-08-06 12:14  SummerMingQAQ  阅读(165)  评论(0)  编辑  收藏  举报