4 、 数论

4.1 阶乘最后非 0 位

//求阶乘最后非零位,复杂度 O(nlogn)
//返回该位,n 以字符串方式传入
#include <string.h>
#define MAXN 10000
int lastdigit(char* buf){
const int mod[20]={1,1,2,6,4,2,2,4,2,8,4,4,8,4,6,8,8,6,8,2};
int len=strlen(buf),a[MAXN],i,c,ret=1;
if (len==1)
return mod[buf[0]-'0'];
for (i=0;i<len;i++)
a[i]=buf[len-1-i]-'0';
for (;len;len-=!a[len-1]){
ret=ret*mod[a[1]%2*10+a[0]]%5;
for (c=0,i=len-1;i>=0;i--)
67
c=c*10+a[i],a[i]=c/5,c%=5;
}
return ret+ret%2*5;
}

4.2 模线性方程组

#ifdef WIN32
typedef __int64 i64;
#else
typedef long long i64;
#endif
//扩展 Euclid 求解 gcd(a,b)=ax+by
int ext_gcd(int a,int b,int& x,int& y){
int t,ret;
if (!b){
x=1,y=0;
return a;
}
ret=ext_gcd(b,a%b,x,y);
t=x,x=y,y=t-a/b*y;
return ret;
}
//计算 m^a, O(loga), 本身没什么用, 注意这个按位处理的方法 :-P
int exponent(int m,int a){
int ret=1;
for (;a;a>>=1,m*=m)
if (a&1)
ret*=m;
return ret;
}
//计算幂取模 a^b mod n, O(logb)
int modular_exponent(int a,int b,int n){ //a^b mod n
int ret=1;
for (;b;b>>=1,a=(int)((i64)a)*a%n)
if (b&1)
ret=(int)((i64)ret)*a%n;
return ret;
}
//求解模线性方程 ax=b (mod n)
//返回解的个数,解保存在 sol[]中
68
//要求 n>0,解的范围 0..n-1
int modular_linear(int a,int b,int n,int* sol){
int d,e,x,y,i;
d=ext_gcd(a,n,x,y);
if (b%d)
return 0;
e=(x*(b/d)%n+n)%n;
for (i=0;i<d;i++)
sol[i]=(e+i*(n/d))%n;
return d;
}
//求解模线性方程组(中国余数定理)
// x = b[0] (mod w[0])
// x = b[1] (mod w[1])
// ...
// x = b[k-1] (mod w[k-1])
//要求 w[i]>0,w[i]与 w[j]互质,解的范围 1..n,n=w[0]*w[1]*...*w[k-1]
int modular_linear_system(int b[],int w[],int k){
int d,x,y,a=0,m,n=1,i;
for (i=0;i<k;i++)
n*=w[i];
for (i=0;i<k;i++){
m=n/w[i];
d=ext_gcd(w[i],m,x,y);
a=(a+y*m*b[i])%n;
}
return (a+n)%n;
}

4.3 素数

//用素数表判定素数,先调用 initprime
int plist[10000],pcount=0;
int prime(int n){
int i;
if ((n!=2&&!(n%2))||(n!=3&&!(n%3))||(n!=5&&!(n%5))||(n!=7&&!(n%7)))
return 0;
for (i=0;plist[i]*plist[i]<=n;i++)
if (!(n%plist[i]))
return 0;
return n>1;
}
69
void initprime(){
int i;
for (plist[pcount++]=2,i=3;i<50000;i++)
if (prime(i))
plist[pcount++]=i;
}
//miller rabin
//判断自然数 n 是否为素数
//time 越高失败概率越低,一般取 10 到 50
#include <stdlib.h>
#ifdef WIN32
typedef __int64 i64;
#else
typedef long long i64;
#endif
int modular_exponent(int a,int b,int n){ //a^b mod n
int ret;
for (;b;b>>=1,a=(int)((i64)a)*a%n)
if (b&1)
ret=(int)((i64)ret)*a%n;
return ret;
}
// Carmicheal number: 561,41041,825265,321197185
int miller_rabin(int n,int time=10){
if (n==1||(n!=2&&!(n%2))||(n!=3&&!(n%3))||(n!=5&&!(n%5))||(n!=7&&!(n%7)))
return 0;
while (time--)
if
(modular_exponent(((rand()&0x7fff<<16)+rand()&0x7fff+rand()&0x7fff)%(n-1)+1,n-1,n)!=1)
return 0;
return 1;
}

4.4 欧拉函数

int gcd(int a,int b){
return b?gcd(b,a%b):a;
}
inline int lcm(int a,int b){
return a/gcd(a,b)*b;
}
//求 1..n-1 中与 n 互质的数的个数
int eular(int n){
int ret=1,i;
for (i=2;i*i<=n;i++)
if (n%i==0){
n/=i,ret*=i-1;
while (n%i==0)
n/=i,ret*=i;
}
if (n>1)
ret*=n-1;
return ret;
}