[Luogu4239] 【模板】多项式求逆（加强版） [多项式求逆][MTT][三模数NTT]

_er

于 2020-08-25 23:53:35 发布

阅读量149

点赞数

本文链接：https://blog.csdn.net/estia_/article/details/89186270

版权

Link
Luogu - https://www.luogu.org/problemnew/show/P4239

笔记本写题看题调题的效率都好低啊啊啊rz

共轭FFT参考 zjt’s blog 和 cmxrynp’s blog

$P = a + b i$
DFT：
$f a + i f b = f (a + b i) = f P$
$fa-ifb=\left[\overline{f(a+bi)}\right]'=fQ$
$fQ[k]=\overline{fP[n-k\pmod{n}]}$

IDFT：鸽

另外单次FFT还可以拆分奇偶。
这样的话MTT（拆系数FFT）就可以被优化到只做 3.5 次 FFT。
听说精度也会高一些（不过巨难调
具体参考 cmxrynp’s blog
注意预处理FFT的时候单位根指数没有乘2，但是在用到的时候是一样的

重要性质：
$\overline{z_1z_2}=\overline{z_1}\cdot\overline{z_2}$
$\overline{f}a=(fa)'$
$F=f^{-1}=\overline{f}/n$
由上有 $Fa=\frac{\overline fa}{n}=\frac{(fa)'}{n}$

myy Paper里面的奇偶优化是直接基于多项式乘法考虑的
~~那那那我手推IDFT好了吧~~
DFT是
$f(x)=f_0(x^2)+xf_1(x^2)$
$f(\omega_n^k)=f_0(\omega_{n/2}^k)+\omega_n^kf_1(\omega_{n/2}^k)$
把 $f_0$ 和 $f_1$ 共轭FFT即可。合并得到 $f(\omega_n^k)$ 。

IDFT
先把 $f a$ 变成 $F a$ （当然你也可以按别人博客说的那样在最后再IDFT？）
$F(\omega_n^k)=F_0(\omega_{n/2}^k)+\omega_n^kF_1(\omega_{n/2}^k)$
$F(\omega_n^{k+n/2})=F_0(\omega_{n/2}^k)-\omega_n^kF_1(\omega_{n/2}^k)$
$F_0(\omega_{n/2}^k)=\frac{F(\omega_n^k)+F(\omega_n^{k+n/2})}{2}$
$\omega_n^kF_1(\omega_{n/2}^k)=\frac{F(\omega_n^k)-F(\omega_n^{k+n/2})}{2}$
$P=F_0+iF_1$
$fP=\frac{F(\omega_n^k)+F(\omega_n^{k+n/2})}{2}+\frac{F(\omega_n^k)-F(\omega_n^{k+n/2})}{2}\cdot i/\omega_n^k=fF_0+ifF_1$
$fP=\frac{F(\omega_n^k)+F(\omega_n^{k+n/2})}{2}+\frac{F(\omega_n^k)-F(\omega_n^{k+n/2})}{2i}\cdot \omega_n^k=fF_0+ifF_1$
因为 $fF_0$ 和 $fF_1$ 都是实数，所以取 $f P$ 的实部和虚部即可。

PS：多项式求逆的时候，求逆是没有必要强制长度为 2 的幂的。。
处理很简单用vector的话就更简单了

#include<cstdio>
#include<cmath>
#include<iostream>
#include<algorithm>
#include<cstdlib>
#include<cstring>
#include<vector>
using namespace std;
char frBB[1<<12],*frS=frBB,*frT=frBB;
#define getchar() (frS==frT&&(frT=(frS=frBB)+fread(frBB,1,1<<12,stdin),frS==frT)?EOF:*frS++)
template<typename T>inline void read(T&x)
{
    x=0;bool w=0;char ch=getchar();
    while(!isdigit(ch))w|=(ch=='-'),ch=getchar();
    while(isdigit(ch))x=x*10+(ch^48),ch=getchar();
    w?(x=-x):0;
}
#define R register
const int MAXN = 3e5 + 5;
const double tpi = acos(-1.0);
const int P = 1e9 + 7;
const long long lP = 1e9 + 7;
int qpowm(long long a, int b)
{
	static long long ret;
	ret = 1;
	while (b)
	{
		if (b & 1) ret = ret * a % lP;
		a = a * a % lP; b >>= 1;
	}
	return ret;
}
namespace Poly
{
    struct Complex
    {
        double r, i;
        inline Complex(const double& Re = 0, const double& Im = 0) { r = Re; i = Im;}
        inline Complex operator + (const Complex& o) const { return Complex(r + o.r, i + o.i);}
        inline Complex operator - (const Complex& o) const { return Complex(r - o.r, i - o.i);}
        inline Complex operator * (const Complex& o) const { return Complex(r * o.r - i * o.i, r * o.i + i * o.r);}
        inline Complex operator * (const double& x) const { return Complex(r * x, i * x);}
        inline Complex operator ~ () const { return Complex(r, -i);}
        inline void operator += (const Complex& o) { r += o.r; i += o.i;}
        inline void operator -= (const Complex& o) { r -= o.r; i -= o.i;}
        inline void operator *= (const Complex& o) { *this = *this * o;}
        inline void operator *= (const double& x) { r = r * x; i = i * x;}
    }A[MAXN], B[MAXN], C[MAXN], D[MAXN], Wn[MAXN];
    int Calc(const int& x)
    {
        register int t = 1;
        while (t < x) t <<= 1;
        return t;
    }
    void PreWn(const int& n)
    {
        for (R int i = 1; i < n; i <<= 1)
        {
            Wn[i] = Complex(1, 0);
            for (R int j = 1; j < i; ++j)
            {
                Wn[i + j] = ((j & 31) == 1) ? Complex(cos(tpi * j / i), sin(tpi * j / i)) : Wn[i + j - 1] * Wn[i + 1];
            }
        }
    }
    void FFT(Complex* F, const int& n)
    {
        for (R int j = 0, i = 0; i < n; ++i)
        {
            if (i > j) swap(F[i], F[j]);
            for (R int k = n >> 1; (j ^= k) < k; k >>= 1);
        }
        R Complex Y;
        for (R int Cur, Len, Mid = 1; Mid < n; Mid <<= 1)
        {
            Len = Mid << 1;
            for (R int Pos = 0; Pos < n; Pos += Len)
            {
                Cur = Pos;
                for (R int Sub = 0; Sub < Mid; ++Sub, ++Cur)
                {
                    Y = Wn[Mid + Sub] * F[Cur + Mid]; // omega_n^k\cdot F(n^{\frac{n}{2}+k})
                    F[Cur + Mid] = F[Cur] - Y;
                    F[Cur] += Y;
                }
            }
        }
    }
    void DFT(Complex* F, int n)
    {
        if (n == 1) return;
        static Complex P[MAXN>>1];
        n >>= 1;
        for (R int i = 0; i < n; ++i) P[i].r = F[i<<1].r, P[i].i = F[i<<1|1].r;
        FFT(P, n);
        R Complex Q, X, Y;
        for (R int i = 0; i < n; ++i)
        {
            Q = ~P[(n-i)&(n-1)];
            X = (P[i] + Q) * .5; // * .5 = / 2
            Y = (P[i] - Q) * Complex(0, -.5) * Wn[n + i]; // * -.5i = / 2i
            F[i] = X + Y;
            F[n + i] = X - Y; 
        }
    }
    void IDFT(Complex* F, int n)
    {
        if (n == 1) return;
        static Complex P[MAXN>>1];
        reverse(F + 1, F + n); n >>= 1;
        for (R int i = 0; i < n; ++i) P[i] = (F[i] + F[n + i]) * .5 + (F[i] - F[n + i]) * Complex(0, .5) * Wn[n + i];
        FFT(P, n);
        R double t = 1. / n;
        for (R int i = 0; i < n; ++i) F[i<<1] = Complex(P[i].r * t, 0), F[i<<1|1] = Complex(P[i].i * t, 0);
    }
    inline long long Estimate(const double& x) { return (long long)(x + .5)%lP;}
    inline long long Adjust(long long x)
    {
        x = x % lP + lP;
        return (x >= lP) ? (x - lP) : x;
    }
    vector<int> Multi(const vector<int>& a, const vector<int>& b)
    {
    	R int l = a.size() + b.size() - 1;
    	vector<int> ans(l); l = Calc(l);
    	fill(A, A + l, Complex(0, 0)); fill(B, B + l, Complex(0, 0));
    	fill(C, C + l, Complex(0, 0)); fill(D, D + l, Complex(0, 0));
        for (R int i = 0; i < a.size(); ++i) A[i].r = a[i] >> 15, B[i].r = a[i] & 32767;
		for (R int i = 0; i < b.size(); ++i) C[i].r = b[i] >> 15, D[i].r = b[i] & 32767;
        DFT(A, l), DFT(B, l), DFT(C, l), DFT(D, l);
        R Complex t;
        for (R int i = 0; i < l; ++i) 
        {
        	t = A[i];
			A[i] = A[i] * D[i] + B[i] * C[i];
			B[i] *= D[i]; C[i] *= t;
        }
        IDFT(A, l); IDFT(B, l); IDFT(C, l);
        for (R int i = 0; i < ans.size(); ++i) ans[i] = Adjust((Estimate(C[i].r) << 30) + (Estimate(A[i].r) << 15) + Estimate(B[i].r));
        return ans;
    }
    int c[MAXN], d[MAXN];
    vector<int> Inv(const vector<int>& a, int n)
    {
    	if (n == -1) n = a.size();
    	if (n == 1) { vector<int> ans; return ans.push_back(qpowm(a[0], P - 2)), ans;}
    	vector<int> inv = Inv(a, n + 1 >> 1), tmp(&a[0], &a[0] + n);
    	tmp = Multi(Multi(tmp, inv), inv), tmp.resize(n);
    	for (R int i = 0; i < inv.size(); ++i) { tmp[i] = (2 * inv[i] - tmp[i]) % P; (tmp[i]<0) ? (tmp[i]+=P) : 0;}
    	for (R int i = inv.size(); i < n; ++i) tmp[i] ? (tmp[i] = P - tmp[i]) : 0;
    	return tmp;
    }
}
int n, l;
vector<int> F;
int main()
{
	read(n); F.resize(n);
	for (R int i = 0; i < n; ++i) read(F[i]);
	Poly::PreWn(Poly::Calc(n<<1));
	F = Poly::Inv(F, -1);
	for (R int i = 0; i < n; ++i) printf("%d ", F[i]);
    return 0;
}

_er

关注

0
点赞
踩
0

收藏

觉得还不错? 一键收藏
0
评论
[Luogu4239] 【模板】多项式求逆（加强版） [多项式求逆][MTT][三模数NTT]

#include<cstdio>#include<iostream>#include<algorithm>#include<cstdlib>#include<cstring>#include<cmath>#include<ctime>#include<cctype>using names...
复制链接

扫一扫