jzptab
Time Limit: 10 Sec Memory Limit: 512 MB[Submit][Status][Discuss]
Description
求
Input
第一行一个 T 表示数据组数
接下来T行 每行两个正整数 表示N、M
Output
T行 每行一个整数 表示第i组数据的结果
Sample Input
1
4 5
4 5
Sample Output
122
HINT
T <= 10000
N, M<=10000000
Solution
我们先根据BZOJ2154运用莫比乌斯反演推到一个式子，然后优化求解：
Code
1 #include<iostream>
2 #include<string>
3 #include<algorithm>
4 #include<cstdio>
5 #include<cstring>
6 #include<cstdlib>
7 #include<cmath>
8 using namespace std;
9 typedef long long s64;
10
11 const int ONE = 10000005;
12 const int MOD = 100000009;
13
14 int T;
15 int n,m;
16 bool isp[ONE];
17 int prime[700005],p_num;
18 int f[ONE];
19 s64 Ans,sum[ONE];
20
21 int get()
22 {
23 int res=1,Q=1; char c;
24 while( (c=getchar())<48 || c>57)
25 if(c=='-')Q=-1;
26 if(Q) res=c-48;
27 while((c=getchar())>=48 && c<=57)
28 res=res*10+c-48;
29 return res*Q;
30 }
31
32 void Getf(int MaxN)
33 {
34 f[1] = 1;
35 for(int i=2; i<=MaxN; i++)
36 {
37 if(!isp[i])
38 prime[++p_num] = i, f[i] = (-(s64)i*i%MOD+i+MOD)%MOD;
39 for(int j=1; j<=p_num, i*prime[j]<=MaxN; j++)
40 {
41 isp[i * prime[j]] = 1;
42 if(i % prime[j] == 0)
43 {
44 f[i * prime[j]] = (s64)f[i] * prime[j] % MOD;
45 break;
46 }
47 f[i * prime[j]] = (s64)f[i] * f[prime[j]] % MOD;
48 }
49 }
50 for(int i=1; i<=MaxN; i++)
51 sum[i] = (sum[i-1] + f[i]) % MOD;
52 }
53
54 s64 Sum(int n,int m)
55 {
56 return ((s64)n*(n+1)/2%MOD) * ((s64)m*(m+1)/2%MOD) % MOD;
57 }
58
59 void Solve()
60 {
61 n=get(); m=get();
62 if(n > m) swap(n,m);
63 Ans = 0;
64 for(int i=1, j=0; i<=n; i=j+1)
65 {
66 j = min(n/(n/i), m/(m/i));
67 Ans += Sum(n/i,m/i) * ((s64)sum[j] - sum[i-1] + MOD) % MOD;
68 Ans %= MOD;
69 }
70 printf("%lld
",Ans);
71 }
72
73 int main()
74 {
75 Getf(ONE-1);
76 T=get();
77 while(T--)
78 Solve();
79 }