int qmod(int a,int b){ int res = 1; while(b){ if(b&1) res=res*a%MOD; b>>=1; a=a*a%MOD; } return res; } int inv(int a){ return qmod(a,MOD-2); }
求解高次同余方程
还是用 babystep_gaintstep算法求解。但是这题并不能用POJ_2417的算法，直接套该 算法，下面简要说明一下不能用的原因。首先我们有必要归纳一下用babystep算法解题 的步骤： (1) 求M = ceil( sqrt(C) ) ； (2) for(i=0;i< M;i++) hash( i , A^i ) ； (3) 求D = A^M%C； (4) r = 1 ; for( i = 0 ; i < M ; i++ ) ex_gcd(r , C , x , y ) ; res = x * B % C ; jj = find( res) 如果找到了这时候的jj，则答案就是i*M+jj，如果没有找到，则res = res * D % C，继续循 环查找，如果最终都没有找到，则输出无解。 在上述的步骤中，如果题目中没有告诉我们 gcd(A , C) = 1，则我们上述的方法是错误的，原因就在于第4步，求res的时候。因为如果 我们无法保证gcd(A , C) = 1 ，也就不能保证gcd(r ,C) = 1（因为D=A^M, r = D^i），所有在 用 扩展欧几里得求出r*x + C*y = gcd(r,C ) 的一个解x0之后，原方程:r*x+C*y = B的解 x = x0 * B / gcd(r,C) + i*C / gcd(r,C) ，但是我们这个时候并不能计算出gcd(r,C)，因为此时 的r本来就是经过取余之后得出的，并不能直接用来求gcd，因此我们上述的普通babystep 算法就会出错了。 这样我们就要换一种处理的方法了，这里介绍一种AC大牛博客上的一种“消因子”的方法， 具体内容请看这里：AC大牛。经过上面的分析我们很清楚接下去的处理应该从哪方面着 手，就是应该从不能求出gcd(r , C)入手。一种思想就是既然无法求， 那我每次只要保证 gcd(r, C) = 1那样就可以想普通babystep一样求解了，既然要保证gcd(r,C) =1 ，而 r = (A^M)^i,因此归根到底还是要求gcd( A , C ) = 1。下面就是从AC大牛博客上参考的“消因子” 法了，每次我们 都消去A,C的一个因子，然后对B,C, D进行如下的处理：B/=tmp;C/=tmp ; D = D* A/tmp%C ，这样经过b轮的消因子之后，gcd(A,C) = 1， 接下去我们就可以用普通 的babystep求解出方程：A^x = B'( mod C' ) 的解 res1， 原方程的解就是 res = res1 + b。 下面给出这种方法正确的简要证明；一开始我们要求的方程是：A^x = B( mod C )，也就是 求一个最小的x，使得A^x + C*y = B，通过消因子， 我们不断在方程两遍消去gcd(A,C)，这 样方程就可以变成 D*A^x1 + C'*y1 = B'，很简单就可以证明上式中 x = x1 + b ; y = y1 的（只要 在方程的两边分别将消去的因子乘回去等式还是保持不变的）。这样我们的问题就转化为了 求x1和y1，即D*A^x1 = B'( mod C' )，此时gcd( A , C') = 1，这样我们就可以用普通的babystep 求出上述式子的解x1，同时也就求出了x，这样本题就解决了。 但是上述的方法还是有一个bug的，也就是说，我们用babystep求出的x1>=0，所以上述的 方法只能求出x >= b的解，这样我们自然就会想到如果有一个解x < b怎么办，上述方法就会出 先错误了，因此我们这里还需要改进。考虑b的最大值是多少，考虑每次我们消去的因子数都 最小也就是2，这样我们就可以得到b的最大值就是log(C)，这样我们只要保证每次log(C)之内的 解都特判一下， 就不会出现我们刚才的问题了， 所以我们要在进行上述处理之前进行一次for 循环 ，特判0 - log(C)直接的x是否能成为解，接下去再用上述的“消因子”算法。 最后不得不佩服发明这种算法的人的神奇，将O(C)复杂度的判断，分两级判断将复杂度降低 到O( sqrt(C) )，所以就是为什么叫" babystep_gaintstep "了， 哈哈。
babystep算法模板
# define CC(m ,what) memset(m , what , sizeof(m)) LL A, B ,C ; const int NN = 99991 ; bool has[NN] ; int idx[NN] , val[NN] ;
void insert_(int id , LL vv) { LL v = vv % NN ; while( has[v] && val[v]!=vv) { v++ ; if(v == NN) v-=NN ; } if( !has[v] ) { has[v] = 1; val[v] = vv ; idx[v] = id ; } } int findi(LL vv) { LL v = vv % NN ; while( has[v] && val[v]!=vv) { v++ ; if(v == NN) v-=NN ; } if( !has[v] ) return -1; return idx[v] ; } void ex_gcd(LL a , LL b , LL& x , LL& y) { if(b == 0) { x = 1 ; y = 0 ; return ; } ex_gcd(b , a%b , x, y) ; LL t = x ; x = y; y = t - a/b*y ; } LL gcd(LL a,LL b) { while( a%b != 0) { LL c = a ; a = b ; b = c % b ; } return b ; } LL baby_step(LL A, LL B , LL C) { LL ans = 1 ; for(LL i=0; i<=50; i++) { if(ans == B) return i ; ans = ans * A % C ; } LL tmp , d = 0 ; LL D = 1 % C ; while( (tmp=gcd(A,C)) != 1 ) { if(B % tmp) return -1 ; d++ ; B/=tmp ; C/=tmp ; D = D*A/tmp%C ; } CC( has , 0) ; CC( idx, -1) ; CC(val , -1) ; LL M = ceil( sqrt(C*1.0) ) ; LL rr = 1 ; for(int i=0; i<M; i++) { insert_(i, rr) ; rr = rr * A % C ; } LL x, y ; for(int i=0; i<M; i++) { ex_gcd(D, C , x, y) ; LL r = x * B % C; r = (r % C + C) % C ; int jj = findi( r ) ; if(jj != -1) { return LL(i)*M + LL(jj) + d ; } D = D * rr % C ; } return -1 ; }
素数测试
线性筛法打素数表
int prime[20000],kp=0; int Is_or[65536]; void Prime() { int n =65536; //2~n之间的素数 kp=0; memset(Is_or,1,sizeof(Is_or)); Is_or[0]=Is_or[1]=0; for(int i=2;i<n;i++) { if(Is_or[i]) prime[kp++]=i;
于是我们就得到了一个定理的直接应用，对于待验证的数p，我们不断取a∈［1，p-1]且a∈Z，验证a^(p-1) mod p是否等于1，不是则p果断不是素数，共取s次。其中a^(p-1) mod p可以通过把p-1写成二进制，由(a*b)mod c=(a mod c)*b mod c，可以在t=log(p-1)的时间内计算出解，如考虑整数相乘的复杂度，则一次计算的总复杂度为log³(p-1)。这个方法叫快速幂取模。
为了提高算法的准确性，我们又有一个可以利用的定理。 定理二：对于0 < x < p，x^2 mod p =1 => x=1或p-1。
我们令p-1=(2^t)*u，即p-1为u二进制表示后面跟t个0。我们先计算出x[0]=a^u mod p ，再平方t次并在每一次模p，每一次的结果记为x[i]，最后也可以计算出a^(p-1) mod p。若发现x[i]=1而x[i-1]不等于1也不等于p-1，则发现p果断不是素数。
可以证明，使用以上两个定理以后，检验s次出错的概率至多为2^(-s)，所以这个算法是很可靠的。
需要注意的是，为了防止溢出（特别大的数据），a*b mod c 也应用类似快速幂取模的方法计算。当然，数据不是很大就可以免了。
下面是我的程序。
//**************************************************************** // Miller_Rabin 算法进行素数测试 //速度快，而且可以判断 <2^63的数 //**************************************************************** const int S=20;//随机算法判定次数，S越大，判错概率越小 //计算 (a*b)%c. a,b都是long long的数，直接相乘可能溢出的 // a,b,c <2^63 long long mult_mod(long long a,long long b,long long c) { a%=c; b%=c; long long ret=0; while(b) { if(b&1){ret+=a;ret%=c;} a<<=1; if(a>=c)a%=c; b>>=1; } return ret; } //计算 x^n %c long long pow_mod(long long x,long long n,long long mod)//x^n%c { if(n==1)return x%mod; x%=mod; long long tmp=x; long long ret=1; while(n) { if(n&1) ret=mult_mod(ret,tmp,mod); tmp=mult_mod(tmp,tmp,mod); n>>=1; } return ret; }
//以a为基,n-1=x*2^t a^(n-1)=1(mod n) 验证n是不是合数 //一定是合数返回true,不一定返回false bool check(long long a,long long n,long long x,long long t) { long long ret=pow_mod(a,x,n); long long last=ret; for(int i=1;i<=t;i++) { ret=mult_mod(ret,ret,n); if(ret==1&&last!=1&&last!=n-1) return true;//合数 last=ret; } if(ret!=1) return true; return false; }
bool Miller_Rabin(long long n) { if(n<2)return false; if(n==2)return true; if((n&1)==0) return false;//偶数 long long x=n-1; long long t=0; while((x&1)==0){x>>=1;t++;} for(int i=0;i<S;i++) { long long a=rand()%(n-1)+1;//rand()需要stdlib.h头文件 if(check(a,n,x,t)) return false;//合数 } return true; }
pollard_rho(longlong质因子分解)
//************************************************ //pollard_rho 算法进行质因数分解 //************************************************ long long factor[100];//质因数分解结果（刚返回时是无序的） int tol;//质因数的个数。数组小标从0开始
long long gcd(long long a,long long b) { if(a==0)return 1;//??????? if(a<0) return gcd(-a,b); while(b) { long long t=a%b; a=b; b=t; } return a; }
long long Pollard_rho(long long x,long long c) { long long i=1,k=2; long long x0=rand()%x; long long y=x0; while(1) { i++; x0=(mult_mod(x0,x0,x)+c)%x; long long d=gcd(y-x0,x); if(d!=1&&d!=x) return d; if(y==x0) return x; if(i==k){y=x0;k+=k;} } } //对n进行素因子分解 void findfac(long long n) { if(Miller_Rabin(n))//素数 { factor[tol++]=n; //值得注意的是 这里的factor并不是有序的!!!!! return; } long long p=n; while(p>=n)p=Pollard_rho(p,rand()%(n-1)+1); findfac(p); findfac(n/p); }
const float EPS = 1e-5; int sqrt(double x) { if(x == 0) return 0; double result = x; /*Use double to avoid possible overflow*/ double lastValue; do{ lastValue = result; result = result / 2.0f + x / 2.0f / result; }while(abs(result - lastValue) > EPS); return (double)result; }
更牛逼的一种开跟方式 快的一笔
/** 来自雷神之锤III的源代码中q_math.c的文件中。 */ float Q_rsqrt( float number ) { long i; float x2, y; const float threehalfs = 1.5F; x2 = number * 0.5F; y = number; i = * ( long * ) &y; // evil floating point bit level hacking i = 0x5f3759df - ( i >> 1 ); // what the fuck? y = * ( float * ) &i; y = y * ( threehalfs - ( x2 * y * y ) ); // 1st iteration // y = y * ( threehalfs - ( x2 * y * y ) ); // 2nd iteration, this can be removed #ifndef Q3_VM # ifdef __linux__ assert( !isnan(y) ); // bk010122 - FPE? #endif #endif return y; } //整理得到 int sqrt(float x) { if(x == 0) return 0; float result = x; float xhalf = 0.5f*result; int i = *(int*)&result; i = 0x5f375a86- (i>>1); // what the fuck? result = *(float*)&i; result = result*(1.5f-xhalf*result*result); // Newton step, repeating increases accuracy result = result*(1.5f-xhalf*result*result); return 1.0f/result; }
$O(n)$预处理逆元
方法一 i的逆元 inv[1] = 1; for (int i = 2; i<MAXN; i++) inv[i] = inv[MOD%i]*(MOD-MOD/i)%MOD; 方法二 inv{(n-i)!} = inv(n!)*n //阶乘逆元 Fac[0] = 1; for (int i = 1; i < N; i++) Fac[i] = (Fac[i-1] * i) % MOD; Inv[N-1] = pow_mod(Fac[N-1], MOD-2);//Fac[N]^{MOD-2} for (int i = N - 2; i >= 0; i--) Inv[i] = Inv[i+1] * (i + 1) % MOD; 方法三 费马小定理 fac[0]=1; for(int i=1;i<N;i++)fac[i]=fac[i-1]*i%MOD; for(int i=1;i<N;i++)inv[i]=qmod(fac[i],MOD-2);
