转载请注明出处:http://blog.csdn.net/crescent__moon/article/details/18653111
数论
高效素数打表:
void init_prim(int n)//prime存的是下标,visit存的是数。visit[5]==true。
{
memset(visit,true,sizeof(visit));
for(int i=2;i<=n;++i)
{
if(visit[i]==true)
{
num++;
prime[num]=i;
}
for(int j=1;((j<=num)&&(i*prime[j]<=n));++j)
{
visit[i*prime[j]]=false;
if(i%prime[j]==0) break;
}
}
}
算数基本定理:
任何一个大于1的自然数
任何一个大于1的自然数N,都可以唯一分解成有限个质数的乘积N=p1^a1*p2^a2...*pn^an,这里p1<p2<p3<..pn均为质数,其诸指数ai是正整数。
#include<iostream>
#include<cstdio>
#include<cmath>
#include<cstring>
using namespace std;
#define N 500010
int num; //素数的个数
int p[500010];
void fun() //素数打表
{
num = 0;
memset(p, 0, sizeof(p));
int res = sqrt(N * 1.0);
for(int i = 2; i < res; ++i)
{
if(!p[i])
{
for(int j = i * i; j < N; j += i)
p[j] = 1;
}
}
for(int i = 2; i < N; ++i) //存素数表
if(!p[i])
p[num++] = i;
}
int judge(int n) //求n的因子之和
{
int sum, ans, sumall = 1;
for(int i = 0; i < num; ++i)
{
sum = 0;
while(n % p[i] == 0 && n != 1)
{
sum++;
n /= p[i];
}
ans = 0;
for(int j = 0; j <= sum; ++j)
{
ans += pow(p[i] * 1.0, j);
}
sumall *= ans;
if(n == 1)
break;
}
return sumall;
}
int judge2(int n) //判断n的因子个数
{
int sum, sumall = 1;
for(int i = 0; i < num; ++i)
{
sum = 0;
while(n % p[i] == 0 && n != 1)
{
sum++;
n /= p[i];
}
sumall *= (sum + 1);
if(n == 1)
break;
}
return sumall; //因子个数
}
int main()
{
int ncase;
int n;
fun();
scanf("%d", &ncase);
while(ncase--)
{
scanf("%d", &n);
printf("%d的因子个数为:\n%d\n", n, judge2(n));
printf("%d的因子之和为:\n%d\n", n, judge(n));
}
return 0;
}
long long get(long long n) //S[i][0]代表第i个素数,S[i][1]代表第i个素数的个数
{
len=0;
for(long long i=2;i*i<=n;i++)
{
if(n%i==0)
{
s[len][0]=i;s[len][1]=0;
do{n/=i;s[len][1]++;}while(n%i==0);
len++;
}
}
if(n>1){s[len][0]=n;s[len++][1]=1;}
}
设正整数n有素因子分解n=p1^a1*p2^a2*...*ps^as,那么
n的所有正因子之和为S=[ p1 ^ ( a1 + 1 ) - 1 ] / (p1-1 ) ] * [ p2 ^ ( a2 + 1 ) - 1 ] / (p2-1 ) ] * ...* [ ps ^ ( as + 1 ) - 1 ] / (ps-1 ) ]
S也可以表示为S=( 1 + p1 + p1^2 + ... + p1 ^ a1 ) * ( 1 + p2 + p2^2 + ... + p2 ^ a2 ) * ... * ( 1 + ps + ps^2 + ... + ps ^ as )
n的所有正因子个数为T=(a1+1)*(a2+1)*...*(as+1)
同余定理:
定义1 用给定的正整数m分别除整数a、b,如果所得的余数相等,则称a、b对模m同余,记作a≡b(mod m)
定理1 整数a,b对模m同余的充要条件是 a-b能被m整除(即m|a-b)。
推论 a≡b(mod m)的充要条件是a=mt+b(t为整数)。
定理2 同余关系具有反身性、对称性与传递性,即
1)a≡a (mod m);
2)若a≡b (mod m), 则b≡a (mod m);
3)若a≡b (mod m), b≡c (mod m),则a≡c (mod m).
定理3 若a≡b(mod m), c≡d (mod m),则
1)a+c≡b+d (mod m);
2)a-c≡b-d (mod m);
3)ac≡bd (mod m).
推论 若a≡b(mod m),n为自然数,则an≡bn (mod m)。
定理4 若ca≡cb(mod m), (c,m)=d, 且a,b为整数,则a≡b(mod m/d).
推论若ca=cb(mod m), (c,m)=1,且a,b为整数,则a≡b(mod m).
定理5 若a≡b (mod m),a≡b (mod n),则a≡b(mod [m,n]).
推论若a≡b(mod mi), i=1,2,…,n,则a≡b (mod [m1,m2,..,mn]).
碾转相除法(欧几里得算法):
设a=q*b+r,其中a,b,q,r都是整数,则gcd(a,b)=gcd(b,r)
欧几里得算法实现:
int gcd(int a,int b)
{
if(!b)return a;
return gcd(b,a%b);
}
扩展欧几里得算法:
设a和b不全为0,则存在整数x和y,使得:gcd(a,b)=x*a+y*b
扩展欧几里得算法的应用主要在以下三方面:
1:求解不定方程;
2:求解模的逆元;
3:求解同余方程;
对于一个二元方程Ax+By=C (A,B已知),可以用扩展欧几里得算法求得其中一个解x(或者y),不过所求的解并不是最小的,也不是最大的,只不过是方程的一个解,令d=exgcd(A,B),该方程有解的充要条件为 d | C,即C% d==0。
方程Ax+By=C的最小解:x=x*(C/d)
方程Ax+By=C的最小整数解:x=(x%(B/d)+B/d)%(B/d)
如果是模线性方程ax=b (mod n),令d=exgcd(a,n),该方程有解的充要条件为 d | b ,即 b% d==0
方程ax=b(mod n)的最小解 :x=(x*(b/d))%n
方程ax=b(mod n)的最整数小解: x=(x%(n/d)+n/d)%(n/d)
解一元线性同余方程组:
假设同余方程组为:
x≡r1(mod a1)
x≡r2(mod a2)
...............
x≡rn(mod an)
下面给出此方程组的最小非负整数解代码:
long long solve()
{
bool flag=1;
scanf("%I64d%I64d",&a1,&r1);
for(long long i=1; i<n; i++)
{
scanf("%I64d%I64d",&a2,&r2);
a=a1;
b=a2;
c=r2-r1;
long long d=exgcd(a,x,b,y);
if(c%d!=0)
{
flag=0;
}
long long t=b/d;
x=(x*(c/d)%t+t)%t;
r1=a1*x+r1;
a1=a1*(a2/d);
}
if(!flag)return -1;//若返回-1,说明方程组无解
return r1;
}
用中国剩余定理来解:
x≡a1(mod m1)
x≡a2(mod m2)
.....................
x≡ak(mod mk)
int China(int r)
{
M=1;
int x,y,d,ans=0,Mi;
for(int i=1;i<=r;i++)
M*=m[i];
for(int i=1;i<=r;i++)
{
Mi=M/m[i];
int d=exgcd(Mi,x,m[i],y);
ans=(ans+Mi*x*a[i])%M;
}
if(ans<0)ans+=M;
return ans;
}
下面是两种扩展欧几里得算法的代码:
第一种是迭代的扩展欧几里得算法:
int exgcd(int A,int &x,int B,int &y)
{
int x1,y1,x0,y0;
x0=1;y0=0;
x1=0;y1=1;
int r=(A%B+B)%B;
int q=(A-r)/B;
x=0;y=1;
while(r)
{
x=x0-q*x1;
y=y0-q*y1;
x0=x1;
y0=y1;
x1=x;y1=y;
A=B;B=r;r=A%B;
q=(A-r)/B;
}
return B;
}
第二种是递归的欧几里得算法:
long long extended_euclidean(long long n, long long m, long long &x, long long &y) {
if (m == 0) {
x = 1; y = 0; return n;
}
long long g = extended_euclidean(m, n % m, x, y);
long long t = x - n / m * y;
x = y;
y = t;
return g;
}
POJ 1061 青蛙的约会
POJ 2115 C Looooops
POJ 2891 Strange Way to Express Integers
解高次同余方程:
解A^x≡B(mod C)的解的模板。
#define Maxn 65535
struct hash
{
int a,b,next;
}Hash[Maxn*2];
int flg[Maxn+66];
int top,idx;
void ins(int a,int b)
{
int k=b&Maxn;
if(flg[k]!=idx)
{
flg[k]=idx;
Hash[k].next=-1;
Hash[k].a=a;
Hash[k].b=b;
return;
}
while(Hash[k].next!=-1)
{
if(Hash[k].b==b)return;
k=Hash[k].next;
}
Hash[k].next=++top;
Hash[top].next=-1;
Hash[top].a=a;
Hash[top].b=b;
}
int find(int b)
{
int k=b&Maxn;
if(flg[k]!=idx)return -1;
while(k!=-1)
{
if(Hash[k].b==b)return Hash[k].a;
k=Hash[k].next;
}
return -1;
}
int gcd(int a,int b)
{
return b==0?a:gcd(b,a%b);
}
int ex_gcd(int a,int b,int& x,int& y)
{
int t,ret;
if(!b){x=1,y=0;return a;}
ret=ex_gcd(b,a%b,x,y);
t=x,x=y,y=t-a/b*y;
return ret;
}
int Inval(int a,int b,int n)
{
int x,y,e;
ex_gcd(a,n,x,y);
e=(long long)x*b%n;
return e<0?e+n:e;
}
int pow_mod(long long a,int b,int c)
{
long long ret=1%c;a%=c;
while(b)
{
if(b&1){ret=ret*a%c;}
a=a*a%c;
b>>=1;
}
return ret;
}
int BabyStep(int A,int B,int C)
{
top=Maxn;++idx;
long long buf=1%C,D=buf,K;
int i,d=0,tmp;
for(i=0;i<100;buf=buf*A%C,++i)
if(buf==B)return i;
while((tmp=gcd(A,C))!=1)
{
if(B%tmp)return -1;
++d;
C/=tmp;
B/=tmp;
D=D*A/tmp%C;
}
int M=(int)ceil(sqrt(C*1.0));
for(buf=1%C,i=0;i<=M;buf=buf*A%C,++i)ins(i,buf);
for(i=0,K=pow_mod((long long)A,M,C);i<=M;D=D*K%C,++i)
{
tmp=Inval((int)D,B,C);int w;
if(tmp>=0&&(w=find(tmp))!=-1)
return i*M+w+d;
}
return -1;
}
POJ 3243 Clever Y
POJ 2417 Discrete Logging
HDU 2815 Mod Tree
佩尔方程:
形如x^2-d*y^2=1(d>1且d不为完全平方数)的不定方程称为佩尔方程。
若佩尔方程x^2-d*y^2=1的最小特解(最小正整数解)是(x1,y1),那么可有迭代公式:
xn=xn-1*x1+d*yn-1*y1
yn=xn-1*y1+yn-1*x1
求出所有的正整数解(xk,yk)。
求解最小正整数解代码:
void search(int d)
{
int x,y;//x,y为最小正整数解
y=1;
while(1)
{
x=sqrt(d*y*y+1);
if(x*x-d*y*y==1)break;
y++;
}
}
毕达哥拉斯三元组:
毕达哥拉斯三元组:通俗点来说就是满足勾股定理的三个数x,y,z,即满足x^2+y^2=z^2的三个正整数三元组
定理:如果正整数x,y,z构成一个毕达哥拉斯三元组且y为偶数,当且仅当存在互素的正整数m,n(m>n),其中m为
偶数n为奇数,或者n为偶数m为奇数,并且满足
x=m^2-n^2
y=2mn
z=m^2+n^2
枚举m,n求解的代码:
for(int n=1; n<=tmp; n++)
for(int m=n+1; m<=tmp; m++)
{
if(n%2!=m%2)
{
if(gcd(m,n)==1)
{
x=m*m-n*n;
y=2*m*n;
z=m*m+n*n;
}
}
}
整数快速幂a^b%M:
long long multi(long long a,long long b,long long m)//a*b%m
{
long long ret=0;
while(b>0)
{
if(b&1)ret=(ret+a)%m;
b>>=1;
a=(a<<1)%m;
}
return ret;
}
long long quick_mod(long long a,long long b,long long m)//a^b%m
{
long long ans=1;
a%=m;
while(b)
{
if(b&1)
{
ans=multi(ans,a,m);
b--;
}
b>>=1;
a=multi(a,a,m);
}
return ans;
}
POJ 1995 Raising Modulo Numbers
矩阵快速幂A^k%M:
struct Matrax
{
int m[maxn][maxn];
}a,per;
void init()//建立矩阵
{
for(int i=0;i<n;i++)
for(int j=0;j<n;j++)
{
scanf("%d",&a.m[i][j]);
a.m[i][j]%=M;
per.m[i][j]=(i==j);
}
}
Matrax multi(Matrax a,Matrax b)//矩阵相乘
{
Matrax c;
for(int i=0;i<n;i++)
for(int j=0;j<n;j++)
{
c.m[i][j]=0;
for(int k=0;k<n;k++)
c.m[i][j]+=a.m[i][k]*b.m[k][j];
c.m[i][j]%=M;
}
return c;
}
Matrax power(int k)//矩阵快速幂
{
Matrax c,p,ans=per;
p=a;
while(k)
{
if(k&1){ans=multi(ans,p);k--;}
else {k/=2;p=multi(p,p);}
}
return ans;
}
POJ 3233 Matrix Power Series
Miller-Rabin素数测试:
#define Times 11
#define inf ((long long)1<<61)
#define C 201
long long jl[501],mini;//jl里面存的是大数的所有质因子
int ct;
long long gcd(long long a,long long b)//最大公约数
{
return b==0?a:gcd(b,a%b);
}
long long random(long long n)//生成随机数
{
return (long long)((double)rand()/RAND_MAX*n+0.5);
}
long long multi(long long a,long long b,long long m)//a*b%m
{
long long ret=0;
while(b>0)
{
if(b&1)ret=(ret+a)%m;
b>>=1;
a=(a<<1)%m;
}
return ret;
}
long long quick_mod(long long a,long long b,long long m)//a^b%m
{
long long ans=1;
a%=m;
while(b)
{
if(b&1)
{
ans=multi(ans,a,m);
b--;
}
b/=2;
a=multi(a,a,m);
}
return ans;
}
bool Witness(long long a,long long n)
{
long long m=n-1;
int j=0;
while(!(m&1))
{
j++;
m>>=1;
}
long long x=quick_mod(a,m,n);
if(x==1||x==n-1)return false;
while(j--)
{
x=x*x%n;
if(x==n-1)return false;
}
return true;
}
bool miller_rabin(long long n)//素数测试
{
if(n<2)return false;
if(n==2)return true;
if(!(n&1))return false;
for(int i=1; i<=Times; i++)
{
long long a=random(n-2)+1;
if(Witness(a,n))return false;
}
return true;
}
Pollard rho整数分解:
long long jl[501],mini=inf;//jl里面存的是大数的所有质因子,mini为最小的质因数
int ct=0;
long long pollard_rho(long long n,int c)//整数n分解,c一般为201
{
long long x,y,d,i=1,k=2;
x=random(n-1)+1;
y=x;
while(1)
{
i++;
x=(multi(x,x,n)+c)%n;
d=gcd(y-x,n);
if(1<d&&d<n)return d;
if(y==x)return n;
if(i==k)
{
y=x;
k<<=1;
}
}
}
void find(long long n,int k)//k=201
{
if(n==1)return;
if(miller_rabin(n))
{
jl[++ct]=n;
if(n<mini)mini=n;
return;
}
long long p=n;
while(p>=n)p=pollard_rho(p,k--);
find(p,k);
find(n/p,k);
}
欧拉函数:
欧拉函数就是从1到n-1跟n互质的数的个数
int eular(int n)//求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;//返回n的欧拉函数
}
void getphi()
{
memset(phi, 0, sizeof (phi));//phi存的就是欧拉函数
phi[1] = 1;
for (int i = 2; i < N; ++i)
{
if (!phi[i])
{
for (int j = i; j < N; j += i)
{
if (!phi[j])
phi[j] = j;
phi[j] = phi[j] / i * (i - 1);
}
}
}
}
POJ 1284 Primitive Roots
POJ 2773 Happy 2006
前十个完美数:
1:6
2:28
3:496
4:8128
5:33550336
6:8589869056
7:137438691328
8:2305843008139952128
9:2658455991569831744654692615953842176
10:191561942608236107294793378084303638130997321548169216
组合数学
若干排列组合的公式:
A(n,r)=n!/(n-r)!
C(n,r)=n!/((n-r)!*r!)
C(n,k)=C(n-1,k)+C(n-1,k-1)
n-1
C(n,m)=∑C(i,m-1)=C(n-1,m-1)+C(n-2,m-1)+...+C(2,m-1)+C(1,m-1)
i=1
字典序号与全排列的关系点击打开链接
字典序法中,对于数字1、2、3......n的排列,不同排列的先后关系是从左到右逐个比较对应的数字的先后来决定的。例如对于5个数字的排列12354和12345,排列12345在前,排列12354在后。按照这样的规定,5个数字的所有的排列中最前面的是12345,最后面的是54321。
生成下一个序列的字典序:
设P是1~n的一个全排列:p=p1p2......pn=p1p2......pj-1pjpj+1......pk-1pkpk+1......pn
1)从排列的右端开始,找出第一个比右边字符小的字符的序号j(j从左端开始计算),即 j=max{i|pi<pi+1}
2)在pj的右边的字符中,找出所有比pj大的字符中最小的字符pk,即 k=max{i|pi>pj}(右边的字符从右至左是递增的,因此k是所有大于pj的字符中序号最大者)
3)对换pj,pk
4)再将pj+1......pk-1pkpk+1pn倒转得到排列p''=p1p2.....pj-1pjpn.....pk+1pkpk-1.....pj+1,这就是排列p的下一个下一个排列。
bool nextpermutation(char *s)//s为序列的数组
{
int i=len-1,j;//len为长度
while(i>0&&s[i-1]>=s[i])i--;//从右往左找第一个比右边小的序号j
if(!i)return false;//如果找不到,则是最后一个排列
int mp=i;
for(j=i+1; j<len; j++)//在pj的右边的字符中,找出所有比pj大的字符中最小的字符pk
{
if(s[i-1]>=s[j])continue;
if(s[j]<s[mp])mp=j;
}
swap(s[mp],s[i-1]);//交换pj,pk
sort(s+i,s+len);//再将pj+1......pk-1pkpk+1pn倒转得到排列p''=p1p2.....pj-1pjpn.....pk+1pkpk-1.....pj+1
return len;
}
当然在STL里面有next_permutation这个函数,它可以一步就求出下一个字典序。
next_permutation(str.begin(),str.end());//str为string型的
已知排列求序号:
int perm2num(int n,int *p)//n为元素的个数,p为数组存的元素
{
int i,j,num=0,k=1;
for(i=n-2; i>=0; k*=n-(i--))//每一轮后k=(n-i)!,下标从0开始
for(j=i+1; j<n; j++)
if(p[j]<p[i])num+=k;//是否有逆序,如有,统计
return num;//返回的num即为序号(从0号位开始的)
}
已知序号求排列:
void num2perm(int n, int *p,int num)//n是元素个数,排列序号是num,对应排列是数组p
{
int i,j; //求逆序数数组
for(i=n-1; i>=0; i--)p[i]=num%(n-i),num/=n-i;
for (i=n-1; i; i--)
for (j=i-1; j>=0; j--)
if (p[j]<=p[i]) p[i]++; //根据逆序数数组进行调整
}
POJ 1146 ID Codes
POJ 1833 排列
POJ 1716 排列2
卡特兰数:
POJ 2084 Game of Connections
Polya定理总结:
设G是一个集合,*是G上的二元运算,如果(G,*)满足下面的条件:
封闭性:对于任何a,b∈G,有a*b∈G;
结合律:对任何a,b,c∈G有(a*b)*c=a*(b*c);
单位元:存在e∈G,使得对所有的a∈G,都有a*e=e*a=a;
逆元:对于每个元素a∈G,存在x∈G,使得a*x=x*a=e,这个时候记x为a-1,称为a的逆元,那么则称(G,*)为一个群。
定义:设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)。
Polya定理:设G={π1,π2,π3........πn}是X={a1,a2,a3.......an}上一个置换群,用m中颜色对X中的元素进行涂色,那么不同的涂色方案数为:1/|G|*(mC(π1)+mC(π2)+mC(π3)+...+mC(πk)). 其中C(πk)为置换πk的循环节的个数。
POJ 1286 Necklace of Beads
POJ 2409 Let it Beads
哈达玛矩阵:
哈达玛(Hadamard)矩阵是由+1和-1元素构成的且满足Hn*Hn’=nI(这里Hn’为Hn的转置,I为单位方阵)n阶方阵。
性质1:Hn正交方阵,所谓正交矩阵指它的任意两行(或两列)都是正交的。
性质2:任意一行(列)的所有元素的平方和等于方阵的阶数。即:设A为n阶由+1和-1元素构成的方阵,若AA‘=nI(这里A’为A的转置,I为单位方阵)
性质3:Hadamard矩阵的阶数都是2或者是4的倍数。
求哈达玛矩阵:
#include<stdio.h>
int H[2][2]= {{1,1},{1,-1}};
long long dfs(long long x,long long y)
{
if(x<2&&y<2)return H[x][y];
int i=0,j=0;
if(x>=2)
{
i=1;
while(i*2<=x)
i*=2;
x-=i;
}
if(y>=2)
{
j=1;
while(j*2<=y)
j*=2;
y-=j;
}
if(i==j)return -dfs(x,y);//位于4区
else if(i>j)return dfs(x,y+j);//位于3区
else return dfs(x+i,y);//位于2区
}
int main()
{
long long x,y,i,j;
scanf("%lld%lld",&x,&y);
for(i=0; i<x; i++)
{
for(j=0; j<y-1; j++)
printf("%d ",dfs(i,j));
printf("%d\n",dfs(i,j));
}
return 0;
}
POJ 2961 Sylvester construction
拉丁方阵:
拉丁方 阵:用从1开始的n个连续正整数排成n行n列的方阵,如果每一行和每一列都没有重复的数,就称为一个n阶拉丁方。
获取n阶最小拉丁方阵:
#include<stdio.h>
#include<string.h>
#include<stdlib.h>
#define M 107
int g[M][M];
int n;//n代表n阶拉丁方
bool mx[M][M],my[M][M];//mx[x][i]:第x行填过i my[y][i]:第y行填过i
void Run(int x,int y)
{
int i,j;
if(x>n)
{
for(i=1; i<=n; i++)
{
for(j=1; j<=n; j++)
printf("%4d",g[i][j]);
printf("\n");
}
printf("\n");
exit(0);
}
for(i=1; i<=n; i++)
if(!mx[x][i]&&!my[y][i])
{
mx[x][i]=1;
my[y][i]=1;
g[x][y]=i;
Run(x+y/n,y%n+1);
mx[x][i]=0;
my[y][i]=0;
g[x][y]=0;
}
}
int main()
{
while(scanf("%d",&n)!=EOF)
{
memset(g,0,sizeof(g));
memset(mx,false,sizeof(mx));
memset(my,false,sizeof(my));
Run(1,1);
}
return 0;
}