bzoj 3283 扩展BSGS + 快速阶乘

T2  扩展BSGS

T3 快速阶乘

给定整数n，质数p和正整数c，求整数s和b，满足n! / p^b = s mod p^c

考虑每次取出floor(n/p)个p因子，然后将问题转化为子问题。

/**************************************************************
    Problem: 3283
    User: idy002
    Language: C++
    Result: Accepted
    Time:1704 ms
    Memory:12380 kb
****************************************************************/
 
#include <cstdio>
#include <cstring>
#include <cmath>
 
typedef long long dnt;
void exgcd( dnt a, dnt b, dnt &d, dnt &x, dnt &y ) {
    if( b==0 ) {
        x=1, y=0, d=a;
    } else {
        exgcd(b,a%b,d,y,x);
        y-=a/b*x;
    }
}
dnt inv( int a, int n ) {
    dnt d, x, y;
    exgcd(a,n,d,x,y);
    return d==1 ? (x%n+n)%n : -1;
}
dnt mpow( dnt a, dnt b, dnt c ) {
    dnt rt;
    for( rt=1; b; b>>=1,a=a*a%c )
        if( b&1 ) rt=rt*a%c;
    return rt;
}
 
namespace Task1 {
    void sov( dnt a, dnt b, dnt c ) {
        printf( "%lld
", mpow(a,b,c) );
    }
};
namespace Task2 {
    const int mod = 38281;
    const int len = mod<<3;
    struct Hash {
        int head[mod], val[len], rat[len], next[len], ntot;
        void init() { 
            ntot=0; 
            memset( head, 0, sizeof(head) );
        }
        void insert( int v, int r ) {
            int k = v % mod;
            ntot++;
            next[ntot] = head[k];
            val[ntot] = v;
            rat[ntot] = r;
            head[k] = ntot;
        }
        int query( int v ) {
            int k = v % mod;
            for( int t=head[k]; t; t=next[t] ) 
                if( val[t]==v ) return rat[t];
            return -1;
        }
    }hash;
     
    dnt gcd( dnt a, dnt b ) {
        return b ? gcd(b,a%b) : a;
    }
    int bsgs( dnt s, dnt a, dnt b, dnt c ) {
        hash.init();
        int m = ceil(sqrt(c));
        for( int i=0; i<m; i++ ) {
            if( s==b ) return i;
            hash.insert( s, i );
            s = s*a % c;
        }
        dnt am = 1;
        for( int i=0; i<m; i++ )
            am = am*a % c;
        am = inv(am,c);
        b = b*am % c;
        for( int i=m; i<c; i+=m ) {
            int j = hash.query( b );
            if( j!=-1 ) return i+j;
            b = b*am % c;
        }
        return -1;
    }
    int exbsgs( dnt a, dnt b, dnt c ) {
        dnt s = 1;
        for( int i=0; i<32; i++ ) {
            if( s==b ) return i;
            s = s*a % c;
        }
        dnt cd;
        s = 1;
        int rt = 0;
        while( (cd=gcd(a,c))!=1 ) {
            rt++;
            s*=a/cd;
            if( b%cd ) return -1;
            b/=cd;
            c/=cd;
            s%=c;
        }
        int p = bsgs(s,a,b,c);
        if( p==-1 ) return -1;
        return rt + p;
    }
    void sov( int a, int b, int c ) {
        int res = exbsgs(a,b,c);
        if( res==-1 ) printf( "Math Error
" );
        else printf( "%d
", res );
    }
};
namespace Task3 {
    struct Pair {
        int s, k;
        Pair( int s, int k ):s(s),k(k){}
    };
    dnt aa[50], mm[50];
    dnt pres[1001000];
    dnt pp[50], cc[50], ppp[50], tot;
 
    dnt china( int n, dnt *a, dnt *m ) {
        int M=1;
        for( int i=0; i<n; i++ ) 
            M *= m[i];
        int rt = 0;
        for( int i=0; i<n; i++ ) {
            dnt Mi = M/m[i];
            rt = (rt+Mi*inv(Mi,m[i])*a[i]) % M;
        }
        return rt;
    }
    void init( int p, int pp ) {
        pres[0] = 1;
        for( int i=1; i<=pp; i++ ) {
            if( i%p==0 ) {
                pres[i] = pres[i-1];
            } else {
                pres[i] = pres[i-1]*i % pp;
            }
        }
    }
    Pair split( int n, int p, int c, int pp ) {
        int b = n/p;
        if( b==0 ) {
            return Pair( pres[n], 0 );
        } else {
            Pair pr = split( b, p, c, pp );
            return Pair( (pr.s*pres[n%pp]%pp) * mpow(pres[pp],n/pp,pp) % pp, pr.k+b );
        }
    }
    void sov( int m, int n, int c ) {
        tot = 0;
        for( int i=2; i*i<=c; i++ ) {
            if( c%i==0 ) {
                pp[tot] = i;
                cc[tot] = 0;
                ppp[tot] = 1;
                while( c%i==0 ) {
                    cc[tot]++;
                    ppp[tot] *= pp[tot];
                    c/=i;
                }
                tot++;
            }
        }
        if( c!=1 ) {
            pp[tot] = c;
            cc[tot] = 1;
            ppp[tot] = c;
            tot++;
            c = 1;
        }
        for( int i=0; i<tot; i++ ) {
            init(pp[i],ppp[i]);
            Pair pn = split( n, pp[i], cc[i], ppp[i] );
            Pair pa = split( m, pp[i], cc[i], ppp[i] );
            Pair pb = split( n-m, pp[i], cc[i], ppp[i] );
            if( pn.k-pa.k-pb.k >= cc[i] ) {
                aa[i] = 0;
                mm[i] = ppp[i];
            } else {
                aa[i] = pn.s * (inv(pa.s,ppp[i])*inv(pb.s,ppp[i])%ppp[i]) % ppp[i];
                for( int j=0; j<pn.k-pa.k-pb.k; j++ )
                    aa[i] = (dnt) aa[i]*pp[i] % ppp[i];
                mm[i] = ppp[i];
            }
        }
/*
        fprintf( stderr, "tot=%d
", tot );
        for( int i=0; i<tot; i++ )
            fprintf( stderr, "%d %d
", aa[i], mm[i] );
*/
        printf( "%lld
", china(tot,aa,mm) );
    }
};
 
int main() {
    int n;
    scanf( "%d", &n );
    for( int i=1,opt,y,z,p; i<=n; i++ ) {
        scanf( "%d%d%d%d", &opt, &y, &z, &p );
        if( opt==1 )
            Task1::sov( y, z, p );
        else if( opt==2 ) 
            Task2::sov( y, z, p );
        else
            Task3::sov( y, z, p );
    }
}