多项式全家福（缺插值和点值）

写法：

vector写有什么好处？

1.分治NTT的时候不用处理麻烦的下标。

2.传进函数的时候不用传指针而不会爆栈。

3.copy，swap，reverse，clear，resize等vector的函数可以帮助您快速地处理。

vector写的时候注意什么？

由于本人非常非常地懒，所以全部都是用vector写的。

根据测试，在开了 O2优化的情况下，是没有什么区别的。

但是不开O2，你懂的……

还是想用vector怎么办呢？
又经过一轮测试，只要在dft的时候把vector里的东西copy到一个数组，用这个数组做，再copy回去，是差不多的。

---

Dft：

不说了太累了。

我们设一次dft的常数为1，以便后面分析复杂度。

---

求逆：

设对多项式A求逆。

倍增法。

设有两个多项式$B、B0$

$B0\ast A=1(mod{x}^{n/2})$
$B\ast A=1(mod{x}^n)$

其实可以用$B0$去推$B$，推导如下：

首先B0的初值即n=1时$B0=a[0{]}^{-1}$

$B0\ast A=1(mod{x}^{n/2})$
$B0\ast A-1=0(mod{x}^{n/2})$
等式两边同时平方：
$(B0\ast A-1{)}^2=0(mod{x}^n)$
$B{0}^2\ast A^2-2B0\ast A+1=0(mod{x}^n)$
$2B0\ast A-B{0}^2\ast A^2=1(mod{x}^n)$
按A整理：
$A\ast (2B0-B{0}^2\ast A)=1(mod{x}^n)$
所以$B=2B0-B{0}^2\ast A(mod{x}^n)$

我们来思考一下怎么做会快一点：
对B0求dft，对A求dft，然后算一算，再IDFT回去。
那么需要三次，但是每层的次数界必须开到4*n

复杂度为$T(n)=nlogn+T(n/2)\approx nlogn$
因为$n+n/2+n/4+\dots <2n$

所以复杂度大概为6次dft

---

开二次根：

方法类似于求逆，也是用倍增的思想。

已知$B{0}^2=A(mod{x}^{n/2})$，推$B^2=A(mod{x}^n)$

B0的初值可能需要二次剩余，一般多项式开根的题的常数项是个常数，可以直接预处理。

二次剩余戳这儿.

$B{0}^2=B^2(mod{x}^{n/2})$
$B{0}^2-B^2=0(mod{x}^{n/2})$
$(B{0}^2-B^2{)}^2=0(mod{x}^n)$
$B{0}^4-2B{0}^2\ast B^2+B^4=0(mod{x}^n)$
$B{0}^4+2B{0}^2\ast B^2+B^4=4B{0}^2\ast B^2(mod{x}^n)$
$(B{0}^2+B^2{)}^2=(2B0B{)}^2(mod{x}^n)$
$B{0}^2+B^2=2B0B(mod{x}^n)$
$B{0}^2+A=2B0B(mod{x}^n)$
$B=B0/2+A/(2B0)(mod{x}^n)$

每层做的时候注意次数界是2n

每层要$dft$两次和求逆一次，复杂度大概是2*(3+6)=18次dft？

---

取模：

给出$A(x),B(x),deg(A)>=deg(B)$
求满足$A(x)=B(x)\ast C(x)+D(x)(deg(D)<deg(B))$的$C、D$

设$n=deg(A),m=deg(B)$
$A(x)=B(x)\ast C(x)+D(x)$
$A(\frac{1}{x})=B(\frac{1}{x})\ast C(\frac{1}{x})+D(\frac{1}{x})$
$x^n\ast A(\frac{1}{x})=x^m\ast B(\frac{1}{x})\ast x^{n-m}\ast C(\frac{1}{x})+x^n\ast D(\frac{1}{x})$

我们来理一下次数的关系：
$x^n\ast A(\frac{1}{x})\in [0..n]$
$x^m\ast B(\frac{1}{x})\ast x^{n-m}\ast C(\frac{1}{x})\in [0..n]$
$x^n\ast D(\frac{1}{x})\in [n-m+1..n]$

对等式两边同时$mod{x}^{n-m+1}$，则$x^n\ast D(\frac{1}{x})$没了。

$x^n\ast A(\frac{1}{x})=x^m\ast B(\frac{1}{x})\ast x^{n-m}\ast C(\frac{1}{x})(mod{x}^{n-m+1})$

设$A^T$就是把A反过来。

那么：
$A^T=B^T\ast C^T(mod{x}^{n-m+1})$
$C=(A^T/B^{T}mod{x}^{n-m+1}{)}^T$

坑点：必须要先$mod{x}^{n-m+1}$再倒置。

求C要一次求逆和一次NTT，需要6+3=9次dft。

求了C后，$D=A-B\ast C$，需要一次NTT，即3次dft。

---

对数：

前置技能：

1.复合函数求导：

$F(G(x){)}^{\prime }=F^{\prime }(G(x))\ast G^{\prime }(x)$

2.ln函数求导：

$ln(x{)}^{\prime }=\frac{1}{x}$

---

$ln(A)$
$=\int ln(x{)}^{\prime }\mathrm{d}x$
$=\int l{n}^{\prime }(x)\ast x^{\prime }\mathrm{d}x$
$=\int \frac{x^{\prime }}{x}\mathrm{d}x$

只需要求逆一次和NTT一次，9次dft

---

指数：

前置技能：

牛顿迭代：

https://blog.csdn.net/ccnt_2012/article/details/81837154

可以发现这就是个求函数零点的玩意儿，虽然有些局限性，但是面对exp函数这么优美的图像肯定是没有问题的。

---

现在来思考多项式牛顿迭代，即求$f(B)=0$，其中$B$是一个多项式。

思路是倍增，即知道了前$n/2$项的答案，去推前$n$项的答案。

我们已知前$n/2$项为多项式B0，想要求出前$n$项B。

可以直接由牛顿迭代得出：

$B=B0-\frac{f(B0)}{f^{\prime }(B0)}(mod{x}^n)$

下面介绍一下用泰勒展开推的方法：

$f(B)={\sum }_{i=0}^{\infty }{f}^{(i)}(B0)\ast (B-B0{)}^i/i!$

注意当i>1时，每一项最低次项都是$x^n$

即：
$f(B)=f(B0)+f^{\prime }(B0)\ast (B-B0)(mod{x}^n)$
$\because f(B)=0(mod{x}^n)$
$\therefore f(B0)+f^{\prime }(B0)\ast (B-B0)=0(mod{x}^n)$
$B=B0-\frac{f(B0)}{f^{\prime }(B0)}(mod{x}^n)$

---

会了牛顿迭代之后思考如何把$exp(A)$换出牛顿迭代的形式：
设$exp(A)=B$
等式两边同时取$ln$，得$A=ln(B)$
那么$f(B)=ln(B)-A$，现在就是求$f$的零点。

套多项式牛顿迭代式得：
$f(B)=B0-\frac{f(B0)}{f^{\prime }(B0)}(mod{x}^n)$
把$f(B)=ln(B)-A$代入，得：
$f(B)=B0-(ln(B0)-A)/(ln(B0)-A{)}^{\prime }(mod{x}^n)$
$f(B)=B0-(ln(B0)-A)/(ln(B0{)}^{\prime }-A^{\prime })(mod{x}^n)$
注意这里的函数的自变量是一个多项式，因变量也是。
$A$的话其实算作常数，求导之后是没有的。
因为自变量是多项式，所以$ln(B0{)}^{\prime }$不算复合函数求导，只是对$ln$求导，$ln(B0{)}^{\prime }=1/B0$。
$f(B)=B0-(ln(B0)-A)\ast B0(mod{x}^n)$
$f(B)=B0\ast (1-ln(B0)+A)(mod{x}^n)$

每层需要用到$ln$和一次乘法，复杂度为$(9+3)\ast 2=24$次dft

模板：

#include<ctime>
#include<vector>
#include<cstdio>
#include<algorithm>
#define ll long long
#define pp printf
#define fo(i, x, y) for(int i = x, B = y; i <= B; i ++)
#define ff(i, x, y) for(int i = x, B = y; i <  B; i ++)
#define fd(i, x, y) for(int i = x, B = y; i >= B; i --)
#define pb push_back
#define si size()
using namespace std;

typedef vector<ll> V;

const int N = (1 << 19) + 5;

const int mo = 998244353;

const int ni2 = (mo + 1) / 2;

ll ksm(ll x, ll y) {
    ll s = 1;
    for(; y; y /= 2, x = x * x % mo)
        if(y & 1) s = s * x % mo;
    return s;
}

int r[N]; ll aa[N];

void dft(V &b, int f) {
	int n = b.si;
    ff(i, 0, n) aa[i] = b[i];
    ff(i, 0, n) {
        r[i] = r[i / 2] / 2 + (i & 1) * (n / 2);
        if(i < r[i]) swap(aa[i], aa[r[i]]);
    }
    for(int h = 1; h < n; h *= 2) {
        ll wn = ksm(ksm(3, (mo - 1) / 2 / h), f == -1 ? mo - 2 : 1);
        for(int j = 0; j < n; j += 2 * h) {
            ll A, W = 1, *l = aa + j, *r = aa + j + h;
            ff(i, 0, h) {
                A = *r * W, *r = (*l - A) % mo, *l = (*l + A) % mo;
                W = W * wn % mo; l ++; r ++;
            }
        }
    }
    if(f == -1) {
        ll v = ksm(n, mo - 2);
        ff(i, 0, n) aa[i] = (aa[i] + mo) * v % mo;
    }
    ff(i, 0, n) b[i] = aa[i];
}

V operator * (V a, V b) {
    int z = a.si + b.si - 1;
    int n = 1; while(n < z) n *= 2;
    a.resize(n); b.resize(n);
    dft(a, 1); dft(b, 1);
    ff(i, 0, n) a[i] = a[i] * b[i] % mo;
    dft(a, -1); a.resize(z); return a;
}

V operator + (V a, V b) {
	int sz = max(a.si, b.si); a.resize(sz); b.resize(sz);
	ff(i, 0, sz) a[i] = (a[i] + b[i] + mo) % mo;
	return a;
}

V operator - (V a, V b) {
	int sz = max(a.si, b.si); a.resize(sz); b.resize(sz);
	ff(i, 0, sz) a[i] = (a[i] - b[i] + mo) % mo;
	return a;
}

V qni(V a) {
    int n0 = 1; while(n0 < a.si) n0 *= 2;
    V b; b.clear(); b.pb(ksm(a[0], mo - 2));
    for(int n = 2; n <= n0; n *= 2) {
        b.resize(n * 2); dft(b, 1);
        V c = a; c.resize(n); c.resize(n * 2); dft(c, 1);
        ff(i, 0, n * 2) b[i] = (b[i] * 2 - b[i] * b[i] % mo * c[i]) % mo;
        dft(b, -1); b.resize(n);
    }
    b.resize(a.si); return b;
}

void fan(V &a) {
	reverse(a.begin(), a.end());
}

V div(V a, V b) {
	int n = a.size() - b.size() + 1;
    fan(a); fan(b);
    b.resize(a.size()); b = qni(b);
    a = a * b; a.resize(n);
    fan(a);
    return a;
} 

V qmo(V a, V b) {
	return a - (b * (div(a, b)));
}

ll w;
struct P {
	ll x, y;
	P(){}
	P(ll _x, ll _y) {x = _x, y = _y;}
};
P operator *(P a, P b) { return P((a.x * b.x + a.y * b.y % mo * w) % mo, (a.x * b.y + a.y * b.x) % mo);}
int pd(ll x) {
	return ksm(x, (mo - 1) / 2) == 1;
}
P ksm(P x, ll y) {
	P s = P(1, 0);
	for(; y; y /= 2, x = x * x)
		if(y & 1) s = s * x;
	return s;
}
ll Sqrt(ll n) {
	if(!n) return 0;
	while(1) {
		ll a = rand() % mo;
		w = (a * a - n + mo) % mo;
		if(pd(w)) continue;
		return ksm(P(a, 1), (mo + 1) / 2).x;
	}
}

V Sqrt(V a) {
	int n0 = 1; while(n0 < a.si) n0 *= 2;
	V b; b.clear(); b.push_back(Sqrt(a[0])); b[0] = 1;
	for(int n = 2; n <= n0; n *= 2) {
		V c = a; c.resize(n);
		V d = b; d.resize(n);
		ff(i, 0, d.si) d[i] = d[i] * 2 % mo;
		ff(i, 0, b.si) b[i] = b[i] * ni2 % mo;
		c = c * qni(d);
		b = b + c;
	}
	b.resize(n0); return b;
}

V qd(V a) {
    a[0] = 0;
    ff(i, 1, a.si) a[i - 1] = a[i] * i % mo;
    return a;
}

V jf(V a) {
    a.pb(0);
    fd(i, a.si, 1) a[i] = a[i - 1] * ksm(i, mo - 2) % mo;
    a[0] = 0;
    return a;
}

V ln(V a) {
    int sa = a.si;
    V b = a; b = qni(b); a = qd(a);
    a = a * b; a = jf(a); a.resize(sa);
    return a;
}

V exp(V a) {
    int n0 = 1; while(n0 < a.si) n0 *= 2;
    V b; b.clear(); b.pb(1);
    for(int n = 2; n <= n0; n *= 2) {
		V c = b; c.resize(n); c = ln(c);
		V d = a; d.resize(n);
		c = c - d;
		c.resize(2 * n); dft(c, 1);
		b.resize(2 * n); dft(b, 1);
		ff(i, 0, 2 * n) b[i] = (b[i] - c[i] * b[i]) % mo;
		dft(b, -1); b.resize(n);
	}
	b.resize(a.si);
	return b;
}

int n, m;

V a, b;

int main() {
	srand(time (0));
}