组合数取模深度解析

来源：互联网发布：手机怎么注册淘宝店铺编辑：程序博客网时间：2024/04/30 05:32

组合数取模深度解析

分类：数论2013-02-13 15:47 63人阅读评论(0) 收藏举报

今天我们来讨论如何求：C(n,m)%p

首先如果n,m,p的范围均为[1,1000]，那么我们完全可以用两个循环搞定，C(n,m)=(C(n-1,m)%p+C(n-1,m-1)%p)%p。

但是如果n，m，p的范围均为[1,10^9]，另外p为素数，则完全可以用Lucas定理搞定。

如果n，m，p范围还是[1,10^9]，另外p如果为合数，则我们可以先把p分解：p=p1^a1*p2^a2*...*pk^ak。

问题主要转化为如何求C(n,m)%(p^c) 其中p为素数。求完直接合并一个模方程即可。p^c 的规模大约是10^5。c 不是1，lucas阻止不了它。

n,m太大，因子分解也阻止不了它。

下面介绍一种好方法:

假设 p = 3, c = 2，也就是mod 9，假设n = 19，n! = 1 * 2 * 3 * 4 * 5 * 6 * 7 * 8 * …… * 19

要是可以快速得到 n! 中除掉3 以后 mod 9的结果，那么多好呀！看3多讨厌，直接砍fac( int n) ：

n! = [ 1 * 2 * 4 * 5 * 7 * 8 * … * 16 * 17 * 19 ] * (3 * 6 * 9 * 12 * 15 * 18)=[ 1 * 2 * 4 * 5 * 7 * 8 * … * 16 * 17* 19 ] * 3^6*( 1 * 2 * 3 * 4 * 5 * 6)

然后发现后面的一坨实际上是 fac( n / p) ，再看前半部分，尼玛是以 p^c 为周期的啊！！！

[1 * 2 * 4 * 5 * 7 * 8 ] = [10 * 11 * 13 * 14 * 16 *17 ] = (mod 9)，于是说白了，对于前面的部分，由于周期，都是浮云了，下面是孤立出来的19

可以知道孤立出来的长度不超过 p^c ,于是暴力啊，暴力啊！于是完美解决n! 中和 p无关的项 mod p^c的值！！！

接下来是分母部分，一模一样，无非多了一个求逆元(因为都和p没关系了，逆元必然存在)

我们来分析一下，这样的复杂度是如何的呢，每次递归，规模变为原来的 1/p，logp N的啊！！！

当然是层数= =于是问题完美解决

[cpp] view plaincopy
#include <iostream>  
#define LL long long  
  
LL p,k,r,x,y;  
  
void extend_Euclid(LL a,LL b)  
{  
    if(b==0)  
    {  
        x=1;  
        y=0;  
        return;  
    }  
    extend_Euclid(b,a%b);  
    LL temp=x;  
    x=y;  
    y=temp-(a/b)*y;  
}  
  
LL quick_mod(LL a,LL b,LL m)  
{  
    LL ans=1;  
    a%=m;  
    while(b)  
    {  
  
        if(b&1)  
        {  
            ans=ans*a%m;  
            b--;  
        }  
        b>>=1;  
        a=a*a%m;  
    }  
    return ans;  
}  
  
LL multi(LL p,LL k)  
{  
    LL i,ret=1;  
    for(i=1;i<=k;i++)  
       ret*=p;  
    return ret;  
}  
  
LL sum(LL n)  
{  
    LL ans=0;  
    while(n)  
    {  
        ans+=n/p;  
        n/=p;  
    }  
    return ans;  
}  
  
LL Solve(LL n)  
{  
    LL i;  
    if(n==0) return 1;  
    LL ans1=1;  
    for(i=1;i<=r;i++)  
    {  
        if(i%p!=0)  
           ans1=ans1*i%r;  
    }  
    LL ret=quick_mod(ans1,n/r,r);  
    LL x=n-n%r+1;  
    LL ans2=1;  
    for(i=x;i<=n;i++)  
    {  
        if(i%p!=0)  
           ans2=ans2*i%r;  
    }  
    ret=ret*ans2%r*Solve(n/p)%r;  
    return ret;  
}  
  
int main()  
{  
    LL n,m;  
    LL ans,ret;  
    LL x1,x2,x3,x4,x5;  
    while(std::cin>>n>>m>>p>>k)  
    {  
        r=multi(p,k);  
        x1=sum(n);  
        x2=sum(n-m);  
        x3=sum(m);  
        x4=x1-x2-x3;  
        ans=multi(p,x4);  
        x1=Solve(n);  
        x2=Solve(n-m);  
        x3=Solve(m);  
        x4=x2*x3%r;  
        extend_Euclid(x4,r);  
        x=(x%r+r)%r;  
        x5=x1*x%r;  
        ret=ans*x5%r;  
        std::cout<<ret<<std::endl;  
    }  
    return 0;  
}  

Petr在CodeForces上给世界级选手出的10道难题2 ：

输出C(n,k)的最后10位，0<=k<=n<=10^18,n>=1

HDU3802深度解析

分类：数论2013-02-24 18:15 66人阅读评论(0) 收藏举报

题目：Ipad,IPhone AC大神的题目

对于这个题我们首先得知道二次剩余的一个结论，就是如果p为奇素数，a为p的二次剩余，那么有a^((p-1)/2)=1(mod p)，

如果a为p的二次非剩余，则a^((p-1)/2)=-1(mod p)，如果为非剩余则答案直接为0，否则前面那两坨为4，只需考虑后面的，后面的用矩阵构造。

在降幂的时候注意必须是a和b都为p的二次剩余时才可以，否则就错误。还有在p已经为素数的情况下，对于a^bmod p 那么当b大于p-1时

a^(b%(p-1))%p与a^(b%(p-1)+p-1)%p是一样的结果。

[cpp] view plaincopy
#include <iostream>  
#define N 2  
  
typedef struct  
{  
    long long m[N][N];  
}Matrix;  
  
Matrix I={1,0,0,1};  
  
Matrix multi(Matrix a,Matrix b,long long p)  
{  
    Matrix c;  
    int i,j,k;  
    for(i=0;i<N;i++)  
    {  
        for(j=0;j<N;j++)  
        {  
            c.m[i][j]=0;  
            for(k=0;k<N;k++)  
            {  
                c.m[i][j]+=a.m[i][k]*b.m[k][j]%p;  
            }  
            c.m[i][j]%=p;  
        }  
    }  
    return c;  
}  
  
Matrix multi_mod(Matrix a,long long k,long long m)  
{  
    Matrix ans=I,p=a;  
    while(k)  
    {  
        if(k&1)  
        {  
            ans=multi(ans,p,m);  
            k--;  
        }  
        k>>=1;  
        p=multi(p,p,m);  
    }  
    return ans;  
}  
  
long long quick_mod(long long a,long long b,long long m)  
{  
    long long ans=1;  
    a%=m;  
    while(b)  
    {  
        if(b&1)  
        {  
            ans=ans*a%m;  
            b--;  
        }  
        b>>=1;  
        a=a*a%m;  
    }  
    return ans;  
}  
  
int main()  
{  
    int t;  
    Matrix res;  
    long long a,b,n,p;  
    long long x,y,z,ret;  
    std::cin>>t;  
    while(t--)  
    {  
        Matrix A={0,1,1,1};  
        std::cin>>a>>b>>n>>p;  
        x=quick_mod(a,(p-1)>>1,p);  
        y=quick_mod(b,(p-1)>>1,p);  
        if(x!=1||y!=1)  
        {  
            std::cout<<"0"<<std::endl;  
            continue;  
        }  
        Matrix ans=multi_mod(A,n,p-1);  
        z=ans.m[0][0]+ans.m[0][1];  
  
        z%=(p-1);  
        res.m[0][0]=2*(a+b)%p;  
        res.m[0][1]=-(a-b)*(a-b)%p;  
        res.m[1][0]=1;  
        res.m[1][1]=0;  
        ans=multi_mod(res,z-1,p);  
        ret=(ans.m[0][0]*((2*a+2*b)%p)%p+ans.m[0][1]*2%p)%p;  
        ret=(ret+p)%p;  
        ret=(ret*4)%p;  
        std::cout<<ret<<std::endl;  
    }  
    return 0;  
} 