Millar-Rabin + Pollard-rho + 矩阵乘法 一道比较基础的递推题...... 相关题目:「数学」约数个数和

题目描述

小 F 在数学课上与同学闲聊。

给你一个正整数,计算他的约数个数。

「你们竞赛就学这个啊?太简单了吧。」

「……」

给你一个正整数,计算他的所有约数的约数个数和。

「我想想……嗯,还是不难的。诶,你们竞赛还能报名吗?」

「……」

给你一个正整数,计算他的所有约数的所有约数的约数个数和。

「哎呀反正你们电脑总是能爆算出来的嘛,快给我说在哪报名。」

「……」

给你一个正整数,计算他的所有约数的所有约数的所有约数的约数个数和。

「有完没完了?」

「……」

被嘲讽的小 F 将这个问题交给了你,请展示你的爆算实力。

给你一个正整数 NN ,请计算 NN(( 所有约数的 )×K)\times K 约数个数和。

答案可能很大,请输出对 998244353998244353 取模的结果。

输入输出格式

输入格式: 输入一行两个整数 N,K(1N1018,0K1018)N, K(1 \leq N \leq 10 ^ {18}, 0 \leq K \leq 10^{18})

输出格式:

输出一行一个非负整数,表示所求的答案对 998244353998244353 取模的结果。

输入样例#1:

4 0

输出样例#1:

3

输入样例#2:

4 1

输出样例#2:

6

输入样例#3:

4 2

输出样例#3:

10

说明

样例 1,2,31, 2, 3 解释:

4, 0:44,\ 0:4 的约数  1 2 4\ 1\ 2\ 4 4, 1:44,\ 1:4 的所有约数的约数  (1) (1 2) (1 2 4)\ (1)\ (1\ 2)\ (1\ 2\ 4) 4, 2:44,\ 2:4 的所有约数的所有约数的约数 ((1)) ((1) (1 2)) ((1) (1 2) (1 2 4))((1))\ ((1)\ (1\ 2))\ ((1)\ (1\ 2)\ (1\ 2\ 4))

子任务

子任务 1(11pts):N,K1041(11 \mathrm{pts}) : N, K \leq 10 ^ 4 子任务 2(31pts):N1042(31 \mathrm{pts}) : N \leq 10 ^ 4 子任务 3(41pts):N1093(41 \mathrm{pts}) : N \leq 10 ^ 9 子任务 4(67pts):1N1018,0K10184(67 \mathrm{pts}) : 1 \leq N \leq 10 ^ {18}, 0 \leq K \leq 10^{18}

Solution

讲一下出题人的原意吧。

Subtask #1:暴力 Subtask #2:找递推式 Subtask #3:矩阵乘法 Subtask #4:Millar-Rabin + Pollard-rho + 矩阵乘法

然后发现递推出来的东西和组合数有关,所以直接求组合数也是可以的。

正解: 首先将 NN 分解因数,注意到 。对于朴素的 n\sqrt{n} 的分解方法估计时间复杂度承受不起,我们可以套用 Pollard-rho 进行分解。

将分解出来的每个 pkp^{k} 放在第 kk 位。构造递推矩阵:

(1000011000111001111011111) \begin{pmatrix} 1 & 0 & 0 & 0 & 0\\ 1 & 1 & 0 & 0 & 0\\ 1 & 1 & 1 & 0 & 0\\ 1 & 1 & 1 & 1 & 0\\ 1 & 1 & 1 & 1 & 1 \end{pmatrix}

然后统计 pp 对答案的贡献,对于所有的 pp 的贡献,相乘即为答案。

Code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
#include <iostream>
#include <cstdlib>
#include <climits>
#include <cmath>
#include <algorithm>
#include <cstring>
using namespace std;
#define re register int
#define rep(i,a,b) for(re i=(a);i<=(b);++i)
typedef long long int ll;
const int p = 998244353;

inline ll mul(ll a, ll b, ll p) {
ll d = (ll)((long double)a*b/p);
ll res = (a*b-d*p)%p;
if(res < 0) res += p;
return res;
}
ll qpow(ll a, ll b, ll p) {
ll res = 1;
while(b) {
if(b&1) res = mul(res, a, p);
a = mul(a, a, p);
b >>= 1;
}
return res;
}

namespace Millar_Rabin {
const int p[] = {0,2,3,5,7,11,13,17,19,23,29,31,37,41};
const int num = 13;

bool witness(int pri, ll z, int tim, long long p) {
ll res = qpow(pri, z, p);
while(tim--) {
if(res==1) return 0;
if(res==p-1) return 0;
res = mul(res, res, p);
}
return 1;
}
bool main(long long n) {
int tim=0; ll z = n-1;
rep(i,1,num)
if(n==p[i]) return 1;
else if(n%p[i]==0) return 0;
while(!(z&1))
z >>= 1, ++tim;
rep(i,1,num)
if(witness(p[i], z, tim, n)) return 0;
return 1;
}
}

ll q[62]; int len=0;
namespace Pollard_rho {
inline ll f(ll x, ll a, ll p) {return (mul(x,x,p)+a)%p;}

void sol(long long n) {
if(Millar_Rabin::main(n)) {
q[++len] = n;
return ;
}
ll a, b, c, g;
while(1) {
c = rand()%n;
a = b = rand()%n;
b = f(b, c, n);
while(a!=b) {
g = __gcd(abs(a-b), n);
if(g > 1) {
sol(g), sol(n/g);
return ;
}
a = f(a, c, n);
b = f(f(b, c, n), c, n);
}
}
}
void main(long long n) {
sol(n); sort(q+1, q+len+1);
}
}

class Matrix{
public: int a[60][60];
Matrix() {memset(a, 0, sizeof a);}
inline void init() {
rep(i,0,59) a[i][i] = 1;
}
inline long long count() {
long long res=0;
rep(i,0,59) {
res+=a[0][i];
if(res >= p) res-=p;
}
return res;
}
inline Matrix operator * (const Matrix& rhs)const {
Matrix Mat;
rep(i,0,59)
rep(j,0,59)
rep(k,0,59)
(Mat.a[i][j]+=1ll*a[i][k]*rhs.a[k][j]%p)%=p;
return Mat;
}
inline void operator *= (const Matrix& rhs) {
*this = *this * rhs;
}
};
Matrix quick_pow(Matrix a, long long b) {
Matrix Mat; Mat.init();
while(b) {
if(b&1) Mat *= a;
a *= a, b >>= 1;
}
return Mat;
}
long long n, k;

inline long long change(int cnt, ll k) {
Matrix res; res.a[0][cnt] = 1;
Matrix skd;
rep(i,0,59)
rep(j,i,59) skd.a[j][i] = 1;
res *= quick_pow(skd, k);

return res.count();
}

int main() {
srand(1008611);
cin >> n >> k;
if(n==1) {
cout << 1 << endl;
return 0;
}
Pollard_rho::main(n);
int cnt=1; long long ans=1;
rep(i,2,len+1)
if(q[i]==q[i-1]) ++cnt;
else (ans*=change(cnt, k+1))%=p, cnt=1;
cout << ans << endl;
return 0;
}