数学模版

一.数论部分

1.素数测试

#include <iostream>  
#include <cmath>  
#include <algorithm>  
#include <cstdlib>  
#include <cstdio>  
#include <map>  
using namespace std;  
#define Times 10  
typedef long long LL;  
map<LL,int>mp;  
LL random(LL n)  
{  
    return LL ((double)rand()/RAND_MAX*n+0.5);  
}  
//避免出现了相乘取余前造成的溢出  
LL multi(LL a,LL b,LL mod)  
{  
    LL ans=0;  
    while(b)  
    {  
        if(b&1)  
        {  
            b--;  
            ans=(ans+a)%mod;  
        }  
        else  
        {  
            b/=2;  
            a=(a+a)%mod;  
        }  
    }  
    return ans;  
}  
LL Pow(LL a,LL b,LL mod)  
{  
    LL ans=1;  
    while(b)  
    {  
        if(b&1)  
        {  
            b--;  
            ans=multi(ans,a,mod);  
        }  
        else  
        {  
            b/=2;  
            a=multi(a,a,mod);  
        }  
    }  
    return ans;  
}  
// 二次探测找到一个不是x=1或者x=p-1的解就一定不是素数  
bool witness(LL a,LL n)  
{  
    LL d=n-1;  
    while(!(d&1))//幂是偶数时有探测的必要  
        d>>=1;  
    LL t=Pow(a,d,n);  
    while(d!=n-1 && t!=1 && t!=n-1)//如果存在一个非1或非p-1的解，就不是素数  
    {  
        t=multi(t,t,n);  
        d<<=1;  
    }  
    return t==n-1 || d&1 ;  
}  
bool miller_rabin(LL n)  
{  
    if(n==2)  
        return true;  
    if(n<2 || !(n&1))  
        return false;  
    for(int i=1;i<=Times;i++)  
    {  
        LL a=random(n-2)+1;  
        if(!witness(a,n))  
            return false;  
    }  
    return true;  
}  
LL gcd(LL a,LL b)  
{  
    if(b==0)  
        return a;  
    return gcd(b,a%b);  
}  
LL pollard_rho(LL n,int c)  
{  
    LL x,y,d,i=1,k=2;  
    x=random(n-2)+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(LL n,LL c)  
{  
    if(n==1)  
        return ;  
    if(miller_rabin(n))  
    {  
        mp[n]++;  
        return ;  
    }  
    LL p=n;  
    while(p>=n)  
        p=pollard_rho(p,c--);  
    find(p,c);  
    find(n/p,c);  
}  

2.素数筛法等

分段筛素数：  
  
100W 差不多接近8W素数  
  
void make_prime()    
{    
    m.clear();    
    unsigned __int64 tmp;    
    memset(cnt,true,sizeof(cnt));    
    for(int i=1;i<=n_prime;i++)    
    {    
        tmp=l/num[i];    
        while(tmp*num[i]<l||tmp<=1)    
            tmp++;    
        for(unsigned __int64 j=tmp*num[i];j<=r;j+=num[i])    
        {    
            if(j>=l&&j<=r)    
                cnt[j-l]=0;    
        }    
    }    
    if(l==1)    
        cnt[0]=0;    
    for(unsigned __int64 i=0;i+l<=r;i++)    
        if(cnt[i]&&i+l<=r&&i>=0)    
            m[i+l]++;    
}

求因子和

void factor_sum()  
{  
    for(int i=1;i<=maxn;i++)  
        prime[i]=1;  
    for(int i=2;2*i<=maxn;i++)  
            for(int j=2*i;j<=maxn;j+=i)  
                prime[j]+=i;  
    prime[1]=0;  
}  
另外一种求法：  
(pn^(xn+1)-1)/(pn-1)的累乘

3.欧拉函数

phi[i]=i*(1-1/p1)*(1-1/p2)*(1-1/p3)...

单独求欧拉值：  
long long euler(long long x)  
{  
    if(x==1)  
        return 0;  
    long long i,res=x;  
    for(i=2;i*i<=x;i++)  
        if(x%i==0)  
        {  
            res=res/i*(i-1);  
            while(x%i==0)  
                x/=i;  
        }  
    if(x>1)  
        res=res/x*(x-1);  
    return res;  
}  

递推求欧拉值：  
void Phi()  
{  
    for(int i=1;i<=maxn;i++)  
        phi[i]=i;  
    for(int i=2;i<=maxn;i+=2)  
        phi[i]/=2;  
    for(int i=3;i<=maxn;i+=2)  
        if(phi[i]==i)  
            for(int j=i;j<=maxn;j+=i)  
                phi[j]=phi[j]/i*(i-1);  
}

欧拉公式的引伸：  
    小于或等于n的数中，与n互质的数的总和为：φ(x) * x / 2。(n>1)

线性求欧拉值：  
const int maxn=1000001;  
long long euler[maxn];  
long long prime[maxn/3];  
bool isprime[maxn];  
void solve()  
{  
    int cnt=0;  
    memset(isprime,1,sizeof(isprime));  
    for(int i=2;i<maxn;i++)  
    {  
        if(isprime[i])  
        {  
            prime[++cnt]=i;  
            euler[i]=i-1;  
        }  
        for(int j=1;j<=cnt&&prime[j]*i<maxn;j++)  
        {  
            isprime[prime[j]*i]=0;  
            if(i%prime[j]!=0)  
                euler[i*prime[j]]=euler[i]*euler[prime[j]];  
            else  
            {  
                euler[i*prime[j]]=euler[i]*prime[j];  
                break;  
            }  
        }  
    }  
    euler[1]=0;  
}  

4.Lucas

对于组合数C(n,m)方便求出来的情况，直接求fac,然后用逆元解  
  
void init()  
{  
    fac[0]=1;  
    for(int i=1;i<=p;i++)  
        fac[i]=fac[i-1]*i%p;  
}  
long long C(long long n,long long m)  
{  
    if(n<m)  
        return 0;  
    long long ans=fac[n]*Pow(fac[n-m]*fac[m]%p,p-2,p)%p;  
    return ans;  
}  
  
不然的话，就只好一步一步逆元的求解组合数  
  
long long C(long long n,long long m)  
{  
    if(n<m)  
        return 0;  
    long long ans=1;  
    for(int i=1;i<=m;i++)  
       ans=(((n-m+i)%p)*Pow(i,p-2,p)%p)*ans%p;  
    return ans;  
}  
long long Lucas(long long n,long long m)  
{  
    if(m==0)  
        return 1;  
    return Lucas(n/p,m/p)*C(n%p,m%p)%p;  
}

二.矩阵

一类是具体题目具体 构造的，类型比较多，也有涉及到数的实部虚部，状态什么的，没什么好写的。

这里翻墙可以看到一篇比较好的资料和习题：http://zobayer.blogspot.com/2010/11/matrix-exponentiation.html

另一种就是对平面内点的变换了，具体看下图：

                           转自 matrix67

具体实现如下：（注意全部操作应为，操作矩阵 * 目标矩阵）

#define PI acos(-1.0)    
struct Matrix    
{    
    double m[4][4];    
}E,matrix[10010];    
void init()    
{    
    for(int i=1;i<=3;i++)    
        for(int j=1;j<=3;j++)    
        {    
            if(i==j)    
                E.m[i][j]=1.0;    
            else    
                E.m[i][j]=0.0;    
        }    
}    
Matrix init_translation(double _x,double _y)    
{    
    Matrix ans=E;    
    ans.m[1][3]=_x;    
    ans.m[2][3]=_y;    
    return ans;    
}    
Matrix init_zoom(double l)    
{    
    Matrix ans=E;    
    ans.m[1][1]=l;    
    ans.m[2][2]=l;    
    return ans;    
}    
Matrix init_up_down()    
{    
    Matrix ans=E;    
    ans.m[2][2]=-1;    
    return ans;    
}    
Matrix init_left_right()    
{    
    Matrix ans=E;    
    ans.m[1][1]=-1;    
    return ans;    
}    
Matrix Multi(Matrix A,Matrix B)    
{    
    Matrix ans;    
    for(int i=1;i<=3;i++)    
        for(int j=1;j<=3;j++)    
        {    
            ans.m[i][j]=0;    
            for(int k=1;k<=3;k++)    
                ans.m[i][j]+=A.m[i][k]*B.m[k][j];    
        }    
    return ans;    
}    
Matrix init_rotate(double r)    
{    
    Matrix ans=E;    
    ans.m[1][1]=cos(r);    
    ans.m[1][2]=-sin(r);    
    ans.m[2][1]=sin(r);    
    ans.m[2][2]=cos(r);    
    return ans;    
}    
Matrix last(Matrix A,Matrix B)    
{    
    Matrix ans;    
    for(int i=1;i<=3;i++)    
        for(int j=1;j<=1;j++)    
        {    
            ans.m[i][j]=0;    
            for(int k=1;k<=3;k++)    
                ans.m[i][j]+=A.m[i][k]*B.m[k][j];    
        }    
    return ans;    
}

三.方程之类

扩展欧几里得：

long long exgcd(long long a,long long b,long long &x,long long &y)  
{  
    if(b==0)  
    {  
        x=1;  
        y=0;  
        return a;  
    }  
    else  
    {  
        long long ans=exgcd(b,a%b,x,y);  
        long long t=x;  
        x=y;  
        y=t-a/b*y;  
        return ans;  
    }  
}  
//开始的时候对于方程 ax+by=c 可以方程两边同时除以gcd（a,b）   
//    c%gcd(a,b)!=0 时无解  
// 得到的解x0,y0实际上为特解，通解是 x=x0+k*b/gcd(a,b) y=y0-k*a/gcd(a,b)  
//若是求x的最小值，首先要还原(因为原方程解的是ax+by=1的解)x*=c,然后是x%=b...  
//求方程的r=gcd（a,b）个解:For（i=1...r） ans[i]=x0+（i-1）*b//这里的b就是前面的b/gcd(a,b)

扩展欧几里得求逆元：

求 (a/b)%mod   
  =((a%mod)*(b^-1%mod))%mod  
    
由exgcd可以得到  
    方程 b*b^-1 = 1(%mod)  
    即 b*b^-1+k*mod = 1  
  so  
    b^-1 = (x%mod+mod)%mod;  
      
特殊情况：  
    当b是质数时，由费马小定理可知a^(p-1)=1(mod p)  
    so  
       a^-1 = a^(p-2) 

一元线性同余方程组：

// 已知 n 和 x=a[i](mod b[i])  
  
//大体思路是两两运算，推出新的方程，然后运算  
//已知输入的线性方程为  
// x=a1(mod b1)  
// x=a2(mod b2)  
// ...  
// x=an(mod bn)  
//先看第一个 x-a1=b1y1  x-a2=b2y2  
// 联立：b1y1-b2y2=a2-a1  
bool solve()  
{  
    a[0]=a[1],b[0]=b[1];  
    for(int i=2;i<=n;i++)  
    {  
        int a0=b[0];  
        int b0=b[i];  
        int c0=a[i]-a[0];  
        int x0,y0;  
        int r=gcd(a0,b0);  
        if(c0%r!=0)  
            return false;  
        else  
        {  
            a0/=r;  
            b0/=r;  
            c0/=r;  
            exgcd(a0,b0,x0,y0);  
            x0*=c0;  
            x0=(x0%b0+b0)%b0;  
            a[0]=b[0]*x0+a[0];  
            b[0]=b[0]*b0;  
        }  
    }  
    return true; // a[0]为小于[b1,b2...bn]的非负整数解  
}

中国剩余定理：  
int CRT(int r)  
{  
    long long M=1;  
    long long i,d,x0,y0,ans=0;  
    for(i=1;i<=r;i++)  
        M*=w[i];//好像很容易溢出的样子，CRT弱爆了  
    for(i=1;i<=r;i++)  
    {  
        d=M/w[i];  
        exgcd(d,w[i],x0,y0);  
        ans=(ans+d*x0*a[i])%M;  
    }  
    while(ans<= 0)  
       ans+=M;  
    return ans;  
}  
//用w[]表示模数（两两互质），a[]表示余数，返回值是以M为模的解  

高次同余方程：

// 已知 A^x=B(mod C)  
  
#include <iostream>    
#include <cstdio>    
#include <cmath>    
#define maxn 65535    
using namespace std;    
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?a:gcd(b,a%b);    
}    
int exgcd(int a,int b,int &x,int &y)    
{    
    int t,ans;    
    if(!b)    
    {    
        x=1;    
        y=0;    
        return a;    
    }    
    ans=exgcd(b,a%b,x,y);    
    t=x;    
    x=y;    
    y=t-a/b*y;    
    return ans;    
}    
int Inval(int a,int b,int n)    
{    
    int x,y,e;    
    exgcd(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 ans=1;    
    a%=c;    
    while(b)    
    {    
        if(b&1)    
        {    
            ans=ans*a%c;    
        }    
        a=a*a%c;    
        b>>=1;    
    }    
    return ans;    
}    
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((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=find(tmp);    
        if(tmp>=0&&w!=-1)    
          return i*M+w+d;    
    }    
    return -1;    
}    
int main()    
{    
    int A,B,C;    
    while(scanf("%d%d%d",&A,&C,&B)!=EOF,A||B||C)    
    {    
        B%=C;    
        int tmp=BabyStep(A,B,C);    
        if(tmp<0)    
            printf("No Solution\n");    
        else    
            printf("%d\n",tmp);    
    }    
    return 0;    
}

佩尔方程：

已知方程 x^2 - d*y^2 = 1  
   
 可知迭代方程：  
     x(n) = x(n-1)*x1 + d*y(n-1)*y1  
     y(n) = x(n-1)*y1 + y(n-1)*x1  
       
 用矩阵快速幂求解：  
     |  x(n) |   =  | x1  dy1 |^(n-1) |  x1 |  
     |  y(n) |   =  | y1   x1 |       |　y1 |

高斯模板  // czyuan大神

int equ,var;/// equ个方程，var个变量    
int a[maxn][maxn];    
int x[maxn];///解    
bool free_x[maxn];///不确定变元    
int free_num;    
void debug()    
{    
    for(int i=0;i<equ;i++)    
    {    
        for(int j=0;j<var+1;j++)    
            printf("%d ",a[i][j]);    
        printf("\n");    
    }    
    printf("\n");    
}    
int gcd(int a,int b)    
{    
    if(b==0)    
        return a;    
    return gcd(b,a%b);    
}    
int lcm(int a,int b)    
{    
    return a/gcd(a,b)*b;    
}    
/// -2表示浮点解 -1表示无解 0表示唯一解 大于0表示自由变元个数    
int gauss()    
{    
    int i,j,k;    
    int max_r;    
    int col=0;    
    int ta,tb;    
    int LCM;    
    int tmp;    
    int free_x_num;    
    int free_index;    
    for(k=0;k<equ&&col<var;k++,col++)    
    {    
        max_r=k;    
        for(i=k+1;i<equ;i++)    
        {    
            if(abs(a[i][col])>abs(a[max_r][col]))    
                max_r=i;    
        }    
        if(max_r!=k)    
        {    
            for(j=k;j<var+1;j++)    
                swap(a[k][j],a[max_r][j]);    
        }    
        if(a[k][col]==0)    
        {    
            k--;    
            continue;    
        }    
        for(i=k+1;i<equ;i++)    
        {    
            if(a[i][col]!=0)    
            {    
                LCM=lcm(abs(a[i][col]),abs(a[k][col]));    
                ta=LCM/abs(a[i][col]);    
                tb=LCM/abs(a[k][col]);    
                if(a[i][col]*a[k][col]<0)    
                    tb=-tb;    
                for(j=col;j<var+1;j++)    
                    a[i][j]=a[i][j]*ta-a[k][j]*tb;    
            }    
        }    
    }    
    for(i=k;i<equ;i++)    
        if(a[i][col]!=0)    
            return -1; /// 无解    
    if(k<var) /// 无穷解    
    {    
        for(i=k-1;i>=0;i--)    
        {    
            free_x_num=0;    
            for(j=0;j<var;j++)    
            {    
                if(a[i][j]!=0&&free_x[j])    
                {    
                    free_x_num++;    
                    free_index=j;    
                }    
            }    
            if(free_x_num>1)    
                continue;    
            tmp=a[i][var];    
            for(j=0;j<var;j++)    
            {    
                if(a[i][j]!=0&&j!=free_index)    
                    tmp-=a[i][j]*x[j];    
            }    
            x[free_index]=tmp/a[i][free_index];    
            free_x[free_index]=0;    
        }    
        return var-k;    
    }    
    for(i=var-1;i>=0;i--)    
    {    
        tmp=a[i][var];    
        for(j=i+1;j<var;j++)    
            if(a[i][j]!=0)    
                tmp-=a[i][j]*x[j];    
        if(tmp%a[i][i]!=0)    
            return -2;    
        x[i]=tmp/a[i][i];    
    }    
    return 0;    
}

java版的高斯消元（分数） hdu2449.java

四.不是很重要的东西

[1].n!的素因子分解中素数p的幂为[n/p]+[n/p^2]+[n/p^3]+...

[2].C(2n,n)中素数p的幂为（[2n/p]-2[n/p]）+([2n/p^2]+2[n/p^2])+...+([2n/p^t]-2[n/p^t])  
      t=[logp(2n)]

[3].梅森素数：  
Lucas-Lehmer判定法:(注意期间平方运算可能会溢出)  
对于第p个梅森数Mp=2^p-1,R1=4,对于k>=2,根据Rk=(R(k-1))^2-2(mod Mp)得到Rk序列，则Mp是素数，当且仅当R(p-1)=0（mod Mp）  
bool Lucas-Lehmer(int p)  
{  
    long long ans=1;// 开始写的是 long long ans=(1<<p)  
    ans<<=p;//然后试出来，那个1是默认的int，必须改成 ans=((long long)1<<p);  
    ans--;    
    printf("%lld",ans);  
    long long R[65]={0,4};  
    for(int i=2;i<=p-1;i++)  
    {  
        long long tmp=multi(R[i-1],R[i-1],ans);  
        R[i]=(tmp-2)%ans;// multi写一个a*b%m就好  
    }  
    if(p==2) //特判  
        return true;  
    else  
    {  
        if(R[p-1]==0)  
            return true;  
        return false;  
    }  
}

[4]毕达哥拉斯三元组  
x=m*m-n*n;   y=2*m*n;  z=m*m+n*n (m，n不同为奇数和偶数)

[5]佩尔方程  
X^2 - dY^2 =1 令最小解为X1,Y1  
Xn=Xn-1*X1+d*Yn-1*Y1  
Yn=Xn-1*Y1+Yn-1*X1

[6]回文数个数  
求10^n范围内回文数的个数：  
if(n%2==0)           
   f(n)=2*(10^(n/2)-1);  
else  
   f(n)=11*10^((n-1)/2)-2;

[7]大数数字前n项  
可以对10取对数，得到的值减去floor部分，再将剩下的部分还原

[8]勒让德两平方数之和定理：R（N）=4*D1（N）-4*D3（N）  
D1(N)=(整除N 且满足d=1(mod 4)的正约数d 的个数）,  
D3(N)=(整除N 且满足d=3(mod4)的正约数d 的个数)。

[9]指数循环节  
a^b mod c = a^(b%phi(c))%c (b>phi(c))

[10]Catalan数 C（2n,n）/(n+1)  
高精度样例可以参照：hdu1131

[11] sigma(gcd(x,y))
sigma(gcd(i,n))=sigma(j*euler(n/j)+n/j*euler(j))  (n%j==0)

[12] 五边形数定理  参照：http://zh.wikipedia.org/wiki/%E4%BA%94%E9%82%8A%E5%BD%A2%E6%95%B8%E5%AE%9A%E7%90%86

微积分啥的，遇到了，拿微积分表当模版吧

傅里叶变换，不会，这模版估计也不会用，跳过

莫比乌斯反演和polya计数等还得学，或者直接丢掉

其他遇到了，再学吧，再补充。