快速沃尔什变换

新的卷积

卷积的东西都很好玩呀。有时间一定好好学泛函.jpg

前面提到了,离散傅里叶变换的多项式卷积是对于每个 k 获得 i+j=k\sum a_i \times a_j,而在数论中比较多的狄利克雷卷积是对于每个 k 获得 i \times j = k\sum a_i \times a_j。在实际题目中,我们可能会遇到这样的一类卷积:

对于每个 k,获得 i \otimes j = k\sum a_i \times a_j。其中 \otimes 是任意按位逻辑运算,如 \&, |, \oplus 等。

很显然,我们如果按照最朴素的枚举方法:先枚举 i,再获得可行的 k,复杂度是 O(n^2) 朝上走的。如果 n \ge 10^5 不就完蛋了吗。按照这个架势,我们应该是要找复杂度为 O(n\log n)O(n \log \log n)O(n \log^2 n) 的算法。

很幸运的是,我们可以通过分治来完成这样的一件事情。当然了,由于位运算的值域是 2^n 级别的,我们也要像 FFT 那样将数列长度补齐。

快速沃尔什变换

注意到,位运算是有位独立性的。每一位之间并不会互相影响。那么,我们每次只考虑某一位?

我们令 A_{xy} = (A_{0y}, A_{1y})B_{xy} = (B_{0y}, B_{1y}),行向量分块,其中 y 是任意01序列,那么我们想要
C_{xy} = A_{xy} * B_{xy} = \left(\sum_{a \otimes b = 0} A_{ay} * B_{by}\ , \sum_{a \otimes b = 1} A_{ay} * B_{by} \right)
由于逻辑运算是确定的,我们可以很容易的计算出 \otimes 的真值表,并且后续又变成了规模更小的卷积,就是一个分治的过程啦。

假如构造了一个变换 \mathfrak{T}(A),满足 \mathfrak{T} (C) = \mathfrak{T} (A) \times \mathfrak{T} (B),并且可以很短时间内求得 \mathfrak{T}(A)\mathfrak{T}^{-1}(A),就可以像FFT一样卷积了。

Fast Walsh-Hadamard Transform,这个神奇的变换,我们分别称变换与逆变换为 FWT(A)UFWT(A),那么

对于与运算 \& 来说,有

FWT(A) = \left(\ FWT(A_0) + FWT(A_1)\ ,\ FWT(A_1)\ \right)
UFWT(A) = \left(\ UFWT(A_0) – UFWT(A_1)\ ,\ UFWT(A_1)\ \right)

并且

FWT(A)[i] = \sum_{i \& j = i} A[j]

void FWT(lld X[], int l, int r, int dft = 1)
{
    if (l == r) return;
    int m = (l+r)>>1;
    FWT(X, l, m, dft);
    FWT(X, m+1, r, dft);
    for (int i = 0; i <= m-l; i++)
        X[l+i] += X[m+1+i] * dft;
}

对于或运算 | 来说,有

FWT(A) = \left(\ FWT(A_0)\ ,\ FWT(A_0) + FWT(A_1)\ \right)
UFWT(A) = \left(\ UFWT(A_0)\ ,\ UFWT(A_1) – UFWT(A_0)\ \right)

并且

FWT(A)[i] = \sum_{i | j = i} A[j]

void FWT(lld X[], int l, int r, int dft = 1)
{
    if (l == r) return;
    int m = (l+r)>>1;
    FWT(X, l, m, dft);
    FWT(X, m+1, r, dft);
    for (int i = 0; i <= m-l; i++)
        X[m+1+i] += X[l+i] * dft;
}

异或运算 \oplus
FWT(A) = \left(\ FWT(A_0) + FWT(A_1)\ ,\ FWT(A_0) – FWT(A_1)\ \right)
UFWT(A) = \left(\ UFWT(\frac{A_0 + A_1}2)\ ,\ UFWT(\frac{A_0-A_1}{2})\ \right)

暂时没有找到比较明显的规律

void FWT(lld X[], int l, int r, int dft)
{
    if (l == r) return;
    int m = (l+r)>>1;
    FWT(X, l, m, dft);
    FWT(X, m+1, r, dft);

    for (int i = 0; i <= m-l; i++)
    {
        lld x = X[l+i], y = X[m+1+i];
        X[l+i] = dft == 1 ? x+y : (x+y)/2;
        X[m+1+i] = dft == 1 ? x-y : (x-y)/2;
    }
}

例子

链接:https://ac.nowcoder.com/acm/contest/295/H
来源:牛客网

Niuniu likes playing games. He has n piles of stones. The i-th pile has ai stones. He wants to play with his good friend, UinUin. Niuniu can choose some piles out of the n piles. They will play with the chosen piles of stones. UinUin takes the first move. They take turns removing at least one stone from one chosen pile. The player who removes the last stone from the chosen piles wins the game. Niuniu wants to choose the maximum number of piles so that he can make sure he wins the game. Can you help Niuniu choose the piles? 

给出 n 个数字,记为 a_0, a_1, a_2, …, a_{n-1},且 a_i \le 5\times10^5,求 max|A|, \oplus_{a\in A} a = 0

不妨记 \oplus_{i=0}^{n-1}a_i = sum,那么我们构造一个最小的堆,使得其异或和为 sum 即可。

我们只要记录可以被构造出来的堆值,就可以在第一次这个堆被构造出来时停止。由于 x \oplus x = 0,所以我们直接用已经构造出来的集合和原始给出的集合进行逻辑运算卷积。因为不需要记录方案数,所以直接置1即可。

B_{k}^n = \left[\left(\sum_{i \oplus j = k} A_{i} B_{j}^{n-1}\right) > 0 \right],相当于记录了,有 B 的堆和 A 原来堆,异或后能构造的值的集合。显然这个数字不会超过 2^{19},所以做19次卷积即可。

#include <cstring>
#include <cstdio>
typedef long long lld;
const int MAX = 524287;
lld a[MAX+1], b[MAX+1] = {1};

void FWT(lld X[], int l, int r, int dft)
{
    if (l == r) return;
    int m = (l+r)>>1;
    FWT(X, l, m, dft);
    FWT(X, m+1, r, dft);

    for (int i = 0; i <= m-l; i++)
    {
        lld x = X[l+i], y = X[m+1+i];
        X[l+i] = dft == 1 ? x+y : (x+y)/2;
        X[m+1+i] = dft == 1 ? x-y : (x-y)/2;
    }
}

int main()
{
    int n, tmp, sum = 0; scanf("%d", &n);
    for (int i = 0; i < n; i++)
    {
        scanf("%d", &tmp);
        a[tmp]++;
        sum ^= tmp;
    }

    FWT(a, 0, MAX, 1);
    int max_ans = 0;
    while (!b[sum])
    {
        if (max_ans++ == 19) break;
        FWT(b, 0, MAX, 1);
        for (int j = 0; j <= MAX; j++)
            b[j] = a[j] * b[j];
        FWT(b, 0, MAX, -1);
        for (int j = 0; j <= MAX; j++)
            b[j] = !!b[j];
    }

    printf("%d\n", n - max_ans);
    return 0;
}

Polya置换

置换

M 是一个非空的有限集合,M 的一个一对一变换称为一个置换。

M 的元素为 a_1, a_2, \cdots, a_n,则 M 的置换 \sigma 可以简记为

\sigma = \left(
\begin{array}{cccc}
a_1 & a_2 & \cdots & a_n\\
b_1 & b_2 & \cdots & b_n
\end{array}
\right)

其中 b_i = \sigma(a_i)。如果

\sigma = \left(
\begin{array}{cccc}
a_1 & a_2 & \cdots & a_{n-1} & a_n\\
a_2 & a_3 & \cdots & a_n & a_1
\end{array}
\right)

则简记 \sigma = (a_1 a_2 \cdots a_n)

\sigma\tau 是两个不相杂的轮换,则其乘法适合交换律,即 \sigma\tau = \tau\sigma。可以证明,任意置换 \sigma 恰有一法写成不相交轮换的乘积。

可以分解为奇数个对换之积的置换称为奇置换,可以分解为偶数个对换之积的置换称为偶置换。

Burnside 引理

Gn 元集 S = \lbrace 1, 2, \cdots, n \rbrace 上的置换群 G = \lbrace \sigma_1, \sigma_2, \cdots, \sigma_n \rbrace,把每个 \sigma_i 都表示成不相杂的轮换的乘积,令 \lambda_k(\sigma_i) 表示置换 \sigma_ik 阶轮换的个数,则 GS 上诱导出的等价关系将 S 划分为不同等价类的个数为

L = \frac{1}{|G|} \sum_{j=1}^r \lambda_1(\sigma_j)

其中 \lambda_1(\sigma) 为置换 \sigma 中的不动点(即1阶轮换)的个数。

Polya 定理

Hn 个对象上的一个置换群,用 m 种颜色给这 n 个对象涂色,一个对象涂任意种颜色,则在置换群 H 的作用下不等价的方案数为

L = \frac{1}{|H|} \sum_{\tau \in H} m^{\lambda(\tau)}

其中 \lambda(\tau) 表示置换 \tau 表示为不相杂轮换乘积的形式中轮换的个数。

例如,对于置换群H’ = \lbrace (1)(2)(3), (1 2 3), (1 3 2), (1)(2 3), (2)(1 3), (3)(1 2) \rbrace,不等价的方案数为L’ = \frac{1}{6} \left[3^3 + 2 \times 3^1 + 3 \times 3^2 \right]。记颜色为 \lbrace r, g, b \rbrace,则其生成函数为L(r,g,b) = \frac{1}{6} \left[(r+g+b)^3 + 2 \times (r^3+b^3+g^3) + 3 \times (r+g+b)(r^2+b^2+g^2) \right],其中 r^xg^yb^z 前的系数表示将对应颜色染对应个数时的方案数。

组合数学知识计算。

常见模型举例

圆环的旋转与反射

例题:POJ 2409 / POJ 1286

在轴对称的作用下。奇数阶环的对称轴一定通过其中的一个珠子,那么会形成 (1)(2, n)(3, n-1)\cdots 这样的置换,总共 n 个,不相杂轮换 \frac{n+1}{2} 个。偶数阶环的对称轴可能经过珠子,也可能不经过珠子,那么会形成 (1,n)(2,n-1)(3,n-2)\cdots 型置换 \frac{n}{2} 个、(1)(2,n)(3,n-1)\cdots(\frac{n}{2}+1) 型置换 \frac{n}{2} 个。

在中心对称的作用下,考虑珠子旋转 k\in[0,n-1] 个单位时的情况,不相杂轮换有 gcd(k,n) 个。对于同一个 k,又可以由欧拉函数求出相似者,于是枚举 n 所有的因子即可,对于每个因子 d 都有 \varphi(\frac{n}{d}) 个共轭的置换。

while (!factors.empty()) {
    int d = factors.back();
    sum += phi[n/d] * fpow(c, d);
    factors.pop_back();
}

if (n & 1) {
    sum += n * fpow(c, (n+1)>>1);
} else {
    sum += (n>>1) * fpow(c, n>>1);
    sum += (n>>1) * fpow(c, 1+(n>>1));
}

sum /= n*2;

正方形的旋转与反射

例题:HDU 1812

对于 N\times N 的正方形,有1个 n^2 的恒等置换,有1个 \lceil \frac{n^2}{2} \rceil 的 180° 旋转,有2个 \lceil \frac{n^2}{4} \rceil 的 90° 旋转,有2个 \lceil \frac{n}{2} \rceil n 的纵向横向反射,有2个 \frac{n(n+1)}{2} 的斜反射。

BigInteger sum = BigInteger.valueOf(0);
sum = sum.add(colors.pow(n*n));
sum = sum.add(colors.pow((n*n+1)/2));
sum = sum.add(colors.pow((n*n+3)/4).shiftLeft(1));
sum = sum.add(colors.pow((n+1)/2*n).shiftLeft(1));
sum = sum.add(colors.pow(n*(n+1)/2).shiftLeft(1));
System.out.println(sum.shiftRight(3));

正方体的旋转

1x静止、4x2x体对角线、6x对棱中心连线、3x3x对立面中心连线

lld solve2(int i = 0, int s = 6) {
    if (i >= 6) return 1;
    if (cnt[i] % 2) return 0;
    return C[s][cnt[i]/2] * solve2(i+1, s-(cnt[i]/2));
} // 提取 (a2+b2+c2+d2+e2+f2)^6 的某项系数

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];
}

Xamarin.Forms 安卓沉浸式状态栏

一开始用MasterDetailPage做DrawerLayout的时候发现没办法透过状态栏,记录踩坑全过程。

values/style.xmlMainTheme.Base 里面加入:

<item name="android:windowTranslucentStatus">true</item>

然后状态栏是透过去了,但是Toolbar也跟着上去了……

values/style.xmlMainTheme 里面加入:

<item name="actionBarSize">@dimen/action_bar_default_height_material_overlay</item>

action_bar_default_height_material_overlay竖屏76.0dip,横屏68.0dip

layout/Toolbar.axml里加入两个属性

<android.support.v7.widget.Toolbar
...
android:layout_height="?attr/actionBarSize"
android:minHeight="?attr/actionBarSize"
android:paddingTop="12dp"
app:titleMarginTop="24dp"
...
/>

大致是完成了。

后期发现ActionMode的状态栏会变宽,后来这样解决

values/style.xml里面加入:

<style name="MainTheme.ActionMode" parent="Widget.AppCompat.ActionMode">
<item name="height">@dimen/action_mode_default_height_material_overlay</item>
<item name="actionBarSize">@dimen/action_mode_default_height_material_overlay</item>
</style>

values/style.xmlMainTheme.Base里面加入:

<item name="actionModeStyle">@style/MainTheme.ActionMode</item>

action_mode_default_height_material_overlay竖屏56.0dip,横屏48.0dip

业界毒瘤(误)OpenLitespeed

:confused: 前段时间发现伪静态设置一直都不对,后来发现OLS根本没有读取伪静态规则,后台INFO级别没有记录任何。先挂在这里。

但是突然又莫名其妙可以伪静态了,可能是要把日志级别调整为9。后来发现要restart才行,reload还是没有用的。

–2017/7/17更新–

再带点今晚读OLS的admin网站源码的感受吧。

OLS的配置文件是一种很奇特的格式,

继续阅读“业界毒瘤(误)OpenLitespeed”