50pts 拆系数FFT求助

P4245 【模板】任意模数多项式乘法

Albert_Wei @ 2023-04-27 12:31:19

rt

#include <bits/stdc++.h>
#define LL long long
using namespace std;

namespace Polynomial{
 const double Pi = acos(-1), eps = 0.5;
 const int MAXN = 4000005, p = 998244353, g = 3, invg = 332748118;
 struct comp{
  double x, y;
  comp(double a = 0, double b = 0) {x = a, y = b;}
 };
 comp operator + (comp x, comp y) {return comp(x.x + y.x, x.y + y.y);}
 comp operator - (comp x, comp y) {return comp(x.x - y.x, x.y - y.y);}
 comp operator * (comp x, comp y) {return comp(x.x * y.x - x.y * y.y, x.x * y.y + x.y * y.x);}
 comp operator / (comp x, comp y) {
  double a = x.x, b = x.y, c = y.x, d = y.y;
  return comp((a * c + b * d) / (c * c + d * d), (b * c - a * d) / (c * c + d * d));
 }
 int num[MAXN];
 comp A[MAXN], B[MAXN], C[MAXN], Wn[MAXN];
 inline void FFT(int Lim, comp *A, int f) {
  for (int i = 0; i < Lim; i++)
   if (i < num[i]) swap(A[i], A[num[i]]);
  for (int i = 1; i < Lim; i <<= 1) {
   comp w = Wn[i]; w.y *= f;
   for (int j = 0; j < Lim; j += (i << 1)) {
    comp w0 = comp(1, 0);
    for (int k = j; k < i + j; k++, w0 = w * w0) {
     comp f = A[k], g = A[i + k] * w0;
     A[k] = f + g, A[i + k] = f - g;
    }
   }
  }
 }
 inline void modMul(LL *a, LL *b, int n, int m, int mod) {
  const int t = 32768;
  for (int i = 0; i <= n; i++) A[i] = comp(a[i] / t, a[i] % t);
  for (int i = 0; i <= m; i++) B[i] = b[i] / t, C[i] = b[i] % t;
  int len = 1; while (len <= n + m) len <<= 1;
  for (int i = 0; i < len; i++) Wn[i] = comp(cos(Pi / i), sin(Pi / i));
  for (int i = 0; i < len; i++) num[i] = (num[i >> 1] >> 1) | (i & 1) * (len >> 1);
  FFT(len, A, 1), FFT(len, B, 1), FFT(len, C, 1);
  for (int i = 0; i < len; i++) B[i] = A[i] * B[i], C[i] = A[i] * C[i];
  FFT(len, B, -1), FFT(len, C, -1);
  for (int i = 0; i < len; i++) {
   a[i] = ((LL)(B[i].x / len + eps) * (t * t % mod)) % mod;
   a[i] = (a[i] + (((LL)(B[i].y / len + eps) + (LL)(C[i].x / len + eps)) * t)) % mod;
   a[i] = ((LL)(C[i].y / len + eps) + a[i]) % mod;
  }
 }
}
using namespace Polynomial; 

int n, m, k;
LL a[MAXN], b[MAXN];
signed main() {
 cin >> n >> m >> k;
 for (int i = 0; i <= n; i++) cin >> a[i];
 for (int i = 0; i <= m; i++) cin >> b[i];
 modMul(a, b, n, m, k);
 for (int i = 0; i <= n + m; i++) cout << a[i] << " ";
 return 0;
}

by 飞雨烟雁 @ 2023-04-27 13:22:18

@Unbeliveble

应该是精度问题和取模问题。前者可暂时用 define double long double 解决,后者只需把 modMul 函数最后 a[i] 赋值多加几个取模就好了,具体如下:

a[i] = ((LL)(B[i].x / len + eps) * (t * t % mod)) % mod;

改为:

a[i] = ((LL)(B[i].x / len + eps) % mod * (t * t % mod)) % mod;


by Albert_Wei @ 2023-04-27 18:18:58

@飞雨烟雁 对不起,刚刚出去了亿会。A了,thx


by bloodstalk @ 2023-06-22 19:22:49

@飞雨烟雁 感谢提供思路


|