图个安心的数论类简单模板

Posted

tags:

篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了图个安心的数论类简单模板相关的知识,希望对你有一定的参考价值。

/***                              ****
    数论模板
    1.EXGCD
    2.CRT (互质 与 非互质)
    3.逆元(线性预处理 欧拉(费小)定理)
    4.筛 (线性筛与数论函数 状压筛)
    5.大素数计数 (Meisell-Lehmer O2/3 O3/4 )
    6.SG函数 mex方法
    7.Gauss消元
    8.FFT NTT FWT
    9.二次剩余  (二次域方法)
    10.素性测试 (Miller_rabin)
    11.LUCAS 大组合取模 
/**                               ****/
//1.EXGCD
//
//用于解 ax + by = c mod n x与y的值
//前提: (a,b) = d | c
//先求 ax‘+by‘=(a,b) 得特解x‘ y‘
//其中一组解
//x0 = x‘·c/d
//y0 = y‘·c/d
//通解:
//x = (x‘+ k·b/d)c/d
//y = (y‘- k·a/d)c/d
//最后最小 x=(x‘ mod b/d) c/d
LL exgcd(LL a, LL b, LL &x, LL &y)
{
	LL d = a;
	if(b == 0)
		x = 1, y = 0;
	else
	{
		d = exgcd(b, a % b, y, x);
		y -= (a / b) * x;
	}
	return d;
}
/********************************************************/
//2.CRT
//
//用于求同余方程的 x ≡ r[i] mod m[i] x的值
//可用于化简 模数较大但可分解成多个小模数的公式,通过列出同余方程得到最终的X值最后取模
//
//(1)模数互质
//设M = PI(i:1->n){ mod[i] } 即模数之积
//定义M[i] = M/mod[i] inv[i] = mod[i]^-1
//则X=(r[1]·M[1]·inv[1]) mod M   实则M除以mod[i]的平方.
//param1:同余方程数量
//param2:余数
//param3:模数
LL CRT(LL n, LL rem[], LL mod[])
{
	LL M = 1, x, y;
	for(int i = 0; i < n; i++)
		M *= mod[i];
	LL res = 0;
	for(int i = 0; i < n; i++)
	{
		LL t = M / mod[i];
		exgcd(t, mod[i], x, y);//EXGCD求mod[i]逆元
		res = (res + t * rem[i] * x) % M;
	}
	return (res % M + M) % M;
}
//(2)模数非互质 (合并方程)
//利用EXGCD性质 m1·x + m2·y = r2 - r1 求出特解x‘ y‘
//再求出x最小解 x = x‘ * (c/g) mod (b/d)
//X ≡ (m1 * x + r1) mod (LCM(m1, m2))
//当前余数 模数 下个传入的余数 模数
LL crtCD(LL &rem, LL &mod, LL newRem, LL newMod)
{
	LL c = newRem - rem;// c= r2 -r1
	LL g = __gcd(mod, newMod);
	if(c % g!=0)
		return 0;
	LL x, y;//特解x‘ y‘
	exgcd(mod, newMod, x, y);
	LL t = newMod / g;// b/d
	LL minx = (c / g * x % t + t) % t;//x = x‘ * (c/g) mod (b/d) 最小x
	rem = minx * mod + rem;// newrem = mod * x + rem,
	mod = mod / g * newMod; //newmod = lcm(mod, newMod)
	return 1;
}
/********************************************************/
//3.逆元
//
//大多线性预处理用于组合数
void init(){
    inv[1] = inv[0] = 1;
    for(int i = 1; i < N; i++)
        inv[i] = (mod - mod / i) * inv[mod % i];
}
//单次求解
//等价方程 a/b % m = a mod (m·b)/b
LL Inv(LL a, LL b)
{
    LL d = __gcd(a, b);
    if(d != 1)
        return -1;
    LL x, y;
    extend_Euclid(a, b, x, y);
    return (x % b + b) % b;
}
//费小定理 略
/********************************************************/
//4.筛
//
//(1)线性筛
//常用于筛出数论函数
LL pri[N];
LL mu[N];
LL phi[N];
int c = 0;
bool vis[N];
void prime()
{
	MMF(vis);
	MMF(sum);
	mu[1] = 1;
	for(int i = 2; i < N; i++) {
		if(!vis[i])
			pri[c++] = i,
			mu[i] = -1,
			phi[i] = i - 1;
		for(int j = 0; j < c && i * pri[j] < N; j++)  {
			vis[i * pri[j]] = 1;
			if(i % pri[j] == 0){//倍数i=kp
			    phi[i* pri[j]] = 1;
				mu[i * pri[j]] = 0;
				break;
			}
			else
                phi[i* pri[j]] = phi[i] * (pri[j] - 1);
                mu[i * pri[j]] = -mu[i];
		}
	}
}
//(2)状压素数筛 不知道拿来干嘛用。
int pri[5800000];
unsigned int sum[5800000];
int vis[N/32+10];
int c = 0;
void prime()
{
    pri[0] = 2;
    sum[0] = 2;
    c = 1;
    for(int i = 3; i < 100000007; i+=2)
    {
        if(!(vis[i/32] & (1 << ((i/2)%16))) )
        {
            pri[c++] = i;
            for(int j = 3*i; j < 100000007; j+=2*i)
                vis[j/32] |= (1 << ((j/2)%16));
        }
    }
}
/********************************************************/
//5.大素数计数
//
//用于论文题
//(1) Meisell-Lehmer
const int N = 5e6 + 2;
bool np[N];
int prime[N], pi[N];
int getprime()
{
    int cnt = 0;
    np[0] = np[1] = true;
    pi[0] = pi[1] = 0;
    for(int i = 2; i < N; ++i)
    {
        if(!np[i]) prime[++cnt] = i;
        pi[i] = cnt;
        for(int j = 1; j <= cnt && i * prime[j] < N; ++j)
        {
            np[i * prime[j]] = true;
            if(i % prime[j] == 0)   break;
        }
    }
    return cnt;
}
const int M = 7;
const int PM = 2 * 3 * 5 * 7 * 11 * 13 * 17;
int phi[PM + 1][M + 1], sz[M + 1];
void init()
{
    getprime();
    sz[0] = 1;
    for(int i = 0; i <= PM; ++i)  phi[i][0] = i;
    for(int i = 1; i <= M; ++i)
    {
        sz[i] = prime[i] * sz[i - 1];
        for(int j = 1; j <= PM; ++j) phi[j][i] = phi[j][i - 1] - phi[j / prime[i]][i - 1];
    }
}
int sqrt2(LL x)
{
    LL r = (LL)sqrt(x - 0.1);
    while(r * r <= x)   ++r;
    return int(r - 1);
}
int sqrt3(LL x)
{
    LL r = (LL)cbrt(x - 0.1);
    while(r * r * r <= x)   ++r;
    return int(r - 1);
}
LL getphi(LL x, int s)
{
    if(s == 0)  return x;
    if(s <= M)  return phi[x % sz[s]][s] + (x / sz[s]) * phi[sz[s]][s];
    if(x <= prime[s]*prime[s])   return pi[x] - s + 1;
    if(x <= prime[s]*prime[s]*prime[s] && x < N)
    {
        int s2x = pi[sqrt2(x)];
        LL ans = pi[x] - (s2x + s - 2) * (s2x - s + 1) / 2;
        for(int i = s + 1; i <= s2x; ++i) ans += pi[x / prime[i]];
        return ans;
    }
    return getphi(x, s - 1) - getphi(x / prime[s], s - 1);
}
LL getpi(LL x)
{
    if(x < N)   return pi[x];
    LL ans = getphi(x, pi[sqrt3(x)]) + pi[sqrt3(x)] - 1;
    for(int i = pi[sqrt3(x)] + 1, ed = pi[sqrt2(x)]; i <= ed; ++i) ans -= getpi(x / prime[i]) - i + 1;
    return ans;
}
LL lehmer_pi(LL x)//答案
{
    if(x < N)   return pi[x];
    int a = (int)lehmer_pi(sqrt2(sqrt2(x)));int bit[64];
int a[N], cnt;
void gauss()
{
	cnt = 1;
	for(int k = 30; k >= 0; k--)
	{
		for(int i = cnt; i <= n; i++)
		{
			if(a[i] & (1 << k))
			{
				for(int j = 1; j <= n; j++)
				{
					if(j != i && (a[j] & (1 << k)))
						a[j] ^= a[i];
				}
				swap(a[i], a[cnt]);//因为存的是从高到低 放到前面
				cnt++;
				break;
			}
		}
	}
	cnt--;
	/*for(int i = 1; i <= n; i++)
		printf("%d", a[i]);
	cout << endl;
	cout << "~" << cnt << endl;*/
}
    int b = (int)lehmer_pi(sqrt2(x));
    int c = (int)lehmer_pi(sqrt3(x));
    LL sum = getphi(x, a) +(LL)(b + a - 2) * (b - a + 1) / 2;
    for (int i = a + 1; i <= b; i++)
    {
        LL w = x / prime[i];
        sum -= lehmer_pi(w);
        if (i > c) continue;
        LL lim = lehmer_pi(sqrt2(w));
        for (int j = i; j <= lim; j++) sum -= lehmer_pi(w / prime[j]) - (j - 1);
    }
    return sum;
}
//O 3/4 略
/********************************************************/
//6.SG
//
//用于组合游戏SG局面值求解
int sg[N];
int dfs(int x)
{
	if(x < 0)
		return 0;
	if(sg[x] != -1)
		return sg[x];
	int vis[2010];//注意函数内声明
	MMF(vis);//注意清空
	for(int i = 1; i <= x; i++)
		vis[dfs(i - 3) ^ dfs(x - i - 2)] = 1;
	for(int i = 0; ;i++)//后继局面最小未出现
		if(!vis[i])
			return sg[x] = i;
	return sg[x];
}
/********************************************************/
//7.Gauss
//
//线代用于求多元方程组解
//求线性基等
int bit[64];
int a[N], cnt;
void gauss()
{
	cnt = 1;
	for(int k = 30; k >= 0; k--)
	{
		for(int i = cnt; i <= n; i++)
		{
			if(a[i] & (1 << k))//消元部分
			{
				for(int j = 1; j <= n; j++)
				{
					if(j != i && (a[j] & (1 << k)))
						a[j] ^= a[i];
				}
				swap(a[i], a[cnt]);//因为存的是从高到低 放到前面
				cnt++;
				break;
			}
		}
	}
	cnt--;
	/*for(int i = 1; i <= n; i++)
		printf("%d", a[i]);
	cout << endl;*/
}
/********************************************************/
//8.快速傅里叶XX变化
//
//高精度 快速符号卷积运算
//(1) FFT
const double eps = 1e-8;
const double pi = acos(-1.0);

struct Complex
{
	double a, b;
	Complex(){}
	Complex(double _a, double _b):a(_a),b(_b){}
	Complex(double _a):a(_a),b(0.0){}
	inline Complex operator +(const Complex &z)const{
		return Complex(a + z.a, b + z.b);
	}
	inline Complex operator -(const Complex &z)const{
		return Complex(a - z.a, b - z.b);
	}
	inline Complex operator *(const Complex &z)const{
		return Complex(a * z.a - b * z.b, a * z.b + b * z.a);
	}
	inline Complex operator /(const Complex &z)const{
		double m = z.a * z.a + z.b * z.b;
		return Complex((a * z.a + b * z.b) / m, (z.a * b - z.b * a) / m);
	}
}A[N], B[N];
int L;
int rev[N];
int ans[N];
char a[N], b[N];

void init()
{
	MMF(A);
	MMF(B);
	MMF(ans);
	L = 0;//?
}

int reverse(int x,int r)  //蝴蝶操作
{
    int ans=0;
    for(int i=0; i<r; i++){
        if(x&(1<<i)){
            ans+=1<<(r-i-1);
        }
    }
    return ans;
}

void rader(LL f[], int len)//雷德二进制反转
{
	for(int i = 1, j = len >> 1; i < len - 1; i++)
	{
		if(i < j) swap(f[i], f[j]);
		int k = len >> 1;
		while(j >= k)
		{
			j -= k;
			k >>= 1;
		}
		if(j < k) j += k;
	}
}
void FFT(Complex c[], int nlen, int on) //主要部分 on为正逆开关
{
	Complex wn, w, x, y;
	for(int i = 0; i < nlen; i++)
		if(i < rev[i]) swap(c[i], c[rev[i]]);
	for(int i = 1; i < nlen; i <<= 1)
	{
		wn = Complex(cos(pi/i), sin(pi/i) * on);
		for(int p = i << 1, j = 0; j < nlen; j+= p)
		{
			w = Complex(1, 0);
			for(int k = 0; k < i; k++, w = w * wn)
			{
				x = c[j + k];
				y = w * c[j + k + i];
				c[j + k] = x + y;
				c[j + k + i] = x - y;
			}
		}
	}
	if(on == -1)
		for(int i = 0; i < nlen; i++)
			c[i].a /= (double)nlen;
}

void work(Complex a[], Complex b[], int len)
{
	FFT(a, len, 1); //
	FFT(b, len, 1); //
	for(int i = 0; i < len; i++)
		a[i] = a[i] * b[i];
	FFT(a, len, -1);
}
int main()
{
	while(~scanf("%s%s", a, b))
	{
		init();
		int alen = strlen(a);
		int blen = strlen(b);
		int len = alen + blen;
		//getlen
		int n;
		for(n = 1; n < len - 1; n <<= 1, L++);
		//rev
		for(int i = 0; i < n; i++)
			rev[i] = (rev[i>>1] >> 1) | ((i & 1) << (L - 1));
		//deal
		for(int i = 0; i < alen; i++)
			A[i] = a[alen - i - 1] - ‘0‘;
		for(int i = 0; i < blen; i++)
			B[i] = b[blen - i - 1] - ‘0‘;
		work(A, B, n);
		///
		for(int i = 0; i <= len; i++)
			ans[i] = (int)(A[i].a + 0.5);
		for(int i = 0; i <= len; i++)
			ans[i + 1] += ans[i] / 10, ans[i] %= 10;
		while(ans[len] == 0 && len)
			len--;
		for(int i = len; ~i; i--)
			printf("%c", ans[i] + ‘0‘);
		printf("\n");
	}
    return 0;
}
//(2)NTT
//适用于有取模运算 模数为费马素数 p = c·2^k + 1
//HDU 6061
const LL N = 1e5 + 10;
const LL mod = 998244353;
const int g = 3;//原根注意
const int maxlen = 1 << 18;


LL wn[maxlen], fac[N], inv[N];
LL fpow(LL a, LL n)
{
	LL ans = 1;
	while(n)
	{
		if(n & 1)
			ans = (ans * a % mod + mod) %mod;
		a = (a * a + mod) % mod;
		n >>= 1;
	}
	return ans;
}

void init()
{
	wn[0] = 1;
	wn[1] = fpow(g, ((mod - 1)>>18));
	for(int i = 2; i < maxlen; i++)
		wn[i] = wn[i - 1] * wn[1] % mod;
	fac[1] = inv[1] = 1;
	fac[0] = inv[0] = 1;
	for(int i = 2; i < N; i++)
	{
		fac[i] = fac[i - 1] * i % mod;
		inv[i] = (mod - mod / i) * inv[mod % i] % mod;
	}
	for(int i = 1; i < N; i++)
		(inv[i] *= inv[i - 1]) %= mod;
}

void rader(LL f[], int len)
{
	for(int i = 1, j = len >> 1; i < len - 1; i++)
	{
		if(i < j) swap(f[i], f[j]);
		int k = len >> 1;
		while(j >= k)
		{
			j -= k;
			k >>= 1;
		}
		if(j < k) j += k;
	}
}
void ntt(LL f[], int len, int on)
{
	/*for(int i = 0, j = 0; i < len;i++)
	{
		if(i > j) swap(f[i], f[j]);
		for(int l = len >> 1; (j^=l) < l; l>>=1);
	}*/
	rader(f, len);
	for(int i = 1, d = 1; d < len; i++, d <<= 1)
	{
		//LL wnn = fpow(g, (mod-1)/(d<<1));
		for(int j = 0; j < len; j += (d << 1))
		{
			//LL w = 1;
			for(int k = 0; k < d; k++)
			{
				LL t = wn[(maxlen >> i) * k] * f[j + k + d] % mod;
				//LL t = w*f[j+k+d]%mod;
				//w = w*wnn % mod;
				f[j + k + d] = ((f[j + k] - t) % mod + mod) % mod;
				f[j + k] = ((f[j + k] + t) % mod + mod) % mod;
			}
		}
	}
	if(on == -1)
	{
		reverse(f + 1, f + len);
		LL inv2 = fpow(len, mod - 2);
		for(int i = 0; i < len; i++)
			f[i] = f[i] % mod * inv2 % mod;
	}
}

void work(LL a[], LL b[], int len)
{
	ntt(a, len, 1);
	ntt(b, len, 1);
	for(int i = 0; i < len; i++)
		a[i] = (a[i] * b[i] % mod + mod) % mod;
	ntt(a, len, -1);
}

LL A[maxlen], B[maxlen], Suma, c[N];
int n, m;
int main()
{
	init();
	while(~scanf("%d", &n))
	{
		for(int i = 0; i <= n; i++) scanf("%lld", c + i);
		scanf("%d", &m);
		Suma = 0;
		LL t;
		for(int i = 0; i < m; i++)
			scanf("%lld", &t), Suma -= t;
		Suma = (Suma + mod) % mod;
		while(Suma < 0)
			Suma += mod;
		if(Suma == 0)
		{
			for(int i = 0; i <= n; i++)
				printf("%lld ", c[i]);
			printf("\n");
			continue;
		}
		//getLength
		int len = 1;
		while(len <= 2*n) len <<= 1;
		//
		LL ae = 1;
		for(int i = 0; i < len; i++)
		{
			if(i <= n)
			{
				B[i] = ae * inv[i] % mod;
				A[i] =(fac[n - i] * c[n - i]) % mod;
			}
			else A[i] = B[i] = 0;
			ae = (ae * Suma) % mod;
			//cout << A[i] << "~"<< B[i] << endl;
		}
		work(A, B, len);
		for(int i = 0; i <= n; i++)
			A[i] = A[i] * inv[n - i] % mod;
		for(int i = 0; i <= n; i++)
			printf("%lld ", A[n - i]);
		printf("\n");
	}
}
//(3)FWT 暂略
//用于快速计算位间操作符号的卷积
/********************************************************/
//9.二次剩余解
//
//常用于数论意义下的化简式子
//可用于复杂方程中替代有多重开根的值(一般为平方)
//因此需要求多重剩余的解x
//
//二次复数域方法求解

LL mod;
//二次域乘法????

//O1乘法取模黑科技
LL mul(LL x,LL y)
{
    return (x * y-(LL)(x /(long double)mod * y + 1e-3) * mod + mod) % mod;
}

LL fpow(LL a, LL n)
{
	LL ans = 1;
	while(n)
	{
		if(n & 1)
			//ans = (a * ans + mod)% mod;
			ans = mul(ans, a);
		n >>= 1;
		//a = (a * a % mod + mod) % mod;
		a = mul(a, a);
	}
	return ans;
}


LL QuadraticResidue()
{
	LL rem = (-3 % mod + mod) % mod;
	if(rem == 0)//特判mod==3
		return 0;
	if(mod == 2)//特判非奇素数
		return 1;
	if(fpow(rem, (mod - 1)/2) + 1 == mod)//欧拉判别条件 非剩余
		return -1;
	LL a = -1, b;
	while(1)
	{
		a = (1 + rand() % (mod - 1));//注意模板是错的 原来a可能取到0
		b = (a * a - rem) % mod;
		while(b < 0)
			b += mod;
		if(fpow(a, (mod-1)/2) + 1 == mod)//找一个非剩余的b 推导上一阶剩余
			break;
	}
	a = fpow(b, (mod + 1)/2);
	return a;
}

LL QuadraticResidue2()//求解二次剩余根 x
{
	LL rem = ((-3) % mod + mod) % mod;
	if(rem == 0)//特判mod==3
		return 0;
	if(mod == 2)//特判非奇素数
		return rem % mod;
	if(fpow(rem, (mod - 1)/2) != 1)//欧拉判别条件 非剩余
		return -1;
	LL x = -1, b;
	LL i, k;
	if(mod % 4 == 3)
		x = fpow(rem, (mod + 1)/4);
	else
    {
        for(b = 1; fpow(b,(mod - 1)/2) == 1;
        b = 1 + rand() % (mod - 1));

        i = (mod - 1) / 2;
        k = 0;
        do
        {
            i /= 2,k /= 2;
            if(( mul(fpow(rem, i),fpow(b, k)) + 1)% mod==0)
                k += (mod - 1)/2;
        }
        while(i % 2 == 0);
        x = mul(fpow(rem,(i+1)/2),fpow(b,k/2));
    }
    if(x * 2 > mod)
        x = mod - x;
	return x;
}

vectora;
int main()
{
	int T;
	cin >> T;
	while(T--)
	{
		a.clear();
		LL n, p;
		scanf("%lld%lld", &n, &mod);
		for(int i = 0; i < n; i++)
		{
			LL t;
			scanf("%lld", &t);
			if(t > 0)//注意有0的...
				a.PB(t);
		}
		LL cnt = a.size();
		sort(a.begin(), a.end());

		///////////
		LL ans = 0;
		if(mod == 2)//特殊情况
			ans = cnt * (cnt - 1) / 2LL;
		else
		{
			LL t = QuadraticResidue2();
			if(t == -1)
			{
				printf("0\n");
				continue;
			}
			LL inv = (mod + 1) >> 1;
			LL x = mul((mod + t - 1LL)%mod, inv);
			LL y = mul((mod - t - 1LL)%mod, inv);
			for(int i = 0; i < cnt; i++)
			{
				LL tmp = mul(x , a[i]);
				ans += upper_bound(a.begin(), a.begin() + i, tmp)
				- lower_bound(a.begin(), a.begin() + i, tmp);
			}
			if(x != y)//两个解
			{
				for(int i = 0; i < cnt; i++)
				{
					LL tmp = mul(y, a[i]);
					ans += upper_bound(a.begin(), a.begin() + i, tmp)
					- lower_bound(a.begin(), a.begin() + i, tmp);
				}
			}
		}
		printf("%lld\n", ans);
	}
    return 0;
}
/********************************************************/
//10.Miller-rabin 素性测试
//
//40个素数足以判别出64位整数
bool Miller_Rabbin(int n,int a)//米勒拉宾素数测试
{
    int r = 0, s = n - 1;
    if(!(n % a))
        return false;
    while(!(s & 1))
    {
        s>>=1;
        r++;
    }
    LL k = fpow(a, s, n);
    if(k == 1)
        return true;
    for(int j = 0; j < r; j++, k = k * k % n)
        if(k == n-1)
            return true;
    return false;
}
bool IsPrime(int n)//判断是否是素数
{
    int tab[14] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41};
    for(int i = 0; i < 13; i++)
    {
        if(n == tab[i])
            return true;
        if(!Miller_Rabbin(n, tab[i]))
            return false;
    }
    return true;
}
/********************************************************/
//11.Lucas
//
//
LL C(LL n, LL m, LL mod)
{
	if(m > n)
        return 0;
    LL ans = 0;
    ans = ((fac[n] * inv[m] % mod)* inv[n - m]) % mod;
    return ans;
}
LL lucas(LL n, LL m, LL mod)
{
	if(m == 0)
		return 1;
	return C(n % mod, m % mod, mod) * lucas(n / mod, m / mod, mod) % mod;
}


以上是关于图个安心的数论类简单模板的主要内容,如果未能解决你的问题,请参考以下文章

[数论+模板] 快速幂及快速幂求逆元算法模板(模板)

数论——快速幂,模运算及快速幂求逆元

数论Lucas定理

算法分级

简单数论(poj 1152)

P5091 模板扩展欧拉定理 题解