• 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;    
    }
  • 相关阅读:
    servlet上传图片 服务器路径(转)
    图片和提交servlet的相对和绝对路径
    Intel 的面试经历中国研究院
    CentOS-6.5-x86_64 最小化安装,已安装包的总数,这些包?
    西门子PLC学习笔记8-(计时器)
    这个周末我太累了
    windows7股票的,win8残疾人,安装Han澳大利亚sinoxn个时间,sinox它支持大多数windows软体
    net.sf.json 迄今 时刻 格式 办法
    ar命令提取.a时刻,一个错误 is a fat file (use libtool(1) or lipo(1) and ar(1) on it)
    POJ 2187: Beauty Contest(旋转卡)
  • 原文地址:https://www.cnblogs.com/invoid/p/5645307.html
Copyright © 2020-2023  润新知