最新文章

快速数论变换(NTT)原理及实现

快速数论变换(NTT)原理及实现

概述 上次写了一篇狗屎文章快速傅里叶变换(FFT)原理与实现 | KSkun’ 

[洛谷1919]【模板】A*B Problem升级版(FFT快速傅里叶) 题解

[洛谷1919]【模板】A*B Problem升级版(FFT快速傅里叶) 题解

题目地址:洛谷:【P1919】【模板】A*B Problem升级版(FFT快速傅里叶) & 

快速傅里叶变换(FFT)原理与实现

快速傅里叶变换(FFT)原理与实现

这篇文章写得很烂,而且大部分来自于算导,如有感兴趣者可以直接阅读算导或参考资料QAQ!本文省略大量证明,请各位谨慎食用QAQ!

概述

快速傅里叶变换(FFT)是一种快速计算序列的离散傅里叶变换(DFT)及其逆变换(IDFT)的算法,其复杂度为O(n \log n),比朴素的DFT计算更优。OI中常用它来计算多项式乘法等需要应用卷积的地方。
除了本文章中介绍的数学知识,你还需要知道:

多项式的表达法

系数表达

一个n次多项式可以写成A(x) = \sum_{i=0}^{n} a_ix^i的形式,每一项对应有系数,我们把这些系数组成的向量a = (a_0, a_1, \ldots, a_n)叫做这个多项式的系数表达。

点值表达

一个n次多项式A(x)的点值表达是一个由n+1个点值对组成的集合\{ (x_0, y_0), (x_1, y_1), \ldots, (x_n, y_n) \},其中,x值互不相同,且y_i = A(x_i)
对于n+1个点值能表达唯一的n次多项式的证明,参见[学习笔记] 多项式与快速傅里叶变换(FFT)基础 – rvalue – 博客园
有了点值表达,我们来想想多项式相乘的结果:A(x) \cdot B(x)应该怎么计算。假设A多项式的次数是n,B的次数是m,那么得到的结果一定是一个n+m次多项式,因此我们需要n+m+1个点值才能确定它。因此我们需要把A和B的点值表示扩展为n+m+1个点值,然后对它们的点值表示中的y值对应相乘,结果就是新多项式的点值表示。

单位根的使用与多项式乘法的步骤

我们知道,朴素进行多项式乘法需要O(n^2)时间。如果我们能快速进行系数转点值(DFT)及其逆运算,通过点值表达做乘法只需要O(n)的时间,就能够实现加速这一过程的目的。我们可以通过选取精妙的点值来实现这一加速。
还记得复数中的单位根吗?它有一个折半引理,也就证明了 (\omega_n^k)^2 = \omega_{\frac{n}{2}}^k 。这个结论会在FFT中用到。这决定了我们可以分治实现O(n \log n)的系数转点值及其逆运算。这种加速后的算法被称为FFT。
那么我们计算多项式乘法步骤就是:

  1. 系数转点值,一次FFT,O(n \log n)
  2. 计算乘积的点值,O(n)
  3. 点值转系数,一次FFT的逆运算,O(n \log n)

离散傅里叶变换(DFT)

离散傅里叶变换是一种完成系数表达转换为点值表达的变换。这一变换是对多项式A在n个n次单位根处求点值,也就是说,是在完成下面这样一个事情:
多项式的系数表达给出为a = (a_0, a_1, \ldots, a_n),求y = (y_0, y_1, \ldots, y_n),其中y_i = A(\omega_{n+1}^i)
但是这样看起来其实也是一个O(n^2)的算法啊,教练加速去哪里了?往下看↓

快速傅里叶变换(FFT)

下面的叙述中,n代表多项式的次数界,次数界是指严格大于一个多项式的次数的正整数。为了方便,我们默认n是2的某次幂。
上面我们提到了要把这个问题分治解决,那么就要构造出适合分治解决的形式。首先是怎么构造分治的分支。我们考虑把系数按奇偶性分开,构造下面这样两个多项式
\begin{aligned} A_1(x) &= a_0 + a_2x + a_4x^2 + \cdots \\ A_2(x) &= a_1 + a_3x + a_5x^2 + \cdots \end{aligned}
那么就有A(x) = A_1(x^2) + A_2(x^2)x
那么我们考虑那个折半引理的式子 (\omega_n^k)^2 = \omega_{\frac{n}{2}}^k ,要求的一堆平方实际上变成了\frac{n}{2}个数字,那么我们就分别计算这\frac{n}{2}个数处的点值,然后合并到当前点值上就可以了。这样我们成功实现了缩小问题规模的目的。
具体而言,我们需要对\frac{n-1}{2}之前和之后的点值分别讨论
y_k = A(\omega_n^k) = A_1(\omega_n^{2k}) + \omega_n^kA_2(\omega_n^{2k}) = y^{[1]}_k + \omega_n^ky^{[2]}_k \ (k = 0, 1, \ldots, \frac{n}{2} - 1)
y_{k+\frac{n}{2}} = A(\omega_n^{k+\frac{n}{2}}) = A_1(\omega_n^{2k+n}) + \omega_n^{k+\frac{n}{2}}A_2(\omega_n^{2k+n}) = A_1(\omega_n^{2k}) - \omega_n^kA_2(\omega_n^{2k}) = y^{[1]}_k - \omega_n^ky^{[2]}_k \ (k = 0, 1, \ldots, \frac{n}{2} - 1)
对于递归边界,因为\omega_1^0 = 1,所以y_0 = a_0\omega_1^0 = a_0,直接不做处理即可。

逆变换:点值转系数(IDFT)

我们发现我们可以用下面这张矩阵来表示DFT的过程
\begin{bmatrix} y_0 \\ y_1 \\ \vdots \\ y_{n-1} \end{bmatrix} = \begin{bmatrix} 1 & 1 & 1 & \cdots & 1 \\ 1 & \omega_n & \omega_n^2 & \cdots & \omega_n^{n-1} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \omega_n^{n-1} & \omega_n^{2(n-1)} & \cdots & \omega_n^{(n-1)^2} \end{bmatrix} \begin{bmatrix} a_0 \\ a_1 \\ \vdots \\ a_{n-1} \end{bmatrix}
其实中间那个大块头矩阵就是一个范德蒙矩阵V_n,实际上我们还有一个定理来求这个矩阵的逆矩阵: (V_n^{-1})_{j, k} = \frac{\omega_n^{-kj}}{n} 。那么我们就可以把逆矩阵往点值上一乘,点值就乖乖地变回系数啦~具体而言,变化是这样的a_j = \frac{1}{n} \sum_{k=0}^{n-1} y_k\omega_n^{-kj}。我们发现这个玩意和FFT在做的事情很相似,差就差在一个除以n还有单位根指数上面的负号。负号我们可以通过函数参数传进去,除以n可以在函数外部做,那么直接用FFT的过程不就好了。

FFT的优化方式

蝴蝶操作

我们发现上面的代码中\omega_n^ky_k^{[2]}的计算似乎是重复的,我们可以临时存起来而不是计算两次。

迭代实现

我们来观察FFT的调用树:
FFT调用树
(图片来自《算法导论(第3版)》,参见参考资料)
我们发现我们可以从最下层逐对合并结果,而需要把最下层的顺序处理出来。通过观察我们发现神奇的事:最下层第i个位置的下标刚好是它的二进制倒过来。那么交换一个位置和它二进制倒过来位置的元素即可。二进制倒过来的数可以预处理好。

不要用STL

STL比手写complex慢一些。

代码

本代码可以通过【P3803】【模板】多项式乘法(FFT) – 洛谷一题。

// Code by KSkun, 2018/3
#include <cstdio>
#include <cmath>

#include <algorithm>

typedef long long LL;

inline char fgc() {
    static char buf[100000], *p1 = buf, *p2 = buf;
    return p1 == p2 && (p2 = (p1 = buf) + fread(buf, 1, 100000, stdin), p1 == p2) ? EOF : *p1++;
}

inline LL readint() {
    register LL res = 0, neg = 1;
    char c = fgc();
    while(c < '0' || c > '9') {
        if(c == '-') neg = -1;
        c = fgc();
    }
    while(c >= '0' && c <= '9') {
        res = res * 10 + c - '0';
        c = fgc();
    }
    return res * neg;
}

const int MAXN = 1 << 22;
const double PI = acos(-1);

struct Complex {
    double real, imag;
    Complex(double real = 0, double imag = 0): real(real), imag(imag) {}
    inline Complex operator+(const Complex &rhs) const {
        return Complex(real + rhs.real, imag + rhs.imag);
    }
    inline Complex operator-(const Complex &rhs) const {
        return Complex(real - rhs.real, imag - rhs.imag);
    }
    inline Complex operator*(const Complex &rhs) const {
        return Complex(real * rhs.real - imag * rhs.imag, real * rhs.imag + imag * rhs.real);
    }
    inline Complex& operator*=(const Complex &rhs) {
        Complex old = *this;
        real = old.real * rhs.real - old.imag * rhs.imag;
        imag = old.real * rhs.imag + old.imag * rhs.real;
        return *this;
    }
};

int n, m, len, rev[MAXN];
Complex a[MAXN], b[MAXN];

inline void fft(Complex *arr, int f) {
    for(int i = 0; i < n; i++) {
        if(i < rev[i]) std::swap(arr[i], arr[rev[i]]);
    }
    for(int i = 1; i < n; i <<= 1) {
        Complex wn(cos(PI / i), f * sin(PI / i));
        for(int j = 0; j < n; j += i << 1) {
            Complex w(1, 0);
            for(int k = 0; k < i; k++) {
                Complex x = arr[j + k], y = w * arr[j + k + i];
                arr[j + k] = x + y;
                arr[j + k + i] = x - y;
                w *= wn;
            }
        }
    }
}

int main() {
    n = readint(); m = readint();
    for(int i = 0; i <= n; i++) {
        a[i].real = readint();
    }
    for(int i = 0; i <= m; i++) {
        b[i].real = readint();
    }
    m += n;
    for(n = 1; n <= m; n <<= 1) len++;
    for(int i = 0; i < n; i++) {
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (len - 1));
    }
    fft(a, 1);
    fft(b, 1);
    for(int i = 0; i <= n; i++) {
        a[i] *= b[i];
    }
    fft(a, -1);
    for(int i = 0; i <= m; i++) {
        printf("%d ", int(a[i].real / n + 0.5));
    }
    return 0;
}

参考资料

数学笔记:矩阵、行列式

数学笔记:矩阵、行列式

矩阵(Matrix) 一个的矩阵是一个由行列元素排列成的矩形阵列。元素可以是数字、符号或数 

数学笔记:复数

数学笔记:复数

虚数单位(Imaginary unit) 虚数单位定义为二次方程式的两个解答中的一个解答。 

[CF364D]Ghd 题解

[CF364D]Ghd 题解

题目地址:Codeforces:Problem – 364D – Codeforces、洛谷:【CF364D】Ghd – 洛谷

题目描述

John Doe offered his sister Jane Doe find the gcd of some set of numbers a.
Gcd is a positive integer g, such that all number from the set are evenly divisible by g and there isn’t such g’ (g’ > g), that all numbers of the set are evenly divisible by g’.
Unfortunately Jane couldn’t cope with the task and John offered her to find the ghd of the same subset of numbers.
Ghd is a positive integer g, such that at least half of numbers from the set are evenly divisible by g and there isn’t such g’ (g’ > g) that at least half of the numbers from the set are evenly divisible by g’.
Jane coped with the task for two hours. Please try it, too.
有一个数列,求最大的数,使得该数是数列中一半以上的数的因数。

输入输出格式

输入格式:
The first line contains an integer n (1 ≤ n ≤ 10^6) showing how many numbers are in set a. The second line contains space-separated integers a1, a2, …, an (1 ≤ ai ≤ 10^12). Please note, that given set can contain equal numbers.
Please, do not write the %lld specifier to read or write 64-bit integers in С++. It is preferred to use the %I64d specifier.

输出格式:
Print a single integer g — the Ghd of set a.

输入输出样例

输入样例#1:

6
6 2 3 4 5 6

输出样例#1:

3

输入样例#2:

5
5 5 6 10 15

输出样例#2:

5

题解

参考资料:CFR 364 D. Ghd ( Random, Math ) – 0w1
随机化的题目,直接求显然不好办,n的数据范围看着就吓人,直观感觉像O(n \log n)
我们知道这个数字一定是这个数列中某个数的因数,因此我们需要选择一个数求它与其他数的最大公因数。这些公因数中的一个就是答案。求一遍公因数是O(n \log n)的。
对于求出来的公因数,我们去从大到小找一个会成为超过一半数的因数的数字。具体做法是,选择一个因数,去找比它大的因数,如果它能整除大因数,说明大因数对应的数字也可以被这个数整除,应当把加到这个数的计数上。这个过程直观看上去是O(n^2)的,但是实际上我们不会对比当前最优解还小的因数计算,并且当计数超过n的一半的时候就可以停止计数并计入答案,在这样的处理下我们可以让这个计数的时间变得可接受。
答案肯定是这些数字中超过一半数的因数,因此每一次随机找到正确解的概率超过0.5。我们考虑重复找解的过程多次,这里我使用10次(超过10次似乎会TLE 18)。

代码

// Code by KSkun, 2018/3
#include <cstdio>
#include <ctime>
#include <cstdlib>

#include <map>

typedef long long LL;

inline char fgc() {
    static char buf[100000], *p1 = buf, *p2 = buf;
    return p1 == p2 && (p2 = (p1 = buf) + fread(buf, 1, 100000, stdin), p1 == p2) ? EOF 
        : *p1++;
}

inline LL readint() {
    register LL res = 0, neg = 1;
    char c = fgc();
    while(c < '0' || c > '9') {
        if(c == '-') neg = -1;
        c = fgc();
    }
    while(c >= '0' && c <= '9') {
        res = res * 10 + c - '0';
        c = fgc();
    }
    return res * neg;
}

const int MAXN = 1000005;

inline LL gcd(LL a, LL b) {
    LL t;
    while(b) {
        t = a % b;
        a = b;
        b = t;
    }
    return a;
}

int n;
LL a[MAXN];

int main() {
    srand(time(NULL));
    n = readint();
    for(int i = 1; i <= n; i++) {
        a[i] = readint();
    }
    LL ans = 1;
    for(int rep = 1; rep <= 10; rep++) { // bad rand
        int rnd = (rand() * RAND_MAX + rand()) % n + 1;
        std::map<LL, int> fact;
        for(int i = 1; i <= n; i++) {
            LL t = gcd(a[i], a[rnd]);
            if(!fact.count(t)) fact[t] = 1;
            else fact[t]++;
        }
        std::map<LL, int>::iterator it = fact.end();
        do {
            it--;
            if((*it).first <= ans) continue;
            int cnt = 0;
            for(std::map<LL, int>::iterator it1 = it; 
                it1 != fact.end() && cnt << 1 < n; it1++) {
                if(!((*it1).first % (*it).first)) {
                    cnt += (*it1).second;
                }
            }
            if(cnt << 1 >= n) ans = (*it).first;
        } while(it != fact.begin());
    }
    printf("%I64d", ans);
    return 0;
}
[WF2015]Cutting Cheese 题解

[WF2015]Cutting Cheese 题解

题目地址:UVa:UVa Online Judge、洛谷:【UVA1712】Cutting 

[CF16E]Fish 题解

[CF16E]Fish 题解

题目地址:Codeforces:Problem – 16E – C 

[ZJOI2015]地震后的幻想乡 题解

[ZJOI2015]地震后的幻想乡 题解

题目地址:洛谷:【P3343】[ZJOI2015]地震后的幻想乡 – 洛谷、BZOJ:Problem 3925. — [Zjoi2015]地震后的幻想乡

题目描述

傲娇少女幽香是一个很萌很萌的妹子,而且她非常非常地有爱心,很喜欢为幻想乡的人们做一些自己力所能及的事情来帮助他们。 这不,幻想乡突然发生了地震,所有的道路都崩塌了。现在的首要任务是尽快让幻想乡的交通体系重新建立起来。
幻想乡一共有n个地方,那么最快的方法当然是修复n-1条道路将这n个地方都连接起来。 幻想乡这n个地方本来是连通的,一共有m条边。现在这m条边由于地震的关系,全部都毁坏掉了。每条边都有一个修复它需要花费的时间,第i条边所需要的时间为ei。地震发生以后,由于幽香是一位人生经验丰富,见得多了的长者,她根据以前的经验,知道每次地震以后,每个ei会是一个0到1之间均匀分布的随机实数。并且所有ei都是完全独立的。
现在幽香要出发去帮忙修复道路了,她可以使用一个神奇的大魔法,能够选择需要的那n-1条边,同时开始修复,那么修复完成的时间就是这n-1条边的ei的最大值。当然幽香会先使用一个更加神奇的大魔法来观察出每条边ei的值,然后再选择完成时间最小的方案。 幽香在走之前,她想知道修复完成的时间的期望是多少呢?

输入输出格式

输入格式:
第一行两个数n,m,表示地方的数量和边的数量。其中点从1到n标号。 接下来m行,每行两个数a,b,表示点a和点b之间原来有一条边。 这个图不会有重边和自环。

输出格式:
一行输出答案,四舍五入保留6位小数。

输入输出样例

输入样例#1:

5 4
1 2
1 5
4 3
5 3

输出样例#1:

0.800000

说明

提示:
(以下内容与题意无关,对于解题也不是必要的。)
对于n个[0,1]之间的随机变量x1,x2,…,xn,第k小的那个的期望值是k/(n+1)。
样例解释:
对于第一个样例,由于只有4条边,幽香显然只能选择这4条,那么答案就是4条边的ei中最大的数的期望,由提示中的内容,可知答案为0.8。
数据范围:
对于所有数据:n<=10, m<=n(n-1)/2, n,m>=1。
对于15%的数据:n<=3。
另有15%的数据:n<=10, m=n。
另有10%的数据:n<=10, m=n(n-1)/2。
另有20%的数据:n<=5。
另有20%的数据:n<=8。

题解

题目要求最小生成树上的最大边权Y的期望,因此我们只需要关心这个最大边权在所有边权中的排名,我们发现有下面的式子
\mathrm{E}(Y) = \sum_{i=1}^n \mathrm{E}(i)\mathrm{P}(i) = \sum_{i=1}^n \frac{i}{m+1} \cdot \mathrm{P}(i) = \frac{\sum_{i=1}^n i\mathrm{P}(i)}{m+1} = \frac{\mathrm{E}(Z)}{m+1}
最后我们发现要求的就是最小生成树上最大边权排名的期望值Z。考虑怎么求这个Z,设L为同定义的随机变量,有
\mathrm{E}(Z) = \sum_{i=1}^m i\mathrm{P}(L=i) = \sum_{i=1}^m \mathrm{P}(L \geq i)
如果说最大边排名大于等于i不好求,我们就正难则反,求排名小于i的边集无法构成生成树的概率,由于所有的边权都是随机分布的变量,这个概率可以转化为从边集中随机选择i条边无法构成生成树的概率,我们考虑计算出选择i条边无法构成生成树的方案数,进而计算概率。
通过观察题目数据范围,我们发现n的范围很适合状压。考虑状压DP,用f[S][i]表示在S点集构成的子图中,选i条边使得该点集不连通的方案数,这个状态似乎没法转移,但是我们考虑与其意义互补的量:g[S][i]表示在S点集构成的子图中,选i条边使得该点集连通的方案数,显然有
f[S][i] + g[S][i] = \mathrm{C}_{ecnt[S]}^i
ecnt表示该点集构成子图内的边数。如果知道了其中一个,我们就可以知道另外一个。那么怎么转移成了问题,我们考虑固定S中的某个点,对于S这个点集的一个包含定点的真子集T,只要使T连通但\complement_S T与T之间没有连边就能保证S不连通了,而枚举每一个T可以实现对S不连通情况的遍历。最后我们的转移式就是
f[S][i] = \sum_{T \subsetneqq S, U \in T} \sum_{j=0}^{ecnt[T]} g[T][j] \cdot \mathrm{C}_{ecnt[\complement_S T]}^{i-j}
最后,我们获得了 \mathrm{E}(Z) = \sum_{i=0}^{n-1} f[\mathbb{U}][i] ,直接计算答案即可。
注意爆int。

附加内容:证明提示结论

需要证明的内容:对于n个[0,1]之间的随机变量x1,x2,…,xn,第k小的那个的期望值是k/(n+1)。
我们首先令x_1为第k小值,我们现在来求它的期望。
对于这种情况,我们显然有其中k-1个值比它小,n-k个值比它大,我们分别计算概率,乘起来,就得到了概率密度,对其积分就是x_1为第k小值时x_1的期望值,即\mathrm{E}(x_1, x_1是第k小数)
\begin{aligned} & \int_0^1 x_1 \cdot x_1^{k-1} (1-x_1)^{n-k} \mathrm{d}x_1 \\ = & \int_0^1 x_1^k (1-x_1)^{n-k} \mathrm{d}x_1 \\ = & \mathrm{B}(k+1, n-k+1) \\ = & \frac{k!(n-k)!}{(n+1)!} \end{aligned}
我们知道了这个,但是要求的是E(第k小数),即\mathrm{E}(x_1|x_1是第k小数)。我们记x_1是第k小数这一事件为A,根据下面的关系
\mathrm{E}(x_1|A) = \frac{\mathrm{E}(x_1, A)}{\mathrm{P}(A)}
我们又容易知道
\mathrm{P}(A) = \frac{1}{n \cdot \mathrm{C}_{n-1}^{k-1}}
直接就可以算出来了,结果就是
\mathrm{E}(x_1|A) = \frac{k}{n+1}

代码

// Code by KSkun, 2018/3
#include <cstdio>

typedef long long LL;

inline char fgc() {
    static char buf[100000], *p1 = buf, *p2 = buf;
    return p1 == p2 && (p2 = (p1 = buf) + fread(buf, 1, 100000, stdin), p1 == p2) ? EOF : *p1++;
}

inline LL readint() {
    register LL res = 0, neg = 1;
    char c = fgc();
    while(c < '0' || c > '9') {
        if(c == '-') neg = -1;
        c = fgc();
    }
    while(c >= '0' && c <= '9') {
        res = res * 10 + c - '0';
        c = fgc();
    }
    return res * neg;
}

const int MAXN = 13, MAXM = 50;

int n, m, ut, vt;
LL gra[MAXN], cnt[1 << MAXN], ecnt[1 << MAXN], f[1 << MAXN][MAXM], g[1 << MAXN][MAXM];

LL C[MAXM][MAXM];

inline void calc() {
    for(int i = 0; i <= m; i++) {
        C[i][0] = 1;
        for(int j = 1; j <= i; j++) {
            C[i][j] = C[i - 1][j - 1] + C[i - 1][j];
        }
    }
}

int main() {
    n = readint(); m = readint();
    calc();
    for(int i = 0; i < m; i++) {
        ut = readint(); vt = readint();
        gra[ut] |= (1 << (vt - 1));
        gra[vt] |= (1 << (ut - 1));
    }
    for(int i = 0; i < 1 << n; i++) {
        cnt[i] = cnt[i >> 1] + (i & 1);
    }
    for(int i = 0; i < 1 << n; i++) {
        for(int j = 1; j <= n; j++) {
            if(i & (1 << (j - 1))) {
                ecnt[i] += cnt[gra[j] & i];
            }
        }
        ecnt[i] >>= 1;
    }
    for(int i = 0; i < 1 << n; i++) {
        if(cnt[i] == 1) {
            g[i][0] = 1;
            continue;
        }
        int t = i & (-i);
        for(int j = (i - 1) & i; j; j = (j - 1) & i) {
            if(j & t) {
                for(int k1 = 0; k1 <= ecnt[j]; k1++) {
                    for(int k2 = 0; k2 <= ecnt[i ^ j]; k2++) {
                        f[i][k1 + k2] += g[j][k1] * C[ecnt[i ^ j]][k2];
                    }
                }
            }
        }
        for(int j = 0; j <= ecnt[i]; j++) {
            g[i][j] = C[ecnt[i]][j] - f[i][j];
        }
    }
    double ans = 0;
    for(int i = 0; i <= m; i++) {
        ans += double(f[(1 << n) - 1][i]) / C[m][i];
    }
    printf("%.6lf", ans / (m + 1));
    return 0;
}
数学笔记:概率、期望

数学笔记:概率、期望

期望(Expectation)的定义 如果X是在概率空间中的随机变量,那么它的期望值的定义