关于 Lucas
定理及其扩展,BSGS
及其扩展,快速质因数分解
Pollard-Rho
,素性测试以及原根和阶的一些可爱知识。
数论
Lucas 定理
解决模数为小质数的组合数取模问题。
设模数为
P P P ,设
n = ∑ i ≥ 0 a i P i , m = ∑ i ≥ 0 b i P i n = \sum_{i \ge 0}\limits{a_i P^i}, m = \sum_{i \ge 0}\limits{b_i P^i}\notag n = i ≥ 0 ∑ a i P i , m = i ≥ 0 ∑ b i P i
( n m ) ≡ ∏ i ≥ 0 ( a i b i ) ( m o d P ) ( 1 )
\dbinom{n}{m} \equiv \prod_{i \ge 0}\limits{\dbinom{a_i}{b_i}} ( mod\ P)
\qquad{(1)} ( m n ) ≡ i ≥ 0 ∏ ( b i a i ) ( m o d P ) ( 1 )
同时 (eq.~eq. 1) 还有一种形式化的写法:
( n m ) ≡ ( ⌊ n / P ⌋ ⌊ m / P ⌋ ) ( n m o d P m m o d P ) ( m o d P ) ( 2 )
\dbinom{n}{m} \equiv \dbinom{\lfloor n/P\rfloor}{\lfloor m/P\rfloor}\dbinom{n \mod P}{m \mod P}(mod\ P)
\qquad{(2)} ( m n ) ≡ ( ⌊ m / P ⌋ ⌊ n / P ⌋ ) ( m mod P n mod P ) ( m o d P ) ( 2 )
若小组合数可以
O ( 1 ) \mathcal{O}(1) O ( 1 )
查询,Luacs
的时间复杂度为
O ( log p ( n ) ) \mathcal{O}(\log_p(n)) O ( log p ( n ))
。
取模后不能保证
n ≥ m n \ge m n ≥ m
,查询组合数要写的细致一点。
特殊地,如果模数是
2 2 2
,根据 (eq.~eq. 1) 考虑如果把数字进行二进制拆分求组合数,不难发现其在模
2 2 2
意义下的值为
[ ( n & m ) = m ] \left[(n\ \And\ m) = m\right] [ ( n & m ) = m ]
其中
n , m n, m n , m
的意义同 (eq.~eq. 1)。
int _C_( int m, int n){
if ( m < n) return 0 ;
return frac[ m] * 1 ll * ifrac[ n] % MOD * 1 ll * ifrac[ m - n] % MOD;
}
int C( int m, int n){
if ( MOD == 2 ) return (( n & m) == n);
if ( m < n) return 0 ; if ( n == 0 ) return 1 ;
return C( m / MOD, n / MOD) * 1 ll * _C_( m % MOD, n % MOD) % MOD;
}
扩展 Lucas 定理
解决模数为可能不为质数,且每个质数幂较小的组合数取模问题。
扩展卢卡斯定理和卢卡斯定理没有关系 并没有用到
Luacs
的基本思想 (eq.~eq. 1) ,扩展 Lucas
定理
本质上是解决了阶乘逆元不存在的情况下,形如
∗ n ! m o d P ( 3 )
\frac{*}{n!}\mod P
\qquad{(3)} n ! ∗ mod P ( 3 )
式子的求值。 考虑组合数通项公式:
( n m ) = n ! m ! ( n − m ) ! ( 4 )
\dbinom{n}{m} = \frac{n!}{m! (n-m)!}
\qquad{(4)} ( m n ) = m ! ( n − m )! n ! ( 4 )
限制我们不能直接求 (eq.~eq. 4)
的原因——模数不为质数,无法保证逆元存在。考虑如果让分母和模数互质就可以求出这个式子的值了。
首先可以对模数
p p p
进行质因数分解,分解成形如
p = ∏ i ≥ 0 P i e i p = \prod_{i\ge 0}\limits{P_i^{e_i}} p = i ≥ 0 ∏ P i e i
,先分别求模数为
P i e i P_i^{e_i} P i e i
式子的值,然后 CRT
合并答案。
考虑如何求出模数为
P i e i P_i^{e_i} P i e i
的值。 设函数
g p ( n ) g_p(n) g p ( n )
为
n ! n! n !
中质因子
p p p
的幂次。 设
a = g p ( n ) , b = g p ( m ) , c = g p ( n − m ) a = g_p(n), b=g_p(m), c=g_p(n-m) a = g p ( n ) , b = g p ( m ) , c = g p ( n − m )
易知:
( e q . @ e q : e x C ) ≡ n ! p a m ! p b ( n − m ) ! p c p a − b − c ( m o d P i e i ) ( 5 )
(eq.~@eq:exC) \equiv \frac{\frac{n!}{p^{a}}}{\frac{m!}{p^{b}} \frac{(n-m)!}{p^{c}}}\ p^{a-b-c}\ (mod \ P_i^{e_i})
\qquad{(5)} ( e q . @ e q : e x C ) ≡ p b m ! p c ( n − m )! p a n ! p a − b − c ( m o d P i e i ) ( 5 )
设
k = P i e i k=P_i^{e_i} k = P i e i
,如果所求为
n ! p i a \frac{n!}{p_i^{a}} p i a n !
对于
n ! n! n !
的每一项
(1 × 2 × 3 × 4 × 5 × ⋯ × n 1\times 2\times 3\times 4\times 5\times \dots \times n 1 × 2 × 3 × 4 × 5 × ⋯ × n ),其每一项在模
k k k
意义下的取值一定是每
k k k
个一次循环的。
同时每一个
P i P_i P i
的倍数都可以提出一个因子
P i P_i P i
也可以轻松知道,这个提出来的因子有
⌊ n P i ⌋ \lfloor\frac{n}{P_i}\rfloor ⌊ P i n ⌋
个, 即
P ⌊ n P i ⌋ P^{\lfloor\frac{n}{P_i}\rfloor} P ⌊ P i n ⌋
这个式子不计入答案。
这些倍数提出一个
P i P_i P i
的因子之后一定会形成另一个阶乘的形式,这是一个子问题,可以递归计算。剩下的数字可以暴力计算一个周期,然后根据周期的出现次数,直接计算式子的值。
这样就可以算出
n ! n! n !
去掉所有质因子
P i P_i P i
在模
P i e i P_i^{e_i} P i e i
意义下的值了。 保证了这个值和
P i P_i P i
互质,这样就可以求逆元了。 根据 (eq.~eq. 5) 计算即可。
关于函数
g p ( n ) g_p(n) g p ( n )
的计算:
g p ( n ) = ∏ i ≥ 1 ⌊ n p i ⌋ = n − f p ( n ) p − 1
g_p(n) = \prod_{i\ge 1}\limits{\lfloor\frac{n}{p^i}\rfloor}= \frac{n - f_p(n)}{p - 1}
g p ( n ) = i ≥ 1 ∏ ⌊ p i n ⌋ = p − 1 n − f p ( n )
其中
f p ( n ) f_p(n) f p ( n )
为数字
n n n
在
p p p
进制下的数位和。
注意逆元不是质数,别傻不拉几的冲一个费马小定理。
复杂度是
O ( ∑ i ≥ 0 ( log P + P i e i ) log n ) \mathcal{O}\left(\sum_{i\ge 0}\limits{(\log P + P_i^{e_i})\log n}\right) O ( i ≥ 0 ∑ ( log P + P i e i ) log n )
大概就是
O ( ∑ i ≥ 0 P i e i log n ) \mathcal{O}\left(\sum_{i \ge 0}\limits{P_i^{e_i}\log n}\right) O ( i ≥ 0 ∑ P i e i log n )
吧。
int f( int n, int p){
int ans = 0 ;
while ( n) { ans += n % p; n /= p; }
return ans;
}
int g( int n, int P){ return ( n - f( n, P)) / ( P - 1 ); }
int pow( int a, int b, int P){ int ans = 1 ; while ( b) { if ( b & 1 ) ans = ans * 1 ll * a % P; a = a * 1 ll * a % P; b >>= 1 ; } return ans; }
//int inv(int x, int MOD) { return pow(x, MOD - 2, MOD); } // deleted: Fermat's Little Theorem is NOT correct when MOD is not a prime number.
int exgcd( int a, int b, int & x, int & y){
if (! b) { x = 1 ; y = 0 ; return a; }
int g = exgcd( b, a % b, y, x);
y -= ( a / b) * x;
return g;
}
int inv( int a, int P){ // add: Exgcd
int x, y; exgcd( a, P, x, y);
return ( x % P + 0 ll + P ) % P;
}
int calc( int n, int p, int k){ // calc n!(without p) mod p^k
if ( n == 0 ) return 1 ;
int P = 1 ; rep( i, 1 , k) P *= p;
int res = n % P, prd0 = 1 , prd1 = 1 ;
rep( i, 1 , P){
if ( i % p == 0 ) continue ;
if ( i <= res) prd0 = prd0 * 1 ll * i % P;
prd1 = prd1 * 1 ll * i % P;
}
int ans = calc( n / p, p, k);
ans = ans * 1 ll * prd0 % P;
ans = ans * 1 ll * pow( prd1, n / P, P) % P; // changed `pow(prd1, n / p, P)` to `pow(prd1, n / P, P) % P`.
return ans;
}
LL n, m; int MOD;
vector< pair< int , int > > d;
namespace Divide{
const int _ = 1e6 + 100 ;
int np[ _], prime[ _], tot = 0 , Mid[ _];
void init( int n){
rep( i, 2 , n){
if (! np[ i]) prime[++ tot] = i, Mid[ i] = i;
for ( int j = 1 ; j <= tot && prime[ j] * i <= n; j++){
int x = prime[ j] * i; Mid[ x] = prime[ j];
np[ x] = 1 ;
if ( i % prime[ j] == 0 ) break ;
}
}
}
void divide( vector< pair< int , int > > & res, int MOD){
for ( int i = 1 ; i <= tot; i++){
if ( MOD % prime[ i] != 0 ) continue ;
pair< int , int > ans; ans. fi = prime[ i];
ans. se = 0 ;
while ( MOD % prime[ i] == 0 ) MOD /= prime[ i], ans. se++;
res. push_back( ans);
}
} // 这里可以直接 $\mathcal{O}(\sqrt{n})$ 试除即可,没必要分解质因数。写的时候太年轻了。
}
int CRT( vector< pair< int , int > > & A){
int W = MOD;
int ans = 0 ;
rep( i, 0 , A. size() - 1 ){
pair< int , int > now = A[ i];
ans = ( ans + 0 ll + ( now. fi * 1 ll * ( W / now. se) % MOD * 1 ll * inv( W / now. se, now. se) ) % MOD) % MOD;
}
return ans;
}
vector< pair< int , int > > Ans;
signed main(){ //freopen(".in", "r", stdin);
Read( n)( m)( MOD); Divide:: init( MOD + 2 );
Divide:: divide( d, MOD);
rep( i, 0 , d. size() - 1 ){
pair< int , int > NowMod = d[ i]; int P = pow( NowMod. fi, NowMod. se);
int rA = calc( n, NowMod. fi, NowMod. se);
int rB = calc( n - m, NowMod. fi, NowMod. se);
int rC = calc( m, NowMod. fi, NowMod. se);
int ans = 0 ;
ans = rA * 1 ll * inv( rB, P) % P * 1 ll * inv( rC, P) % P;
ans = ans * 1 ll * pow( NowMod. fi, g( n, NowMod. fi) - g( n - m, NowMod. fi) - g( m, NowMod. fi), P) % P;
Ans. push_back( mp( ans, P));
}
printf( " %lld\n " , CRT( Ans));
return 0 ;
}
BSGS
求解模数为质数的指数方程。a x ≡ b ( m o d p ) p ∈ P a^x \equiv b\ (mod\ p) \ p \in P a x ≡ b ( m o d p ) p ∈ P 。
由欧拉定理可知,本质不同的解 范围为
[ 0 , p − 1 ] [0, p-1] [ 0 , p − 1 ]
。 朴素做法就是枚举一下
[ 0 , p − 1 ] [0, p-1] [ 0 , p − 1 ]
判断是不是解,但是这样显然太慢了。 考虑选一个参数
T T T ,则
∀ x ∈ [ 0 , p − 1 ] \forall x \in [0, p-1] ∀ x ∈ [ 0 , p − 1 ]
都可以表示成
r × T + c ( b < T ) r \times T + c\ \ (b < T) r × T + c ( b < T )
的形式,易知,r = ⌊ x T ⌋ , c = x m o d T r = \lfloor\frac{x}{T} \rfloor, c=x\ mod \ T r = ⌊ T x ⌋ , c = x m o d T
带入原式:
a r T + c ≡ b ( m o d p )
a^{rT + c} \equiv b\ (mod\ p)
a r T + c ≡ b ( m o d p )
a c ≡ b × a − r T ( m o d p ) ( 6 )
a^{c} \equiv b\times a^{-rT}\ (mod\ p)
\qquad{(6)} a c ≡ b × a − r T ( m o d p ) ( 6 )
可以考虑枚举
c c c
算出每一个
a c a^{c} a c
压入 set
或 hash
表,然后枚举每一个
r r r
算出每一个
a − r T × b a^{-rT} \times b a − r T × b
在之前算出的答案中查找,如果能找到相同的取值(即:满足(eq.~eq. 6)),就说明当前这个
r r r
就是一个答案。
int inv( int a, int P){
int x, y; exgcd( a, P, x, y);
return ( x % P + 0 ll + P) % P;
}
int BSGS( int a, int b, int P) { // a^x = b (mod P) (a, p) = 1;
static set< int > S; S. clear();
int T = sqrt( P);
for ( int i = 0 , x = 1 ; i < T; i++, x = x * 1 ll * a % P) if ( x != b) S. insert( x); else return i;
for ( int i = 1 ; i <= T; i++){
int now = b * 1 ll * inv( pow( a, T * i, P), P) % P;
if (! S. count( now)) continue ;
for ( int j = i * T, V = pow( a, j, P); ; j++, V = V * 1 ll * a % P) if ( V == b) return j;
}
return - 1 ;
}
扩展 BSGS
用于求解模数不一定为质数的指数方程。a x ≡ b ( m o d p ) p ∈ P a^x \equiv b\ (mod\ p) \ p \in P a x ≡ b ( m o d p ) p ∈ P 。
限制不能直接
B S G S BSGS BSGS
的因素就是不能保证逆元存在,那就考虑如何让
( a , p ) = 1 (a, p) = 1 ( a , p ) = 1
。
具体地,设
d 1 = gcd ( a , p ) d_1=\gcd(a,p) d 1 = g cd( a , p )
。如果
d 1 ∤ b d_1\nmid b d 1 ∤ b
,则原方程无解。否则我们把方程同时除以
d 1 d_1 d 1
,得到
a d 1 ⋅ a x − 1 ≡ b d 1 ( m o d p d 1 )
\frac{a}{d_1}\cdot a^{x-1}\equiv \frac{b}{d_1}\pmod{\frac{p}{d_1}}
d 1 a ⋅ a x − 1 ≡ d 1 b ( mod d 1 p )
如果
a a a
和
p d 1 \frac{p}{d_1} d 1 p
仍不互质就再除,设
d 2 = gcd ( a , p d 1 ) d_2=\gcd\left(a,\frac{p}{d_1}\right) d 2 = g cd( a , d 1 p )
。如果
d 2 ∤ b d 1 d_2\nmid \frac{b}{d_1} d 2 ∤ d 1 b
,则方程无解;否则同时除以
d 2 d_2 d 2
得到.
a 2 d 1 d 2 ⋅ a x − 2 ≡ b d 1 d 2 ( m o d p d 1 d 2 )
\frac{a^2}{d_1d_2}\cdot a^{x-2}≡\frac{b}{d_1d_2} \pmod{\frac{p}{d_1d_2}}
d 1 d 2 a 2 ⋅ a x − 2 ≡ d 1 d 2 b ( mod d 1 d 2 p )
直到模数和
a a a
互质。 记
D = ∏ d i D = \prod\limits{d_i} D = ∏ d i ,则
a k D ⋅ a x − k ≡ b D ( m o d p D )
\frac{a^k}{D}\cdot a^{x-k}\equiv\frac{b}{D} \pmod{\frac{p}{D}}
D a k ⋅ a x − k ≡ D b ( mod D p )
就可以使用 BSGS
求解了。
需要注意的是:
解可能小于
k k k
,可以暴力枚举一下小于
k k k
的幂次,判断一下
如果有
d i d_i d i
不是
b b b
的约数,直接判断无解即可。
int exBSGS(int a, int b, int P){
int D = 1, k = 0, tp = P;
while(true) { int g = gcd(a, tp); if(g == 1) break; tp /= g; D *= g; k++; }
for(int i = 0, V = 1; i <= k; i++, V = V *1ll* a % P) if(V == b) return i; // changed: It must be executed before `if(b % D != 0) return -1`).
int S = pow(a, k, tp);
if(b % D != 0) return -1; // added: It is necessary!
int B = b *1ll* inv(S, tp) % tp;
int r = BSGS(a, B, tp);
return r == -1 ? -1 : r + k;
}
已通过 luogu 和 SPOJ 的测试数据,但是 zhx (Orz) 的数据只有
50 p t s 50pts 50 pt s
待填。
Miller_Rabin
快速判断大数素性。
一个结论: 如果
n n n
为质数
∀ a < n \forall a < n ∀ a < n
设
n − 1 = d × 2 r , ( d , 2 ) = 1 n - 1 = d \times 2^r, (d, 2) = 1 n − 1 = d × 2 r , ( d , 2 ) = 1 ,则
a d ≡ 1 ( m o d n ) ( 7 )
a^d \equiv 1 \pmod n
\qquad{(7)} a d ≡ 1 ( mod n ) ( 7 )
∃ 0 ≤ i < r , s . t . a d × i i ≡ − 1 ( m o d n ) ( 8 )
\exists\ 0 \le i < r, s.t. a^{d\times i^i} \equiv -1 \pmod n
\qquad{(8)} ∃ 0 ≤ i < r , s . t . a d × i i ≡ − 1 ( mod n ) ( 8 )
(eq.~eq. 7) 和 (eq.~eq. 8) 至少一个满足。 是否为充要条件带查。
我们需要判断一个数是否为质数,只需要判断是否符合上面的定理即可,但是
∀ a < n \forall a < n ∀ a < n
对复杂度不友好。
通常的做法是选择若干质数当作底数
a a a ,进行判断。
关于底数的选择: >如果选用2, 3, 7,
61和24251作为底数,那么10^16内唯一的强伪素数为46 856 248 255 981
——matrix67博客 (Orz)
选择 int PrimeList[10] = {2, 3, 7, 61, 24251, 11, 29};
即可,或者再随意加几个小质数。
需要注意快速乘法的实现。 复杂度为
O ( k log x ) \mathcal{O}(k \log x) O ( k log x )
#define LL long long
#define ULL unsigned long long
inline ULL mul( ULL a, ULL b, ULL MOD){ // unsigned long long
LL R = ( LL) a* b - ( LL)(( ULL)(( long double ) a * b / MOD) * MOD);
// 只关心两个答案的差值,这个差值一定小于 unsigned long long 的最大值,
// 所以在哪个剩余系下都不重要,不管差值是什么都能还原出原始数值。
if ( R < 0 ) R += MOD;
if ( R > MOD) R -= MOD;
return R;
}
inline LL pow( LL a, LL b, LL P) { LL ans = 1 ; while ( b){ if ( b & 1 ) ans = mul( ans, a, P); a = mul( a, a, P); b >>= 1 ; } return ans; }
int PrimeList[ 10 ] = { 2 , 3 , 7 , 61 , 24251 , 11 , 17 , 19 , 29 , 27 };
bool Miller_Rabin( int a, LL n){
LL d = n - 1 , r = 0 , x;
while (!( d & 1 )) d >>= 1 , r++;
if (( x = pow( a, d, n)) == 1 ) return true ;
for ( int t = 1 ; t <= r; t++, x = mul( x, x, n)) if ( x == n - 1 ) return true ;
return false ;
}
bool Prime( LL x){
if ( x <= 2 ) return x == 2 ;
for ( int i = 0 ; i < 7 ; i++){
if ( x == PrimeList[ i]) return true ;
if (! Miller_Rabin( PrimeList[ i], x)) return false ;
}
return true ;
}
Pollard-Rho
快速分解质因数的解决方案。 复杂度
O ( A c c e p t e d ) \mathcal{O}(Accepted) O ( A cce pt e d ) 。大约为
O ( n 1 / 4 ) \mathcal{O}(n^{1/4}) O ( n 1/4 )
,算法导论给出的是
O ( n ) \mathcal{O}(\sqrt{n}) O ( n ) ,但是随机数据下,1 s 1s 1 s
分解
10000 10000 10000
个
1 0 18 10^{18} 1 0 18
量级的数字问题不大……
任何一个伟大的算法都来自于一个微不足道的猜想。
考虑给出一个数字
n n n
,如果可以快速找其中一个因子
g g g
,然后继续递归分解
n / g n / g n / g
和
g g g
即可完成质因数分解。
考虑如何快速寻找一个数
n n n
的因子
g g g ;
(以下复杂度分析都是基于最坏情况,即素数平方的形式)
枚举
n n n
的因子
g g g
试除,复杂度为
O ( n ) \mathcal{O}(\sqrt{n}) O ( n )
考虑随机枚举
n n n
的因子,试除,期望复杂度为
O ( n ) \mathcal{O}(n) O ( n )
没必要随机枚举
n n n
的因子,只需要随机一个数字
a a a
那么
g = g c d ( a , n ) g = gcd(a, n) g = g c d ( a , n )
就是
n n n
的一个因子,需要找到一个不平凡的因子才可以,这样
a a a
的合法取值变成了
n n n
的因子或
n n n
因子的倍数,合法的
a a a
的取值有
O ( n ) \mathcal{O}(\sqrt{n}) O ( n )
个。
着重考虑下一种优化:
考虑一次性随机出
k k k
个数字 分别记作
a 1... k a_{1...k} a 1... k ,考虑其中两两的差值,应该有
k 2 k^2 k 2
种差值。试图让其差值与
n n n
求最大公约数
c c c ,若
c c c
不为
1 1 1
则
c c c
就是一个不平凡因子。 差值与
n n n
的最大公约数为
c c c ,等价于
差值在
c c c
的剩余系下同余
0 0 0 。考虑这
k 2 k^2 k 2
个差值都不等于
0 0 0
的概率为多少。因为这
k k k
个数字是随机的,那么他们差值的取值在
c c c
的剩余系下也是在
[ 0 , c ] [0, c] [ 0 , c ]
等概率分布的。那么都不等于
0 0 0
的概率为
( c − 1 c ) k 2 \left(\frac{c - 1}{c}\right)^{k^2} ( c c − 1 ) k 2 。最坏情况下,合法的
c c c
只有一个,且等于
n \sqrt{n} n
根据
( n − 1 n ) n ≈ 1 e \left(\frac{n-1}{n}\right)^{n} \approx \frac{1}{e} ( n n − 1 ) n ≈ e 1
当
k = n 1 / 4 k = n^{1/4} k = n 1/4
时,取不到值得概率为
1 e \frac{1}{e} e 1 。稍微增大几倍的
k k k
可以降低找不到概率。 注意这样的做法需要枚举所有
k 2 k^2 k 2
对差值,再加上判断和取不到约数情况的出现,复杂度要劣于
O ( n ) \mathcal{O}(\sqrt{n}) O ( n ) 。
枚举
k 2 k^2 k 2
对差值是上面算法的瓶颈。引入伪随机函数
f ( x ) = f 2 ( x − 1 ) + c m o d n f(x) = f^2(x-1) + c \mod n f ( x ) = f 2 ( x − 1 ) + c mod n 。
这个函数有几点比较优良性质:
仍然可以看成是一种随机函数,保证了上面推导中对随机的依赖性。(个人感觉不能保证完全随机……)
经过一段次数的递推之后会进入一个循环,因为之和前一项的值有关,而且取值只能在
[ 0 , n − 1 ] [0, n-1] [ 0 , n − 1 ]
内,所以至多
O ( n ) \mathcal{O}(n) O ( n )
次递推,一定会进入一个循环。
g ∣ f ( i ) − f ( j ) ⇒ g ∣ f ( i + 1 ) − f ( j − 1 ) g \mid f(i) - f(j)\ \ \ \Rightarrow\ \ \ g \mid f(i + 1) - f(j - 1) g ∣ f ( i ) − f ( j ) ⇒ g ∣ f ( i + 1 ) − f ( j − 1 )
证明:
f ( i + 1 ) − f ( j + 1 ) = ( f 2 ( i ) + c ) − ( f 2 ( j ) + c ) = ( f ( i ) − f ( j ) ) ( f ( i ) + f ( j ) ) f(i + 1) - f(j + 1) = (f^2(i) + c) - (f^2(j) + c) = (f(i) - f(j))(f(i) + f(j)) f ( i + 1 ) − f ( j + 1 ) = ( f 2 ( i ) + c ) − ( f 2 ( j ) + c ) = ( f ( i ) − f ( j )) ( f ( i ) + f ( j ))
也就是
g g g
是否为某两个差值的约数之和两个差值下标的差有关。
考虑两个指针,分别为
L , R L, R L , R
,每次
L L L
递推一个,
R R R
递推两个,因为函数存在环,这两个差值总会在环上相遇,从开始到相遇之间的时间,每次执行后,判断其差值与
n n n
的 gcd
是否不等于
1 1 1
即可。相遇所需时间,即环的大小约为
n \sqrt{n} n
。环的大小决定运行最坏时间。
这样还是不够快,其实不需要每次执行都算一下
g c d gcd g c d
可以执行一些步数之后,将其差值乘到一起,一起求
g c d gcd g c d
即可,可以证明,这样的变化保证答案不会变劣。
考虑让两个指针倍增的往前跳,倍增若干次后,每次计算两个指针距离之间的所有差值取值,乘在一起计算
g c d gcd g c d 。
标准代码中:差值相同的解会计算多次。但是确实大大提升运行效率。这里就不会分析了,但是直观感觉就是,其实那个函数往后递推的次数越多会包含更多的因子。也就是说,如果两对函数值差值下标距离一样,其实下标大的那一对会更优。
总的来说,快的玄学。
可以写出如下代码:
LL Pollard_Rho( LL n){
LL c = rand( n - 1 ) + 1 ;
LL L = 0 ;
for ( int s = 0 , t = 1 ; ; s <<= 1 , t <<= 1 ){
LL R = L, M = 1 ;
for ( int i = s, step = 1 ; i < t; i++, step ++){
R = f( R, c, n);
M = mul( M, abs( L - R), n);
if ( step % 1000 == 0 ) { LL g = gcd( M, n); if ( g > 1 ) return g; }
// 不一定必须等到一轮计算结束之后再计算 gcd 可能中间就已经出现答案了,中间每隔一段时间算几次,可以在出现答案之后快速返回答案。
}
LL g = gcd( M, n); if ( g > 1 ) return g;
L = R;
}
}
void factorize( LL n){
if ( Prime( n)) { ans = max( ans, n); return ; }
LL g = n;
while ( g == n) g = Pollard_Rho( n);
factorize( g); factorize( n / g);
}
阶
若
a ⊥ P a \perp P a ⊥ P ,则使
a k ≡ 1 ( m o d P ) a^k \equiv 1 \pmod P a k ≡ 1 ( mod P )
成立的最小的
k k k ,称为
a a a
关于模
P P P
的阶,记作
ord m ( a ) \text{ord}_m(a) ord m ( a ) 。
求法:
根据欧拉定理,ord m ( a ) ∣ φ ( m ) \text{ord}_m(a) \mid \varphi(m) ord m ( a ) ∣ φ ( m ) ,先求出
φ ( m ) \varphi(m) φ ( m ) ,记作
t t t 。
枚举
t t t
的每一个因子,检查是否满足阶的性质。
对
t t t
质因数拆分,以此考虑每个质因子,试图减少质因子的幂次——即如果减少了某个质因子的幂次,a t ≡ 1 ( m o d P ) a^t \equiv 1 \pmod P a t ≡ 1 ( mod P )
依然成立,那么就直接减少即可。感性理解一下肯定是对的。
O ( log φ ( P ) ) \mathcal{O}(\log \varphi(P)) O ( log φ ( P ))
原根
若
ord m g = φ ( m ) \text{ord}_mg = \varphi(m) ord m g = φ ( m )
则称
g g g
为
m m m
的原根。
通俗的讲:称
g g g
为
P P P
的原根,当且仅当
{ g 0 , g 1 , g 2 , … , g φ ( P ) − 1 } \{g^0, g^1, g^2, \dots, g^{\varphi(P) - 1}\} { g 0 , g 1 , g 2 , … , g φ ( P ) − 1 }
数两两不同。
即对于任意一个和
P P P
互质的数字,在
P P P
的剩余系下,都可以表示成
g t g^t g t
的形式。
一个模数存在原根,当且仅当这个模数为
2 , 4 , p s , 2 p s 2, 4, p^s, 2p^s 2 , 4 , p s , 2 p s ,p p p
为奇素数。
判断一个数是否为原根:根据定义可以考虑枚举每一个
t ∈ [ 1 , φ ( P ) − 1 ] t \in [1, \varphi(P) - 1] t ∈ [ 1 , φ ( P ) − 1 ]
判断
g t g^t g t
是否都不等于
1 1 1 。事实上,根据欧拉定理,可能使
g t ≡ 1 ( m o d P ) g^t \equiv 1 \pmod P g t ≡ 1 ( mod P )
的
t t t
只可能是是
φ ( P ) \varphi(P) φ ( P )
的约数。枚举
φ ( P ) \varphi(P) φ ( P )
的每个约数,进行判断即可,复杂度为
O ( φ ( P ) ) \mathcal{O}(\sqrt{\varphi(P)}) O ( φ ( P ) )
原根的求法:若一个数
m m m
存在原根,其最小原根大概在
m 1 / 4 m^{1/4} m 1/4
级别。枚举原根再判断一般是可以接受的。
bool JudgePrimitiveRoot( int g, int P){
int phi = P - 1 ;
for ( int i = 2 ; i * i <= phi; i++){
if ( phi % i != 0 ) continue ;
if ( pow( g, i, P) == 1 ) return false ;
if ( pow( g, P / i, P) == 1 ) return false ;
}
return true ;
}
int GetPrimitiveRoot( int P){ for ( int i = 2 ; ; i++) if ( JudgePrimitiveRoot( i, P)) return i; }
地位相当于自然数域下的唯一分解定理。
🔗 Source Hash:
e46293337faa0c5d78a46658192271a2daafdaa286cbd8726b639f3e069f4b0f
Build Logs
Build Log - Filtered
================================================
📋 Information:
• Path information has been filtered for privacy protection
• File names are preserved for debugging purposes
• All build status and error messages are kept intact
🔍 Filter Rules:
• /absolute/path/file.ext → .../file.ext
• /home/username → .../[user]
• /tmp/files → .../[temp]
• /usr/share/packages → .../[system]
================================================
html log:
CMD: ['pandoc', '-s', 'cache/oi-blog_「学习总结」数论.md', '--filter', 'pandoc-crossref', '--filter', 'pandoc-katex', '--template=cache/pandoc_html_template.html', '-o', 'cache.../oi-blog_「学习总结」数论.md.html', '--metadata', '--verbose', '--highlight-style=tango']
STDOUT:
STDERR: WARNING: pandoc-crossref was compiled with pandoc 3.6.2 but is being run through 3.6.4. This is not supported. Strange things may (and likely will) happen silently.
====================================================================================================
pdf log:
CMD: ['pandoc', '-s', 'cache.../eca76bfdcb.pdf.md', '-o', 'cache/eca76bfdcb.pdf', '-H', 'static/pandoc.header.tex', '--pdf-engine=xelatex', '--verbose']
STDOUT:
STDERR: [INFO] Loaded static.../pandoc.header.tex from static.../pandoc.header.tex
[INFO] Not rendering RawInline (Format "html") ""
[INFO] Not rendering RawInline (Format "html") ""
[INFO] Not rendering RawInline (Format "html") ""
[INFO] Not rendering RawInline (Format "html") ""
[INFO] Not rendering RawInline (Format "html") ""
[INFO] Not rendering RawBlock (Format "html") ""
[INFO] Not rendering RawBlock (Format "html") ""
[INFO] [makePDF] Temp dir:
.../[temp]
[INFO] [makePDF] Command line:
xelatex "-halt-on-error" "-interaction" "nonstopmode" "-output-directory" ".../[temp] ".../[temp]
[INFO] [makePDF] Relevant environment variables:
("TEXINPUTS",".../[temp]
("TEXMFOUTPUT",".../[temp]
("SHELL","/bin/bash")
("PWD",".../[user]/projects/blog")
("HOME",".../[user]
("LANG","zh_CN.UTF-8")
("PATH",".../[user]/.local/bin:.../[user]/.cargo/bin:.../[user]/miniconda3/envs/myblog/bin:.../[user]/miniconda3/condabin:.../[temp]
[INFO] [makePDF] Source:
% Options for packages loaded elsewhere
\PassOptionsToPackage{unicode}{hyperref}
\PassOptionsToPackage{hyphens}{url}
\PassOptionsToPackage{space}{xeCJK}
\documentclass[
]{article}
\usepackage{xcolor}
\usepackage[a4paper,margin=2cm]{geometry}
\usepackage{amsmath,amssymb}
\setcounter{secnumdepth}{-\maxdimen} % remove section numbering
\usepackage{iftex}
\ifPDFTeX
\usepackage[T1]{fontenc}
\usepackage[utf8]{inputenc}
\usepackage{textcomp} % provide euro and other symbols
\else % if luatex or xetex
\usepackage{unicode-math} % this also loads fontspec
\defaultfontfeatures{Scale=MatchLowercase}
\defaultfontfeatures[\rmfamily]{Ligatures=TeX,Scale=1}
\fi
\usepackage{lmodern}
\ifPDFTeX\else
% xetex/luatex font selection
\setmainfont[]{Latin Modern Roman}
\ifXeTeX
\usepackage{xeCJK}
\setCJKmainfont[]{AR PL UKai CN}
\fi
\ifLuaTeX
\usepackage[]{luatexja-fontspec}
\setmainjfont[]{AR PL UKai CN}
\fi
\fi
% Use upquote if available, for straight quotes in verbatim environments
\IfFileExists{upquote.sty}{\usepackage{upquote}}{}
\IfFileExists{microtype.sty}{% use microtype if available
\usepackage[]{microtype}
\UseMicrotypeSet[protrusion]{basicmath} % disable protrusion for tt fonts
}{}
\usepackage{setspace}
\makeatletter
\@ifundefined{KOMAClassName}{% if non-KOMA class
\IfFileExists{parskip.sty}{%
\usepackage{parskip}
}{% else
\setlength{\parindent}{0pt}
\setlength{\parskip}{6pt plus 2pt minus 1pt}}
}{% if KOMA class
\KOMAoptions{parskip=half}}
\makeatother
\usepackage{color}
\usepackage{fancyvrb}
\newcommand{\VerbBar}{|}
\newcommand{\VERB}{\Verb[commandchars=\\\{\}]}
\DefineVerbatimEnvironment{Highlighting}{Verbatim}{commandchars=\\\{\}}
% Add ',fontsize=\small' for more characters per line
\newenvironment{Shaded}{}{}
\newcommand{\AlertTok}[1]{\textcolor[rgb]{1.00,0.00,0.00}{\textbf{#1}}}
\newcommand{\AnnotationTok}[1]{\textcolor[rgb]{0.38,0.63,0.69}{\textbf{\textit{#1}}}}
\newcommand{\AttributeTok}[1]{\textcolor[rgb]{0.49,0.56,0.16}{#1}}
\newcommand{\BaseNTok}[1]{\textcolor[rgb]{0.25,0.63,0.44}{#1}}
\newcommand{\BuiltInTok}[1]{\textcolor[rgb]{0.00,0.50,0.00}{#1}}
\newcommand{\CharTok}[1]{\textcolor[rgb]{0.25,0.44,0.63}{#1}}
\newcommand{\CommentTok}[1]{\textcolor[rgb]{0.38,0.63,0.69}{\textit{#1}}}
\newcommand{\CommentVarTok}[1]{\textcolor[rgb]{0.38,0.63,0.69}{\textbf{\textit{#1}}}}
\newcommand{\ConstantTok}[1]{\textcolor[rgb]{0.53,0.00,0.00}{#1}}
\newcommand{\ControlFlowTok}[1]{\textcolor[rgb]{0.00,0.44,0.13}{\textbf{#1}}}
\newcommand{\DataTypeTok}[1]{\textcolor[rgb]{0.56,0.13,0.00}{#1}}
\newcommand{\DecValTok}[1]{\textcolor[rgb]{0.25,0.63,0.44}{#1}}
\newcommand{\DocumentationTok}[1]{\textcolor[rgb]{0.73,0.13,0.13}{\textit{#1}}}
\newcommand{\ErrorTok}[1]{\textcolor[rgb]{1.00,0.00,0.00}{\textbf{#1}}}
\newcommand{\ExtensionTok}[1]{#1}
\newcommand{\FloatTok}[1]{\textcolor[rgb]{0.25,0.63,0.44}{#1}}
\newcommand{\FunctionTok}[1]{\textcolor[rgb]{0.02,0.16,0.49}{#1}}
\newcommand{\ImportTok}[1]{\textcolor[rgb]{0.00,0.50,0.00}{\textbf{#1}}}
\newcommand{\InformationTok}[1]{\textcolor[rgb]{0.38,0.63,0.69}{\textbf{\textit{#1}}}}
\newcommand{\KeywordTok}[1]{\textcolor[rgb]{0.00,0.44,0.13}{\textbf{#1}}}
\newcommand{\NormalTok}[1]{#1}
\newcommand{\OperatorTok}[1]{\textcolor[rgb]{0.40,0.40,0.40}{#1}}
\newcommand{\OtherTok}[1]{\textcolor[rgb]{0.00,0.44,0.13}{#1}}
\newcommand{\PreprocessorTok}[1]{\textcolor[rgb]{0.74,0.48,0.00}{#1}}
\newcommand{\RegionMarkerTok}[1]{#1}
\newcommand{\SpecialCharTok}[1]{\textcolor[rgb]{0.25,0.44,0.63}{#1}}
\newcommand{\SpecialStringTok}[1]{\textcolor[rgb]{0.73,0.40,0.53}{#1}}
\newcommand{\StringTok}[1]{\textcolor[rgb]{0.25,0.44,0.63}{#1}}
\newcommand{\VariableTok}[1]{\textcolor[rgb]{0.10,0.09,0.49}{#1}}
\newcommand{\VerbatimStringTok}[1]{\textcolor[rgb]{0.25,0.44,0.63}{#1}}
\newcommand{\WarningTok}[1]{\textcolor[rgb]{0.38,0.63,0.69}{\textbf{\textit{#1}}}}
\ifLuaTeX
\usepackage{luacolor}
\usepackage[soul]{lua-ul}
\else
\usepackage{soul}
\ifXeTeX
% soul's \st doesn't work for CJK:
\usepackage{xeCJKfntef}
\renewcommand{\st}[1]{\sout{#1}}
\fi
\fi
\setlength{\emergencystretch}{3em} % prevent overfull lines
\providecommand{\tightlist}{%
\setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}
% \usepackage{xeCJK}
% \setCJKmainfont{AR PL UKai CN}
% \usepackage{unicode-math}
\setmathfont{Latin Modern Math}
\usepackage{bookmark}
\IfFileExists{xurl.sty}{\usepackage{xurl}}{} % add URL line breaks if available
\urlstyle{same}
\hypersetup{
pdftitle={「学习总结」数论},
pdfauthor={Jiayi Su (ShuYuMo)},
hidelinks,
pdfcreator={LaTeX via pandoc}}
\title{「学习总结」数论}
\author{Jiayi Su (ShuYuMo)}
\date{2020-12-20 23:41:12}
\begin{document}
\maketitle
\setstretch{1.3}
关于 \texttt{Lucas} 定理及其扩展,\texttt{BSGS}及其扩展,快速质因数分解
\texttt{Pollard-Rho},素性测试以及原根和阶的一些可爱知识。
\section{数论}\label{ux6570ux8bba}
\subsection{Lucas 定理}\label{lucas-ux5b9aux7406}
解决模数为小质数的组合数取模问题。
设模数为 \(P\),设
\[n = \sum_{i \ge 0}\limits{a_i P^i}, m = \sum_{i \ge 0}\limits{b_i P^i}\notag\]
\[
\dbinom{n}{m} \equiv \prod_{i \ge 0}\limits{\dbinom{a_i}{b_i}} ( mod\ P)
\] \{\#eq:Lucas\}
同时 (eq.\textasciitilde{}@eq:Lucas) 还有一种形式化的写法:
\[
\dbinom{n}{m} \equiv \dbinom{\lfloor n/P\rfloor}{\lfloor m/P\rfloor}\dbinom{n \mod P}{m \mod P}(mod\ P)
\] \{\#eq:Luacs\_ex\}
若小组合数可以 \(\mathcal{O}(1)\) 查询,\texttt{Luacs} 的时间复杂度为
\(\mathcal{O}(\log_p(n))\) 。
取模后不能保证 \(n \ge m\) ,查询组合数要写的细致一点。
特殊地,如果模数是 \(2\) ,根据 (eq.\textasciitilde{}@eq:Lucas)
考虑如果把数字进行二进制拆分求组合数,不难发现其在模 \(2\) 意义下的值为
\(\left[(n\ \And\ m) = m\right]\) 其中 \(n, m\) 的意义同
(eq.\textasciitilde{}@eq:Lucas)。
\begin{Shaded}
\begin{Highlighting}[]
\DataTypeTok{int}\NormalTok{ \_C\_}\OperatorTok{(}\DataTypeTok{int}\NormalTok{ m}\OperatorTok{,} \DataTypeTok{int}\NormalTok{ n}\OperatorTok{)\{}
\ControlFlowTok{if}\OperatorTok{(}\NormalTok{m }\OperatorTok{\textless{}}\NormalTok{ n}\OperatorTok{)} \ControlFlowTok{return} \DecValTok{0}\OperatorTok{;}
\ControlFlowTok{return}\NormalTok{ frac}\OperatorTok{[}\NormalTok{m}\OperatorTok{]} \OperatorTok{*}\DecValTok{1}\BuiltInTok{ll}\OperatorTok{*}\NormalTok{ ifrac}\OperatorTok{[}\NormalTok{n}\OperatorTok{]} \OperatorTok{\%}\NormalTok{ MOD }\OperatorTok{*}\DecValTok{1}\BuiltInTok{ll}\OperatorTok{*}\NormalTok{ ifrac}\OperatorTok{[}\NormalTok{m }\OperatorTok{{-}}\NormalTok{ n}\OperatorTok{]} \OperatorTok{\%}\NormalTok{ MOD}\OperatorTok{;}
\OperatorTok{\}}
\DataTypeTok{int}\NormalTok{ C}\OperatorTok{(}\DataTypeTok{int}\NormalTok{ m}\OperatorTok{,} \DataTypeTok{int}\NormalTok{ n}\OperatorTok{)\{}
\ControlFlowTok{if}\OperatorTok{(}\NormalTok{MOD }\OperatorTok{==} \DecValTok{2}\OperatorTok{)} \ControlFlowTok{return} \OperatorTok{((}\NormalTok{n }\OperatorTok{\&}\NormalTok{ m}\OperatorTok{)} \OperatorTok{==}\NormalTok{ n}\OperatorTok{);}
\ControlFlowTok{if}\OperatorTok{(}\NormalTok{m }\OperatorTok{\textless{}}\NormalTok{ n}\OperatorTok{)} \ControlFlowTok{return} \DecValTok{0}\OperatorTok{;} \ControlFlowTok{if}\OperatorTok{(}\NormalTok{n }\OperatorTok{==} \DecValTok{0}\OperatorTok{)} \ControlFlowTok{return} \DecValTok{1}\OperatorTok{;}
\ControlFlowTok{return}\NormalTok{ C}\OperatorTok{(}\NormalTok{m }\OperatorTok{/}\NormalTok{ MOD}\OperatorTok{,}\NormalTok{ n }\OperatorTok{/}\NormalTok{ MOD}\OperatorTok{)} \OperatorTok{*}\DecValTok{1}\BuiltInTok{ll}\OperatorTok{*}\NormalTok{ \_C\_}\OperatorTok{(}\NormalTok{m }\OperatorTok{\%}\NormalTok{ MOD}\OperatorTok{,}\NormalTok{ n }\OperatorTok{\%}\NormalTok{ MOD}\OperatorTok{)} \OperatorTok{\%}\NormalTok{ MOD}\OperatorTok{;}
\OperatorTok{\}}
\end{Highlighting}
\end{Shaded}
\subsection{扩展 Lucas 定理}\label{ux6269ux5c55-lucas-ux5b9aux7406}
解决模数为可能不为质数,且每个质数幂较小的组合数取模问题。
扩展卢卡斯定理和卢卡斯定理没有关系 并没有用到 \texttt{Luacs} 的基本思想
(eq.\textasciitilde{}@eq:Lucas) ,扩展 \texttt{Lucas} 定理
本质上是解决了阶乘逆元不存在的情况下,形如
\[
\frac{*}{n!}\mod P
\] \{\#eq:exl\}
式子的求值。 考虑组合数通项公式:
\[
\dbinom{n}{m} = \frac{n!}{m! (n-m)!}
\] \{\#eq:exC\}
限制我们不能直接求 (eq.\textasciitilde{}@eq:exC)
的原因——模数不为质数,无法保证逆元存在。考虑如果让分母和模数互质就可以求出这个式子的值了。
首先可以对模数 \(p\) 进行质因数分解,分解成形如
\(p = \prod_{i\ge 0}\limits{P_i^{e_i}}\) ,先分别求模数为 \(P_i^{e_i}\)
式子的值,然后 \texttt{CRT} 合并答案。
考虑如何求出模数为 \(P_i^{e_i}\) 的值。 设函数 \(g_p(n)\) 为 \(n!\)
中质因子 \(p\) 的幂次。 设 \(a = g_p(n), b=g_p(m), c=g_p(n-m)\) 易知:
\[
(eq.~@eq:exC) \equiv \frac{\frac{n!}{p^{a}}}{\frac{m!}{p^{b}} \frac{(n-m)!}{p^{c}}}\ p^{a-b-c}\ (mod \ P_i^{e_i})
\] \{\#eq:exLucas\}
设 \(k=P_i^{e_i}\) ,如果所求为 \(\frac{n!}{p_i^{a}}\) 对于 \(n!\)
的每一项
(\(1\times 2\times 3\times 4\times 5\times \dots \times n\)),其每一项在模
\(k\) 意义下的取值一定是每 \(k\) 个一次循环的。
同时每一个 \(P_i\) 的倍数都可以提出一个因子 \(P_i\)
也可以轻松知道,这个提出来的因子有 \(\lfloor\frac{n}{P_i}\rfloor\) 个,
即 \(P^{\lfloor\frac{n}{P_i}\rfloor}\) 这个式子不计入答案。
这些倍数提出一个 \(P_i\)
的因子之后一定会形成另一个阶乘的形式,这是一个子问题,可以递归计算。剩下的数字可以暴力计算一个周期,然后根据周期的出现次数,直接计算式子的值。
这样就可以算出 \(n!\) 去掉所有质因子 \(P_i\) 在模 \(P_i^{e_i}\)
意义下的值了。 保证了这个值和 \(P_i\) 互质,这样就可以求逆元了。 根据
(eq.\textasciitilde{}@eq:exLucas) 计算即可。
关于函数 \(g_p(n)\) 的计算:
\[
g_p(n) = \prod_{i\ge 1}\limits{\lfloor\frac{n}{p^i}\rfloor}= \frac{n - f_p(n)}{p - 1}
\]
其中 \(f_p(n)\) 为数字 \(n\) 在 \(p\) 进制下的数位和。
注意逆元不是质数,别\st{傻不拉几的}冲一个费马小定理。
复杂度是
\(\mathcal{O}\left(\sum_{i\ge 0}\limits{(\log P + P_i^{e_i})\log n}\right)\)
大概就是
\(\mathcal{O}\left(\sum_{i \ge 0}\limits{P_i^{e_i}\log n}\right)\) 吧。
\begin{Shaded}
\begin{Highlighting}[]
\DataTypeTok{int}\NormalTok{ f}\OperatorTok{(}\DataTypeTok{int}\NormalTok{ n}\OperatorTok{,} \DataTypeTok{int}\NormalTok{ p}\OperatorTok{)\{}
\DataTypeTok{int}\NormalTok{ ans }\OperatorTok{=} \DecValTok{0}\OperatorTok{;}
\ControlFlowTok{while}\OperatorTok{(}\NormalTok{n}\OperatorTok{)} \OperatorTok{\{}\NormalTok{ ans }\OperatorTok{+=}\NormalTok{ n }\OperatorTok{\%}\NormalTok{ p}\OperatorTok{;}\NormalTok{ n }\OperatorTok{/=}\NormalTok{ p}\OperatorTok{;} \OperatorTok{\}}
\ControlFlowTok{return}\NormalTok{ ans}\OperatorTok{;}
\OperatorTok{\}}
\DataTypeTok{int}\NormalTok{ g}\OperatorTok{(}\DataTypeTok{int}\NormalTok{ n}\OperatorTok{,} \DataTypeTok{int}\NormalTok{ P}\OperatorTok{)\{} \ControlFlowTok{return} \OperatorTok{(}\NormalTok{n }\OperatorTok{{-}}\NormalTok{ f}\OperatorTok{(}\NormalTok{n}\OperatorTok{,}\NormalTok{ P}\OperatorTok{))} \OperatorTok{/} \OperatorTok{(}\NormalTok{P }\OperatorTok{{-}} \DecValTok{1}\OperatorTok{);} \OperatorTok{\}}
\DataTypeTok{int}\NormalTok{ pow}\OperatorTok{(}\DataTypeTok{int}\NormalTok{ a}\OperatorTok{,} \DataTypeTok{int}\NormalTok{ b}\OperatorTok{,} \DataTypeTok{int}\NormalTok{ P}\OperatorTok{)\{} \DataTypeTok{int}\NormalTok{ ans }\OperatorTok{=} \DecValTok{1}\OperatorTok{;} \ControlFlowTok{while}\OperatorTok{(}\NormalTok{b}\OperatorTok{)} \OperatorTok{\{} \ControlFlowTok{if}\OperatorTok{(}\NormalTok{b }\OperatorTok{\&} \DecValTok{1}\OperatorTok{)}\NormalTok{ ans }\OperatorTok{=}\NormalTok{ ans }\OperatorTok{*}\DecValTok{1}\BuiltInTok{ll}\OperatorTok{*}\NormalTok{ a }\OperatorTok{\%}\NormalTok{ P}\OperatorTok{;}\NormalTok{ a }\OperatorTok{=}\NormalTok{ a }\OperatorTok{*}\DecValTok{1}\BuiltInTok{ll}\OperatorTok{*}\NormalTok{ a }\OperatorTok{\%}\NormalTok{ P}\OperatorTok{;}\NormalTok{ b }\OperatorTok{\textgreater{}\textgreater{}=} \DecValTok{1}\OperatorTok{;} \OperatorTok{\}} \ControlFlowTok{return}\NormalTok{ ans}\OperatorTok{;} \OperatorTok{\}}
\CommentTok{//int inv(int x, int MOD) \{ return pow(x, MOD {-} 2, MOD); \} // deleted: Fermat\textquotesingle{}s Little Theorem is NOT correct when MOD is not a prime number.}
\DataTypeTok{int}\NormalTok{ exgcd}\OperatorTok{(}\DataTypeTok{int}\NormalTok{ a}\OperatorTok{,} \DataTypeTok{int}\NormalTok{ b}\OperatorTok{,} \DataTypeTok{int} \OperatorTok{\&}\NormalTok{x}\OperatorTok{,} \DataTypeTok{int} \OperatorTok{\&}\NormalTok{y}\OperatorTok{)\{}
\ControlFlowTok{if}\OperatorTok{(!}\NormalTok{b}\OperatorTok{)} \OperatorTok{\{}\NormalTok{ x }\OperatorTok{=} \DecValTok{1}\OperatorTok{;}\NormalTok{ y }\OperatorTok{=} \DecValTok{0}\OperatorTok{;} \ControlFlowTok{return}\NormalTok{ a}\OperatorTok{;} \OperatorTok{\}}
\DataTypeTok{int}\NormalTok{ g }\OperatorTok{=}\NormalTok{ exgcd}\OperatorTok{(}\NormalTok{b}\OperatorTok{,}\NormalTok{ a }\OperatorTok{\%}\NormalTok{ b}\OperatorTok{,}\NormalTok{ y}\OperatorTok{,}\NormalTok{ x}\OperatorTok{);}
\NormalTok{ y }\OperatorTok{{-}=} \OperatorTok{(}\NormalTok{a }\OperatorTok{/}\NormalTok{ b}\OperatorTok{)} \OperatorTok{*}\NormalTok{ x}\OperatorTok{;}
\ControlFlowTok{return}\NormalTok{ g}\OperatorTok{;}
\OperatorTok{\}}
\DataTypeTok{int}\NormalTok{ inv}\OperatorTok{(}\DataTypeTok{int}\NormalTok{ a}\OperatorTok{,} \DataTypeTok{int}\NormalTok{ P}\OperatorTok{)\{} \CommentTok{// add: Exgcd}
\DataTypeTok{int}\NormalTok{ x}\OperatorTok{,}\NormalTok{ y}\OperatorTok{;}\NormalTok{ exgcd}\OperatorTok{(}\NormalTok{a}\OperatorTok{,}\NormalTok{ P}\OperatorTok{,}\NormalTok{ x}\OperatorTok{,}\NormalTok{ y}\OperatorTok{);}
\ControlFlowTok{return} \OperatorTok{(}\NormalTok{x }\OperatorTok{\%}\NormalTok{ P }\OperatorTok{+}\DecValTok{0}\BuiltInTok{ll}\OperatorTok{+}\NormalTok{ P }\OperatorTok{)} \OperatorTok{\%}\NormalTok{ P}\OperatorTok{;}
\OperatorTok{\}}
\DataTypeTok{int}\NormalTok{ calc}\OperatorTok{(}\DataTypeTok{int}\NormalTok{ n}\OperatorTok{,} \DataTypeTok{int}\NormalTok{ p}\OperatorTok{,} \DataTypeTok{int}\NormalTok{ k}\OperatorTok{)\{} \CommentTok{// calc n!(without p) mod p\^{}k}
\ControlFlowTok{if}\OperatorTok{(}\NormalTok{n }\OperatorTok{==} \DecValTok{0}\OperatorTok{)} \ControlFlowTok{return} \DecValTok{1}\OperatorTok{;}
\DataTypeTok{int}\NormalTok{ P }\OperatorTok{=} \DecValTok{1}\OperatorTok{;}\NormalTok{ rep}\OperatorTok{(}\NormalTok{i}\OperatorTok{,} \DecValTok{1}\OperatorTok{,}\NormalTok{ k}\OperatorTok{)}\NormalTok{ P }\OperatorTok{*=}\NormalTok{ p}\OperatorTok{;}
\DataTypeTok{int}\NormalTok{ res }\OperatorTok{=}\NormalTok{ n }\OperatorTok{\%}\NormalTok{ P}\OperatorTok{,}\NormalTok{ prd0 }\OperatorTok{=} \DecValTok{1}\OperatorTok{,}\NormalTok{ prd1 }\OperatorTok{=} \DecValTok{1}\OperatorTok{;}
\NormalTok{ rep}\OperatorTok{(}\NormalTok{i}\OperatorTok{,} \DecValTok{1}\OperatorTok{,}\NormalTok{ P}\OperatorTok{)\{}
\ControlFlowTok{if}\OperatorTok{(}\NormalTok{i }\OperatorTok{\%}\NormalTok{ p }\OperatorTok{==} \DecValTok{0}\OperatorTok{)} \ControlFlowTok{continue}\OperatorTok{;}
\ControlFlowTok{if}\OperatorTok{(}\NormalTok{i }\OperatorTok{\textless{}=}\NormalTok{ res}\OperatorTok{)}\NormalTok{ prd0 }\OperatorTok{=}\NormalTok{ prd0 }\OperatorTok{*}\DecValTok{1}\BuiltInTok{ll}\OperatorTok{*}\NormalTok{ i }\OperatorTok{\%}\NormalTok{ P}\OperatorTok{;}
\NormalTok{ prd1 }\OperatorTok{=}\NormalTok{ prd1 }\OperatorTok{*}\DecValTok{1}\BuiltInTok{ll}\OperatorTok{*}\NormalTok{ i }\OperatorTok{\%}\NormalTok{ P}\OperatorTok{;}
\OperatorTok{\}}
\DataTypeTok{int}\NormalTok{ ans }\OperatorTok{=}\NormalTok{ calc}\OperatorTok{(}\NormalTok{n }\OperatorTok{/}\NormalTok{ p}\OperatorTok{,}\NormalTok{ p}\OperatorTok{,}\NormalTok{ k}\OperatorTok{);}
\NormalTok{ ans }\OperatorTok{=}\NormalTok{ ans }\OperatorTok{*}\DecValTok{1}\BuiltInTok{ll}\OperatorTok{*}\NormalTok{ prd0 }\OperatorTok{\%}\NormalTok{ P}\OperatorTok{;}
\NormalTok{ ans }\OperatorTok{=}\NormalTok{ ans }\OperatorTok{*}\DecValTok{1}\BuiltInTok{ll}\OperatorTok{*}\NormalTok{ pow}\OperatorTok{(}\NormalTok{prd1}\OperatorTok{,}\NormalTok{ n }\OperatorTok{/}\NormalTok{ P}\OperatorTok{,}\NormalTok{ P}\OperatorTok{)} \OperatorTok{\%}\NormalTok{ P}\OperatorTok{;} \CommentTok{// changed \textasciigrave{}pow(prd1, n / p, P)\textasciigrave{} to \textasciigrave{}pow(prd1, n / P, P) \% P\textasciigrave{}.}
\ControlFlowTok{return}\NormalTok{ ans}\OperatorTok{;}
\OperatorTok{\}}
\NormalTok{LL n}\OperatorTok{,}\NormalTok{ m}\OperatorTok{;} \DataTypeTok{int}\NormalTok{ MOD}\OperatorTok{;}
\NormalTok{vector}\OperatorTok{\textless{}}\NormalTok{pair}\OperatorTok{\textless{}}\DataTypeTok{int}\OperatorTok{,} \DataTypeTok{int} \OperatorTok{\textgreater{}} \OperatorTok{\textgreater{}}\NormalTok{ d}\OperatorTok{;}
\KeywordTok{namespace}\NormalTok{ Divide}\OperatorTok{\{}
\AttributeTok{const} \DataTypeTok{int}\NormalTok{ \_ }\OperatorTok{=} \FloatTok{1e6} \OperatorTok{+} \DecValTok{100}\OperatorTok{;}
\DataTypeTok{int}\NormalTok{ np}\OperatorTok{[}\NormalTok{\_}\OperatorTok{],}\NormalTok{ prime}\OperatorTok{[}\NormalTok{\_}\OperatorTok{],}\NormalTok{ tot }\OperatorTok{=} \DecValTok{0}\OperatorTok{,}\NormalTok{ Mid}\OperatorTok{[}\NormalTok{\_}\OperatorTok{];}
\DataTypeTok{void}\NormalTok{ init}\OperatorTok{(}\DataTypeTok{int}\NormalTok{ n}\OperatorTok{)\{}
\NormalTok{ rep}\OperatorTok{(}\NormalTok{i}\OperatorTok{,} \DecValTok{2}\OperatorTok{,}\NormalTok{ n}\OperatorTok{)\{}
\ControlFlowTok{if}\OperatorTok{(!}\NormalTok{np}\OperatorTok{[}\NormalTok{i}\OperatorTok{])}\NormalTok{ prime}\OperatorTok{[++}\NormalTok{tot}\OperatorTok{]} \OperatorTok{=}\NormalTok{ i}\OperatorTok{,}\NormalTok{ Mid}\OperatorTok{[}\NormalTok{i}\OperatorTok{]} \OperatorTok{=}\NormalTok{ i}\OperatorTok{;}
\ControlFlowTok{for}\OperatorTok{(}\DataTypeTok{int}\NormalTok{ j }\OperatorTok{=} \DecValTok{1}\OperatorTok{;}\NormalTok{ j }\OperatorTok{\textless{}=}\NormalTok{ tot }\OperatorTok{\&\&}\NormalTok{ prime}\OperatorTok{[}\NormalTok{j}\OperatorTok{]} \OperatorTok{*}\NormalTok{ i }\OperatorTok{\textless{}=}\NormalTok{ n}\OperatorTok{;}\NormalTok{ j}\OperatorTok{++)\{}
\DataTypeTok{int}\NormalTok{ x }\OperatorTok{=}\NormalTok{ prime}\OperatorTok{[}\NormalTok{j}\OperatorTok{]} \OperatorTok{*}\NormalTok{ i}\OperatorTok{;}\NormalTok{ Mid}\OperatorTok{[}\NormalTok{x}\OperatorTok{]} \OperatorTok{=}\NormalTok{ prime}\OperatorTok{[}\NormalTok{j}\OperatorTok{];}
\NormalTok{ np}\OperatorTok{[}\NormalTok{x}\OperatorTok{]} \OperatorTok{=} \DecValTok{1}\OperatorTok{;}
\ControlFlowTok{if}\OperatorTok{(}\NormalTok{i }\OperatorTok{\%}\NormalTok{ prime}\OperatorTok{[}\NormalTok{j}\OperatorTok{]} \OperatorTok{==} \DecValTok{0}\OperatorTok{)} \ControlFlowTok{break}\OperatorTok{;}
\OperatorTok{\}}
\OperatorTok{\}}
\OperatorTok{\}}
\DataTypeTok{void}\NormalTok{ divide}\OperatorTok{(}\NormalTok{vector}\OperatorTok{\textless{}}\NormalTok{pair}\OperatorTok{\textless{}}\DataTypeTok{int}\OperatorTok{,} \DataTypeTok{int} \OperatorTok{\textgreater{}} \OperatorTok{\textgreater{}} \OperatorTok{\&}\NormalTok{ res}\OperatorTok{,} \DataTypeTok{int}\NormalTok{ MOD}\OperatorTok{)\{}
\ControlFlowTok{for}\OperatorTok{(}\DataTypeTok{int}\NormalTok{ i }\OperatorTok{=} \DecValTok{1}\OperatorTok{;}\NormalTok{ i }\OperatorTok{\textless{}=}\NormalTok{ tot}\OperatorTok{;}\NormalTok{ i}\OperatorTok{++)\{}
\ControlFlowTok{if}\OperatorTok{(}\NormalTok{MOD }\OperatorTok{\%}\NormalTok{ prime}\OperatorTok{[}\NormalTok{i}\OperatorTok{]} \OperatorTok{!=} \DecValTok{0}\OperatorTok{)} \ControlFlowTok{continue}\OperatorTok{;}
\NormalTok{ pair}\OperatorTok{\textless{}}\DataTypeTok{int}\OperatorTok{,} \DataTypeTok{int} \OperatorTok{\textgreater{}}\NormalTok{ ans}\OperatorTok{;}\NormalTok{ ans}\OperatorTok{.}\NormalTok{fi }\OperatorTok{=}\NormalTok{ prime}\OperatorTok{[}\NormalTok{i}\OperatorTok{];}
\NormalTok{ ans}\OperatorTok{.}\NormalTok{se }\OperatorTok{=} \DecValTok{0}\OperatorTok{;}
\ControlFlowTok{while}\OperatorTok{(}\NormalTok{MOD }\OperatorTok{\%}\NormalTok{ prime}\OperatorTok{[}\NormalTok{i}\OperatorTok{]} \OperatorTok{==} \DecValTok{0}\OperatorTok{)}\NormalTok{ MOD }\OperatorTok{/=}\NormalTok{ prime}\OperatorTok{[}\NormalTok{i}\OperatorTok{],}\NormalTok{ ans}\OperatorTok{.}\NormalTok{se}\OperatorTok{++;}
\NormalTok{ res}\OperatorTok{.}\NormalTok{push\_back}\OperatorTok{(}\NormalTok{ans}\OperatorTok{);}
\OperatorTok{\}}
\OperatorTok{\}}\CommentTok{// 这里可以直接 $\textbackslash{}mathcal\{O\}(\textbackslash{}sqrt\{n\})$ 试除即可,没必要分解质因数。写的时候太年轻了。}
\OperatorTok{\}}
\DataTypeTok{int}\NormalTok{ CRT}\OperatorTok{(}\NormalTok{vector}\OperatorTok{\textless{}}\NormalTok{pair}\OperatorTok{\textless{}}\DataTypeTok{int}\OperatorTok{,} \DataTypeTok{int} \OperatorTok{\textgreater{}} \OperatorTok{\textgreater{}} \OperatorTok{\&}\NormalTok{ A}\OperatorTok{)\{}
\DataTypeTok{int}\NormalTok{ W }\OperatorTok{=}\NormalTok{ MOD}\OperatorTok{;}
\DataTypeTok{int}\NormalTok{ ans }\OperatorTok{=} \DecValTok{0}\OperatorTok{;}
\NormalTok{ rep}\OperatorTok{(}\NormalTok{i}\OperatorTok{,} \DecValTok{0}\OperatorTok{,}\NormalTok{ A}\OperatorTok{.}\NormalTok{size}\OperatorTok{()} \OperatorTok{{-}} \DecValTok{1}\OperatorTok{)\{}
\NormalTok{ pair}\OperatorTok{\textless{}}\DataTypeTok{int}\OperatorTok{,} \DataTypeTok{int} \OperatorTok{\textgreater{}}\NormalTok{ now }\OperatorTok{=}\NormalTok{ A}\OperatorTok{[}\NormalTok{i}\OperatorTok{];}
\NormalTok{ ans }\OperatorTok{=} \OperatorTok{(}\NormalTok{ans }\OperatorTok{+}\DecValTok{0}\BuiltInTok{ll}\OperatorTok{+} \OperatorTok{(}\NormalTok{ now}\OperatorTok{.}\NormalTok{fi }\OperatorTok{*}\DecValTok{1}\BuiltInTok{ll}\OperatorTok{*} \OperatorTok{(}\NormalTok{W }\OperatorTok{/}\NormalTok{ now}\OperatorTok{.}\NormalTok{se}\OperatorTok{)} \OperatorTok{\%}\NormalTok{ MOD }\OperatorTok{*}\DecValTok{1}\BuiltInTok{ll}\OperatorTok{*}\NormalTok{ inv}\OperatorTok{(}\NormalTok{W }\OperatorTok{/}\NormalTok{ now}\OperatorTok{.}\NormalTok{se}\OperatorTok{,}\NormalTok{ now}\OperatorTok{.}\NormalTok{se}\OperatorTok{)} \OperatorTok{)} \OperatorTok{\%}\NormalTok{ MOD}\OperatorTok{)} \OperatorTok{\%}\NormalTok{ MOD}\OperatorTok{;}
\OperatorTok{\}}
\ControlFlowTok{return}\NormalTok{ ans}\OperatorTok{;}
\OperatorTok{\}}
\NormalTok{vector}\OperatorTok{\textless{}}\NormalTok{pair}\OperatorTok{\textless{}}\DataTypeTok{int}\OperatorTok{,} \DataTypeTok{int} \OperatorTok{\textgreater{}} \OperatorTok{\textgreater{}}\NormalTok{ Ans}\OperatorTok{;}
\DataTypeTok{signed}\NormalTok{ main}\OperatorTok{()\{} \CommentTok{//freopen(".in", "r", stdin);}
\NormalTok{ Read}\OperatorTok{(}\NormalTok{n}\OperatorTok{)(}\NormalTok{m}\OperatorTok{)(}\NormalTok{MOD}\OperatorTok{);}\NormalTok{ Divide}\OperatorTok{::}\NormalTok{init}\OperatorTok{(}\NormalTok{MOD }\OperatorTok{+} \DecValTok{2}\OperatorTok{);}
\NormalTok{ Divide}\OperatorTok{::}\NormalTok{divide}\OperatorTok{(}\NormalTok{d}\OperatorTok{,}\NormalTok{ MOD}\OperatorTok{);}
\NormalTok{ rep}\OperatorTok{(}\NormalTok{i}\OperatorTok{,} \DecValTok{0}\OperatorTok{,}\NormalTok{ d}\OperatorTok{.}\NormalTok{size}\OperatorTok{()} \OperatorTok{{-}} \DecValTok{1}\OperatorTok{)\{}
\NormalTok{ pair}\OperatorTok{\textless{}}\DataTypeTok{int}\OperatorTok{,} \DataTypeTok{int}\OperatorTok{\textgreater{}}\NormalTok{ NowMod }\OperatorTok{=}\NormalTok{ d}\OperatorTok{[}\NormalTok{i}\OperatorTok{];} \DataTypeTok{int}\NormalTok{ P }\OperatorTok{=}\NormalTok{ pow}\OperatorTok{(}\NormalTok{NowMod}\OperatorTok{.}\NormalTok{fi}\OperatorTok{,}\NormalTok{ NowMod}\OperatorTok{.}\NormalTok{se}\OperatorTok{);}
\DataTypeTok{int}\NormalTok{ rA }\OperatorTok{=}\NormalTok{ calc}\OperatorTok{(}\NormalTok{n}\OperatorTok{,}\NormalTok{ NowMod}\OperatorTok{.}\NormalTok{fi}\OperatorTok{,}\NormalTok{ NowMod}\OperatorTok{.}\NormalTok{se}\OperatorTok{);}
\DataTypeTok{int}\NormalTok{ rB }\OperatorTok{=}\NormalTok{ calc}\OperatorTok{(}\NormalTok{n }\OperatorTok{{-}}\NormalTok{ m}\OperatorTok{,}\NormalTok{ NowMod}\OperatorTok{.}\NormalTok{fi}\OperatorTok{,}\NormalTok{ NowMod}\OperatorTok{.}\NormalTok{se}\OperatorTok{);}
\DataTypeTok{int}\NormalTok{ rC }\OperatorTok{=}\NormalTok{ calc}\OperatorTok{(}\NormalTok{m}\OperatorTok{,}\NormalTok{ NowMod}\OperatorTok{.}\NormalTok{fi}\OperatorTok{,}\NormalTok{ NowMod}\OperatorTok{.}\NormalTok{se}\OperatorTok{);}
\DataTypeTok{int}\NormalTok{ ans }\OperatorTok{=} \DecValTok{0}\OperatorTok{;}
\NormalTok{ ans }\OperatorTok{=}\NormalTok{ rA }\OperatorTok{*}\DecValTok{1}\BuiltInTok{ll}\OperatorTok{*}\NormalTok{ inv}\OperatorTok{(}\NormalTok{rB}\OperatorTok{,}\NormalTok{ P}\OperatorTok{)} \OperatorTok{\%}\NormalTok{ P }\OperatorTok{*}\DecValTok{1}\BuiltInTok{ll}\OperatorTok{*}\NormalTok{ inv}\OperatorTok{(}\NormalTok{rC}\OperatorTok{,}\NormalTok{ P}\OperatorTok{)} \OperatorTok{\%}\NormalTok{ P}\OperatorTok{;}
\NormalTok{ ans }\OperatorTok{=}\NormalTok{ ans }\OperatorTok{*}\DecValTok{1}\BuiltInTok{ll}\OperatorTok{*}\NormalTok{ pow}\OperatorTok{(}\NormalTok{NowMod}\OperatorTok{.}\NormalTok{fi}\OperatorTok{,}\NormalTok{ g}\OperatorTok{(}\NormalTok{n}\OperatorTok{,}\NormalTok{ NowMod}\OperatorTok{.}\NormalTok{fi}\OperatorTok{)} \OperatorTok{{-}}\NormalTok{ g}\OperatorTok{(}\NormalTok{n }\OperatorTok{{-}}\NormalTok{ m}\OperatorTok{,}\NormalTok{ NowMod}\OperatorTok{.}\NormalTok{fi}\OperatorTok{)} \OperatorTok{{-}}\NormalTok{ g}\OperatorTok{(}\NormalTok{m}\OperatorTok{,}\NormalTok{ NowMod}\OperatorTok{.}\NormalTok{fi}\OperatorTok{),}\NormalTok{ P}\OperatorTok{)} \OperatorTok{\%}\NormalTok{ P}\OperatorTok{;}
\NormalTok{ Ans}\OperatorTok{.}\NormalTok{push\_back}\OperatorTok{(}\NormalTok{mp}\OperatorTok{(}\NormalTok{ans}\OperatorTok{,}\NormalTok{ P}\OperatorTok{));}
\OperatorTok{\}}
\NormalTok{ printf}\OperatorTok{(}\StringTok{"}\SpecialCharTok{\%lld\textbackslash{}n}\StringTok{"}\OperatorTok{,}\NormalTok{ CRT}\OperatorTok{(}\NormalTok{Ans}\OperatorTok{));}
\ControlFlowTok{return} \DecValTok{0}\OperatorTok{;}
\OperatorTok{\}}
\end{Highlighting}
\end{Shaded}
\subsection{BSGS}\label{bsgs}
求解模数为质数的指数方程。\(a^x \equiv b\ (mod\ p) \ p \in P\)。
由欧拉定理可知,本质不同的解 范围为 \([0, p-1]\) 。 朴素做法就是枚举一下
\([0, p-1]\) 判断是不是解,但是这样显然太慢了。 考虑选一个参数 \(T\),则
\(\forall x \in [0, p-1]\) 都可以表示成 \(r \times T + c\ \ (b < T)\)
的形式,易知,\(r = \lfloor\frac{x}{T} \rfloor, c=x\ mod \ T\)
带入原式:
\[
a^{rT + c} \equiv b\ (mod\ p)
\]
\[
a^{c} \equiv b\times a^{-rT}\ (mod\ p)
\] \{\#eq:BSGS0\}
可以考虑枚举 \(c\) 算出每一个 \(a^{c}\) 压入 \texttt{set} 或
\texttt{hash} 表,然后枚举每一个 \(r\) 算出每一个 \(a^{-rT} \times b\)
在之前算出的答案中查找,如果能找到相同的取值(即:满足(eq.\textasciitilde{}@eq:BSGS0)),就说明当前这个
\(r\) 就是一个答案。
\begin{Shaded}
\begin{Highlighting}[]
\DataTypeTok{int}\NormalTok{ inv}\OperatorTok{(}\DataTypeTok{int}\NormalTok{ a}\OperatorTok{,} \DataTypeTok{int}\NormalTok{ P}\OperatorTok{)\{}
\DataTypeTok{int}\NormalTok{ x}\OperatorTok{,}\NormalTok{ y}\OperatorTok{;}\NormalTok{ exgcd}\OperatorTok{(}\NormalTok{a}\OperatorTok{,}\NormalTok{ P}\OperatorTok{,}\NormalTok{ x}\OperatorTok{,}\NormalTok{ y}\OperatorTok{);}
\ControlFlowTok{return} \OperatorTok{(}\NormalTok{x }\OperatorTok{\%}\NormalTok{ P }\OperatorTok{+}\DecValTok{0}\BuiltInTok{ll}\OperatorTok{+}\NormalTok{ P}\OperatorTok{)} \OperatorTok{\%}\NormalTok{ P}\OperatorTok{;}
\OperatorTok{\}}
\DataTypeTok{int}\NormalTok{ BSGS}\OperatorTok{(}\DataTypeTok{int}\NormalTok{ a}\OperatorTok{,} \DataTypeTok{int}\NormalTok{ b}\OperatorTok{,} \DataTypeTok{int}\NormalTok{ P}\OperatorTok{)} \OperatorTok{\{} \CommentTok{// a\^{}x = b (mod P) (a, p) = 1;}
\AttributeTok{static}\NormalTok{ set}\OperatorTok{\textless{}}\DataTypeTok{int}\OperatorTok{\textgreater{}}\NormalTok{ S}\OperatorTok{;}\NormalTok{ S}\OperatorTok{.}\NormalTok{clear}\OperatorTok{();}
\DataTypeTok{int}\NormalTok{ T }\OperatorTok{=}\NormalTok{ sqrt}\OperatorTok{(}\NormalTok{P}\OperatorTok{);}
\ControlFlowTok{for}\OperatorTok{(}\DataTypeTok{int}\NormalTok{ i }\OperatorTok{=} \DecValTok{0}\OperatorTok{,}\NormalTok{ x }\OperatorTok{=} \DecValTok{1}\OperatorTok{;}\NormalTok{ i }\OperatorTok{\textless{}}\NormalTok{ T}\OperatorTok{;}\NormalTok{ i}\OperatorTok{++,}\NormalTok{ x }\OperatorTok{=}\NormalTok{ x }\OperatorTok{*}\DecValTok{1}\BuiltInTok{ll}\OperatorTok{*}\NormalTok{ a }\OperatorTok{\%}\NormalTok{ P}\OperatorTok{)} \ControlFlowTok{if}\OperatorTok{(}\NormalTok{x }\OperatorTok{!=}\NormalTok{ b}\OperatorTok{)}\NormalTok{ S}\OperatorTok{.}\NormalTok{insert}\OperatorTok{(}\NormalTok{x}\OperatorTok{);} \ControlFlowTok{else} \ControlFlowTok{return}\NormalTok{ i}\OperatorTok{;}
\ControlFlowTok{for}\OperatorTok{(}\DataTypeTok{int}\NormalTok{ i }\OperatorTok{=} \DecValTok{1}\OperatorTok{;}\NormalTok{ i }\OperatorTok{\textless{}=}\NormalTok{ T}\OperatorTok{;}\NormalTok{ i}\OperatorTok{++)\{}
\DataTypeTok{int}\NormalTok{ now }\OperatorTok{=}\NormalTok{ b }\OperatorTok{*}\DecValTok{1}\BuiltInTok{ll}\OperatorTok{*}\NormalTok{ inv}\OperatorTok{(}\NormalTok{pow}\OperatorTok{(}\NormalTok{a}\OperatorTok{,}\NormalTok{ T }\OperatorTok{*}\NormalTok{ i}\OperatorTok{,}\NormalTok{ P}\OperatorTok{),}\NormalTok{ P}\OperatorTok{)} \OperatorTok{\%}\NormalTok{ P}\OperatorTok{;}
\ControlFlowTok{if}\OperatorTok{(!}\NormalTok{S}\OperatorTok{.}\NormalTok{count}\OperatorTok{(}\NormalTok{now}\OperatorTok{))} \ControlFlowTok{continue}\OperatorTok{;}
\ControlFlowTok{for}\OperatorTok{(}\DataTypeTok{int}\NormalTok{ j }\OperatorTok{=}\NormalTok{ i }\OperatorTok{*}\NormalTok{ T}\OperatorTok{,}\NormalTok{ V }\OperatorTok{=}\NormalTok{ pow}\OperatorTok{(}\NormalTok{a}\OperatorTok{,}\NormalTok{ j}\OperatorTok{,}\NormalTok{ P}\OperatorTok{);} \OperatorTok{;}\NormalTok{ j}\OperatorTok{++,}\NormalTok{ V }\OperatorTok{=}\NormalTok{ V }\OperatorTok{*}\DecValTok{1}\BuiltInTok{ll}\OperatorTok{*}\NormalTok{ a }\OperatorTok{\%}\NormalTok{ P}\OperatorTok{)} \ControlFlowTok{if}\OperatorTok{(}\NormalTok{V }\OperatorTok{==}\NormalTok{ b}\OperatorTok{)} \ControlFlowTok{return}\NormalTok{ j}\OperatorTok{;}
\OperatorTok{\}}
\ControlFlowTok{return} \OperatorTok{{-}}\DecValTok{1}\OperatorTok{;}
\OperatorTok{\}}
\end{Highlighting}
\end{Shaded}
\subsection{扩展 BSGS}\label{ux6269ux5c55-bsgs}
用于求解模数不一定为质数的指数方程。\(a^x \equiv b\ (mod\ p) \ p \in P\)。
限制不能直接 \(BSGS\) 的因素就是不能保证逆元存在,那就考虑如何让
\((a, p) = 1\) 。
具体地,设 \(d_1=\gcd(a,p)\) 。如果 \(d_1\nmid b\)
,则原方程无解。否则我们把方程同时除以 \(d_1\) ,得到\\
\[
\frac{a}{d_1}\cdot a^{x-1}\equiv \frac{b}{d_1}\pmod{\frac{p}{d_1}}
\]
如果 \(a\) 和 \(\frac{p}{d_1}\) 仍不互质就再除,设
\(d_2=\gcd\left(a,\frac{p}{d_1}\right)\) 。如果
\(d_2\nmid \frac{b}{d_1}\) ,则方程无解;否则同时除以 \(d_2\) 得到. \[
\frac{a^2}{d_1d_2}\cdot a^{x-2}≡\frac{b}{d_1d_2} \pmod{\frac{p}{d_1d_2}}
\] 直到模数和 \(a\) 互质。 记 \(D = \prod\limits{d_i}\),则
\[
\frac{a^k}{D}\cdot a^{x-k}\equiv\frac{b}{D} \pmod{\frac{p}{D}}
\] 就可以使用 \texttt{BSGS} 求解了。
需要注意的是:
\begin{itemize}
\tightlist
\item
解可能小于 \(k\) ,可以暴力枚举一下小于 \(k\) 的幂次,判断一下
\item
如果有 \(d_i\) 不是 \(b\) 的约数,直接判断无解即可。
\end{itemize}
\begin{verbatim}
int exBSGS(int a, int b, int P){
int D = 1, k = 0, tp = P;
while(true) { int g = gcd(a, tp); if(g == 1) break; tp /= g; D *= g; k++; }
for(int i = 0, V = 1; i <= k; i++, V = V *1ll* a % P) if(V == b) return i; // changed: It must be executed before `if(b % D != 0) return -1`).
int S = pow(a, k, tp);
if(b % D != 0) return -1; // added: It is necessary!
int B = b *1ll* inv(S, tp) % tp;
int r = BSGS(a, B, tp);
return r == -1 ? -1 : r + k;
}
\end{verbatim}
已通过 luogu 和 SPOJ 的测试数据,但是 zhx (Orz) 的数据只有 \(50pts\)
待填。
\subsection{Miller\_Rabin}\label{miller_rabin}
快速判断大数素性。
一个结论: 如果 \(n\) 为质数 \(\forall a < n\) 设
\(n - 1 = d \times 2^r, (d, 2) = 1\),则 \[
a^d \equiv 1 \pmod n
\] \{\#eq:MR0\}
\[
\exists\ 0 \le i < r, s.t. a^{d\times i^i} \equiv -1 \pmod n
\] \{\#eq:MR1\}
(eq.\textasciitilde{}@eq:MR0) 和 (eq.\textasciitilde{}@eq:MR1)
至少一个满足。 是否为充要条件带查。
我们需要判断一个数是否为质数,只需要判断是否符合上面的定理即可,但是
\(\forall a < n\) 对复杂度不友好。
通常的做法是选择若干质数当作底数 \(a\),进行判断。
关于底数的选择: \textgreater 如果选用2, 3, 7,
61和24251作为底数,那么10\^{}16内唯一的强伪素数为46 856 248 255 981
——matrix67博客 (Orz)
选择
\texttt{int\ PrimeList{[}10{]}\ =\ \{2,\ 3,\ 7,\ 61,\ 24251,\ 11,\ 29\};}
即可,或者再随意加几个小质数。
需要注意快速乘法的实现。 复杂度为 \(\mathcal{O}(k \log x)\)
\begin{Shaded}
\begin{Highlighting}[]
\PreprocessorTok{\#define LL }\DataTypeTok{long}\PreprocessorTok{ }\DataTypeTok{long}\PreprocessorTok{ }
\PreprocessorTok{\#define ULL }\DataTypeTok{unsigned}\PreprocessorTok{ }\DataTypeTok{long}\PreprocessorTok{ }\DataTypeTok{long}\PreprocessorTok{ }
\KeywordTok{inline}\NormalTok{ ULL mul}\OperatorTok{(}\NormalTok{ULL a}\OperatorTok{,}\NormalTok{ ULL b}\OperatorTok{,}\NormalTok{ ULL MOD}\OperatorTok{)\{} \CommentTok{// unsigned long long }
\NormalTok{ LL R }\OperatorTok{=} \OperatorTok{(}\NormalTok{LL}\OperatorTok{)}\NormalTok{a}\OperatorTok{*}\NormalTok{b }\OperatorTok{{-}} \OperatorTok{(}\NormalTok{LL}\OperatorTok{)((}\NormalTok{ULL}\OperatorTok{)((}\DataTypeTok{long} \DataTypeTok{double}\OperatorTok{)}\NormalTok{a }\OperatorTok{*}\NormalTok{ b }\OperatorTok{/}\NormalTok{ MOD}\OperatorTok{)} \OperatorTok{*}\NormalTok{ MOD}\OperatorTok{);}
\CommentTok{// 只关心两个答案的差值,这个差值一定小于 unsigned long long 的最大值,}
\CommentTok{// 所以在哪个剩余系下都不重要,不管差值是什么都能还原出原始数值。 }
\ControlFlowTok{if}\OperatorTok{(}\NormalTok{R }\OperatorTok{\textless{}} \DecValTok{0}\OperatorTok{)}\NormalTok{ R }\OperatorTok{+=}\NormalTok{ MOD}\OperatorTok{;}
\ControlFlowTok{if}\OperatorTok{(}\NormalTok{R }\OperatorTok{\textgreater{}}\NormalTok{ MOD}\OperatorTok{)}\NormalTok{ R }\OperatorTok{{-}=}\NormalTok{ MOD}\OperatorTok{;}
\ControlFlowTok{return}\NormalTok{ R}\OperatorTok{;}
\OperatorTok{\}}
\KeywordTok{inline}\NormalTok{ LL pow}\OperatorTok{(}\NormalTok{LL a}\OperatorTok{,}\NormalTok{ LL b}\OperatorTok{,}\NormalTok{ LL P}\OperatorTok{)} \OperatorTok{\{}\NormalTok{ LL ans }\OperatorTok{=} \DecValTok{1}\OperatorTok{;} \ControlFlowTok{while}\OperatorTok{(}\NormalTok{b}\OperatorTok{)\{} \ControlFlowTok{if}\OperatorTok{(}\NormalTok{b }\OperatorTok{\&} \DecValTok{1}\OperatorTok{)}\NormalTok{ ans }\OperatorTok{=}\NormalTok{ mul}\OperatorTok{(}\NormalTok{ans}\OperatorTok{,}\NormalTok{ a}\OperatorTok{,}\NormalTok{ P}\OperatorTok{);}\NormalTok{ a }\OperatorTok{=}\NormalTok{ mul}\OperatorTok{(}\NormalTok{a}\OperatorTok{,}\NormalTok{ a}\OperatorTok{,}\NormalTok{ P}\OperatorTok{);}\NormalTok{ b }\OperatorTok{\textgreater{}\textgreater{}=} \DecValTok{1}\OperatorTok{;} \OperatorTok{\}} \ControlFlowTok{return}\NormalTok{ ans}\OperatorTok{;} \OperatorTok{\}}
\DataTypeTok{int}\NormalTok{ PrimeList}\OperatorTok{[}\DecValTok{10}\OperatorTok{]} \OperatorTok{=} \OperatorTok{\{}\DecValTok{2}\OperatorTok{,} \DecValTok{3}\OperatorTok{,} \DecValTok{7}\OperatorTok{,} \DecValTok{61}\OperatorTok{,} \DecValTok{24251}\OperatorTok{,} \DecValTok{11}\OperatorTok{,} \DecValTok{17}\OperatorTok{,} \DecValTok{19}\OperatorTok{,} \DecValTok{29}\OperatorTok{,} \DecValTok{27}\OperatorTok{\};}
\DataTypeTok{bool}\NormalTok{ Miller\_Rabin}\OperatorTok{(}\DataTypeTok{int}\NormalTok{ a}\OperatorTok{,}\NormalTok{ LL n}\OperatorTok{)\{}
\NormalTok{ LL d }\OperatorTok{=}\NormalTok{ n }\OperatorTok{{-}} \DecValTok{1}\OperatorTok{,}\NormalTok{ r }\OperatorTok{=} \DecValTok{0}\OperatorTok{,}\NormalTok{ x}\OperatorTok{;}
\ControlFlowTok{while}\OperatorTok{(!(}\NormalTok{d }\OperatorTok{\&} \DecValTok{1}\OperatorTok{))}\NormalTok{ d }\OperatorTok{\textgreater{}\textgreater{}=} \DecValTok{1}\OperatorTok{,}\NormalTok{ r}\OperatorTok{++;}
\ControlFlowTok{if}\OperatorTok{((}\NormalTok{x }\OperatorTok{=}\NormalTok{ pow}\OperatorTok{(}\NormalTok{a}\OperatorTok{,}\NormalTok{ d}\OperatorTok{,}\NormalTok{ n}\OperatorTok{))} \OperatorTok{==} \DecValTok{1}\OperatorTok{)} \ControlFlowTok{return} \KeywordTok{true}\OperatorTok{;}
\ControlFlowTok{for}\OperatorTok{(}\DataTypeTok{int}\NormalTok{ t }\OperatorTok{=} \DecValTok{1}\OperatorTok{;}\NormalTok{ t }\OperatorTok{\textless{}=}\NormalTok{ r}\OperatorTok{;}\NormalTok{ t}\OperatorTok{++,}\NormalTok{ x }\OperatorTok{=}\NormalTok{ mul}\OperatorTok{(}\NormalTok{x}\OperatorTok{,}\NormalTok{ x}\OperatorTok{,}\NormalTok{ n}\OperatorTok{))} \ControlFlowTok{if}\OperatorTok{(}\NormalTok{x }\OperatorTok{==}\NormalTok{ n }\OperatorTok{{-}} \DecValTok{1}\OperatorTok{)} \ControlFlowTok{return} \KeywordTok{true}\OperatorTok{;}
\ControlFlowTok{return} \KeywordTok{false}\OperatorTok{;}
\OperatorTok{\}}
\DataTypeTok{bool}\NormalTok{ Prime}\OperatorTok{(}\NormalTok{LL x}\OperatorTok{)\{}
\ControlFlowTok{if}\OperatorTok{(}\NormalTok{x }\OperatorTok{\textless{}=} \DecValTok{2}\OperatorTok{)} \ControlFlowTok{return}\NormalTok{ x }\OperatorTok{==} \DecValTok{2}\OperatorTok{;}
\ControlFlowTok{for}\OperatorTok{(}\DataTypeTok{int}\NormalTok{ i }\OperatorTok{=} \DecValTok{0}\OperatorTok{;}\NormalTok{ i }\OperatorTok{\textless{}} \DecValTok{7}\OperatorTok{;}\NormalTok{ i}\OperatorTok{++)\{}
\ControlFlowTok{if}\OperatorTok{(}\NormalTok{x }\OperatorTok{==}\NormalTok{ PrimeList}\OperatorTok{[}\NormalTok{i}\OperatorTok{])} \ControlFlowTok{return} \KeywordTok{true}\OperatorTok{;}
\ControlFlowTok{if}\OperatorTok{(!}\NormalTok{Miller\_Rabin}\OperatorTok{(}\NormalTok{PrimeList}\OperatorTok{[}\NormalTok{i}\OperatorTok{],}\NormalTok{ x}\OperatorTok{))} \ControlFlowTok{return} \KeywordTok{false}\OperatorTok{;}
\OperatorTok{\}}
\ControlFlowTok{return} \KeywordTok{true}\OperatorTok{;}
\OperatorTok{\}}
\end{Highlighting}
\end{Shaded}
\subsection{Pollard-Rho}\label{pollard-rho}
快速分解质因数的解决方案。 复杂度 \(\mathcal{O}(Accepted)\)。大约为
\(\mathcal{O}(n^{1/4})\) ,算法导论给出的是
\(\mathcal{O}(\sqrt{n})\),但是随机数据下,\(1s\) 分解 \(10000\) 个
\(10^{18}\) 量级的数字问题不大……
\begin{quote}
任何一个伟大的算法都来自于一个微不足道的猜想。
\end{quote}
考虑给出一个数字 \(n\) ,如果可以快速找其中一个因子 \(g\)
,然后继续递归分解 \(n / g\) 和 \(g\) 即可完成质因数分解。
考虑如何快速寻找一个数 \(n\) 的因子 \(g\);
(以下复杂度分析都是基于最坏情况,即素数平方的形式)
\begin{enumerate}
\def\labelenumi{\arabic{enumi}.}
\tightlist
\item
枚举 \(n\) 的因子 \(g\) 试除,复杂度为 \(\mathcal{O}(\sqrt{n})\)
\item
考虑随机枚举 \(n\) 的因子,试除,期望复杂度为 \(\mathcal{O}(n)\)
\item
没必要随机枚举 \(n\) 的因子,只需要随机一个数字 \(a\) 那么
\(g = gcd(a, n)\) 就是 \(n\)
的一个因子,需要找到一个不平凡的因子才可以,这样 \(a\)
的合法取值变成了 \(n\) 的因子或 \(n\) 因子的倍数,合法的 \(a\)
的取值有 \(\mathcal{O}(\sqrt{n})\) 个。
\end{enumerate}
着重考虑下一种优化:
考虑一次性随机出 \(k\) 个数字 分别记作
\(a_{1...k}\),考虑其中两两的差值,应该有 \(k^2\) 种差值。试图让其差值与
\(n\) 求最大公约数 \(c\),若 \(c\) 不为 \(1\) 则 \(c\)
就是一个不平凡因子。 差值与 \(n\) 的最大公约数为 \(c\),等价于 差值在
\(c\) 的剩余系下同余 \(0\)。考虑这 \(k^2\) 个差值都不等于 \(0\)
的概率为多少。因为这 \(k\) 个数字是随机的,那么他们差值的取值在 \(c\)
的剩余系下也是在 \([0, c]\) 等概率分布的。那么都不等于 \(0\) 的概率为
\(\left(\frac{c - 1}{c}\right)^{k^2}\)。最坏情况下,合法的 \(c\)
只有一个,且等于 \(\sqrt{n}\) 根据
\[\left(\frac{n-1}{n}\right)^{n} \approx \frac{1}{e}\]
当 \(k = n^{1/4}\) 时,取不到值得概率为 \(\frac{1}{e}\)。稍微增大几倍的
\(k\) 可以降低找不到概率。 注意这样的做法需要枚举所有 \(k^2\)
对差值,再加上判断和取不到约数情况的出现,复杂度要劣于
\(\mathcal{O}(\sqrt{n})\)。 枚举 \(k^2\)
对差值是上面算法的瓶颈。引入伪随机函数 \(f(x) = f^2(x-1) + c \mod n\)。
这个函数有几点比较优良性质:
\begin{enumerate}
\def\labelenumi{\arabic{enumi}.}
\tightlist
\item
仍然可以看成是一种随机函数,保证了上面推导中对随机的依赖性。(个人感觉不能保证完全随机……)
\item
经过一段次数的递推之后会进入一个循环,因为之和前一项的值有关,而且取值只能在
\([0, n-1]\) 内,所以至多 \(\mathcal{O}(n)\)
次递推,一定会进入一个循环。
\item
\[g \mid f(i) - f(j)\ \ \ \Rightarrow\ \ \ g \mid f(i + 1) - f(j - 1)\]
证明:
\[f(i + 1) - f(j + 1) = (f^2(i) + c) - (f^2(j) + c) = (f(i) - f(j))(f(i) + f(j))\]
也就是 \(g\) 是否为某两个差值的约数之和两个差值下标的差有关。
\end{enumerate}
考虑两个指针,分别为 \(L, R\) ,每次 \(L\) 递推一个, \(R\)
递推两个,因为函数存在环,这两个差值总会在环上相遇,从开始到相遇之间的时间,每次执行后,判断其差值与
\(n\) 的 \texttt{gcd} 是否不等于 \(1\)
即可。相遇所需时间,即环的大小约为 \(\sqrt{n}\)
。环的大小决定运行最坏时间。
这样还是不够快,其实不需要每次执行都算一下 \(gcd\)
可以执行一些步数之后,将其差值乘到一起,一起求 \(gcd\)
即可,可以证明,这样的变化保证答案不会变劣。
考虑让两个指针倍增的往前跳,倍增若干次后,每次计算两个指针距离之间的所有差值取值,乘在一起计算
\(gcd\)。
标准代码中:差值相同的解会计算多次。但是确实大大提升运行效率。这里就不会分析了,但是直观感觉就是,其实那个函数往后递推的次数越多会包含更多的因子。也就是说,如果两对函数值差值下标距离一样,其实下标大的那一对会更优。
总的来说,快的玄学。
可以写出如下代码:
\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{LL Pollard\_Rho}\OperatorTok{(}\NormalTok{LL n}\OperatorTok{)\{}
\NormalTok{ LL c }\OperatorTok{=}\NormalTok{ rand}\OperatorTok{(}\NormalTok{n }\OperatorTok{{-}} \DecValTok{1}\OperatorTok{)} \OperatorTok{+} \DecValTok{1}\OperatorTok{;}
\NormalTok{ LL L }\OperatorTok{=} \DecValTok{0}\OperatorTok{;}
\ControlFlowTok{for}\OperatorTok{(}\DataTypeTok{int}\NormalTok{ s }\OperatorTok{=} \DecValTok{0}\OperatorTok{,}\NormalTok{ t }\OperatorTok{=} \DecValTok{1}\OperatorTok{;} \OperatorTok{;}\NormalTok{ s }\OperatorTok{\textless{}\textless{}=} \DecValTok{1}\OperatorTok{,}\NormalTok{ t }\OperatorTok{\textless{}\textless{}=} \DecValTok{1}\OperatorTok{)\{}
\NormalTok{ LL R }\OperatorTok{=}\NormalTok{ L}\OperatorTok{,}\NormalTok{ M }\OperatorTok{=} \DecValTok{1}\OperatorTok{;}
\ControlFlowTok{for}\OperatorTok{(}\DataTypeTok{int}\NormalTok{ i }\OperatorTok{=}\NormalTok{ s}\OperatorTok{,}\NormalTok{ step }\OperatorTok{=} \DecValTok{1}\OperatorTok{;}\NormalTok{ i }\OperatorTok{\textless{}}\NormalTok{ t}\OperatorTok{;}\NormalTok{ i}\OperatorTok{++,}\NormalTok{ step }\OperatorTok{++)\{}
\NormalTok{ R }\OperatorTok{=}\NormalTok{ f}\OperatorTok{(}\NormalTok{R}\OperatorTok{,}\NormalTok{ c}\OperatorTok{,}\NormalTok{ n}\OperatorTok{);}
\NormalTok{ M }\OperatorTok{=}\NormalTok{ mul}\OperatorTok{(}\NormalTok{M}\OperatorTok{,}\NormalTok{ abs}\OperatorTok{(}\NormalTok{L }\OperatorTok{{-}}\NormalTok{ R}\OperatorTok{),}\NormalTok{ n}\OperatorTok{);}
\ControlFlowTok{if}\OperatorTok{(}\NormalTok{step }\OperatorTok{\%} \DecValTok{1000} \OperatorTok{==} \DecValTok{0}\OperatorTok{)} \OperatorTok{\{}\NormalTok{ LL g }\OperatorTok{=}\NormalTok{ gcd}\OperatorTok{(}\NormalTok{M}\OperatorTok{,}\NormalTok{ n}\OperatorTok{);} \ControlFlowTok{if}\OperatorTok{(}\NormalTok{g }\OperatorTok{\textgreater{}} \DecValTok{1}\OperatorTok{)} \ControlFlowTok{return}\NormalTok{ g}\OperatorTok{;} \OperatorTok{\}}
\CommentTok{// 不一定必须等到一轮计算结束之后再计算 gcd 可能中间就已经出现答案了,中间每隔一段时间算几次,可以在出现答案之后快速返回答案。}
\OperatorTok{\}}
\NormalTok{ LL g }\OperatorTok{=}\NormalTok{ gcd}\OperatorTok{(}\NormalTok{M}\OperatorTok{,}\NormalTok{ n}\OperatorTok{);} \ControlFlowTok{if}\OperatorTok{(}\NormalTok{g }\OperatorTok{\textgreater{}} \DecValTok{1}\OperatorTok{)} \ControlFlowTok{return}\NormalTok{ g}\OperatorTok{;}
\NormalTok{ L }\OperatorTok{=}\NormalTok{ R}\OperatorTok{;}
\OperatorTok{\}}
\OperatorTok{\}}
\DataTypeTok{void}\NormalTok{ factorize}\OperatorTok{(}\NormalTok{LL n}\OperatorTok{)\{}
\ControlFlowTok{if}\OperatorTok{(}\NormalTok{Prime}\OperatorTok{(}\NormalTok{n}\OperatorTok{))} \OperatorTok{\{}\NormalTok{ ans }\OperatorTok{=}\NormalTok{ max}\OperatorTok{(}\NormalTok{ans}\OperatorTok{,}\NormalTok{ n}\OperatorTok{);} \ControlFlowTok{return} \OperatorTok{;} \OperatorTok{\}}
\NormalTok{ LL g }\OperatorTok{=}\NormalTok{ n}\OperatorTok{;}
\ControlFlowTok{while}\OperatorTok{(}\NormalTok{g }\OperatorTok{==}\NormalTok{ n}\OperatorTok{)}\NormalTok{ g }\OperatorTok{=}\NormalTok{ Pollard\_Rho}\OperatorTok{(}\NormalTok{n}\OperatorTok{);}
\NormalTok{ factorize}\OperatorTok{(}\NormalTok{g}\OperatorTok{);}\NormalTok{ factorize}\OperatorTok{(}\NormalTok{n }\OperatorTok{/}\NormalTok{ g}\OperatorTok{);}
\OperatorTok{\}}
\end{Highlighting}
\end{Shaded}
\subsection{阶}\label{ux9636}
若 \(a \perp P\),则使 \(a^k \equiv 1 \pmod P\) 成立的最小的 \(k\),称为
\(a\) 关于模 \(P\) 的阶,记作 \(\text{ord}_m(a)\)。
求法: 根据欧拉定理,\(\text{ord}_m(a) \mid \varphi(m)\),先求出
\(\varphi(m)\),记作 \(t\)。
\begin{itemize}
\tightlist
\item
枚举 \(t\) 的每一个因子,检查是否满足阶的性质。
\item
对 \(t\)
质因数拆分,以此考虑每个质因子,试图减少质因子的幂次——即如果减少了某个质因子的幂次,\(a^t \equiv 1 \pmod P\)
依然成立,那么就直接减少即可。感性理解一下肯定是对的。
\(\mathcal{O}(\log \varphi(P))\)
\end{itemize}
\subsection{原根}\label{ux539fux6839}
若 \(\text{ord}_mg = \varphi(m)\) 则称 \(g\) 为 \(m\) 的原根。
通俗的讲:称 \(g\) 为 \(P\) 的原根,当且仅当
\(\{g^0, g^1, g^2, \dots, g^{\varphi(P) - 1}\}\) 数两两不同。
即对于任意一个和 \(P\) 互质的数字,在 \(P\) 的剩余系下,都可以表示成
\(g^t\) 的形式。
一个模数存在原根,当且仅当这个模数为 \(2, 4, p^s, 2p^s\),\(p\)
为奇素数。
判断一个数是否为原根:根据定义可以考虑枚举每一个
\(t \in [1, \varphi(P) - 1]\) 判断 \(g^t\) 是否都不等于
\(1\)。事实上,根据欧拉定理,可能使 \(g^t \equiv 1 \pmod P\) 的 \(t\)
只可能是是 \(\varphi(P)\) 的约数。枚举 \(\varphi(P)\)
的每个约数,进行判断即可,复杂度为 \(\mathcal{O}(\sqrt{\varphi(P)})\)
原根的求法:若一个数 \(m\) 存在原根,其最小原根大概在 \(m^{1/4}\)
级别。枚举原根再判断一般是可以接受的。
\begin{Shaded}
\begin{Highlighting}[]
\DataTypeTok{bool}\NormalTok{ JudgePrimitiveRoot}\OperatorTok{(}\DataTypeTok{int}\NormalTok{ g}\OperatorTok{,} \DataTypeTok{int}\NormalTok{ P}\OperatorTok{)\{}
\DataTypeTok{int}\NormalTok{ phi }\OperatorTok{=}\NormalTok{ P }\OperatorTok{{-}} \DecValTok{1}\OperatorTok{;}
\ControlFlowTok{for}\OperatorTok{(}\DataTypeTok{int}\NormalTok{ i }\OperatorTok{=} \DecValTok{2}\OperatorTok{;}\NormalTok{ i }\OperatorTok{*}\NormalTok{ i }\OperatorTok{\textless{}=}\NormalTok{ phi}\OperatorTok{;}\NormalTok{ i}\OperatorTok{++)\{}
\ControlFlowTok{if}\OperatorTok{(}\NormalTok{phi }\OperatorTok{\%}\NormalTok{ i }\OperatorTok{!=} \DecValTok{0}\OperatorTok{)} \ControlFlowTok{continue}\OperatorTok{;}
\ControlFlowTok{if}\OperatorTok{(}\NormalTok{pow}\OperatorTok{(}\NormalTok{g}\OperatorTok{,}\NormalTok{ i}\OperatorTok{,}\NormalTok{ P}\OperatorTok{)} \OperatorTok{==} \DecValTok{1}\OperatorTok{)} \ControlFlowTok{return} \KeywordTok{false}\OperatorTok{;}
\ControlFlowTok{if}\OperatorTok{(}\NormalTok{pow}\OperatorTok{(}\NormalTok{g}\OperatorTok{,}\NormalTok{ P }\OperatorTok{/}\NormalTok{ i}\OperatorTok{,}\NormalTok{ P}\OperatorTok{)} \OperatorTok{==} \DecValTok{1}\OperatorTok{)} \ControlFlowTok{return} \KeywordTok{false}\OperatorTok{;}
\OperatorTok{\}}
\ControlFlowTok{return} \KeywordTok{true}\OperatorTok{;}
\OperatorTok{\}}
\DataTypeTok{int}\NormalTok{ GetPrimitiveRoot}\OperatorTok{(}\DataTypeTok{int}\NormalTok{ P}\OperatorTok{)\{} \ControlFlowTok{for}\OperatorTok{(}\DataTypeTok{int}\NormalTok{ i }\OperatorTok{=} \DecValTok{2}\OperatorTok{;} \OperatorTok{;}\NormalTok{ i}\OperatorTok{++)} \ControlFlowTok{if}\OperatorTok{(}\NormalTok{JudgePrimitiveRoot}\OperatorTok{(}\NormalTok{i}\OperatorTok{,}\NormalTok{ P}\OperatorTok{))} \ControlFlowTok{return}\NormalTok{ i}\OperatorTok{;} \OperatorTok{\}}
\end{Highlighting}
\end{Shaded}
地位相当于自然数域下的唯一分解定理。
\end{document}
[INFO] [makePDF] LaTeX run number 1
[INFO] [makePDF] LaTeX output
This is XeTeX, Version 3.141592653-2.6-0.999995 (TeX Live 2023/Debian) (preloaded format=xelatex)
restricted \write18 enabled.
entering extended mode
(.../input.tex
LaTeX2e <2023-11-01> patch level 1
L3 programming layer <2024-01-22>
(.../article.cls
Document Class: article 2023/05/17 v1.4n Standard LaTeX document class
(.../[system]
(.../xcolor.sty
(.../color.cfg)
(.../xetex.def)
(.../[system]
(.../geometry.sty
(.../keyval.sty)
(.../ifvtex.sty
(.../iftex.sty)))
(.../amsmath.sty
For additional information on amsmath, use the `?' option.
(.../amstext.sty
(.../amsgen.sty))
(.../amsbsy.sty)
(.../amsopn.sty))
(.../amssymb.sty
(.../amsfonts.sty))
(.../unicode-math.sty
(.../expl3.sty
(.../l3backend-xetex.def))
(.../unicode-math-xetex.sty
(.../xparse.sty)
(.../l3keys2e.sty)
(.../fontspec.sty
(.../fontspec-xetex.sty
(.../fontenc.sty)
(.../fontspec.cfg)))
(.../fix-cm.sty
(.../ts1enc.def))
(.../unicode-math-table.tex)))
(.../lmodern.sty)
(.../xeCJK.sty
(.../ctexhook.sty)
(.../xtemplate.sty)
(.../xeCJK.cfg))
(.../upquote.sty
(.../textcomp.sty))
(.../microtype.sty
(.../etoolbox.sty)
(.../microtype-xetex.def)
(.../microtype.cfg))
(.../setspace.sty)
(.../parskip.sty
(.../kvoptions.sty
(.../ltxcmds.sty)
(.../kvsetkeys.sty)))
(.../fancyvrb.sty)
(.../soul.sty
(.../soul-ori.sty)
(.../infwarerr.sty)
(.../etexcmds.sty))
(.../xeCJKfntef.sty
(.../ulem.sty))
(.../bookmark.sty
(.../hyperref.sty
(.../kvdefinekeys.sty)
(.../pdfescape.sty
(.../pdftexcmds.sty))
(.../hycolor.sty)
(.../auxhook.sty)
(.../nameref.sty
(.../refcount.sty)
(.../gettitlestring.sty))
(.../pd1enc.def)
(.../intcalc.sty)
(.../puenc.def)
(.../url.sty)
(.../bitset.sty
(.../bigintcalc.sty))
(.../atbegshi-ltx.sty))
(.../hxetex.def
(.../stringenc.sty)
(.../rerunfilecheck.sty
(.../atveryend-ltx.sty)
(.../uniquecounter.sty)))
(.../bkm-dvipdfm.def))
(.../xurl.sty)
No file input.aux.
*geometry* driver: auto-detecting
*geometry* detected driver: xetex
(.../mt-LatinModernRoman.cfg)
Package hyperref Warning: Rerun to get /PageLabels entry.
(.../omllmm.fd)
(.../umsa.fd)
(.../mt-msa.cfg)
(.../umsb.fd)
(.../mt-msb.cfg)
Package xeCJK Warning: Unknown CJK family `\CJKttdefault' is being ignored.
(xeCJK)
(xeCJK) Try to use `\setCJKmonofont[<...>]{<...>}' to define
(xeCJK) it.
LaTeX Font Warning: Font shape `TU/ARPLUKaiCN(0)/b/n' undefined
(Font) using `TU/ARPLUKaiCN(0)/m/n' instead on input line 128.
[1] [2]
LaTeX Font Warning: Font shape `TU/ARPLUKaiCN(0)/m/it' undefined
(Font) using `TU/ARPLUKaiCN(0)/m/n' instead on input line 290.
[3] [4]
Underfull \hbox (badness 10000) in paragraph at lines 371--373
Overfull \hbox (282.80313pt too wide) in paragraph at lines 408--408
[] \TU/lmtt/m/n/10 for(int i = 0, V = 1; i <= k; i++, V = V *1ll* a % P) if(
V == b) return i; // changed: It must be executed before []if(b % D != 0) retur
n -1[]).[]
[5] [6] [7] [8] [9] (.../input.aux)
LaTeX Font Warning: Some font shapes were not available, defaults substituted.
LaTeX Warning: Label(s) may have changed. Rerun to get cross-references right.
)
(see the transcript file for additional information)
Output written on .../input.pdf (9 pages).
Transcript written on .../input.log.
[INFO] [makePDF] Rerun needed
Package hyperref Warning: Rerun to get /PageLabels entry.
LaTeX Warning: Label(s) may have changed. Rerun to get cross-references right.
[INFO] [makePDF] LaTeX run number 2
[INFO] [makePDF] LaTeX output
This is XeTeX, Version 3.141592653-2.6-0.999995 (TeX Live 2023/Debian) (preloaded format=xelatex)
restricted \write18 enabled.
entering extended mode
(.../input.tex
LaTeX2e <2023-11-01> patch level 1
L3 programming layer <2024-01-22>
(.../article.cls
Document Class: article 2023/05/17 v1.4n Standard LaTeX document class
(.../[system]
(.../xcolor.sty
(.../color.cfg)
(.../xetex.def)
(.../[system]
(.../geometry.sty
(.../keyval.sty)
(.../ifvtex.sty
(.../iftex.sty)))
(.../amsmath.sty
For additional information on amsmath, use the `?' option.
(.../amstext.sty
(.../amsgen.sty))
(.../amsbsy.sty)
(.../amsopn.sty))
(.../amssymb.sty
(.../amsfonts.sty))
(.../unicode-math.sty
(.../expl3.sty
(.../l3backend-xetex.def))
(.../unicode-math-xetex.sty
(.../xparse.sty)
(.../l3keys2e.sty)
(.../fontspec.sty
(.../fontspec-xetex.sty
(.../fontenc.sty)
(.../fontspec.cfg)))
(.../fix-cm.sty
(.../ts1enc.def))
(.../unicode-math-table.tex)))
(.../lmodern.sty)
(.../xeCJK.sty
(.../ctexhook.sty)
(.../xtemplate.sty)
(.../xeCJK.cfg))
(.../upquote.sty
(.../textcomp.sty))
(.../microtype.sty
(.../etoolbox.sty)
(.../microtype-xetex.def)
(.../microtype.cfg))
(.../setspace.sty)
(.../parskip.sty
(.../kvoptions.sty
(.../ltxcmds.sty)
(.../kvsetkeys.sty)))
(.../fancyvrb.sty)
(.../soul.sty
(.../soul-ori.sty)
(.../infwarerr.sty)
(.../etexcmds.sty))
(.../xeCJKfntef.sty
(.../ulem.sty))
(.../bookmark.sty
(.../hyperref.sty
(.../kvdefinekeys.sty)
(.../pdfescape.sty
(.../pdftexcmds.sty))
(.../hycolor.sty)
(.../auxhook.sty)
(.../nameref.sty
(.../refcount.sty)
(.../gettitlestring.sty))
(.../pd1enc.def)
(.../intcalc.sty)
(.../puenc.def)
(.../url.sty)
(.../bitset.sty
(.../bigintcalc.sty))
(.../atbegshi-ltx.sty))
(.../hxetex.def
(.../stringenc.sty)
(.../rerunfilecheck.sty
(.../atveryend-ltx.sty)
(.../uniquecounter.sty)))
(.../bkm-dvipdfm.def))
(.../xurl.sty)
(.../input.aux)
*geometry* driver: auto-detecting
*geometry* detected driver: xetex
(.../mt-LatinModernRoman.cfg)
(.../omllmm.fd)
(.../umsa.fd)
(.../mt-msa.cfg)
(.../umsb.fd)
(.../mt-msb.cfg)
Package xeCJK Warning: Unknown CJK family `\CJKttdefault' is being ignored.
(xeCJK)
(xeCJK) Try to use `\setCJKmonofont[<...>]{<...>}' to define
(xeCJK) it.
LaTeX Font Warning: Font shape `TU/ARPLUKaiCN(0)/b/n' undefined
(Font) using `TU/ARPLUKaiCN(0)/m/n' instead on input line 128.
[1] [2]
LaTeX Font Warning: Font shape `TU/ARPLUKaiCN(0)/m/it' undefined
(Font) using `TU/ARPLUKaiCN(0)/m/n' instead on input line 290.
[3] [4]
Underfull \hbox (badness 10000) in paragraph at lines 371--373
Overfull \hbox (282.80313pt too wide) in paragraph at lines 408--408
[] \TU/lmtt/m/n/10 for(int i = 0, V = 1; i <= k; i++, V = V *1ll* a % P) if(
V == b) return i; // changed: It must be executed before []if(b % D != 0) retur
n -1[]).[]
[5] [6] [7] [8] [9] (.../input.aux)
LaTeX Font Warning: Some font shapes were not available, defaults substituted.
)
(see the transcript file for additional information)
Output written on .../input.pdf (9 pages).
Transcript written on .../input.log.