数论小专题总结

（主要是各基本知识点的例题，查缺补漏，希望以后比赛不要有本来应该掌握的知识点，却不知道的情况）

1、一元线性同余方程组的一般性解法（拓展gcd、剩余定理）

/*poj2891*解一次线性同余方程组/
/x%a[i]==r[i]*/
/*通解*/
#include<iostream>
#include<string.h>
#include<stdio.h>
#define maxn  10000+5
#define LL long long
using namespace std;
LL r[maxn],m[maxn];

void Ex_gcd(LL a,LL b,LL &d,LL &x,LL &y){
    if (!b){
        x=1;y=0;d=a;
    }else{
        Ex_gcd(b,a%b,d,y,x);
        y=y-x*(a/b);
    }
}
LL Ex_crt(LL N){
    LL M=m[1],R=r[1],x,y,d;
    for(int i=2;i<=N;i++){
        Ex_gcd(M,m[i],d,x,y);
        if ((r[i]-R)%d!=0) return -1;//无解
        x=(r[i]-R)/d*x%(m[i]/d);
        R=R+x*M;
        M=M/d*m[i];
        R%=M;
    }
    if (R<0) R+=M;
    return R;
}
int n;
int main(){
    while(~scanf("%I64d",&n)){
        for(int i=1;i<=n;i++){
            scanf("%I64d%I64d",&m[i],&r[i]);
        }
        LL ans=Ex_crt(n);
        printf("%I64d
",ans);
    }
    return 0;
}

/*hdu3579解一次线性同余方程组/
/x%a[i]==r[i]*/
/*通解*/
/*注意，题目要求的是最小“正”整数解，如果同余方程组的解为0，则解是A[i]的最小公倍数*/
#include<iostream>
#include<string.h>
#include<algorithm>
#include<stdlib.h>
#include<stdio.h>
#define maxn  100+5
#define typed long long
using namespace std;
typed r[maxn],m[maxn];

typed gcd(typed a,typed b){
    return b==0?a:gcd(b,a%b);
}
void Ex_gcd(typed a,typed b,typed &d,typed &x,typed &y){
    if (!b){
        x=1;y=0;d=a;
    }else{
        Ex_gcd(b,a%b,d,y,x);
        y=y-x*(a/b);
    }
}
typed Ex_crt(typed N){
    typed M=m[1],R=r[1],x,y,d;
    for(int i=2;i<=N;i++){
        Ex_gcd(M,m[i],d,x,y);
        if ((r[i]-R)%d!=0) return -1;//无解
        x=(r[i]-R)/d*x%(m[i]/d);
        R=R+x*M;
        M=M/d*m[i];
        R%=M;
    }
    if (R<0) R+=M;
    return R;
}
typed n,T;
typed lcm(){
    if (n==1) return m[1];
    sort(m+1,m+n+1);
    typed GCD=gcd(m[2],m[1]);
    for(int i=3;i<=n;i++)
    GCD=gcd(m[i],GCD);
    typed ans=1;
    for(int i=1;i<=n;i++) ans=ans*m[i]/GCD;
    return ans;
}
int main(){
    cin>>T;
    for(int cas=1;cas<=T;cas++){
        scanf("%I64d",&n);
        for(int i=1;i<=n;i++){
            scanf("%I64d",&m[i]);
        }
        for(int i=1;i<=n;i++){
            scanf("%I64d",&r[i]);
        }
        printf("Case %d: ",cas);
        typed ans=Ex_crt(n);
        if (ans==0) ans=lcm();
        printf("%lld
",ans);
    }
    return 0;
}

2、高次同余方程babystep算法

/*HDU - 2815
题意比较有趣，
但化成数学模型，
就是求K^D=N mod P;
K,P,N已知，求D的最小整数解。
有个小坑，N<P 才可
所以就是个解高次同余方程的题目
*/
#include<cmath>
#include<iostream>
#include<string.h>
#include<stdio.h>
#define maxn  65535
#define LL long long
using namespace std;
struct hash{
    int a,b,next;
}Hash[maxn*2];
int flg[maxn+70];
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 ext_gcd(int a,int b,int &x,int &y){
    int t,ret;
    if (!b) {
        x=1,y=0;return a;
    }
    else{
        ret=ext_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;
    ext_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=a%c;
    while(b){
        if (b&1)
        ret=ret*a%c;
        a=a*a%c;
        b=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=C/tmp;
        B=B/tmp;
        D=D*A/tmp%C;
    }
    int M=(int)ceil(sqrt((double)C));
    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;
}
LL K,P,N;
int main(){
    while(cin>>K>>P>>N){
        if (N>P){
            cout<<"Orz,I can’t find D!"<<endl;
            continue;
        }
        N=N%P;
        int tmp=BabyStep(K,N,P);//参数是底数，被除数，除数
        if (tmp<0)
            cout<<"Orz,I can’t find D!"<<endl;
        else
            printf("%d
",tmp);
    }
    return 0;
}

/*Poj 3243
就是求K^X=N mod P;
K,P,N已知，求X的最小整数解
解高次同余方程的题目
*/
#include<cmath>
#include<iostream>
#include<string.h>
#include<stdio.h>
#define maxn  65535
#define LL long long
using namespace std;
struct hash{
    int a,b,next;
}Hash[maxn*2];
int flg[maxn+70];
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 ext_gcd(int a,int b,int &x,int &y){
    int t,ret;
    if (!b) {
        x=1,y=0;return a;
    }
    else{
        ret=ext_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;
    ext_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=a%c;
    while(b){
        if (b&1)
        ret=ret*a%c;
        a=a*a%c;
        b=b>>1;
    }
    return ret;
}
int BabyStep(int A,int B,int C){//A：底数，B:被除数，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=C/tmp;
        B=B/tmp;
        D=D*A/tmp%C;
    }
    int M=(int)ceil(sqrt((double)C));
    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;
}
LL K,P,N;
//K^X=N mod P;
int main(){
    while(cin>>K>>P>>N){
        if (K==0 && P==0 && N==0) break;
        int ans=BabyStep(K,N,P);
        if (ans==-1) puts("No Solution");else cout<<ans<<endl;
    }
    return 0;
}

3、中国剩余定理

/*Poj1006中国剩余定理：
给定n组数对，使得:
A[i]==X(mod M[i])

剩余定理描述：
设m1,m2,…,mk是两两互素的正整数，则同余
方程组
          x=b1 mod m1
          x=b2 mod m2
               …
          x=bk mod mk
    有唯一解：x=（b1T1M1+…+bkTkMk  ）mod M
    其中，M=m1×m2×…mk（即lcm）
           Mi=M/mi
           Ti=Mi(-1) mod mi 乘法逆元


pps:
 应用中国剩余定理，可以简化一些复杂的运算。
 如计算  47317 mod 713
     解：因为713=2331，令x=47317，
        则计算x mod 713等价于求解同余方程：
            x=47317 mod 23=1317 mod 23
            x=47317 mod 31=817 mod 31
   即： x=6 mod 23
         x=2 mod 31
   根据中国剩余定理可得解x=374(mod 713)
   因此：47317 mod 713=374(mod 713)

*/
/*对于这道题啊，因为数据量太小，手动计算即可，
M=23*28*33=21252
ans=(b1*924*6+b2*759*19+b3*644*2-d + M)% M;//最小正整数解
*/
#include<iostream>
#include<string>
#include<algorithm>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<queue>
#include<map>
#include<stack>
#include<set>
#include<cstdlib>

using namespace std;

int d,b1,b2,b3;
int cas=0;
int main(){
    while(cin>>b1>>b2>>b3>>d){
        if (d==-1 && b1==-1  && b2==-1 && b3==-1) break;
        cas++;
        int M=23*28*33;
        int ans=(b1*924*6+b2*759*19+b3*644*2-d + M)% M;
        if (ans==0) ans=M;//正整数解！
        printf("Case %d: the next triple peak occurs in %d days.
",cas,ans);
    }
    return 0;
}

 4、二元一次不定方程

/*poj2142
解二元一次不定方程
题目描述：给定质量为a和b的砝码，不限天平左右可以放置砝码，称出一个c的质量的物品就可以了
ax+by=d列出数学方程，x和y的正负表示放置的左右，现在要求
|X|+|Y|尽量小，然后是|aX|+|BY|尽量小

http://116.255.173.144/index.php?c=article&a=read&id=56659
因为欧几里德扩展算法只能求ax + by = gcd(a, b)；所以，先求a, b的最大公约数gcd_ab，对于ax + by = d，令a, b, d同除以gcd_ab，然后可求得(a/gcd_ab)x + (b/gcd_ab)y = 1的最小整数解，然后将求得的这组解乘以d/gcd_ab得x0, y0，就是ax + by = d一组解，由数论知识得ax + by = d的所有解为： x = x0 - bt, y = y0 + at（t为参数），那么，要使|x|+|y|最小，可令x = 0求得一组解，再令y = 0求得一组解，两组比较取较小者就好。

一般方法，拓展欧几里得解出一组特解（x0，y0）
x0=x0*(c/d);
y0=y0*(c/d);
令d=gcd(a,b),a0=a/d,b0=b/d;
则通解是x=x0+a0t，y=y0-b0t（t不所谓正负）

这道题，一般都是分析函数的单调性了（需要中学数学知识）：
假设a0>b0，当|x|+|y|最小时，有-b0<y<b0（可以画出函数图像）
尝试几组解即可
//http://hi.baidu.com/lydrainbowcat/item/07913682130d07dcd1f8cd83
*/
#include <iostream>
#include <string.h>
#include <stdio.h>
#include <stdlib.h>
#include <cmath>
#define LL long long
using namespace std;

LL abs(LL x){
    if (x<0) return -x;else return x;
}
void Ex_gcd(LL a,LL b,LL &x,LL &y,LL &d){//d==gcd(a,b)
    if (b==0) {
        d=a;x=1;y=0;
    }else{
        Ex_gcd(b,a%b,y,x,d);y-=x*(a/b);
    }
}
int main(){
    LL a,b,c,d,x0,y0;
    while(cin>>a>>b>>c){
        if (a==0 && b==0 && c==0) break;
        LL chan=0;
        if (a<b){
            chan=1;
            swap(a,b);
        }
        Ex_gcd(a,b,x0,y0,d);
        x0=x0*c/d;
        y0=y0*c/d;
        //令y0-b0t==0解出一个特解t0，然后在t0周围测试
        LL t0=(y0)*d/a;
        LL a0=a/d,b0=b/d;
        LL max_xy=999999999;
        LL ansx,ansy;
        LL max_xayb;
        for(LL t=t0-10;t<=t0+10;t++){
            LL nx=x0+b0*t;
            LL ny=y0-a0*t;
            if(max_xy>abs(nx)+abs(ny)){
                max_xy=abs(nx)+abs(ny);
                max_xayb=abs(nx*a)+abs(ny*b);//开始时忘记更新
                ansx=nx;ansy=ny;
            }
             else if (max_xy==abs(nx)+abs(ny)){
                if (max_xayb>abs(nx*a)+abs(ny*b)){
                    max_xayb=abs(nx*a)+abs(ny*b);
                    ansx=nx;ansy=ny;
                }
            }
        }
        if (chan) swap(ansx,ansy);
        cout<<abs(ansx)<<" "<<abs(ansy)<<endl;
    }
    return 0;
}

5、n元一次不定方程

/*poj1091跳蚤
很经典的一道题，隐藏在活泼的外表下的是，n元一次不定方程。= =
还有，代数系统，容斥原理的应用 = =！

题目描述:
1、有一只跳蚤在一个横坐标上跳着，它有1+N张卡片，从大到小，最后一张上是M。每次可选择一张，按照上面的数字，卡片可以重复使用，随便向左还是向右跳。
2、它的目标是跳到左边一个单位的位置

分析：
1、数学方程：设每张卡片的值是a1，a2.....an,M (1+n张)
假设相应选择了xi次（正负代表左右）
则 a1x1+a2x2+a3x3.....anxn+Mx(n+1)=1
这个方程有解的充要条件是gcd(a1,a2,a3,....,M)==1
2、现在问题转化成了[1，M]的范围内找出几个数（包括M），使它们的gcd等于1，问这种分组有多少？

容斥原理
1、因为必须包含M这个数，所以一组中的数一定没有公因数（来自M因式分解的因子）
2、容斥：
总排列数M^N-sigm（排除含p1的分组，含p2的....）+sigm（含p1*p2的...含pi*pi+1的）-sigm（p1p2p3.....）....
3、举例：含p1公因子的分组数k=M/p1 是M范围内有因子p1的数的个数.选取其中的N张（有顺序，可重复），总共K^N
4、表示：状态压缩，很久没用了啊

这个博客不错：http://www.cnblogs.com/zhaoguanqin/archive/2012/02/25/2368058.html
*/
#include <iostream>
#include <string.h>
#include <stdio.h>
#include <stdlib.h>
#include <cmath>
#define maxn 3200
#define LL long long
using namespace std;
bool flag[maxn];
LL cnt;
LL prim[maxn];
LL dive[25];
void prim_calc(){//素数打表
    cnt=0;
    memset(flag,0,sizeof(flag));
    for(LL i=2;i<maxn;i++){
        if (!flag[i]) {
            prim[cnt++]=i;
            flag[i]=true;
            for(LL j=2*i;j<maxn;j+=i) flag[j]=false;
        }
    }
}
LL count;
LL divide(LL M){
    count = 0;
    for(LL i=0;i<cnt;i++){
        if (M%prim[i]==0){
            dive[count++]=prim[i];
            while(M%prim[i]==0) M=M/prim[i];
            if (M==1) break;
        }
    }
    if (M>1) dive[count++]=M;
//    for(LL i=0;i<count;i++) cout<<dive[i]<<" ";cout<<endl;
    return count;
}
LL N,M;
int main(){
    prim_calc();
    while(cin>>N>>M){
        LL dig=divide(M);//
        LL ans=pow(M,N);
        LL tmp,flag,i,j;
        for(i=1;i<(LL)(1<<dig);i++)//状态压缩
        {
            tmp=1,flag=0;
            for(j=0;j<dig;j++)
                if(i&((LL)(1<<j)))//判断第几个因子目前被用到
                    flag++,tmp*=dive[j];
            if(flag&1)//容斥原理，奇加偶减
                ans-=pow(M/tmp,N);
            else
                ans+=pow(M/tmp,N);
        }
        cout<<ans<<endl;
    }
    return 0;
}

 6、二元二次不定方程（佩尔方程）

/*HDU3292
佩尔方程（二元二次方程）
形如X^2-d*Y^2=1（d>1）
如果d为不完全数，则有解：否则，无解。
求解佩尔方程的两种方法
1、暴力法（本题）
2、连分数法（zoj3759：http://www.cnblogs.com/little-w/p/3588285.html）
求佩尔方程的解（矩阵快速幂）
|x[k]|=|x1 d*y1|^(k-1) |x1|
|y[k]|=|y1 x1  |       |y1|
思路：暴力或者连分数法求x1+矩阵快速幂
*/
#include <iostream>
#include <string.h>
#include <stdio.h>
#include <cmath>
#define LL long long
#define mod 8191
using namespace std;

int x,y,d,k;
bool isquare(){
    int sq=sqrt(d+0.0);
    if (((int)sq*(int)sq)==d) return true;
    else return false;
}
void search(){
    y=1;
    while(1){
        x=(long long)sqrt(d*y*y+1);
        if (x*x-d*y*y==1)
           break;
        y++;
    }
}
struct Mat
{
    LL mat[3][3];
};
Mat operator*(Mat a,Mat b)
{
    Mat c;
    memset(c.mat,0,sizeof(c.mat));
    for (int i=0;i<2;i++)
        for (int j=0;j<2;j++)
            for (int k=0;k<2;k++)
                if (a.mat[i][k] && b.mat[k][j])
                   c.mat[i][j]=(c.mat[i][j]+a.mat[i][k]*b.mat[k][j])%mod;
    return c;
}
Mat operator+(Mat a,Mat b)
{
    Mat c;
    memset(c.mat,0,sizeof(c.mat));
    for (int i=0;i<2;i++)
        for (int j=0;j<2;j++)
            c.mat[i][j]=(a.mat[i][j]+b.mat[i][j])%mod;
    return c;
}
Mat E;
void builtE()
{
    memset(E.mat,0,sizeof(E.mat));
    for (int i=0;i<2;i++)
    {
        E.mat[i][i]=1;
    }
}
Mat operator^(Mat A,int x)
{
    Mat c;
    builtE();
    c=E;
    for ( ; x; x>>=1)
    {
        if ( x & 1)
            c = c * A;
        A = A * A;
    }
    return c;
}
int main(){
    while(cin>>d>>k){
        if (isquare()) {
            printf("No answers can meet such conditions
");
            continue;
        }
        search();
        Mat A;
        A.mat[0][0]=x,A.mat[0][1]=d*y,A.mat[1][0]=y,A.mat[1][1]=x;
        A=A^(k-1);
        Mat B;
        B.mat[0][0]=x,B.mat[1][0]=y;
        B=A*B;
        cout<<B.mat[0][0]<<endl;
    }
    return 0;
}

 7、本原毕达哥拉斯三元组定理

/*fzu1669
毕达哥拉斯三元组的应用

本原毕达哥拉斯定理：
x=m^2-n^2;
y=2mn;
z=m^2+n^2;
m>n,n奇数m偶数或者n偶数m奇数
注意是本原，即(x,y,z)==1,可通过计算本原加上筛选求得一定范围内的所有(x,y,z)

题意：
a+b+c<=L,
a,b是直角边，c是斜边
给定L ,求满足上式的三角形的个数
思路：通过上面的定理，枚举出m,n(需满足gcd(n,m)==1)，找到本原，筛选可得解
具体实现：假设m>n,则m<sqrt(L),再枚举比m小的n即可，因为L≤2000000，m<1000,单次最多不会超过500000左右
*/
#include <iostream>
#include <string.h>
#include <stdio.h>
#include <cmath>
#define LL long long

using namespace std;
int L;
int gcd(int a,int b){
    if (b==0) return a;
    else return gcd(b,a%b);
}
int main(){
    while(cin>>L){
        int up=(int)(sqrt(L+0.0));
        int ans=0;
        for(int m=1;m<=up;m++){
            for(int n=m-1;n>=1;n--){
                if ((m+n)%2==0 || gcd(m,n)!=1) continue;
                int x=m*m-n*n;
                int y=2*n*m;
                int z=m*m+n*n;
                ans+=L/(x+y+z);
            }
        }
        cout<<ans<<endl;
    }
    return 0;
}

8、欧拉函数的应用 POJ 3696

{参考:http://blog.csdn.net/xieshimao/article/details/6689780}

9、lucas定理 HDU3944