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 方括号表示取整函 数)。

时间: 2024-10-07 05:23:01

ACM数论模板的相关文章

ACM数论-欧几里得与拓展欧几里得

ACM数论--欧几里得与拓展欧几里得 欧几里得算法: 欧几里德算法又称辗转相除法,用于计算两个整数a,b的最大公约数. 基本算法:设a=qb+r,其中a,b,q,r都是整数,则gcd(a,b)=gcd(b,r),即gcd(a,b)=gcd(b,a%b). int gcd(int a,int b) { return b ? gcd(b,a%b) : a; } 扩展欧几里德算法: 基本算法:对于不完全为 0 的非负整数 a,b,gcd(a,b)表示 a,b 的最大公约数,必然存在整数对 x,y ,使

【史前巨坑】数论模板整合

跪了一下午数论 整理了一下数论模板 这是个史前巨坑,有空慢慢填 #include<iostream> #include<cstdio> #include<cstring> #include<cmath> #include<algorithm> #define MAXN 1000000 using namespace std; bool not_prime[MAXN]; int prime_number[MAXN]; int nu; int fac

ACM数论中相关定理(不断更新)

费马小定理是数论中的一个重要定理,其内容为: 假如p是质数,且(a,p)=1,那么 a^(p-1) ≡1(mod p).即:假如a是整数,p是质数,且a,p互质,那么a的(p-1)次方除以p的余数恒等于1. 费马大定理,又被称为“费马最后的定理”,由法国数学家费马提出.它断言当整数n >2时,关于x, y, z的方程 x^n + y^n = z^n 没有正整数解.被提出后,经历多人猜想辩证,历经三百多年的历史,最终在1995年被英国数学家安德鲁·怀尔斯证明. 中国剩余定理的结论: 令任意固定整数

ACM函数模板开源

今天平安夜,首先祝大家平平安安. 和众多Android工程师一样,我想写个好的app. 开发这个小应用,我也只是当时的一时兴起,为了提高我的开发速度,我用到了人家大牛写的仿小米启动页界面. Android手机千千万万,要做到全部兼容是很难的,同学老是叫我更新,我都有点烦了,首先我做这个我没有一点收入可言,我没有在里面放任何广告.下载量高了,我也不会有一点钱,提交到市场那只是方便下载罢了.再者,我是二流学校.搞ACM的算法那肯定是远远不及人家那种一流学校的. 还有就是,我比较懒.这就是重点. 我也

【ACM - 搜索模板】

[广搜模板] #include <iostream> #include <stdio.h> #include <string.h> #include <queue> using namespace std; #define MAXX #define MAXY struct Node{ int x,y; int step; }; int n,m; //边界 int dx[4] = {0,1,0,-1}; int dy[4] = {1,0,-1,0}; int

ACM数论中的常见的模板和结论

1:最大公约数的求法 欧几里得算法实现.递归实现 1 #include<stdio.h> 2 #include<string.h> 3 #include<algorithm> 4 #include<iostream> 5 using namespace std; 6 __int64 gcd(__int64 y,__int64 x) 7 { 8 __int64 ans=0; 9 if(x==0) 10 ans=y; 11 else 12 ans=gcd(x,y

ACM数论之旅10---大组合数-卢卡斯定理(在下卢卡斯,你是我的Master吗?(。-`ω&#180;-) )

记得前几章的组合数吧 我们学了O(n^2)的做法,加上逆元,我们又会了O(n)的做法 现在来了新问题,如果n和m很大呢, 比如求C(n, m) % p  , n<=1e18,m<=1e18,p<=1e5 看到没有,n和m这么大,但是p却很小,我们要利用这个p (数论就是这么无聊的东西,我要是让n=1e100,m=1e100,p=1e100你有本事给我算啊(°□°),还不是一样算不出来) 然后,我们著名的卢卡斯(Lucas)在人群中站了出来(`?д?´)说:“让老子来教你这题” 卢卡斯说:

ACM数论之旅13---母函数(又叫生成函数)(痛并快乐着(╭ ̄3 ̄)╭)

(前排出售零食瓜子) 前言:(????不看会吃亏的) 接下来要讲的算法比较难,一般人听不懂,因为你不能用常人的思想去理解 就像高数(说多了都是泪( >﹏<.)) 你要用常规方法去想肯定很累( ̄▽ ̄)~* 有时候一直半解的多读几遍反而高效 就像矩阵,发明矩阵的人是天才,我们只能在使用矩阵的同时一直感叹:“哇!好神奇!这样也可以!” 而不是先研究为什么他会有这么多神奇的性质 因为,毕竟从发明者角度来说,他就是考虑了这么多才创造出了神奇的矩阵 而如果矩阵的每一个性质你都一点点研究过去,终有一天,你也

PWJ的数论模板整理

一些还没学到,但已经听说的就先copy其他博客的 数论 欧拉降幂 求a1^a2^a3^a4^a5^a6 mod m #include<cstdio> #include<cstring> const int N=1e4+11; typedef long long ll; char s[10]; int n,lens,phi[N]; ll md,a[15]; void init(){ for(int i=1;i<N;i++) phi[i]=i; for(int i=2;i<