[施工中]良心数论.

/*
	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;
}