扩展lucas定理bzoj2142待修改（这是我自己yy的。。。所以看看就好23333）

此题还是很好想的知道威尔逊定理就好了。。。然后自己yy，实在不懂要看证明就看看二潘的初等数论就好啦。。。

#include<iostream>#include<cstdio>#include<algorithm>#include<cmath>using namespace std;typedef long long ll;long long prime[1000], times[1000], value[1000];long long b[1000], m[1000];long long functionnum;long long fact[4000000];long long eulerfunction;long long totprime;long long primetimes[1000], resvalue[1000];long long e;long long thingsamount[200];long long  divide(long long q)//质因数分解{for (ll i = 2; i*i <= q; i++){if (q%i == 0){prime[totprime] = i;value[totprime] = i; times[totprime] = 1; q /= i;while (q%i==0){value[totprime] *= i;times[totprime]++, q /= i;}totprime++;}}if (q>1)value[totprime] = q, prime[totprime] = q, times[totprime] = 1, totprime++;//这句话没加就悲剧了！！！return totprime;}long long gcd(long long a, long long b){if (b == 0)return a;elsereturn gcd(b, a%b);}void extgcd(long long a, long long b, long long &x, long long &y){if (b == 0)x = 1, y = 0;else{extgcd(b, a%b, y, x);y -= a / b*x;}}long long init(long long prime,long long mod,int times,long long *fact)//初始化阶乘模 prime^times的数组函数中mod代表这个值 mod=prime^times;{long long temp = mod / prime;//temp=prime^(times-1)temp = mod - temp;//这就是在求欧拉函数fact[0] = 1;//9 的既约剩余系1 2 4 6 5 7 8 fact[1]=1%9 fact[2]=1*2%9 fact[3]=1*2*4%9........ll top = 1;for (ll i = 1; i <= mod; i++){if (gcd(i, mod) != 1)continue;fact[top] = (fact[top - 1] * i)%mod;top++;if (top == temp+1)break;}//for (int i = 1; i < top; i++)cout << "check:"<<fact[i] << endl;return temp; //返回既约剩余系的大小;}long long getrev(long long value, long long mod){long long revs, temp;extgcd(value, mod, revs, temp);return ((revs%mod) + mod) % mod;}long long mod_fact(ll n, long long mod, long long prime,ll &e){e = 0;if (n == 0)return 1;long long res=mod_fact(n / prime, mod,prime, e);e += n / prime;ll resnum = n - n / prime;//剩下的和mod互素的数且是从小到大有周期的。ll ischange = (resnum / eulerfunction) % 2;//因为p^k既约剩余系的乘积mod p^k是等于-1的威尔逊定理。看二潘数论书有介绍。p151long long other = fact[resnum%eulerfunction];//剩余的哪些数不能构成一个既约剩余系。if (ischange)//这有个bug 估计是oj的数据太low没查出来。。因为2^k的既约剩余系莫它不管怎样都是1只要k>=3就不用变号了。。太懒就不想改了。。。other = (mod - other) % mod;return (other*res) % mod;}//void solve(int thingsnum,int n, int num)//k!=left*(prime^q)此函数就是为了得到left%(prime^t) 和q的  prime^t=要莫的那个数值//{//resvalue[thingsnum] = mod_fact(n, value[num], prime[num], primetimes[thingsnum]);//k!%(prime^q)=left*(prime^r) resvalue[thingsnum]=left%(prime^q);primetimes[thingsnum]=r;//}long long fast_pow(ll times, long long mod,long long prime){long long temp = prime; long long ans = 1;while (times){if (times % 2 == 1)ans = (ans*temp) % mod;//temp = (temp*temp) % mod;//times /= 2;}return ans%mod;//}long long getfunction(ll num){eulerfunction=init(prime[num], value[num], times[num], fact);for (int i = 0; i <= e; i++){resvalue[i] = mod_fact(thingsamount[i], value[num], prime[num], primetimes[i]);}ll timess = primetimes[0];for (int i = 1; i <= e; i++)timess -= primetimes[i];long long primepow = fast_pow(timess, value[num], prime[num]);long long ans = resvalue[0] % value[num];for (int i = 1; i <= e; i++){ans = (ans*getrev(resvalue[i], value[num])) % value[num];}return (ans*primepow) % value[num];}void chinese(long long *b, long long *m,ll tot){long long x = 0; long long mm = 1;for (ll i = 0; i < tot; i++){long long aa = mm; long long bb = b[i] - x; long long kk = gcd(aa, m[i]);//这里一定是有解的所以不用判断;long long tt = bb / kk*getrev(aa / kk, m[i] / kk);tt = (tt % (m[i] / kk) + m[i] / kk) % (m[i] / kk);x = x + mm*tt;mm = mm*m[i] / kk;}x = ((x%mm) + mm) % mm;b[e + 1] = x;}void getanser(long long prime){ll  tot=divide(prime);for (ll i = 0; i < tot; i++)b[i] = getfunction(i),m[i]=value[i];chinese(b, m, tot);printf("%lld\n", b[e + 1]);}int main(){long long p;scanf("%lld", &p);scanf("%lld%lld", &thingsamount[0],&e);ll all = 0;for (int i = 1; i <=e; i++){scanf("%lld", &thingsamount[i]);all += thingsamount[i];}thingsamount[e + 1] = thingsamount[0] - all;e++;if (all <= thingsamount[0])getanser(p);elseprintf("Impossible\n");}

下面这份代码是我简洁后的，代码量小多咯。

注意分清中国剩余定理和线性同余方程组的区别

#include<iostream>#include<cstdio>#include<algorithm>using namespace std;typedef long long ll;ll prime[100], times[200], value[200];ll fac[120500];ll p, n, m;ll thingsamount[100];int divide(ll n){int tot = 0;for (ll i = 2; i*i <= n; i++){if (n%i == 0){prime[tot] = i; times[tot] = 1; value[tot] = i; n /= i;while (n%i == 0)value[tot] *= i, times[tot]++, n /= i;tot++;}}if (n > 1)prime[tot] = n, value[tot] = n, times[tot] = 1,tot++;return tot;}ll pow(ll prime, ll times, ll mod){ll value = 1; ll temp = prime;while (times){if (times & 1)value *= temp,value%=mod;temp *= temp; times >>= 1;temp %= mod;}return value;}ll rev(ll a, ll mod,ll prime){ll ruler = (mod / prime)*(prime - 1);return pow(a, ruler-1, mod);}ll cal(ll n, ll prime, ll t, ll value, ll &e,ll euler){if (n < prime){e = 0;return fac[n];}ll res = cal(n / prime, prime, t, value, e,euler);e += n / prime; n -= n / prime;ll anotherres = n%euler; anotherres = fac[anotherres];bool ischange = (n / euler) % 2;//为什么2会比较特殊?因为x mod m=a1 推出x mod m1=a1但是x mod m1=a1退不出x mod m=a1因为有可能解很多if (ischange && (prime != 2 || (prime == 2 && t < 3)))anotherres = value -anotherres;return (anotherres*res) % value;}ll gcd(ll a, ll b){if (b == 0)return a;gcd(b, a%b);}ll solve(ll prime, ll t, ll value){ll ans1,ans2=1,e1,e2=0,tempe;ll euler = value / prime*(prime - 1);int tot = 1; fac[1] = 1; fac[0] = 1;for (int i = 2; tot < euler; i++)if (gcd(i, value) == 1)fac[++tot] = (fac[tot - 1] * i) % value;ans1 = cal(thingsamount[0], prime, t, value, e1,euler);for (int i = 1; i <= m + 1; i++)ans2 *= cal(thingsamount[i], prime,t, value, tempe, euler),e2+=tempe,ans2%=value;ll ans3 = pow(prime, e1 - e2, value);ans2 = pow(ans2, euler - 1, value);return ((ans2*ans1) % value*ans3) % value;}int main(){scanf("%lld%lld%lld", &p, &n, &m);thingsamount[0] = n; int allr = 0;for (int i = 1; i <= m; i++)scanf("%lld", &thingsamount[i]),allr+=thingsamount[i];allr = n - allr;if (allr < 0){printf("Impossible\n");return 0;}thingsamount[m + 1] = allr;int all = divide(p);ll ans = 0;for (int i = 0; i < all; i++){ll M = p / value[i];ll temp = solve(prime[i], times[i], value[i]);ans += (((temp*M) % p)*(rev(M, value[i], prime[i])%p))%p;}printf("%lld\n", ans%p);return 0;}