FFT & NTT 模板

$1004535809=479\times 2^{21}+1$ 相加不会超过 $\text{INT_MAX}$

$2281701377=17\times 2^{27}+1$ 相乘不会超过 $\text{LONG_LONG_MAX}$。

/* * FFT.cc * * @date 2018.2.13 * @author Kvar_ispw17 * @tag Fast Fourier Transform * @brief This is a Fast Fourier Transform-> * Implement by Kvar_ispw17(enkerewpo@gmail.com) * * Copyright 2018 Kvar_ispw17 __Hyperx All rights reserved. * */ #include <algorithm> #include <stdexcept> #include <iostream> #include <cstring> #include <cstdio> #include <cmath> #define S const int N = 1e5 + 5; // COMMENT THE DEFINATION YOU DON'T NEED BELOW TO RUN THE PROGRAM #ifdef S // NUMBER THEORIETIC TRANSFORM #define NTT_INTEGER_VERSION #else // FAST FOURIER TRANSFORM #define FFT_COMPLEX_VERSION #endif #if defined(FFT_COMPLEX_VERSION) && defined(NTT_INTEGER_VERSION) #undef FFT_COMPLEX_VERSION #undef NTT_INTEGER_VERSION int main() { try { throw std::runtime_error("<<<PATTERN NOT SELECTED>>") } catch (...) { fprintf(stderr, "[WARNING] : PLEASE SELECT A PROGRAM TO RUN INTO.\n"); } return 1; } #endif #ifdef FFT_COMPLEX_VERSION typedef long long ll; const long double PI = acos(-1.0); const long double EPS = 1E-8; class complex { public: long double re, im; complex() {} complex(long double re, long double im) : re(re), im(im) {} complex operator + (const complex &x) { return complex(re + x.re, im + x.im); } complex operator - (const complex &x) { return complex(re - x.re, im - x.im); } complex operator * (const complex &x) { return complex(re * x.re - im * x.im, im * x.re + re * x.im); } complex operator / (const complex &x) { return complex((re * x.re + im * x.im) / (x.re * x.re + x.im * x.im), (im * x.re - re * x.im) / (x.re * x.re + x.im * x.im)); } }; void print(complex x, char c) { printf("\t%.10lf", fabs(x.re) < EPS ? 0.0 : (double)x.re); printf("\t\t(%.10lf)", fabs(x.im) < EPS ? 0.0 : (double)x.im); printf("%c", c); } int n, rev[N]; complex F[N], w[N]; void FFT(complex * F, int n, int offset) { for (int i = 0; i < n; i++) if (rev[i] > i) std::swap(F[i], F[rev[i]]); for (int i = 2; i <= n; i <<= 1) { complex wi(cos(offset * 2 * PI / i), sin(offset * 2 * PI / i)); for (int j = 0; j < n; j += i) { complex w(1, 0); for (int k = j, h = i >> 1; k < j + h; k++) { complex t = w * F[k + h], u = F[k]; F[k] = u + t; F[k + h] = u - t; w = w * wi; } } } if (offset == -1) for (int i = 0; i < n; i++) F[i].re = F[i].re / n; } int main() { try { printf("FFT MODE RUNNING...\n"); scanf("%d", &n); for (int i = 0; i < n; i++) { double x; scanf("%lf", &x); F[i].re = x; } n = 1 << (int)ceil(log2(n)); for (int i = 0; i < n; i++) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << ((ll)log2(n) - 1)); printf("\nAFTER INPUT:\n"); for (int i = 0; i < n; i++) print(F[i], '\n'); FFT(F, n, 1); printf("\nAFTER FFT:\n"); for (int i = 0; i < n; i++) print(F[i], '\n'); FFT(F, n, -1); printf("\nAFTER IFFT:\n"); for (int i = 0; i < n; i++) print(F[i], '\n'); return 0; } catch (...) { return 1; } } #endif #ifdef NTT_INTEGER_VERSION typedef long long ll; ll n, inv_n, F[N], rev[N], q; const ll MOD = 1004535809; // = 479 * 2 ^ 21 + 1 const ll g = 3; ll q_pow(ll a, ll b) { ll ret = 1; while (b) { if (b & 1) ret = ret * a % MOD; a = a * a % MOD; b >>= 1; } return ret; } void NTT(ll F[], ll n, int offset) { for (int i = 0; i < n; i++) if (rev[i] > i) std::swap(F[i], F[rev[i]]); for (int i = 2; i <= n; i <<= 1) { ll wi = q_pow(g, offset == 1 ? (MOD - 1) / i : MOD - 1 - (MOD - 1) / i); for (int j = 0; j < n; j += i) { ll w = 1; for (int k = j, h = i >> 1; k < j + h; k++) { ll t = w * F[k + h], u = F[k]; F[k] = (u + t) % MOD; F[k + h] = ((u - t) % MOD + MOD) % MOD; w = w * wi % MOD; } } } if (offset == -1) for (int i = 0; i < n; i++) F[i] = F[i] * inv_n % MOD; } int main() { try { printf("NTT MODE RUNNING...\n"); scanf("%lld", &n); for (int i = 0; i < n; i++) scanf("%lld", &F[i]); n = 1 << (int)ceil(log2(n)); printf("%lld\n", n); inv_n = q_pow(n, MOD - 2); printf("[DEBUG: inverse of n is %lld]\n", inv_n); for (int i = 0; i < n; i++) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << ((ll)log2(n) - 1)); printf("\nAFTER INPUT:\n"); for (int i = 0; i < n; i++) printf("\t%lld\n", F[i]); NTT(F, n, 1); printf("\nAFTER NTT:\n"); for (int i = 0; i < n; i++) printf("\t%lld\n", F[i]); NTT(F, n, -1); printf("\nAFTER INTT:\n"); for (int i = 0; i < n; i++) printf("\t%lld\n", F[i]); return 0; } catch (...) { return 1; } } #endif /* [TEST CASE] - 8 1 2 3 4 5 6 7 8 - 9 1.35131334 3.56336636 -3.64242251 -8.3575486 7.86623785 -2.4457546 1.35332446 -9.94949959 1.23434524 - */
Code language: PHP (php)

发表回复

您的电子邮箱地址不会被公开。 必填项已用*标注