求助 CRT
  • 板块灌水区
  • 楼主haochengw920
  • 当前回复3
  • 已保存回复3
  • 发布时间2024/5/29 08:53
  • 上次更新2024/5/29 15:17:31
查看原帖
求助 CRT
563650
haochengw920楼主2024/5/29 08:53

rt , 我在写 3 模 ntt 的时候尝试了几种不同写法的 crt ,但都无一例外地输出了一个错误的值。

经过调试,似乎只有 Poly::Get() 函数出了问题,求助我的写法有什么问题,或者该怎么改才能对吗

#include<cstdio>
#include<algorithm>
#define int long long
using namespace std;

inline int fsp(int x, int p, int Mod) {
	int res = 1; x %= Mod;
	for (; p; p >>= 1) {
		if (p & 1) res = res * x % Mod;
		x = x * x % Mod;
	}
	return res;
}

const int MAXN = 400005;

int n, m, Mod;
int a[MAXN], b[MAXN];

namespace Poly {
	const int G = 3;
	const int P[3] = {998244353, 1004535809, 469762049};
	int n, rev[MAXN], inv[3];
	__int128 c[3];
	inline void Init(int lenght) {
		int len = 1, l = -1;
		while (len < lenght) len <<= 1, ++ l;
		n = len;
		for (int i = 0; i < 3; ++ i) inv[i] = fsp(n, P[i] - 2, P[i]);
		c[0] = __int128(P[1] * P[2]) * fsp(P[1] * P[2], P[0] - 2, P[0]);
		c[1] = __int128(P[0] * P[2]) * fsp(P[0] * P[2], P[1] - 2, P[1]);
		c[2] = __int128(P[0] * P[1]) * fsp(P[0] * P[1], P[2] - 2, P[2]);
		for (int i = 1; i <= len; ++ i)
			rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << l);
	}
	inline void NTT(int a[], int sgn, int t) {
		int Mod = P[t];
		for (int i = 1; i <= n; ++ i)
			if (i < rev[i]) swap(a[i], a[rev[i]]);
		for (int flr = 2; flr <= n; flr <<= 1) {
			int Wn = fsp((sgn == 1) ? G : (Mod + 1) / G, (Mod - 1) / flr, Mod);
			for (int mid = flr >> 1, l = 0; l < n; l += flr) {
				int w = 1, x, y;
				for (int i = l; i < l + mid; ++ i, w = w * Wn % Mod)
					x = a[i], y = w * a[i + mid] % Mod, 
					a[i] = (x + y) % Mod, a[i + mid] = (x - y + Mod) % Mod;
			}
		}
		if (!~sgn) 
			for (int i = 0; i < n; ++ i)
				a[i] = a[i] * inv[t] % Mod;
	}
	int t[3][MAXN];
	inline int fmul(int a, int b, int Mod) {
		int res = 0; a %= Mod;
		for (; b; b >>= 1) {
			if (b & 1) res = (res + a) % Mod;
			a = (a << 1) % Mod;
		}
		return res;
	}
	inline int Get(int x1, int x2, int x3, int Mod) {
		int res = fmul(c[0] % Mod, x1, Mod);
		res = (res + fmul(c[1] % Mod, x2, Mod)) % Mod;
		res = (res + fmul(c[2] % Mod, x3, Mod)) % Mod;
		return res;
	}
	inline void Mul(int len, int a[], int b[], int c[]) {
		Init(len);
		for (int i = 0; i < 3; ++ i) {
			NTT(a, 1, i); if (&a != &b) NTT(b, 1, i);
			for (int j = 0; j < n; ++ j)
				t[i][j] = a[j] * b[j] % P[i];
			NTT(a, -1, i); if (&a != &b) NTT(b, -1, i);
			NTT(t[i], -1, i);
		}
//		for (int i = 0; i < n; ++ i) printf ("%lld ", a[i]); puts("");
//		for (int i = 0; i < n; ++ i) printf ("%lld ", b[i]); puts("");
		for (int i = 0; i < n; ++ i) printf ("%lld ", t[0][i]); puts("\n");
		for (int i = 0; i < n; ++ i) printf ("%lld ", t[1][i]); puts("\n");
		for (int i = 0; i < n; ++ i) printf ("%lld ", t[2][i]); puts("\n");
		for (int i = 0; i < n; ++ i)
			c[i] = Get(t[0][i], t[1][i], t[2][i], Mod);
	}
}

signed main()
{
	scanf ("%lld%lld%lld", &n, &m, &Mod);
	++ n, ++ m;
	for (int i = 0; i < n; ++ i) scanf ("%lld", a + i);
	for (int i = 0; i < m; ++ i) scanf ("%lld", b + i);
	Poly::Mul(n + m, a, b, a);
	for (int i = 0; i < n + m - 1; ++ i) printf ("%lld ", a[i]);
	return puts(""), 0;
}
2024/5/29 08:53
加载中...