BZOJ 2142 礼物

Posted on 2017-11-02 16:54  ziliuziliu  阅读(152)  评论(0编辑  收藏  举报

扩展lucas(作为noip模拟赛)

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#define maxn 100050
#define maxm 15
using namespace std;
long long p,n,m,tot=0,w[maxm],prime[maxn],d1[maxm],d2[maxm],top=0,cnt=0,phi[maxn];
long long ans[maxm],num=0;
bool vis[maxn];
void get_table()
{
    for (long long i=2;i<=100000;i++)
    {
        if (!vis[i]) {phi[i]=i-1;vis[i]=true;prime[++tot]=i;}
        for (long long j=1;j<=tot && i*prime[j]<=100000;j++)
        {
            vis[i*prime[j]]=true;
            if (i%prime[j]) {phi[i*prime[j]]=phi[i]*prime[j];continue;}
            else {phi[i*prime[j]]=phi[i]*(prime[j]-1);break;}
        }
    }
}
void divide(long long p)
{
    for (long long i=1;i<=tot;i++)
    {
        if (p%prime[i]==0)
        {
            d1[++cnt]=prime[i];d2[cnt]=1;
            while (p%prime[i]==0) {d2[cnt]*=prime[i];p/=prime[i];}
        }
        if (p==1) break;
    }
}
long long f_pow(long long a,long long b,long long p)
{
    if (!b) return 1;
    long long ans=1;
    while (b)
    {
        if (b&1) {ans*=a;ans%=p;}
        a*=a;a%=p;
        b>>=1;
    }
    return ans;
}
long long get_ans(long long x,long long p,long long pk,long long f)
{
    if (x<=1) return 1;
    num+=x/p*f;long long ret=1;
    if (x>=pk)
    {
        for (long long i=1;i<=pk;i++) 
        {
            if (i%p) {ret*=i;ret%=pk;}
            else continue;
        }
        ret=f_pow(ret,x/pk,pk);
    }
    for (long long i=(x/pk)*pk+1;i<=x;i++)
    {
        if (i%p) {ret*=i%pk;ret%=pk;}
        else continue;
    }
    return ret*get_ans(x/p,p,pk,f)%pk;
}
void exgcd(long long a,long long b,long long &x,long long &y)
{
    if (!b) {x=1;y=0;return;}
    exgcd(b,a%b,x,y);
    long long t=x-(a/b)*y;
    x=y;y=t;return;
}
long long inv(long long x,long long p)
{
    long long r1,r2;
    exgcd(x,p,r1,r2);
    return (r1%p+p)%p;
}
void CRT()
{
    long long ret=0;
    for (long long i=1;i<=cnt;i++)
    {
        ret+=ans[i]%p*(p/d2[i])%p*inv(p/d2[i],d2[i])%p;
        ret%=p;
    }
    printf("%lld\n",ret);
}
int main()
{
    scanf("%lld",&p);
    scanf("%lld%lld",&n,&m);
    for (long long i=1;i<=m;i++) {scanf("%lld",&w[i]);tot+=w[i];}
    if (n-tot<0) {printf("Impossible\n");return 0;}
    if (n-tot) w[++m]=n-tot;
    tot=0;get_table();
    divide(p);
    for (long long i=1;i<=cnt;i++)
    {
        num=0;ans[i]=get_ans(n,d1[i],d2[i],1)%d2[i];
        for (long long j=1;j<=m;j++) 
        {
            ans[i]*=inv(get_ans(w[j],d1[i],d2[i],-1)%d2[i],d2[i]);
            ans[i]%=d2[i];
        }
        if (num) ans[i]*=f_pow(d1[i],num,d2[i]);ans[i]%=d2[i]; 
    }
    CRT();
    return 0;
}