我要成为数论king！！！

唯一分解定理

牛牛的幂运算

/*

设a = p^k1   c = p^k2
（0）因为1的幂次永远是1,会使得代码错误,要特判 
（1）k1==k2，此时a == c,那么必然有b == d。答案等于n*n 
（2）k1<k2,
		(2.1)p<=根号n 爆搜1~根号n 爆搜k1与k2 
			（2.1.1）k1与k2不互素,即gcd(k1,k2)!=1,那么在他们等于gcd(k1,k2)的时候
			答案已经被贡献。
				（2.1.1.1）答案怎么贡献？此时我们已经得到k1和k2了,就看有多少个k1*b==k2*d。n/k1：有多少个 
			（2.1.2）k1与k2互素，直接贡献个数 
		(2.2)p == n,此时k1 == k2 == 1.属于（1） 
（3）k2<k1，答案等于（2）。 找到(2)的答案乘2就行
*/
#include<bits/stdc++.h>
using namespace std;
#define fi first
#define se second
#define U unsigned
#define P std::pair<int,int>
#define LL long long
#define pb push_back
#define MP std::make_pair
#define V std::vector<int>
#define all(x) x.begin(),x.end()
#define CLR(i,a) memset(i,a,sizeof(i))
#define FOR(i,a,b) for(int i = a;i <= b;++i)
#define ROF(i,a,b) for(int i = a;i >= b;--i)
#define DEBUG(x) std::cerr << #x << '=' << x << std::endl
#define INF 0x3f3f3f3f
#define MOD 1000000007

int main(){
	int n;scanf("%d",&n);
	LL ans = 1ll*n*n%MOD;
	ans += 1ll*(n-1)*n%MOD;
	LL ghn = sqrt(n);
	for(LL p = 2;p<=ghn;++p){
		for(LL k1 = 1,a = p;a<=n;a*=p,k1++){
			for(LL k2 = 1;k2<k1;k2++){
				if(__gcd(k1,k2)==1){
					ans = (ans + n/k1*2)%MOD;
				}
			}
		}
	}
	printf("%lld\n",ans);
}

欧拉函数

性质

1. $\phi (n)=n{\prod }_{i=1}^s(1-\frac{1}{p_i})$
2. $n={\sum }_{d|n}\phi (d)$
3. 积性函数：if $gcd(a,b)=1$ then $\phi (a\times b)=\phi (a)\times \phi (b)$

筛法

1. 求单个欧拉函数的值
   （可以用Pollard Rho算法优化）
   time complexity: $O(n)$

int euler_phi(int n) {
  int m = int(sqrt(n + 0.5));
  int ans = n;
  for (int i = 2; i <= m; i++)
    if (n % i == 0) {
      ans = ans / i * (i - 1);
      while (n % i == 0) n /= i;
    }
  if (n > 1) ans = ans / n * (n - 1);
  return ans;
}

2. 筛一堆

(a) time complexity: $O(nloglogn)$，基于埃氏筛，不需要记录质数

// 不用筛质数
void phi_table(int n, int* phi) {
  for (int i = 2; i <= n; i++) phi[i] = 0;
  phi[1] = 1;
  for (int i = 2; i <= n; i++)
    if (!phi[i])
      for (int j = i; j <= n; j += i) {
        if (!phi[j]) phi[j] = j;
        phi[j] = phi[j] / i * (i - 1);
      }
}

(b) time complexity: $O(n)$，基于欧拉筛（线性筛），需要记录质数

void init() {
  phi[1] = 1;
  for (int i = 2; i < MAXN; ++i) {
    if (!vis[i]) {
      phi[i] = i - 1;
      pri[cnt++] = i;
    }
    for (int j = 0; j < cnt; ++j) {
      if (1ll * i * pri[j] >= MAXN) break;
      vis[i * pri[j]] = 1;
      if (i % pri[j]) {
        phi[i * pri[j]] = phi[i] * (pri[j] - 1);
      } else {
        phi[i * pri[j]] = phi[i] * pri[j];
        break;
      }
    }
  }
}

欧拉定理&费马小定理

• 欧拉定理：
  若$\gcd (a,m)=1$，则$a^{\phi (m)}\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}m)$.

• 费马小定理：
  欧拉定理的$m$为质数$p$，则$\phi (p)=p-1$，
  故有$a^{p-1}\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}p)$.

• 扩展欧拉定理：
  $a$为任意整数，$b$和$m$是正整数。
  $$a^b\equiv \{\begin{array}{ll}{a}^{b\mathrm{m}\mathrm{o}\mathrm{d}\phi (m)}& \gcd (a,m)=1\\ a^{b\mathrm{m}\mathrm{o}\mathrm{d}\phi (m)+\phi (m)}& \gcd (a,m)\ne 1\end{array}(\mathrm{m}\mathrm{o}\mathrm{d}m)$$

练习题

TIMUS #1673 “Admission to Exam”[Difficulty: High]
URAL上好难的一道题。
这道题分为两步，第一步，推导出题目要求的是欧拉函数的反函数；第二步，求这个反函数。

题意：$[1,N]$中有$k$个数，其倍数通过范围是$[1,N^2]$的模$N$的剩余类。给定$k$，求最小$N$。

思路参考：cf blog
如果我们选一个$\gcd (x,N)>1$的$x$。设$lcm(x,N)=x\times k$（因为$k<N$，$lcm(x,N)$也一定小于$N^2$），因为$\gcd (x,N)>1$，所以$k<N$。所以直到$x$的第$k$个倍数，我们不可能通过$N$的剩余类（因为$N$的剩余类有$N$个）。然而，在$lcm(x,N)$之后每加多一个$x$，即$x$的第$k+1$个之后的倍数，都与前$k$个同余（因为$lcm(x,N)\mathrm{m}\mathrm{o}\mathrm{d}N=0$）。所以，$x$的倍数不能通过$N$的剩余类。反之，我们可以发现，$\gcd (x,N)=1$的时候，循环节在第$N$个倍数之外，能够通过所有的$N$剩余类。（如果不能通过所有，就说明$\gcd (x,N)\ne 1$，归谬）因此我们要找的是$\gcd (x,N)=1$的$x$的数目。

因此，我们知道phi(N)=k，求这个最小的N，k范围2e9
就是求欧拉函数的反函数。

那么如何求欧拉函数的反函数？当函数值范围有2e9的时候。
思路：
（以下思路中的a指代题目中的k，x指代题目中的n）
$$\phi (x)=a=x\prod (1-\frac{1}{p})=p_1^{a_1}{p}_2^{a_2}...p_m^{a_m}\prod \frac{p_i-1}{p_i}$$$$=p_1^{a_1-1}{p}_2^{a_2-1}...p_m^{a_m-1}\prod (p_i-1)=\prod _i^k{p}_i^{a_i}(p_i-1)$$
代码参考：BZOJ 3643
预处理$\sqrt{a}$的质数，（因为$vis$数组不能开2e9那么大，素数筛$O(n)$要用到这个数组），然后对于可能是$a$的大于$\sqrt{a}$的质因数的数，暴力特判是不是质数。
DFS爆搜质数组合，检查$(p_1-1)\times (p_2-1)\times ..\times (p_m-1)$是否可以被$a$整除（因为$x$是整数，所以构成$a$的质数-1一定可以被$a$整除（见上公式）），并检查商是否是这些质数的倍数。
若是，则找到了x的质数。dfs过程中记录$p_1\times p_2\times ...\times p_m$，把这些数乘上刚才的商，即为x。更新最小x。

但是我们编码的时候更加简化了。遇到质数，如果当前能整除他-1，就说明这个质数有可能是x的质因数。那么我们尝试选取他，进行dfs，如果能再除他本身，就除，再dfs，最后商为1，就说明除尽了这些质数的组合，就可以更新答案。

#include<bits/stdc++.h>

using namespace std;
typedef long long ll;
const int N = 1e5+5;
ll ans=200000000000;
int pri[N],vis[N],cnt,m;

void get_prime(){
    for(int i=2;i<=50000;++i){
        if(!vis[i]) pri[++cnt] = i;
        for(int j=1;j<=cnt&&i*pri[j]<=50000;++j){
            vis[i*pri[j]]=1;
            if(i%pri[j]==0) break;
        }
    }
    // 因为空间问题，只能筛出sqrt n的质数. 然后我们在dfs过程中对huge prime进行特判。
}

bool isprime(int n){
    int k=(int)sqrt(n+0.5);
    for(int i=2;i<=k;++i) if(n%i==0) return 0;
    return 1;
}

void dfs(int k,int rem,ll pro){
    // now at k th prime, the previous product is pro, the remainder is rem
    if(pro>=ans) return;  // if the current product is already larger, cut the edge
    if(rem==1){ans = pro;return;} // if the product is not larger, and it can divide n, update
    if(rem>m&&isprime(rem+1)) ans = min(ans,pro*(rem+1)); // huge prime
    for(int i=k+1;pri[i]-1<=m;++i){
        if(pri[i]-1>rem) break;
        if(rem%(pri[i]-1)==0){
            int x = rem/(pri[i]-1);ll y = pro*pri[i];
            dfs(i,x,y);
            while(x%pri[i]==0){
                x/=pri[i];y*=pri[i];   // multiply the prime
                dfs(i,x,y);
            }
        }
    }
}

int main(){
    int k;scanf("%d",&k);
    m = (int)sqrt(k+0.5);
    get_prime();
    dfs(0,k,1);
    if(ans==200000000000) cout<<0<<endl;
    else cout<<ans<<endl;
}

逆元

• 单个求法：扩欧、费马小定理（要求p是质数），时间复杂度都是$O(logn)$
• 存在条件：(a,m)=1.化成二元一次方程可证。

// 扩欧
void exgcd(int a, int b, int& x, int& y) {
  if (b == 0) {
    x = 1, y = 0;
    return;
  }
  exgcd(b, a % b, y, x);
  y -= a / b * x;
}
// 快速幂（费马小定理）
inline int qpow(long long a, int b) {
  int ans = 1;
  a = (a % p + p) % p;
  for (; b; b >>= 1) {
    if (b & 1) ans = (a * ans) % p;
    a = (a * a) % p;
  }
  return ans;
}

• 批量求法：1~n线性求$O(n)$、n个独立数求$O(n+logp)$

// 线性
inv[1] = 1;
for (int i = 2; i <= n; ++i) {
  inv[i] = (long long)(p - p / i) * inv[p % i] % p;
}

// 独立数
s[0] = 1;
for (int i = 1; i <= n; ++i) s[i] = s[i - 1] * a[i] % p;
sv[n] = qpow(s[n], p - 2);
// 当然这里也可以用 exgcd 来求逆元,视个人喜好而定.
for (int i = n; i >= 1; --i) sv[i - 1] = sv[i] * a[i] % p;
for (int i = 1; i <= n; ++i) inv[i] = sv[i] * s[i - 1] % p;

质因数分解

Pollard Rho

#include <bits/stdc++.h>

using namespace std;

typedef long long ll;
#define lll __int128

int t;
ll max_factor, n;

ll gcd(ll a, ll b) {
  if (b == 0) return a;
  return gcd(b, a % b);
}

ll quick_pow(ll x, ll p, ll mod) {
  ll ans = 1;
  while (p) {
    if (p & 1) ans = (lll)ans * x % mod;
    x = (lll)x * x % mod;
    p >>= 1;
  }
  return ans;
}

bool Miller_Rabin(ll p) {
  if (p < 2) return 0;
  if (p == 2) return 1;
  if (p == 3) return 1;
  ll d = p - 1, r = 0;
  while (!(d & 1)) ++r, d >>= 1;
  for (ll k = 0; k < 10; ++k) {
    ll a = rand() % (p - 2) + 2;
    ll x = quick_pow(a, d, p);
    if (x == 1 || x == p - 1) continue;
    for (int i = 0; i < r - 1; ++i) {
      x = (lll)x * x % p;
      if (x == p - 1) break;
    }
    if (x != p - 1) return 0;
  }
  return 1;
}

ll f(ll x, ll c, ll n) { return ((lll)x * x + c) % n; }

ll Pollard_Rho(ll x) {
  ll s = 0, t = 0;
  ll c = (ll)rand() % (x - 1) + 1;
  int step = 0, goal = 1;
  ll val = 1;
  for (goal = 1;; goal <<= 1, s = t, val = 1) {
    for (step = 1; step <= goal; ++step) {
      t = f(t, c, x);
      val = (lll)val * abs(t - s) % x;
      if ((step % 127) == 0) {
        ll d = gcd(val, x);
        if (d > 1) return d;
      }
    }
    ll d = gcd(val, x);
    if (d > 1) return d;
  }
}

void fac(ll x) {
  if (x <= max_factor || x < 2) return;
  if (Miller_Rabin(x)) {
    max_factor = max(max_factor, x);
    return;
  }
  ll p = x;
  while (p >= x) p = Pollard_Rho(x);
  while ((x % p) == 0) x /= p;
  fac(x), fac(p);
}

素性测试miller-rabin

//
// Created by Artist on 2021/7/29.
//
#include<bits/stdc++.h>

using namespace std;
typedef long long ll;

ll Mul(ll a, ll b, ll P) {
    ll s = a * b - ll((long double) a * b / P + 0.5) * P;
    return s >= 0 ? s : s + P;
}

ll quickPow(ll a, ll b, ll P) {
    ll s = 1;

    for (; b; b >>= 1, a = Mul(a, a, P))
        if (b & 1)
            s = Mul(s, a, P);

    return s % P;
}

const int p[] = {3, 5, 7, 11, 13, 17, 19, 23};

bool millerRabin(ll n) {
    if (n <= 1 || !(n & 1)) return n == 2;
    ll cnt = 0, xs = n - 1;
    while (!(xs & 1)) ++cnt, xs >>= 1;
    ll now, pre;
    for (int i = 0; i < 8; ++i) {
        if (n == p[i]) return true;
        if (n % p[i] == 0) return false;
        now = quickPow(p[i], xs, n);
        for (int j = 1; j <= cnt; ++j) {
            pre = now, now = Mul(now, now, n);
            if (now == 1 && (pre != 1 && pre != n - 1)) return false;
        }
        if (now != 1) return false;
    }
    return true;
}

int main() {
    ll n;
    srand(time(NULL));
    while (~scanf("%lld", &n)) {
        printf("%c\n", millerRabin(n) ? 'Y' : 'N');
    }
}

BSGS算法

基本思想是哈希映射和分块。
板子

//
// Created by artist on 2021/7/29.
//

#include<bits/stdc++.h>

using namespace std;
typedef long long ll;

unordered_map<int,int> hs;

int qpow(int a,int n,int p){
    int ans = 1;
    while(n){
        if(n&1) ans=(int)(1ll*ans*a%p);
        a=(int)(1ll*a*a%p);
        n>>=1;
    }
    return ans;
}

int BSGS(int a,int b,int p){
    b%=p;
    int t=(int)ceil(sqrt(p));
    int aj=1;
    for(int i=0;i<t;++i){
        has[(int)(1ll*aj*b%p)]=i+1;
        aj=(int)(1ll*aj*a%p);
    }
    if(!aj) return !b?1:-1;
    int at=qpow(a,t,p),ati=1;
    for(int i=0;i<=t;++i){
        int j=has[ati]-1;
        if(j>=0&&t*i-j>=0) return t*i-j;
        ati=(int)(1ll*at*ati%p);
    }
    return -1;
}

int main() {
    int p,b,n;scanf("%d%d%d",&p,&b,&n);
    int res = BSGS(b,n,p);
    if(~res) printf("%d\n",res);
    else printf("no solution\n");
}

lucas定理

//
// Created by artist on 2021/8/24.
//


#include<bits/stdc++.h>

using namespace std;
typedef long long ll;
#define mkp make_pair
#define fi first
#define se second
#define pb push_back
#define all(x) x.begin(),x.end()
#define DB1(args...) do { cout << #args << " : "; dbg(args); } while (0)

void dbg() { std::cout << "  #\n"; }

template<typename T, typename...Args>
void dbg(T a, Args...args) {
    std::cout << a << ' ';
    dbg(args...);
}
const int maxn = 2e5+4;
int fac[maxn],invf[maxn];

int qpow(int a,int n,int p){
    int ans=1;
    while(n){
        if(n&1) ans=1ll*ans*a%p;
        a=1ll*a*a%p;
        n>>=1;
    }
    return ans;
}

int C(int n,int m,int p){
    return 1ll*fac[n]*invf[m]%p*invf[n-m]%p;
}

int lucas(int n,int m,int p){
  //  DB1(n,m,p);
    if(m==0) return 1;
    return 1ll*C(n%p,m%p,p)*lucas(n/p,m/p,p)%p;
}

/*
1
100000 100000 99991
 */

signed main() {
    int t;scanf("%d",&t);
    fac[1]=1;fac[0]=1;
    invf[0]=1;
    while(t--){
        int n,m,p;scanf("%d%d%d",&n,&m,&p);
        for(int i=2;i<=min(n+m,p-1);++i) {
            fac[i]=1ll*fac[i-1]*i%p;
        }
        invf[min(n+m,p-1)]=qpow(fac[min(n+m,p-1)],p-2,p);
        for(int i=max(1,min(n+m-1,p-2));i>=1;--i) invf[i]=1ll*invf[i+1]*(i+1)%p;
        printf("%d\n",lucas(n+m,m,p));
    }
}

扩展lucas

套板子，缝合怪
洛谷

//
// Created by artist on 2021/8/24.
//


#include<bits/stdc++.h>

using namespace std;
typedef long long ll;
#define mkp make_pair
#define fi first
#define se second
#define pb push_back
#define all(x) x.begin(),x.end()
#define DB1(args...) do { cout << #args << " : "; dbg(args); } while (0)

void dbg() { std::cout << "  #\n"; }

template<typename T, typename...Args>
void dbg(T a, Args...args) {
    std::cout << a << ' ';
    dbg(args...);
}

ll qpow(ll a,ll n,ll P){
    ll ans = 1;
   // DB1(a,n,P);
    while(n){
        if(n&1) ans=ans*a%P;
        a=a*a%P;
        n>>=1;
    }
    return ans;
}

void exgcd(ll a,ll b,ll &x,ll &y){
 //   DB1(a,b);
    if(!b){
        x=1,y=0;
        return;
    }
    exgcd(b,a%b,y,x);
    y-=a/b*x;
}

ll inverse(ll a,ll P){
    ll x,y;
  //  DB1(a,P);
    exgcd(a,P,x,y);
    return (x%P+P)%P;
}

ll calc(ll n,ll x,ll P){
    if(!n) return 1;
    ll s=1;
    for(ll i=1;i<=P;++i)
        if(i%x) s=s*i%P;
    s=qpow(s,n/P,P);
    for(ll i=n/P*P+1;i<=n;++i)
        if(i%x) s=i%P*s%P;
    return s*calc(n/x,x,P)%P;
}

ll multilucas(ll m,ll n,ll x,ll P){
    int cnt = 0;
    for(ll i=m;i;i/=x) cnt+=i/x;
    for(ll i=n;i;i/=x) cnt-=i/x;
    for(ll i=m-n;i;i/=x) cnt-=i/x;
    return qpow(x,cnt,P)%P*calc(m,x,P)%P*inverse(calc(n,x,P),P)%P*inverse(calc(m-n,x,P),P)%P;
}

ll CRT(int k, ll* a, ll* r) {
    ll n = 1, ans = 0;
    for (int i = 1; i <= k; i++) n = n * r[i];
    for (int i = 1; i <= k; i++) {
        ll m = n / r[i], b, y;
        exgcd(m, r[i], b, y);  // b * m mod r[i] = 1
        ans = (ans + a[i] * m * b % n) % n;
    }
    return (ans % n + n) % n;
}

ll exlucas(ll m,ll n,ll P){
    int cnt = 0;
    ll p[20],a[20];
    for(ll i=2;i*i<=P;++i){
        if(P%i==0){
            p[++cnt]=1;
            while(P%i==0) p[cnt]=p[cnt]*i,P/=i;
            a[cnt]=multilucas(m,n,i,p[cnt]);
        }
    }
    if(P>1) p[++cnt] = P,a[cnt] = multilucas(m,n,P,P);
    return CRT(cnt,a,p);
}

signed main() {
    ll n,m,p;scanf("%lld%lld%lld",&n,&m,&p);
    printf("%lld\n",exlucas(n,m,p));
}