A sequence of integer lbrace a_n brace{an} can be expressed as:
displaystyle a_n = left{ egin{array}{lr} 0, & n=0\ 2, & n=1\ frac{3a_{n-1}-a_{n-2}}{2}+n+1, & n>1 end{array} ight.an=⎩⎨⎧0,2,23an−1−an−2+n+1,n=0n=1n>1
Now there are two integers nn and mm. I'm a pretty girl. I want to find all b_1,b_2,b_3cdots b_pb1,b2,b3⋯bp that 1leq b_i leq n1≤bi≤n and b_ibiis relatively-prime with the integer mm. And then calculate:
displaystyle sum_{i=1}^{p}a_{b_i}i=1∑pabi
But I have no time to solve this problem because I am going to date my boyfriend soon. So can you help me?
Input
Input contains multiple test cases ( about 1500015000 ). Each case contains two integers nn and mm. 1leq n,m leq 10^81≤n,m≤108.
Output
For each test case, print the answer of my question(after mod 1,000,000,0071,000,000,007).
Hint
In the all integers from 11 to 44, 11 and 33 is relatively-prime with the integer 44. So the answer is a_1+a_3=14a1+a3=14.
样例输入
4 4
样例输出
14
SOLUTION:
考虑容斥，用所有的和减去不合法的和
也就是减去ai，其中i和m的gcd不为1
考虑对m进行质因子分解
对于m的每一质因子p我们需要减去小标为p的倍数的a
写出来之后发现可以推出来式子o（1）的进行计算，但是俩个质因子p可能
筛掉一个ai多次，加上容斥就行了
CODE:
1 #include <bits/stdc++.h>
2 using namespace std;
3 typedef long long ll;
4 const ll mod = 1e9+7;
5 const ll inv6 = 166666668;
6 const ll inv2 = 500000004;
7 ll a[10005];
8 ll clac(ll n,ll i){
9 n /= i;
10 return (n % mod * (n + 1) % mod * (2 * n + 1) % mod * inv6 % mod * i % mod * i % mod + n % mod * (n + 1) % mod * inv2 % mod * i % mod) % mod;
11 }
12 int main(){
13 ll n,m;
14 while(~scanf("%lld%lld",&n,&m)){
15 int cnt = 0;
16 for(int i = 2; i * i <= m; i++){
17 if(m % i == 0){
18 a[cnt++] = i;
19 while(m % i == 0)
20 m /= i;
21 }
22 }
23 if(m != 1)
24 a[cnt++] = m;
25 ll ans = clac(n,1);
26 ll ans2 = 0;
27 for(int i = 1; i < (1 << cnt); i++){
28 int flag = 0;
29 ll tmp = 1;
30 for(int j = 0; j < cnt; j++){
31 if(i & (1 << j)){
32 flag++;
33 tmp = tmp * a[j] % mod;
34 }
35 }
36 tmp = clac(n,tmp);
37 if(flag & 1) ans2 = (ans2 % mod + tmp % mod) % mod;
38 else ans2 = (ans2 % mod - tmp % mod + mod) % mod;
39 }
40 printf("%lld
",(ans % mod - ans2 % mod + mod) % mod);
41 }
42 return 0;
43 }