关于 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.