zoukankan      html  css  js  c++  java
  • bzoj1951: [Sdoi2010]古代猪文

    数论综合。费马小定理,lucas定理,中国剩余定理,exgcd,快速幂,乘法逆元。

    首先要计算出n的每个约数,简单的sqrt(n)枚举即可。

    然后计算C(i,m)(m个中挑i个的组合数,ps:因为网上正反俩种都有,所以标注一下。。)

    设s=sum(C(i,m))

    题目要求g^(s)%mod,

    由费马小定理得 g^(s)=g^(s%(mod-1))。所以这里需要特判一下g是否等于mod。

    可是组合数计算因为用到除法,所以模数必须为质数。999911659-1=2×3×4679×35617。

    这四个数设为d[i],对这四个数分别进行计算。

    由因为组合数很大,直接计算会超时。

    需要用到lucas定理。 lucas(n,m)=lucas(n/mod,m/mod)*C(n%mod,m%mod)。这里mod不是999911659,而是那4个数。

    所以就得到4个答案m[i]。

    组合数的计算可以根据费马小定理转换成快速幂计算。

    然后再用中国剩余定理合并。

    这时得到的答案可以写成如下形式:

    s%d[0]=m[0],s%d[1]=m[1];

    改写为 s%a1=b1,s%a2=b2 (1)

    则有 a1*x+b1=a2*y+b2 (x,y分别为俩个未知数,和上面的任何数没有任何关系)

    设a=a1,b=a2,c=b2-b1

    则上式变为 ax+by=c (正负号并没有关系,因为不会用到y)

    计算t=extended_gcd(a,b,x,y)

    如果c%t !=0 的话,说明此方程是无解的。在本题中t=1。(一下都基于本题t=1的情况)

    然后我们计算出的是 ax+by=t。 所以我们将每一项乘以c。

    x=(x*c%b+b)%b。

    当然这样的x有无数组为 x1=x+k*b。

    带入(1)式有  s=a1(x+k*b)+b1。

    s=a1*b*k+a1*x+b1.

    所以有 s%(a1*b)=a1*x+b1。

    所以 b1=b1+a1*x,a1=a1*b。

    在本题中这样操作完以后a1>b1,实际上即使t!=1,b1也会小于1。

    最后合并完返回b1即可。

    EN TARO Tassadar

    #include<cstdio>
    #include<algorithm>
    #include<cstring>
    #include<cmath>
    #define LL long long 
    using namespace std;
    const int mod = 999911659;
    const int d[]={2,3,4679,35617};
    int fac[4][40000],m[4];
    int n,g;
    
    int power(LL a,int e,int mod) {
        LL res=1;
        while(e) {
            if(e&1) res=res*a%mod;
            a=a*a%mod;
            e>>=1;
        }
        return res;
    }
    
    int C(int n,int m,int x) {
        if(n<m) return 0;    
        return fac[x][n]*power(fac[x][m]*fac[x][n-m]%d[x],d[x]-2,d[x])%d[x];
    }
    
    int lucas(int n,int m,int x) {
        if(m==0) return 1;
        return lucas(n/d[x],m/d[x],x)*C(n%d[x],m%d[x],x)%d[x];    
    }
    
    int exgcd(int a,int b,int &x,int &y) {
        if(b==0) {
            x=1; y=0; return a;
        }
        else {
            int t=exgcd(b,a%b,y,x);
            y-=(a/b)*x;
            return t;
        }
    }
    
    int solve() {
        int a1,b1,a2,b2,a,b,c,x,y;    
        a1=d[0]; b1=m[0];
        for(int i=1;i<4;i++) {
            a2=d[i];b2=m[i];
            a=a1,b=a2,c=b2-b1;
            exgcd(a,b,x,y);
            x=(c*x%b+b)%b;//if(c%t) printf("Test
    ");
            b1=b1+a1*x;
            a1=a1*b;
        }
        return b1;
    }
    
    int main() {
        for(int i=0;i<4;i++) {
            fac[i][0]=1;
            for(int j=1;j<=40000;j++) fac[i][j]=fac[i][j-1]*j%d[i];
        }
        scanf("%d%d",&n,&g);
        if(g==mod) printf("0
    ");
        else {
            g=g%mod;    
            int n2=(int)sqrt(n+0.5);
            for(int i=1,j;i<=n2;i++) 
            if(n%i==0) {
                j=n/i;
                for(int x=0;x<4;x++) {
                    m[x]=(m[x]+lucas(n,i,x))%d[x];
                    if(i!=j) m[x]=(m[x]+lucas(n,j,x))%d[x];    
                }
            }
            printf("%d
    ",power(g,solve(),mod));
        }
        return 0;    
    }
  • 相关阅读:
    Running ASP.NET Applications in Debian and Ubuntu using XSP and Mono
    .net extjs 封装
    ext direct spring
    install ubuntu tweak on ubuntu lts 10.04,this software is created by zhouding
    redis cookbook
    aptana eclipse plugin install on sts
    ubuntu open folderpath on terminal
    ubuntu install pae for the 32bit system 4g limited issue
    EXT Designer 正式版延长使用脚本
    用 Vagrant 快速建立開發環境
  • 原文地址:https://www.cnblogs.com/invoid/p/5645307.html
Copyright © 2011-2022 走看看