ACM数论模板
w2w原创
1.欧几里得定理
int gcd(int a, int b){
return b==0?a:gcd(b, a%b);
}
2.扩展欧几里得(求ax+by = gcd(a,b)的特解)
void e_gcd(LL a, LL b, LL &d, LL &x, LL &y){
if(b==0){
x = 1; y = 0; d = a;
}
else{
e_gcd(b, a%b, d, y, x);
y-= x*(a/b);
}
}
3.中国剩余定理
同余方程组
x ≡ a1(mod m1)
x ≡ a2(mod m2)
... ...
x ≡ ak(mod mk)
方程组所有的解的集合就是:
x1 = N1*a1 + N2*a2 + ... + Nk*ak
其中 Ni mod mi = 1,Ni = ai * ti , 可用欧几里得扩展定理求 ti. 其中M = m1*m2*m3····*mn;
//互质版
#include <iostream>
using namespace std;
//参数可为负数的扩展欧几里德定理
void exOJLD(int a, int b, int &x, int &y){
//根据欧几里德定理
if(b == 0){//任意数与0的最大公约数为其本身。
x = 1;
y = 0;
}else{
int x1, y1;
exOJLD(b, a%b, x1, y1);
if(a*b < 0){//异号取反
x = - y1;
y = a/b*y1 - x1;
}else{//同号
x = y1;
y = x1 - a/b* y1;
}
}
}
//剩余定理
int calSYDL(int a[], int m[], int k){
int N[k];//这个可以删除
int mm = 1;//最小公倍数
int result = 0;
for(int i = 0; i < k; i++){
mm *= m[i];
}
for(int j = 0; j < k; j++){
int L, J;
exOJLD(mm/m[j], -m[j], L, J);
N[j] = m[j] * J + 1;//1
N[j] = mm/m[j] * L;//2 【注】1和2这两个值应该是相等的。
result += N[j]*a[j];
}
return (result % mm + mm) % mm;//落在(0, mm)之间,这么写是为了防止result初始为负数,本例中不可能为负可以直接 写成:return result%mm;即可。
}
int main(){
int a[3] = {2, 3, 2};
int m[3] = {3, 5, 7};
cout<<"结果:"<<calSYDL(a, m, 3)<<endl;
}
//不互质版
/**
中国剩余定理(不互质)
*/
#include <iostream>
#include <cstdio>
#include <cstring>
using namespace std;
typedef long long LL;
LL Mod;
LL gcd(LL a, LL b)
{
if(b==0)
return a;
return gcd(b,a%b);
}
LL Extend_Euclid(LL a, LL b, LL&x, LL& y)
{
if(b==0)
{
x=1,y=0;
return a;
}
LL d = Extend_Euclid(b,a%b,x,y);
LL t = x;
x = y;
y = t - a/b*y;
return d;
}
//a在模n乘法下的逆元,没有则返回-1
LL inv(LL a, LL n)
{
LL x,y;
LL t = Extend_Euclid(a,n,x,y);
if(t != 1)
return -1;
return (x%n+n)%n;
}
//将两个方程合并为一个
bool merge(LL a1, LL n1, LL a2, LL n2, LL& a3, LL&
n3)
{
LL d = gcd(n1,n2);
LL c = a2-a1;
if(c%d)
return
false;
c = (c%n2+n2)%n2;
c /= d;
n1 /= d;
n2 /= d;
c *= inv(n1,n2);
c %= n2;
c *= n1*d;
c += a1;
n3 = n1*n2*d;
a3 = (c%n3+n3)%n3;
return true;
}
//求模线性方程组x=ai(mod ni),ni可以不互质
LL China_Reminder2(int len, LL* a, LL* n)
{
LL a1=a[0],n1=n[0];
LL a2,n2;
for(int i = 1; i < len; i++)
{
LL aa,nn;
a2 =
a[i],n2=n[i];
if(!merge(a1,n1,a2,n2,aa,nn))
return -1;
a1 = aa;
n1 = nn;
}
Mod = n1;
return (a1%n1+n1)%n1;
}
LL a[1000],b[1000];
int main()
{
int i;
int k;
while(scanf("%d",&k)!=EOF)
{
for(i = 0; i
< k; i++)
scanf("%I64d %I64d",&a[i],&b[i]);
printf("%I64d\n",China_Reminder2(k,b,a));
}
return 0;
}
4.欧拉函数(求一个数前面的所有与这个数互质的数的个数)
Euler函数表达通式:euler(x)=x(1-1/p1)(1-1/p2)(1-1/p3)(1-1/p4)…(1-1/pn),其中p1,p2……pn为x的所有素因数,x是不为0的整数。euler(1)=1(唯一和1互质的数就是1本身)。
Euler函数有几个性质:
1.如果q,p互质,则Euler(p*q) = Euler(p)*Euler(q);
2.如果 a = p^k,则Euler(a) = p^k - p^k-1;
//直接求解欧拉函数
int euler(int n){ //返回euler(n)
int res=n,a=n;
for(int i=2;i*i<=a;i++){
if(a%i==0){
res=res/i*(i-1);//先进行除法是为了防止中间数据的溢出
while(a%i==0) a/=i;
}
}
if(a>1) res=res/a*(a-1);
return res;
}
//线性筛选欧拉函数O(n)用到了一下性质:
//(1) 若(N%a==0 && (N/a)%a==0) 则有:E(N)=E(N/a)*a;
//(2) 若(N%a==0 && (N/a)%a!=0) 则有:E(N)=E(N/a)*(a-1);
//注意:如果范围过大 可能不适宜开数组来做
int euler[maxN],
vis[maxN], prime[maxN/5], e[maxN], cnt = 0;
void make_euler(){
memset(vis, 0, sizeof(vis));
euler[1] = 1;
for(int i=2; i<maxN ; ++i){
if(vis[i] == 0){
prime[cnt++] = i;
euler[i] = i-1;
}
for(int j=0 ; j<cnt
&& i*prime[j] < maxN; ++j){
vis[i*prime[j]] = 1;
if(
i%prime[j] == 0){
euler[i*prime[j]] = euler[i] *prime[j];
break;
}
else euler[i*prime[j]] = euler[i] *(prime[j]-1);
}
}
}
5.求N以前N的约数个数
约数个数的性质,对于一个数N,N=p1^a1 + p2^a2 + ... + pn^an。其中p1 ,p2, p3... pn是N的质因数,a1 ,a2, a2,...an为相应的指数,则
div_num[N]=(p1+1)*(p2+1)*(p3+1)* ... *(pn+1);
结合这个算法的特点,在程序中如下运用:
对于div_num:
(1)如果i|prime[j] 那么 div_num[i*prime[j]]=div_sum[i]/(e[i]+1)*(e[i]+2)
//最小素因子次数加1
(2)否则
div_num[i*prime[j]]=div_num[i]*div_num[prime[j]]
//满足积性函数条件
对于e:
(1)如果i|pr[j] e[i*pr[j]]=e[i]+1; //最小素因子次数加1
(2)否则 e[i*pr[j]]=1;
//pr[j]为1次
#include<string.h>
#include<iostream>
#define M 100000
using namespace std;
int prime[M/3],e[M],div_num[M]; // e[i]表示第i个素数因子的个数
bool flag[M];
void get_prime()
{
int i,j,k;
memset(flag,false,sizeof(flag));
k=0;
for(i=2;i<M;i++){
if(!flag[i]){
prime[k++]=i;
e[i]=1;
div_num[i]=2; //素数的约数个数为2
}
for(j=0;j<k&&i*prime[j]<M;j++){
flag[i*prime[j]]=true;
if(i%prime[j]==0){
div_num[i*prime[j]]=div_num[i]/(e[i]+1)*(e[i]+2);
e[i*prime[j]]=e[i]+1;
break;
}
else{
div_num[i*prime[j]]=div_num[i]*div_num[prime[j]];
e[i*prime[j]]=1;
}
}
}
}
6.莫比乌斯函数
一个讲得比较清楚的PPT:http://wenku.baidu.com/link?url=UARIPTGHjN78vIzedWT2iwICudBIbsuZ5WMrYwJJjp2P5x7hUvtvSoVKiW7a92GiiF7aCJu1FYid2eB5iM9Wh-hW2Bfd1UfJgrstX7nZnrm
线性筛打表莫比乌斯函数:
int mob[maxN], vis[maxN], prime[maxN], cnt=0;
void make_mobius(){
mob[1] = 1;
memset(vis, 0, sizeof(vis));
for(int i = 2; i<maxN ; ++i){
if(!vis[i]){
mob[i] = -1;
prime[cnt++] = i;
}
for(int j= 0; j<cnt && i*prime[j] < maxN ; ++j){
vis[i*prime[j]] = 1;
if(i%prime[j] == 0){
mob[i*prime[j]] = 0;
break;
}
else mob[i*prime[j]] = -mob[i];
}
}
}
7.卢卡斯定理
//卢卡斯定理 C(m, n)%p
LL Lucas(LL m, LL n, LL p){
LL res = 1;
while(n && m){
LL n1 = n%p, m1 = m%p;
//费马小定理求逆元
res = res* fac[n1]* qsm(fac[m1]*fac[n1-m1]%p, p-2, p)%p;
//扩展欧几里得求逆元
//res = res* fac[n1]* reverse(fac[m1],p)* reverse(fac[n1-m1],p)%p;
n /= p;
m /= p;
}
return (res%p + p)%p;
}
8.卡特兰数
递推公式
9.伯努利数
伯努利数满足条件,且有
那么继续得到
void init_ber(){
ber[0] = 1;
for(int i = 1 ; i<maxn; ++i){
LL ans = 0;
for(int j = 0 ; j<i ; ++j)
ans = (ans + comb[i+1][j]*ber[j])%mod;
ans = -ans*inv(i+1, mod)%mod;
ber[i] = (ans%mod + mod)%mod;
}
}
10.乘法逆元
因为可能会很大,超过int范围,所以在快速幂时要二分乘法。
逆元打表 ( O(n)):
inv[1] = 1;
for(int i=2;i<N;i++)
{
if(i >= MOD) break;
inv[i] = (MOD - MOD / i) * inv[MOD % i]% MOD;
}
11. miller rabin算法 pollard rho算法(概率高效判断素数,求因子)
#include <stdlib.h>
#include <string.h>
#include <algorithm>
#include <iostream>
#include <math.h>
const int Times=10;
const int N=5500;
using namespace std;
typedef long long LL;
LL ct,cnt,c,x,y;
LL fac[N],num[N],arr[N];
LL gcd(LL a,LL b)
{
return b? gcd(b,a%b):a;
}
LL multi(LL a,LL b,LL m)
{
LL ans=0;
while(b)
{
if(b&1)
{
ans=(ans+a)%m;
b--;
}
b>>=1;
a=(a+a)%m;
}
return ans;
}
LL quick_mod(LL a,LL b,LL m)
{
LL ans=1;
a%=m;
while(b)
{
if(b&1)
{
ans=multi(ans,a,m);
b--;
}
b>>=1;
a=multi(a,a,m);
}
return ans;
}
bool Miller_Rabin(LL n)
{
if(n==2) return true;
if(n<2||!(n&1)) return false;
LL a,m=n-1,x,y;
int k=0;
while((m&1)==0)
{
k++;
m>>=1;
}
for(int i=0; i<Times; i++)
{
a=rand()%(n-1)+1;
x=quick_mod(a,m,n);
for(int j=0; j<k; j++)
{
y=multi(x,x,n);
if(y==1&&x!=1&&x!=n-1) return false;
x=y;
}
if(y!=1) return false;
}
return true;
}
LL Pollard_rho(LL n,LL c)
{
LL x,y,d,i=1,k=2;
y=x=rand()%(n-1)+1;
while(true)
{
i++;
x=(multi(x,x,n)+c)%n;
d=gcd((y-x+n)%n,n);
if(1<d&&d<n) return d;
if(y==x) return n;
if(i==k)
{
y=x;
k<<=1;
}
}
}
void find(LL n,int c)
{
if(n==1) return;
if(Miller_Rabin(n))
{
fac[ct++]=n;
return ;
}
LL p=n;
LL k=c;
while(p>=n) p=Pollard_rho(p,c--);
find(p,k);
find(n/p,k);
}
void dfs(LL dept, LL product=1)
{
if(dept==cnt)
{
arr[c++]=product;
return;
}
for(int i=0; i<=num[dept]; i++)
{
dfs(dept+1,product);
product*=fac[dept];
}
}
void Solve(LL n)
{
ct=0;
find(n,120);
sort(fac,fac+ct);
num[0]=1;
int k=1;
for(int i=1; i<ct; i++)
{
if(fac[i]==fac[i-1])
++num[k-1];
else
{
num[k]=1;
fac[k++]=fac[i];
}
}
cnt=k;
}
const int M=1000005;
bool prime[M];
int p[M];
int k1;
void isprime()
{
k1=0;
int i,j;
memset(prime,true,sizeof(prime));
for(i=2;i<M;i++)
{
if(prime[i])
{
p[k1++]=i;
for(j=i+i;j<M;j+=i)
{
prime[j]=false;
}
}
}
}
int main()
{
LL n,t,record,tmp;
isprime();
while(cin>>n)
{
LL ans=-1;
record=0;
while(true)
{
tmp=1;
if(Miller_Rabin(n))
{
Solve(n-1);
c=0;
dfs(0,1);
sort(arr,arr+c);
bool flag=0;
for(int i=0; i<c; i++)
{
if(quick_mod(10,arr[i],n)==1)
{
tmp=arr[i];
break;
}
if(i==c-2) flag=1;
}
if(flag)
{
if(ans<tmp) record=n;
cout<<record<<endl;
break;
}
}
else
{
bool flag=false;
LL tmp1=(LL)sqrt(n*1.0);
if(tmp1*tmp1==n&&Miller_Rabin(tmp1))
{
x=tmp1;
y=2;
flag=true;
}
else
{
LL cnt1=0,rea=n;
for(int i=0;i<k1;i++)
{
if(rea%p[i]==0)
{
x=p[i];
while(rea%p[i]==0)
{
rea/=p[i];
cnt1++;
}
break;
}
}
if(rea==1) flag=true;
y=cnt1;
}
if(flag)
{
Solve(x-1);
c=0;
dfs(0,1);
sort(arr,arr+c);
bool flag=0;
for(int i=0; i<c; i++)
{
if(quick_mod(10,arr[i],x)==1)
{
tmp=1;
for(int j=0; j<y-1; j++)
tmp*=x;
tmp*=(x-1);
break;
}
if(i==c-2) flag=1;
}
if(flag)
{
if(ans<tmp)
{
ans=tmp;
record=n;
}
}
}
}
n--;
}
}
return 0;
}
12.快速乘快速幂
LL qsm(LL a, LL n, LL m){
LL ans = 1;
while(n>0){
if(n&1){
ans = (a*ans)%m;
}
a = (a*a)%m;
n>>=1;
}
return ans;
}
long long q_mul( long long a, long long b, long long mod ) //快速计算 (a*b) % mod
{
long long ans = 0; // 初始化
while(b) //根据b的每一位看加不加当前a
{
if(b & 1) //如果当前位为1
{
b--; //也可不要,方便理解而已
ans =(ans+ a)%mod; //ans+=a
}
b /= 2; //b向前移位
a = (a + a) % mod; //更新a
}
return ans;
}
13.博弈
奇异局势: ak =[k(1+√5)/2],bk= ak + k (k=0,1,2,…,n 方括号表示取整函 数)。