4、数論

17251 ワード

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)
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)
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)
n*=w[i];
for (i=0;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 
#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;
}