ACM常用模板-数论
常用数论模板,代码大多数来自网络,个人整理总结,不定期更新
基本函数及其性质:
最大公约数与最小公倍数:
gcd: //最大公约数
int gcd(int a,int b)
{
if(b==0) return a;
else return gcd(b,a%b);
}
lcm: //最小公倍数
int lcm(int a,int b)
{
return a/gcd(a,b)*b; //防止溢出
}
快速幂:
a,n,p都在long long 以内
LL pow(LL a, LL n, LL p) //快速幂 a^n % p
{
LL ans = 1;
while(n)
{
if(n & 1) ans = ans * a % p; //若不取模就去掉p
a = a * a % p;
n >>= 1;
}
return ans;
}
n超过long long
typedef unsigned long long LL;
string str;
int input[10000005];
int output[10000005];
int len;
int sum=1,d=0,k=0;
void to_bin(string str) //将大整数转换为二进制,转换后为逆序
{
len=str.size();
for(int i=0;i<len;i++)
input[i]=str[i]-'0';
memset(output,0,sizeof(output));
sum=1,d=0,k=0;
while(sum)
{
sum = 0;
for(int i=0;i<len;i++){
d = input[i] / 2;
sum += d;
if(i == (len - 1)){
output[k++] = input[i] % 2;
}
else
input[i+1] += (input[i]%2)*10;
input[i] = d;
}
}
}
LL pow(LL a, int n[], LL p) //快速幂 a^n % p
{
LL ans = 1;
int i=0;
while(i<k)
{
if(n[i] == 1) ans = ans * a % p;
a = a * a % p;
i++;
}
return ans;
}
int main()
{
LL a,b,c;
while(cin>>a>>str>>c)
{
to_bin(str);
cout<<pow(a,output,c)<<endl;
}
}
素数筛:
埃氏素数筛:
vector<int>prime; //其中存储2-n的所有素数
bool is_prime[MAXN] //存储第i项是否为素数
void find_prime(int n)
{
int p=0;
for(int i=2;i<=n;i++)
is_prime[i]=true;
is_prime[0]=is_prime[1]=false;
for(int i=2;i<=n;i++)
{
if(is_prime[i])
{
prime.push_back(i);
for(int j=2*i;j<=n;j+=i)
is_prime[j]=false;
}
}
}
线性素数筛:
#include <cstdio>
#include <bitset>
using namespace std;
bitset<10000010> g;
int n,m,x,p[1000010],top=0;
int main()
{
scanf("%d%d",&n,&m);
g[0]=1; g[1]=1;
for(int i=2;i<=n;++i)
{
if(g[i]==0)
p[++top]=i;
for(int j=1;j<=top&&p[j]*i<=n;++j)
{
g[i*p[j]]=1;
if(i%p[j]==0)
break;
}
}
for(int i=1;i<=m;++i)
{
scanf("%d",&x);
if(g[x])
printf("No\n");
else
printf("Yes\n");
}
return 0;
}
素性测试:
miller_rabin:
概率算法,能处理10^19范围数
#include<stdio.h>
#include<string.h>
#include<stdlib.h>
#include<time.h>
#include<iostream>
#include<algorithm>
using namespace std;
//****************************************************************
// 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;
}
// Miller_Rabin()算法素数判定
//是素数返回true.(可能是伪素数,但概率极小)
//合数返回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 算法进行质因数分解
//************************************************
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;
return;
}
long long p=n;
while(p>=n)p=Pollard_rho(p,rand()%(n-1)+1);
findfac(p);
findfac(n/p);
}
int main()
{
//srand(time(NULL));//需要time.h头文件//POJ上G++不能加这句话
long long n;
while(scanf("%I64d",&n)!=EOF)
{
tol=0;
findfac(n);
for(int i=0;i<tol;i++)printf("%I64d ",factor[i]);
printf("\n");
if(Miller_Rabin(n))printf("Yes\n");
else printf("No\n");
}
return 0;
}
逆序数:
#include <iostream>
#include <string.h>
#include <algorithm>
using namespace std;
const int N=100050;
int n;
long long c[N]; //c[n]表示a[1~n]的和,a数组省略
struct node
{
int val,pos;
}a[100005];
int lowbit(int x) //求2^k
{
return x & -x;
}
long long getsum(int n) //区间查询,求a[1~n]的和
{
long long res = 0;
while(n>0)
{
res+=c[n];
n=n-lowbit(n);
}
return res;
}
int change(int x) //单点更新,将c[x]的值加1
{
while(x<=n)
{
c[x]++;
x+=lowbit(x);
}
}
bool cmp(node a,node b) //包含相同数
{
if(a.val!=b.val)
return a.val>b.val;
return a.pos>b.pos;
}
int main()
{
std::ios::sync_with_stdio(false);
while(cin>>n)
{
memset(c,0,sizeof(c));
for(int i=1;i<=n;i++)
{
cin>>a[i].val;
a[i].pos=i;
}
sort(a+1,a+n+1,cmp);
long long cnt=0;
for(int i=1;i<=n;i++)
{
change(a[i].pos); //修改最大数位置的值
cnt+=getsum(a[i].pos-1); //最大数位置之前的所有位置和,即区间求和,可知比当前数小的数有多少个
}
cout<<cnt<<endl;
}
return 0;
}
欧拉函数:
定义: 在数论中,对正整数n,欧拉函数是小于等于n的数中与n互质的数的数目。并且用符号φ(n)表示一个整数的欧拉函数。
性质:
1.对于质数p, φ(p) = p - 1。注意φ(1)=1
2. 欧拉定理:对于互质的正整数a和n,有aφ(n) ≡ 1 mod n。
3. 若m,n互质,φ(mn)=φ(m)φ(n)
4. 当n为奇数时,φ(2n)=φ(n)
5.当n大于2时,所有φ(n)都是偶数
6.当n大于6时,所有φ(n)都是合数
朴素方式O(sqrt(n)):
//直接求解欧拉函数
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(Nlog(m)):
void Init(){
euler[1]=1;
for(int i=2;i<Max;i++)
euler[i]=i;
for(int i=2;i<Max;i++)
if(euler[i]==i)
for(int j=i;j<Max;j+=i)
euler[j]=euler[j]/i*(i-1);//先进行除法是为了防止中间数据的溢出
}
指数降幂公式:
计算exponial(n):
#include <iostream>
using namespace std;
typedef long long ll;
ll n,m,ans;
ll euler(ll n){ //返回euler(n) ,计算欧拉函数值
ll res=n,a=n;
for(ll 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;
}
ll fast_mod(ll x,ll n,ll Max) //快速幂
{
ll res=1;
while(n>0)
{
if(n & 1)
res=(res*x)%Max;
x=(x*x)%Max;
n >>= 1;
}
return res;
}
ll func(ll n,ll m){ //循环求解
if(m==1) return 0;
if(n<=5){
ll ans=1;
for(int i=1;i<=n;i++){
ans=fast_mod(i,ans,m);
}
return ans;
}else{
ll phi=euler(m);
ll z=func(n-1,phi);
ans=fast_mod(n,phi+z,m);
}
return ans;
}
void solve(){ //计算exponial(n)
scanf("%lld%lld",&n,&m);
printf("%lld\n",func(n,m));
}
int main(){
solve();
return 0;
}
组合计数:
排列数:
若n和r都是整数,且0<=r<=n,有
不尽相异元素全排列:
组合数:
重复组合数:
二项式定理:
常见恒等式:
帕斯卡恒等式:
容斥原理:
组合数取模:
typedef long long LL;
const LL maxn(1000005), mod(1e9 + 7);
LL Jc[maxn];
void calJc() //求maxn以内的数的阶乘
{
Jc[0] = Jc[1] = 1;
for(LL i = 2; i < maxn; i++)
Jc[i] = Jc[i - 1] * i % mod;
}
/*
//拓展欧几里得算法求逆元
void exgcd(LL a, LL b, LL &x, LL &y) //拓展欧几里得算法
{
if(!b) x = 1, y = 0;
else
{
exgcd(b, a % b, y, x);
y -= x * (a / b);
}
}
LL niYuan(LL a, LL b) //求a对b取模的逆元
{
LL x, y;
exgcd(a, b, x, y);
return (x + b) % b;
}
*/
//费马小定理求逆元
LL pow(LL a, LL n, LL p) //快速幂 a^n % p
{
LL ans = 1;
while(n)
{
if(n & 1) ans = ans * a % p;
a = a * a % p;
n >>= 1;
}
return ans;
}
LL niYuan(LL a, LL b) //费马小定理求逆元
{
return pow(a, b - 2, b);
}
LL C(LL a, LL b) //计算C(a, b)
{
return Jc[a] * niYuan(Jc[b], mod) % mod
* niYuan(Jc[a - b], mod) % mod;
}
母函数:
http://www.wutianqi.com/?p=596
用于计算组合问题可能情况数
普通母函数:
普通母函数通常解决类似如下的问题: 给5张1元,4张2元,3张5元,要得到15元,有多少种组合? 某些时候会规定至少使用3张1元、1张2元、0张5元。 某些时候会规定有无数张1元、2元、5元。
const int MAX=10000;
const int INF=0x3f3f3f3f;
//为计算结果,b为中间结果。
// K对应具体问题中物品的种类数。
//v[i]表示该乘积表达式第i个因子的权重,对应于具体问题的每个物品的价值或者权重。
//n1[i]表示该乘积表达式第i个因子的起始系数,对应于具体问题中的每个物品的最少个数,即最少要取多少个。0
//n2[i]表示该乘积表达式第i个因子的终止系数,对应于具体问题中的每个物品的最多个数,即最多要取多少个。INF
//P为可能出现的最大指数(所求情况)
int a[MAX],b[MAX],P;
int k,v[MAX],n1[MAX],n2[MAX];
//初始化a
void cal(int n) //n为种类数
{
memset(a,0,sizeof(a));
a[0]=1;
for (int i=1;i<=n;i++)//循环每个因子
{
memset(b,0,sizeof(b));
for (int j=n1[i];j<=n2[i]&&j*v[i]<=P;j++)//循环每个因子的每一项
for (int k=0;k+j*v[i]<=P;k++)//循环a的每个项
b[k+j*v[i]]+=a[k];//把结果加到对应位
memcpy(a,b,sizeof(b));//b赋值给a
}
}
指数母函数:
用于求多重排列数
如以下问题: 假设有8个元素,其中a1重复3次,a2重复2次,a3重复3次。从中取r个排列,求其排列数。 对于上面的问题“假设有8个元素,其中a1重复3次,a2重复2次,a3重复3次。从中取r个排列,求其排列数。
double c1[N],c2[N];
LL fac[N];
void Fac() //求阶乘
{
fac[0]=1;
for(int i=1;i<=N;i++)
fac[i]=fac[i-1]*i;
}
int a[N]; //1~n每种的个数
void cal(int n,int m) //有n种,取m个
{
memset(c1,0,sizeof(c1));
memset(c2,0,sizeof(c2));
c1[0]=1.0/fac[0];
for(int i=1;i<=n;i++)
{
for(int j=0;j<=m;j++)
for(int k=0;k+j<=m && k<=a[i];k++)
c2[k+j]+=c1[j]/fac[k];
for(int j=0;j<=m;j++)
c1[j]=c2[j],c2[j]=0;
}
}
ans=c1[m]*fac[m]; //取m个时的多重排列数
Polya计数:
定义:设G是有限集X上的置换群,点a,b∈X称为”等价”的,当且仅当,存在π∈G使得π(a)=b,记为a~b,这种等价条件下,X的元素形成的等价类称为G的轨道,它是集X的一个子集,G的任意两个不同的轨道之交是空集,所以置换群G的轨道全体是集合X的一个划分,构成若干个等价类,等价类的个数记为L。
Zk (K不动置换类):设G是1…n的置换群。若K是1…n中某个元素,G中使K保持不变的置换的全体,记以Zk,叫做G中使K保持不动的置换类,简称K不动置换类。
Ek(等价类):设G是1…n的置换群。若K是1…n中某个元素,K在G作用下的轨迹,记作Ek。即K在G的作用下所能变化成的所有元素的集合。.
这个时候有:|Ek|*|Zk|=|G|成立(k=1,2,…..n)。
C(π):对于一个置换π∈G,及a∈X,若π(a)=a,则称a为π的不动点。π的不动点的全体记为C(π)。例如π=(123)(3)(45)(6)(7),X={1,2,3,4,5,6,7};那么C(π)={3,6,7}共3个元素。
Burnside引理:L=1/|G|(Z1+Z2+Z3+Z4+……Zk)=1/|G|(C(π1)+C(π2)+C(π3)+…..+C(πn))(其中k∈X,π∈G)。
Polya定理:设G={π1,π2,π3……..πn}是X={a1,a2,a3…….an}上一个置换群,用m中颜色对X中的元素进行涂色,那么不同的涂色方案数为:1/|G|*(m^C(π1)+m^C(π2)+m^C(π3)+…+m^C(πk)). 其中C(πk)为置换πk的循环节的个数。
#define N 100000
int prime[N];
bool is_prime[N];
int sieve(int n)
{
int p = 0;
for (int i = 0; i <= n; ++i) is_prime[i] = true;
is_prime[0] = is_prime[1] = false;
for (int i = 2; i <= n; ++i) {
if (is_prime[i]) {
prime[p++] = i;
for (int j = 2 * i; j <= n; j += i)
is_prime[j] = false;
}
}
return p;
}
int phi(int n)
{
int rea = n;
for(int i = 0; prime[i] * prime[i] <= n; i++)
{
if(n % prime[i] == 0)
{
rea = rea - rea / prime[i];
while (n % prime[i] == 0) n /= prime[i];
}
}
if(n > 1)
rea = rea - rea / n;
return rea;
}
ll polya(int m, int n)
{
ll sum = 0;
ll i;
for (i = 1; i * i < n; ++i) {
if (n % i == 0) {
sum += phi(i) * pow(m, n / i);
sum += phi(n / i) * pow(m, i);
}
}
if (i * i == n) sum += phi(i) * pow(m, i);
if (n & 1) sum += n * pow(m, n / 2 + 1);
else sum += (pow(m, n / 2) + pow(m, n / 2 + 1)) * n / 2;
return sum / 2 / n;
}
快速傅里叶变换(FFT):
https://blog.csdn.net/ggn_2015/article/details/68922404
用于计算多项式乘法O(nlogn)
#include<cstdio>
#include<cstdlib>
#include<cmath>
#include<algorithm>
#include<cstring>
#include<complex>
#include <iostream>
using namespace std;
typedef complex<double> cd;//复数类的定义
const int maxl=2094153;//nlogn的最大长度(来自leo学长的博客)
const double PI=3.14159265358979;//圆周率,不解释
cd a[maxl],b[maxl];//用于储存变换的中间结果
int rev[maxl];//用于储存二进制反转的结果
void getrev(int bit){
for(int i=0;i<(1<<bit);i++){//高位决定二进制数的大小
rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
}//能保证(x>>1)<x,满足递推性质
}
void fft(cd* a,int n,int dft){//变换主要过程
for(int i=0;i<n;i++){//按照二进制反转
if(i<rev[i])//保证只把前面的数和后面的数交换,(否则数组会被翻回来)
swap(a[i],a[rev[i]]);
}
for(int step=1;step<n;step<<=1){//枚举步长的一半
cd wn=exp(cd(0,dft*PI/step));//计算单位复根
for(int j=0;j<n;j+=step<<1){//对于每一块
cd wnk(1,0);//!!每一块都是一个独立序列,都是以零次方位为起始的
for(int k=j;k<j+step;k++){//蝴蝶操作处理这一块
cd x=a[k];
cd y=wnk*a[k+step];
a[k]=x+y;
a[k+step]=x-y;
wnk*=wn;//计算下一次的复根
}
}
}
if(dft==-1){//如果是反变换,则要将序列除以n
for(int i=0;i<n;i++)
a[i]/=n;
}
}
int output[maxl];
char s1[maxl],s2[maxl];
void getmuti() //计算高精度大数乘法,输入两个数a,b
{
scanf("%s%s",s1,s2);//读入两个数
int l1=strlen(s1),l2=strlen(s2);//就算"次数界"
int bit=1,s=2;//s表示分割之前整块的长度
for(bit=1;(1<<bit)<l1+l2-1;bit++){
s<<=1;//找到第一个二的整数次幂使得其可以容纳这两个数的乘积
}
for(int i=0;i<l1;i++){//第一个数装入a
a[i]=(double)(s1[l1-i-1]-'0');
}
for(int i=0;i<l2;i++){//第二个数装入b
b[i]=(double)(s2[l2-i-1]-'0');
}
getrev(bit);fft(a,s,1);fft(b,s,1);//dft
for(int i=0;i<s;i++)a[i]*=b[i];//对应相乘
fft(a,s,-1);//idft
for(int i=0;i<s;i++){//还原成十进制数
output[i]+=(int)(a[i].real()+0.5);//注意精度误差
output[i+1]+=output[i]/10;
output[i]%=10;
}
int i;
for(i=l1+l2;!output[i]&&i>=0;i--);//去掉前导零
if(i==-1)printf("0");//特判长度为0的情况
for(;i>=0;i--){//输出这个十进制数
printf("%d",output[i]);
}
putchar('\n');
}
void getpoly() //计算多项式乘法
{
int n,m;
scanf("%d%d",&n,&m); //输入量多项式最高项次数
for(int i=0;i<=n;i++) scanf("%lf",&a[i].real()); //读入第一个多项式的系数(a0+a1*x+a2*x^2+a3*x^3+.....+an*x^n)
for(int i=0;i<=m;i++) scanf("%lf",&b[i].real()); //读入第二个多项式的系数(b0+b1*x+b2*x^2+b3*x^3+.....+bn*x^n)
int bit=0,s=1;//s表示分割之前整块的长度
for(s=1;s<=n+m;s<<=1)
bit++;
getrev(bit);fft(a,s,1);fft(b,s,1);//dft
for(int i=0;i<s;i++)a[i]*=b[i];//对应相乘
fft(a,s,-1);//idft
for(int i=0;i<=n+m;i++)
printf("%d ",(int)(a[i].real()+0.5)); //表示乘完多项式各项的系数,(a0+a1*x+a2*x^3...)
}
int main(){
getmuti(); //10*10=100
getpoly(); //(1+2x)*(1+2x+x2)=1+4x+5x2+2x3
return 0;
}