数论：简单的模板

/*
Copyright: xjjppm
Author: xjjppm
Date: 08-08-17 11:36
Description: Number Theory
*/

#include <map>
#include <cmath>
#include <cstdio>
#include <cstring>
#include <algorithm>

inline int input() {
char c=getchar();int x=0,a=1;
for(;c<'0'||c>'9';c=getchar())
if(c=='-') a=-1;
for(;c>='0'&&c<='9';c=getchar())
x=(x<<1)+(x<<3)+c-'0';
return x*a;
}

namespace fast_pow { //快速幂;
int n,k;
void read() {
n=input();
k=input();
return;
}
int poow(int n,int k) {
int ret=1;
for(;k;k>>=1,n*=n)
if(k&1) ret*=n;
return ret;
}
int poowmod(int n,int k,int p) {
int ret=1;
for(;k;k>>=1,n=(n*n)%p)
if(k&1) ret=(ret*n)%p;
return ret;
}
void solve() {
read();
printf("%d",poow(n,k));
return;
}
}

namespace matrix { //矩阵乘法 + 矩阵快速幂;
const int maxn=110;
int n,k;
struct Matrix {
int a[maxn][maxn];
void clear() {
memset(a,0,sizeof(a));
return;
}
Matrix() {
this->clear();
}
Matrix operator * (const Matrix &c)const {
Matrix Ans;
for(int k=0;k<n;k++)
for(int i=0;i<n;i++)
for(int j=0;j<n;j++)
Ans.a[i][j]+=a[i][k]*c.a[k][j];
return Ans;
}
void Print() {
for(int i=0;i<n;i++) {
for(int j=0;j<n;j++)
printf("%d ",a[i][j]);
printf("\n");
}
return;
}
}A;
Matrix poow(Matrix A,int k) {
Matrix ret=A,t=A;
for(k--;k;k>>=1,t=t*t)
if(k&1) ret=ret*t;
return ret;
}
void read() {
n=input();k=input();
for(int i=0;i<n;i++)
for(int j=0;j<n;j++)
A.a[i][j]=input();
return;
}
void solve() {
read();
A=poow(A,k);
A.Print();
return;
}
}

namespace Prime {
//只有 1 和它本身两个约数的数叫素数(质数);
const int maxn=1010;
int prime[maxn],n,tot;
bool vis[maxn];
bool is_prime(int x) {
int k=sqrt(x+0.5);
for(int i=2;i<=k;i++)
if(x%i==0) return false;
return true;
}
void read() {
n=input();
return;
}
void Eratosthenes() { //Eratosthenes筛法
int m=sqrt(n+0.5);
vis[1]=true;
for(int i=2;i<=m;i++)
if(vis[i]==false)
for(int j=i*i;j<=n;j+=i)
vis[j]=true;
return;
}
void GetPrime() { //欧拉筛
for(int i=2;i<=n;i++) {
if(vis[i]==false) prime[++tot]=i;
for(int j=1;j<=tot&&i*prime[j]<=n;j++) {
vis[i*prime[j]]=true;
if(i%prime[j]==0) break;
}
}
return;
}
void solve() {
read();
GetPrime();
for(int i=1;i<=tot;i++)
printf("%d ",prime[i]);
return;
}
}

// a ≡b(mod p) 的意思是 a mod p == b mod p ,a-b=kp(k∈ Z);
namespace __inv {
/*

if a*b ≡1(mod p),就称b是a在mod p意义下的逆元,a是b在mod p意义下的逆元;

a/b mod p=a*inv(b,p);

inv(a,p)：表示 a 在mod p意义下的逆元;

如果 a与p 互质，那么inv(a,p)有解，且inv(a,p)<p;

否则 inv(a,p)无解;

逆元一般有3种求法;

1.线性求逆元;

2.费马小定理/欧拉定理;

3.扩展欧几里得(exgcd);

*/
const int maxn=1010;
int inv[maxn],n,p;
void read() {
n=input();
p=input();
return;
}
//1.线性求逆元;
void Getinv(int n,int p) {  //用来求 1-n mod p 的所有逆元(p是素数);
inv[1]=1;
for(int i=2;i<=n;i++)
inv[i]=(p-p/i)*inv[p%i]%p;
return;
}
void solve() {
read();
Getinv(n,p);
for(int i=1;i<=n;i++)
printf("%d ",inv[i]);
return;
}
}

namespace __phi {
const int maxn=1010;
using namespace Prime;
int phi[maxn];
int getphi(int x) { //求phi[x];
int ret=1;
for(int i=1;prime[i]*prime[i]<=x;i++)
if(x%prime[i]==0) {
ret*=prime[i]-1;
x/=prime[i];
while(x%prime[i]==0) ret*=prime[i],x/=prime[i];
}
if(x>1) ret*=x-1;
return ret;
}
/*

在欧拉筛的基础上求phi[1...n]，利用了phi的性质;

如果一个数 x 是素数，则有：

1. phi[x]=x-1,因为如果x是素数，那么<=x的数中只有 x 不与 x 互质;

2.if(i%x==0) 则 phi[i*x]=phi[i]*x;

3.if(i%x!=0) phi[i*x]=phi[i]*phi[x];(phi是积性函数)

由1.得 phi[x]=x-1 即 phi[i*x]=phi[i]*(x-1);

*/
void GetPhi() {
phi[1]=1;
for(int i=2;i<=n;i++) {
if(vis[i]==false) prime[++tot]=i,phi[i]=i-1;
for(int j=1;j<=tot&&i*prime[j]<=n;j++) {
vis[i*prime[j]]=true;
if(i%prime[j]==0) {
phi[i*prime[j]]=phi[i]*prime[j];
break;
} else phi[i*prime[j]]=phi[i]*(prime[j]-1);
}
}
return;
}
void solve() {
read();
GetPhi();
return;
}
/*

欧拉定理：a^phi[p] ≡1(mod p),对于任意互质的 a、p 恒成立;

由上式可得 a*a^(phi[p]-1) ≡1(mod p) 即 inv(a,p)=a^(phi[p]-1);

*/
int inv(int a,int p) {
return fast_pow::poow(a,getphi(p)-1);
}
}

namespace Fermat {
/*

费马小定理(Fermat): a^(p-1) ≡1(mod p),对于任意互质的 a、p 恒成立;

费马小定理是欧拉定理中当 p 是素数时的推论，常用来降低题目难度;

由 Fermat 可得 a* a^(p-2) ≡1(mod p)  即 inv(a,p)=a^ (p-2)，但前提是 p 是素数;

*/
int a,p;
void read() {
a=input();
p=input();
return;
}
int inv(int a,int p) {
return fast_pow::poow(a,p-2);
}
void solve() {
read();
printf("%d",inv(a,p));
return;
}
}

namespace __gcd {
int n,m,a,b,x,y,p;
void read() {
n=input();
m=input();
return;
}
int gcd(int a,int b) {
return b==0?a:gcd(b,a%b);
}
/*

exgcd: 扩展欧几里得,用来求 ax+by=gcd(a,b);

exgcd的应用主要体现在 3 个方面：

1.求解不定方程;(ax+by=c)

2.求线性同余方程;(ax ≡b(mod p))

3.求逆元;(ax ≡1(mod p))

*/
int exgcd(int a,int b,int &x,int &y) { //return 的值是 gcd(a,b);
if(b==0) {
x=1;
y=0;
return a;
}
int ans=exgcd(b,a%b,x,y),t=x;
x=y;
y=t-a/b*y;
return ans;
}
/*

1.求解不定方程;(ax+by=c)

只有 c%gcd(a,b)==0 时，方程有整数解;

设 g=gcd(a,b),则 c=k*g;

首先用 exgcd 求出 ax+by=g 的解(x,y)

同乘 k 得到 a*x*k+b*x*k=k*g ，显然解为(x*k,y*k);

*/
bool ind(int a,int b,int c,int &x,int &y) {
int g=exgcd(a,b,x,y),k=c/g;
if(c%g!=0) return false;
x*=k;y*=k;
return true;
}
/*

2.求线性同余方程;(ax ≡b(mod p))

只有 b%gcd(a,p)==0 时，方程有解，且有 gcd(a,p) 个解，解的间隔为 p/gcd(a,p);

化简同余方程可得 ax-b=kp

移项得到 ax-kp=b 设 y=-k 得到 ax+py=b，显然和 1.中的式子形式一样;

*/
bool modline(int a,int b,int p) {
int g=exgcd(a,p,x,y),k=p/g;
if(b%g!=0) return false;
int x0=x*(b/g)%p;
printf("%d ",x0);
for(int i=1;i<k;i++)
printf("%d ",(x0+k*i)%p);
return true;
}
/*

3.求逆元;(ax ≡1(mod p))

把 2.中的 b变成 1 就ok;

*/
int inv(int a,int p) {
exgcd(a,p,x,y);
return ((x%=p)+p)%p;
}
void solve() {
read();
printf("%d",gcd(n,m));
return;
}
}

namespace Mobius {
/*

Mobius: 莫比乌斯函数，一般用μ表示;

莫比乌斯函数的定义：

1. μ(1) = 1;

2.当 n 存在平方因子时，μ(n) = 0;

3.当 n 是素数或奇数个不同素数之积时，μ(n) = -1;

4.当 n 是偶数个不同素数之积时，μ(n) = 1;

*/
using namespace Prime;
const int maxn=1010;
int mu[maxn];
void Getmu() { //在欧拉筛的基础上筛 mu;
mu[1]=1;
for(int i=2;i<=n;i++) {
if(vis[i]==false) prime[++tot]=i,mu[i]=-1;
for(int j=1;j<=tot&&i*prime[j]<=n;j++) {
vis[i*prime[j]]=true;
if(i%prime[j]==0) {
mu[i*prime[j]]=0;
break;
}
mu[i*prime[j]]=-mu[i];
}
}
return;
}
void solve() {
read();
Getmu();
return;
}
}

namespace Baby_Step_Gaint_Step {
/*

BSGS算法: 求解形如 x^y ≡z(mod p)得式子中 y 的取值的算法;

先设 y = k*m+i，其中 m 是常量，一般取 ceil(sqrt(p+0.5));

代入得到 x^(k*m+i) ≡z(mod p)

x^i ≡z*x^(-k*m) (mod p)

BaByStep: 预处理出 x^i

GaintStep: 枚举 k ，查找是否存在 z*x^(-k*m) 满足条件;

*/
std::map<int,int>mp;
int x,z,p;
void read() {
x=input();
z=input();
p=input();
mp.clear();
return;
}
int bsgs(int x,int z,int p) {
x%=p;
if(x==0) return z==0?1:-1;
int m=ceil(sqrt(p+0.5)),now=1;
mp[1]=m+1;
for(int i=1;i<m;i++) {
now=(now*x)%p;
if(!mp[now]) mp[now]=i;
}
int inv=1,t=fast_pow::poowmod(x,p-m-1,p);
for(int k=0;k<m;k++,inv=(inv*t)%p) {
int i=mp[(z*inv)%p];
if(i) return k*m+(i==m+1?0:i);
}
return -1;
}
void solve() {
read();
printf("%d",bsgs(x,z,p));
return;
}
}

int main() {
/*

一些还不会的数论：

1.狄利克雷卷积

2.莫比乌斯反演

3.扩展 BSGS

4.Lucas 定理，扩展 Lucas

5.Lagrange 插值

6.CRT

*/
printf("The Number Theory is over\n");
return 0;
}