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

最新推荐文章于 2021-09-13 17:14:16 发布

guoshiyuan484

最新推荐文章于 2021-09-13 17:14:16 发布

阅读量355

点赞数

分类专栏：数论

本文链接：https://blog.csdn.net/guoshiyuan484/article/details/76595584

版权

数论专栏收录该内容

6 篇文章 0 订阅

订阅专栏

此题还是很好想的知道威尔逊定理就好了。。。然后自己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;
	else
		return 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的威尔逊定理。看二潘数论书有介绍。p151
	long 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);
	else
		printf("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;
}