CDQ分治做FFT

CDQ分治+FFT

要计算的式子一般有这样的特性:自己的后项由前项推出。时间 O(n\log^2n)

f_n = \sum_{0 \le i \le n} a_i f_{n-i}

// dp[i] = sigma(a[j] * dp[i-j]) (j < i);
num_t x[MAXN], y[MAXN], a[MAXN], dp[MAXN];

void cdqFFT(int L, int R)
{
    if (L == R) return; // (assert that a[0] == 0)
    int mid = (L + R) >> 1;
    cdqFFT(L, mid);

    int len = R - L + 1;
    trans.init(log2_ceil2(len));

    // x = a0, a1, a2, ..., a[len], 0, 0, ...
    // y = dp[L], dp[L+1], dp[L+2], ..., dp[mid], 0, 0, ...
    for (int i = 0; i < trans.N; i++)
    {
        x[i] = i < len ? a[i] : 0;
        y[i] = (i <= mid-L) ? dp[i+L] : 0;
    }

    trans.DFT(x); trans.DFT(y);
    for (int i = 0; i < trans.N; ++i)
        x[i] = x[i] * y[i];
    trans.IDFT(x);

    for (int i = mid + 1; i <= R; ++i)
        dp[i] += x[i-L];
    cdqFFT(mid + 1, R);
}

A_0 = 1; A_{n+1} = \frac{1}{6} \left( \sum_{i+j+k=n} A_iA_jA_k + 3\sum_{i+2j=n} A_iA_j + 2\sum_{3i=n} A_i \right)

num_t tp[MAXN], tp2[MAXN], tp3[MAXN], A[MAXN];

inline void cdqFFT(int l, int r)
{
    if (l == r)
    {
        if (l == 0)
        {
            A[l] = 1;
        }
        else
        {
            A[l] = (l % 3 == 1) ? (A[l] + A[l/3] * 2) % mod : A[l]; // for the last item.
            A[l] = A[l] * inv6 % mod;
        }
        return; // Last modify of A[i], so we can divide 6 here.
    }

    int mid = (l + r) >> 1, len = r - l;
    cdqFFT(l, mid);

    trans.init(log2_ceil2(len * 2 + mid - l));

    // tp[]  = A[l], A[l+1], A[l+2], ..., A[mid], 0, 0, ...
    // tp2[] = A[0], A[1], A[2], ..., A[len-1], 0, 0, ...
    // tp3[] = A[0], 0, A[1], 0, A[2], 0, ..., at most A[len/2].
    for (int i = 0; i < trans.N; i++)
    {
        tp[i] = (i <= mid - l) ? A[i + l] : 0;
        tp2[i] = (i < len) ? A[i] : 0;
        tp3[i] = ((i & 1) == 0 && i <= len) ? A[i / 2] : 0;
    }

    trans.DFT(tp); trans.DFT(tp2); trans.DFT(tp3);

    for (int i = 0; i < trans.N; i++)
        tp2[i] = (tp2[i] * tp2[i] % mod) * tp[i] % mod;
    trans.IDFT(tp2); // tp2-new = A[i]A[j]A[k], i+j+k=n

    if (!l)
    {
        for (int i = mid + 1; i <= r; i++)
            A[i] = (A[i] + tp2[i - l - 1]) % mod;
    }
    else
    {
        for (int i = mid + 1; i <= r; i++)
            A[i] = (A[i] + tp2[i - l - 1] * 3) % mod;
    }

    for (int i = 0; i < trans.N; i++)
        tp[i] = tp[i] * tp3[i] % mod; // tp-new = A[i]A[j], i+2j=n
    trans.IDFT(tp);

    for (int i = mid + 1; i <= r; i++)
        A[i] = (A[i] + tp[i - l - 1] * 3) % mod;

    cdqFFT(mid + 1, r);
}

cdqFFT(L,R) 可以将L~R内的所有项数值更新完成。

cdqFFT(L,mid) 后,左边一半处理完成,然后利用左边一半的数值更新右边一半。

以上述第一个卷积为例,卷积完后,有

x_s = \sum_{i+j=s} a_i f_{j+L}

然后向后项更新时,会有

f_k += \sum_{i+j=k} a_i f_j (k \in [mid+1,R], j \in [L,mid])

此时,f_k 已经内含了 \sum_{i=0}^{mid} a_i f_{k-i},只需要计算 \sum_{i=mid+1}^{k} a_i f_{k-i}

以上述第二个卷积为例,

A_s += \sum_{i+j+k=s} A_i A_j A_k (s \in [mid+1,R], i \in [L,mid])

此时,A_s 已经内含了 \sum_{i=0}^{mid} A_i A_j A_k,此时被卷积的 j+k \in [0,\frac{len}{2}],而 \frac{len}{2} \le mid,没有利用到未计算完成的值。

在进行 cdqFFT(mid+1, R) 时,f_k 会被更新到新的 mid 值,直到最后一次被更新。

时间复杂度 T(n) = 2T(\frac{n}{2}) + O(nlogn), T(1) = O(1), T(n) = O(nlog^2n)

多项式计算

多项式求导

\left(\sum_{k=0}^\infty a_k x^k\right)’ = \sum_{k=1}^\infty a_k k x^{k-1}
void diff(int a[], OUT int b[], int len)
{
    for (int i = 1; i < len; i++)
        b[i-1] = 1LL * i * a[i] % mod;
    b[len] = b[len-1] = 0;
}

多项式积分

\int\left(\sum_{k=0}^\infty a_k x^k\right)dx = \sum_{k=0}^\infty \frac{a_k}{k+1} x^{k+1} + C
void integral(int a[], OUT int b[], int len)
{
    for (int i = 1; i < len; i++)
        b[i] = 1LL * a[i-1] * inv[i] % mod;
    b[0] = 0;
}

多项式求逆

当常数项与模互质时有逆元。

B_{t+1}(x) = 2B_t(x) – A(x)B_t^2(x)
num_t A[_], B[_];

void get_inv(int a[], OUT b[], int len)
{
    if (len == 1)
    {
        b[0] = fpow(a[0], mod-2);
        return;
    }

    get_inv(a, b, len>>1);
    for (int i = 0; i < len; i++)
        A[i] = a[i], B[i] = b[i];
    trans.init(log2_floor(len)+1);
    trans.DFT(A); trans.DFT(B);
    for (int i = 0; i < trans.N; i++)
        A[i] = 1LL * A[i] * B[i] % mod * B[i] % mod;
    trans.IDFT(A);
    for (int i = 0; i < len; i++)
        b[i] = (2LL * b[i] % mod - A[i] + mod) % mod;
    for (int i = 0; i < trans.N; i++)
        A[i] = B[i] = 0;
}

牛顿迭代

B_{t+1}(x) = B_t(x) – \frac{F(B_t(x))}{F'(B_t(x))}

多项式开方

B_{t+1}(x) = \frac{1}{2} \left( \frac{A(x)}{B_t(x)} + B_t(x) \right)
num_t C[_], D[_], inv2 = (mod-1)/2;

void get_sqrt(int a[], OUT int b[], int len)
{
    if (len == 1)
    {
        b[0] = sqrt(a[0]);
        return;
    }

    get_sqrt(a, b, len>>1);
    for (int i = 0; i < len; i++)
        C[i] = a[i];
    get_inv(b, D, len);
    trans.init(log2_floor(len)+1);
    trans.DFT(C); trans.DFT(D);
    for (int i = 0; i < trans.N; i++)
        D[i] = 1LL * D[i] * C[i] % mod;
    trans.IDFT(D);
    for (int i = 0; i < len; i++)
        b[i] = 1LL * (b[i] + D[i]) % mod * inv2 % mod;
    for (int i = 0; i < trans.N; i++)
        C[i] = D[i] = 0;
}

拆系数FFT

为了防止 A(x) \cdot B(x) 爆long long,所以设一个常数 M \approx \sqrt{mod}

A(x) = \frac{F(x)}{M}, B(x) = F(x)\% M, C(x) = \frac{G(x)}{M}, D(x) = G(x) \% M

于是卷积就变成了这样的七次:

M^2 A(x) C(x) + M (A(x)D(x) + B(x)C(x)) + B(x)D(x)

快速傅立叶变换

FFT原理

对于 N-1 次多项式,将系数表达式

f(x) = \sum_{0 \le i < N} a_i x^i

和点值表达式

f(x) = \sum_{0 \le i < N} y_i \frac{ \prod_{j \neq i}(x-x_j) }{ \prod_{j \neq i}(x_i-x_j) }, \ y_i = \sum_{0 \le n < N} c_i x_i^n

之间互相转化。其中,点取值为N次单位根。

考虑 N=2^K,方便分治计算呀。

根据一系列推导,可以得知 DFT 计算是这样的:

\left(
\begin{array}{cccccc}
1 &1 &1 &1 &\cdots &1 \\
1 &\omega &\omega^2 &\omega^3 &\cdots &\omega^{n}\\
1 &\omega^2 &\omega^4 &\omega^6 &\cdots &\omega^{2n}\\
1 &\omega^3 &\omega^6 &\omega^9 &\cdots &\omega^{3n}\\
\vdots &\vdots &\vdots &\vdots &\ddots &\vdots\\
1 &\omega^{n}&\omega^{2n}&\omega^{3n}&\cdots&\omega^{n^2}
\end{array}
\right)
\left(
\begin{array}{c}
a_0\\a_1\\a_2\\a_3\\\vdots\\a_{n}
\end{array}
\right)
=
\left(
\begin{array}{c}
y_0\\y_1\\y_2\\y_3\\\vdots\\y_{n}
\end{array}
\right)

而IDFT的计算相当于在右边左乘它的逆矩阵,可以证明,对矩阵每一项取倒数并除以N,就是它的逆矩阵。

X_k = \sum_{0 \le n < N} x_n \omega^{nk}
x_n = \frac{1}{N} \sum_{0 \le k < N} X_k \omega^{-kn}

NTT原理

考虑上述 N次单位根群 \lbrace 1, \omega, \omega^2, \omega^3, \cdots \rbrace 到某个数论乘法群的同构。

我们令 P = 2^K * Q + 1,则寻找它的原根 g,有 2^K 阶循环群 \lbrace 1, g^Q, g^{2Q}, g^{3Q}, \cdots \rbrace,易知他们同构。

由于FFT涉及到的是浮点数运算,可能会增大误差,所以考虑使用整数。毕竟整数乘法更加精确。

在这里,利用 g^Q 替代 \omega 即可。

一般用NTT的题目中,会取 p = 998244353,因为 p = 119 * 2^{23} + 1。此时 g = 3

X_k = \sum_{0 \le n < N} x_n g^{Qnk} \ (mod\ P)
x_n = \frac{1}{N} \sum_{0 \le k < N} X_k g^{-Qkn} \ (mod\ P)

变换动机

多项式乘法

多项式 A(x) = \sum a_n x^n B(x) = \sum b_n x^n,求 A(x)B(x)

朴素的多项式乘法,计算 c_n = \sum a_k b_{n-k} ,时间复杂度约为 O(n^2)

由变换计算出的点值表示法,如果点是相同的,那么

(A \cdot B)(x_i) = A(x_i) \cdot B(x_i)

可以在 O(n) 的时间内得到,还原 (A \cdot B)为系数表示就完成了多项式乘法。

在实现FFT时,通常可以通过分治达到 O(n logn) 的时间复杂度,所以通过 FFT 的多项式乘法的总时间复杂度为 O(n logn)。完成了一次优化。

卷积

定义对于两个函数的运算
f(x) \circledast g(x) = \int_0^x f(u) * g(x-u) du
则上述 c_n 符合了卷积的定义。

考虑连续傅立叶变换 F(\omega) = \mathfrak{F}[f(t)] = \int_{-\infty}^{\infty}f(t)e^{-i\omega t}dt 下象函数和象原函数有这样的关系,

F(\omega) * G(\omega) = \mathfrak{F}[f(t) \circledast g(t)]

暗示DFT是具有循环卷积性质的。

由于卷积可以表达成为多项式乘法的项系数,所以可以利用上述变换达到目的。

变换细节

分治求点值

分治地求当 x = \omega^n 时多项式的值。

\begin{aligned} F(x) &= a_0 + a_1 x + a_2 x^2 + a_3 x^3 + a_4 x^4 + a_5 x^5 + a_6 x^6 + a_7 x^7 \\\\ &= (a_0 + a_2 x^2 + a_4 x^4 + a_6 x^6) + x (a_1 + a_3 x^2 + a_5 x^4 + a_7 x^6) \\\\ &= G(x^2) + x * H(x^2) \end{aligned}

把当前单位复根的平方分别以DFT的方式带入G函数和H函数求值

DFT(F(x))_k = DFT(G(x^2))_k + \omega^k * DFT(H(x^2))_k

将递归消成迭代

从分治的角度来说,会逐步拆分成如下序列:

\lbrace x_0, x_1, x_2, x_3, x_4, x_5, x_6, x_7 \rbrace
\lbrace x_0, x_2, x_4, x_6 \rbrace, \lbrace x_1, x_3, x_5, x_7 \rbrace
\lbrace x_0, x_4 \rbrace, \lbrace x_2, x_6 \rbrace, \lbrace x_1, x_5 \rbrace, \lbrace x_3, x_7 \rbrace
\lbrace x_0 \rbrace, \lbrace x_4 \rbrace, \lbrace x_2 \rbrace, \lbrace x_6 \rbrace, \lbrace x_1 \rbrace, \lbrace x_5 \rbrace, \lbrace x_3 \rbrace, \lbrace x_7 \rbrace

将下标化为二进制,会发现拆分后的序列下标,恰好为长度为3的二进制数的反转。容易证明这个结论的普遍性。

模板

/*******************************************
* 快速傅立叶变换:时间 O(nlogn)
* 多项式系数表示与点值表示转换
* - Fast Fourier Transform 虚数
* - Number Theory Transform 数论数
*******************************************/

const size_t max_len = 17;
#define _rev(i) (rev[i>>1]>>1)|((i&1)<<len-1)

#ifdef NTT // 整数运算
const int mod = 998244353, g = 3;
typedef long long lld;
struct mint { 整数取模运算+-*/ };
typedef mint num_t;

inline num_t create(int m)
{
   int k = (mod-1)/m;
   if (k < 0) k += mod-1;
   return pow(g, k);
}

#else // 浮点数运算
const double PI = acos(-1.0);
typedef complex<double> num_t;

inline num_t create(int m)
{
   return num_t(cos(2*PI/m), sin(2*PI/m));
}

#endif

class Fourier
{
public:
   int len, N;

private:
   num_t wmk[1 << max_len];
   int rev[1 << max_len];

   void dft(num_t a[], int DFT)
   {
      for (int i = 0; i < N; i++)
         if (i < rev[i])
            swap(a[i], a[rev[i]]);

      for (int m = 2; m <= N; m <<= 1)
      {
         int m2 = m >> 1;
         num_t wm = create(DFT * m);
         wmk[0] = 1;
         for (int j = 1; j < m2; j++)
            wmk[j] = wmk[j-1] * wm;
         num_t t, u;

         for (int k = 0; k < N; k += m)
         {
            for (int j = 0; j < m2; j++)
            {
               t = wmk[j] * a[k+j+m2];
               u = a[k+j];
               a[k+j] = u + t;
               a[k+j+m2] = u - t;
           }
         }
      }

      if (DFT == -1)
         for (int i = 0; i < N; i++)
            a[i] = a[i] / N;
   }

public:
   void init(int _len)
   {
      len = _len, N = 1 << _len;
      for (int i = 0; i < N; i++)
         rev[i] = _rev(i);
   }

   void DFT(num_t a[]) { dft(a, 1); }
   void IDFT(num_t a[]) { dft(a, -1); }
};

莫比乌斯反演

Overview – \mu(n)

莫比乌斯函数:先根据算术基本定理,令 n = {p_1}^{e_1} * {p_2}^{e_2} * {p_3}^{e_3} * {…} * {p_k}^{e_k},则有

\mu(n) = \begin{cases} \begin{array}{lcl} 0 & & when\,\exists\,{e_i}>1 \\\\ (-1)^k & & else \\ \end{array} \end{cases}

例如,\mu(1)=1\mu(2)=-1\mu(3)=-1\mu(2 * 2)=0\mu(2 * 3)=1\mu(2 * 3 * 5)=-1\mu(2 * 3 * 5 * 7)=1……

一个验证较为稳定的莫比乌斯函数+欧拉函数+素数线性筛。

const size_t MAXN = 1e6+5;
char miu[MAXN], isnp[MAXN];
int phi[MAXN], pri[MAXN], p = 0, sum[MAXN];

void init_prime()
{
    miu[1] = 1; phi[1] = 1;

    for (int i = 2; i < MAXN; i++)
    {
        if (!isnp[i])
        {
            phi[i] = i-1;
            miu[i] = -1;
            pri[p++] = i;
        }

        for (int j = 0; j < p; j++)
        {
            int k = pri[j] * i;
            if (k >= MAXN) break;
            isnp[k] = 1;
            if (i % pri[j] == 0)
            {
                miu[k] = 0;
                phi[k] = phi[i]*pri[j];
                break;
            }
            else
            {
                phi[k] = phi[i]*(pri[j]-1);
                miu[k] = -miu[i];
            }
        }
    }

    sum[0] = 0; // 分块求和预处理 miu - 前缀和
    for (int i = 1; i < MAXN; i++)
        sum[i] = sum[i-1] + miu[i];

} // phi前缀和有n^2/3杜教筛

Overview – 反演

F(n)=\sum_{d|n}f(d) \Rightarrow f(n)=\sum_{d|n}\mu(d)F(\frac{n}{d})
F(n)=\sum_{n|d}f(d) \Rightarrow f(n)=\sum_{n|d}\mu(\frac{d}{n})F(d)

这个反演和偏序关系的覆盖、整除关系的覆盖有关。组合数学和数论的结晶。反正我还是有那么一点点懵逼。

一个验证较为稳定的模板:

f(x)x \in [1,b],y \in [1,d],gcd(x,y)==n 的对数,使用了分块求和,加速计算用的。很多题基本上用这个模板都可以解决,少数题略改也可以通过。

lld f(int b, int d, int n)
{
    if (!b || !d || !n) return 0;
    lld ret = 0;
    int last = -1;
    b /= n, d /= n;
    if (b > d) swap(b, d);
    for (int i = 1; i <= b; i = last + 1)
    {
        last = min(b / (b/i), d / (d/i));
        ret += 1LL * (sum[last] - sum[i-1]) * (b/i) * (d/i);
    }
    return ret;
}

A – Visible Lattice Points

a,b,c\in[0,n],使gcd(a,b,c)=1成立的个数。

双区间是\left \lfloor n/i \right \rfloor * \left \lfloor n/i \right \rfloor,三区间就是\left \lfloor n/i \right \rfloor * \left \lfloor n/i \right \rfloor * \left \lfloor n/i \right \rfloor呗。

不过,这个是 [1,n] 的情况,要转为0起始的话,考虑(0,b,c)\ (a,0,c)\ (a,b,0)\ (0,0,0),所以最后的反演公式是

f(n)=\sum_{d|n} \mu(d) \left \lfloor \frac{n}{i} \right \rfloor\left \lfloor \frac{n}{i} \right \rfloor(\left \lfloor \frac{n}{i} \right \rfloor+3)

B – GCD

考虑一下重叠区间内的元素个数,直接除以2好啦。

f(b,d,k)-\lfloor\frac{f(t,t,k)}{2}\rfloor\ , t=min(b,d)

C – Mophues

按n的素因子个数分类 \mu(n),前缀和处理。

D – GCD of Sequence

好像是容斥吧。有待爆肝。

E – Gcd

先说莫比乌斯反演的做法。直接枚举所有小于他们的素数然后求和即可。

有一种欧拉函数前缀和的做法。考虑 gcd(x,y)=p \Leftrightarrow gcd(\frac{x}{p},\frac{y}{p})=1,则使用欧拉函数的定义吧。

递推求欧拉函数表应该也挺快的,不过注意(1,x) (x,1) (1,1)不在此中。

F – 能量采集

考虑 cnt[gcd(a,b)=i],意味着 (a,b), (0,0) 点之间有 i-1 棵植物,花费 2i-1。则

\sum_{i=1}^{\infty}(2 i-1) * f(n,m,i)

注意 i 到某个节点之后就为0了,可以结束计算。

G – Problem b

由于区间不一定是1开始的,考虑这样的一个容斥:

[a,b] \times [c,d]=[1,b] \times [1,d]-[1,a-1] \times [1,d]-[1,b] \times [1,c-1]+[1,a-1] \times [1,c-1]

大衍求一术

int dyqys(int a, int n)
{
    int k[3], r[3], q;
    k[1] = 0, k[2] = 1;
    r[1] = n, r[2] = a;

    while (r[2] != 1)
    {
        r[0] = r[1], r[1] = r[2];
        k[0] = k[1], k[1] = k[2];
        if (r[1] == 0) return -1;
        q = r[0] / r[1];
        r[2] = r[0] % r[1];
        k[2] = (k[0] - q * k[1]) % n;
        if (k[2] < 0) k[2] += n;
    }

    return k[2];
}