这篇笔记写于2019/2/15，原供稿班级公众号，在某人决定博客重更之后重新整理发上来~2019/12/8上传

[TOC]

为什么会有这篇文章？

最近网上也都在吹爆《流浪地球》这部科幻片啊……

然后我自己前几天也慕名去电影院看了一下，感觉特效之类的着实过瘾。但是作为一部定位为硬科幻的科幻片，除去画面以外，科学的严谨性还是很重要的。但是这部片子里面的一些科学设定在一个KSP玩家看来似乎……

还是有些槽点呢……当然瑕不掩瑜。这篇文章就稍微指出一点我发现的《流浪地球》当中的一些Bug，这些Bug主要是轨道力学上的，也就当给同学们做一个科普吧。

自己也是半瓶水，所以，欢迎指正。

定性吐槽

电影里面的运动学乍一看是非常符合直觉的。比如说：

• 地球的发动机一关，整个地球马上就往木星上靠过去了，轨道一下子拐了一个弯
• 空间站发动机一开，整个空间站马上就直线飞离地球了，轨道上看得一清二楚
• 要引爆空间站，就直接对准目标直线冲刺过去了

似乎没啥问题咯？

Uhmmmmmm……

这么说的话以下二位的棺材板就要压不住了：

（牛顿）

（开普勒）

为啥嘞？

开普勒第一定律

回忆一下我们第一学期学的牛顿第一定律：

牛顿第一定律：假若某物体所受合外力为零，则该物体的运动速度不变，即保持静止或做匀速直线运动。

如果所受合力不是零，而是一直受到一个定质点的引力作用，物体的运动状态又会是怎样呢？

开普勒结合牛顿的万有引力定律对牛顿第一定律进行了拓展，得出了开普勒第一定律：

开普勒第一定律：每一个行星都沿各自的椭圆轨道环绕太阳，而太阳则处在椭圆的一个焦点中。

开普勒定律刻画的行星在只受太阳引力作用下的稳定运动轨道，但是真的只有椭圆这一种吗？

其实定律实际推导出的轨道是圆锥曲线，即圆，椭圆，抛物线或双曲线中的任意一种，只是对于普通行星来说，后面两种轨道意味着飞掠而不是围绕，不符合实际，因此开普勒在定律表述里面便只说“椭圆”。

那知道这个定律有啥用呢？

发动机不需要一直点火！

首先，如果地球的行星发动机的推力够足，理论上是用不着像电影里面一直点火的，我们只需要在地球轨道附近把地球加速到足够高速然后熄火，那么根据开普勒定律，地球自己会沿一条圆锥曲线被甩离太阳，并在木星附近被木星捕获，沿以木星为焦点的抛物线或双曲线掠过木星进行引力弹弓加速，然后离开木星被甩出太阳系！这段过程中理想情况下地球只受星体引力影响，不需要发动机一直点火，除了进行一定的轨道修正！‘

这个操作有没有实例呢？

且看旅行者二号：

它在被火箭加速到足够高的速度之后就放飞自我了，之后走的轨道是什么样的呢？

（绿色为木星轨道）

在离开地球的十多年里，旅行者2号几乎没有启用过自己的轨道发动机，包括在飞掠火星木星等星体进行引力加速的时候也是如此。如果它是靠电影中的一直加速才飞到那么远的话，试问这么小的机体燃料放哪里呢？

所以看到地球上的行星发动机一直燃烧了十七年，Emmm...编剧开心就好。（当然不排除是加速度过低导致需要持续加速的情况，这种情况下轨道确实比较难算，出偏差也是正常的，因此这其实不算一个特别大的Bug）

熄火/点火之后的轨道不科学！

看过电影的同学应该记得，在行星发动机一半熄火的时候，地球马上就拐了一个弯，从飞掠木星变成栽进木星了，就好像有什么东西让它大幅度减速了一样。编剧在这里似乎混淆了速度和加速度的概念，就算地球发动机全部停机，地球几乎只受木星引力作用，地球也会沿着熄火前最后一刻所处的抛物线轨道继续运行：

（图找的不是熄火的那个场景，只做演示用）

也就是说，地球会继续按照“原轨道”的红线运行，而不是像影片中的那样当场轨道一拐拐到蓝色轨道上去了。

也就是说，其实发生影片中所谓“危机”的概率会小非常多。

什么？你说抛物线和双曲线本来就有一个拐点？那和影片当中的形状和位置都对不上好吗?

同样的，空间站准备开溜时候的轨道也注定不是一条直线，它同时受着地球和木星的引力影响，做的是三体运动，轨迹怎么会那么漂亮呢？

因为这样，空间站要开到地木之间引爆，也绝非一个“前进三！”的直线冲刺能够与解决的，复杂着呢！

定量吐槽

小马在看完电影的当晚估算了一下：

然后第二天就有同学告诉我我仿佛漏算了什么……确实是这样的，那这个定量的估算到底是怎么来的呢？

Tsiolkovsky公式

在1903年，俄罗斯科学家齐奥尔科夫斯基推导出了著名的火箭公式：

（详细推导见这里）

其中表示火箭工质的喷出速度，表示火箭的初始质量，表示火箭的燃烧后的末质量，表示火箭在发动机加速过程当中所获得的速度增量。

什么叫工质喷出速度呢？“工质”就是火箭发动机喷嘴喷出的物质的代称，也是火箭发动机推力的来源：发动机给工质一个力将工质高速推离火箭，反作用力同时给火箭一个巨大的动力。

如果假设地球的质量是，石头的质量是，那么在我们的情况下，公式可以写作：

我们可以非常简单地估算出（或者借助网络查阅到）的值，因此只要我们知道了和，我们就可以由上式反解。从而算出作为燃料的石头的重量！

计算

直接逃逸

我们首先计算如果不借助木星的引力弹弓效应需要的最小大概是多少：

首先我们计算在地球轨道上相对太阳的逃逸速度，取太阳的标准重力参数，取地日平均距离，若假设轨道偏心率极小可以忽略不计，则根据机械能守恒逃逸速度可以计算如下：

我们当然不用取，因为地球绕太阳公转本来就有一定的线速度，地球的轨道速度如下： 如果行星发动机顺着地球的公转方向加速，则就够了。

借助木星

我们接下来估算一下如果借助木星的引力弹弓效应加速大约需要多少。

因为没有给定地木之间的位置关系，直接算引力弹弓难度极大不会，但是我们可以估算一下的下界：我们假设地球只需要在地球轨道和木星轨道之间执行一次霍曼转移来进行交会。即让地球的椭圆轨道与木星轨道内切。基于Vis-viva公式，取木星轨道半径为近日点的，则在地球轨道上加速所需要的为：

（假设地木轨道共面，不然还要修正轨道倾角）

也就是说，就算地球不是想借助木星加速，只是想碰到木星的轨道，都需要这么多的，要借助木星加速一般轨道会更大一点，因此我们估算出来的是的下界。

发动机种类

我们已经算出了的大小，接下来让我们估一下的大小，我们同时会计算出来所需要燃料的质量。

在下文中，取地球质量为，岩石密度为。

目前科技

目前人类能够造出的最高（比冲最高）的化学燃料发动机大概是液氢液氧发动机，。（根据资料比冲最高的发动机其实是欧空局的DS4G（电推），，可惜推力只有，还贼费电。）

因此，在直接逃逸的情况下我们有：

取相关数据带入，得……

Ummm……不用算了，要把地球挖空的。

在借助木星引力加速的情况下其实也好不到哪里去，计算可以得到。

这个好一点，思考一下人类靠的地球可以干嘛吧，这还没算上后续的加速和减速，估计流浪到最后变成小行星了。细思极恐

以上计算其实反映出一点：因为逆用公式算出的只是燃油质量和火箭质量的比值，这个比值无论大到地球小到火箭都适用，也就是说大家看到巨大的火箭一般都有以上的质量是燃料，而剩下的质量包括了外壳，发动机，管线，最后才是载荷。

外星科技

既然影片当中的发动机叫“重聚变发动机”（先不纠结为啥重元素可以不吸能聚变），那推进方式肯定不是垃圾的氢氧机可以比拟的。这样吧，假设工质可以加速到光速（顺便压住了爱因斯坦的棺材板），取。会怎么样呢？（关于相对论效应：齐奥尔科夫斯基公式只有在被推进物接近光速（Relativistic Rocket）的时候才需要修正，对于工质速度接近光速的情况则不用。）

在直接加速的情况下，我们有。似乎还可以接受哦？

体积是多少呢？

也就是说所需要的燃料可以装满一个边长的立方体！作为参考，国际空间站轨道高度为左右。

如果嫌这个不够直接，取陆地地壳平均厚度，则需要挖掉整个面积为大小的土地，而全球陆地面积才只有！这意味着：

• 大概要把全球陆地都下挖左右，比环球金融中心高一点。
• 把边长为的一块正方形地壳挖光，作为参考，这约占我国国土的四分之一。

利用木星加速需要多少燃料呢？带入可得，可以装满边长为的立方体。需要把全球陆地都下挖左右，如果打算把地壳挖光，则需要挖去的地壳面积为我国国土的五分之一左右。

Emm……Emm……

质能转化角度

换种思想，假设人类已经解锁了一个无敌恐怖的科技点，即把燃料的所有静能量转化为地球的机械能，那会怎么样呢？

根据著名的质能方程式，我们可以列式：

如果直接逃逸，解得燃料质量为。获取这么多燃料需要把整个上海市区的地壳挖空。

如果借助木星。解得燃料质量为。获取这么多燃料也需要把整个上海中心城区的地壳挖空。

这样开起来稍微可行一点了吧……似乎开掘一万座地下城挖剩的石料就够用了？

但是这个视角的问题在于我们的假设，如果静能量全部转化成机械能，则：

• 燃料必须在相互反应当中全部湮灭，作为对比，核聚变和核裂变当中湮灭的仅仅是反应物和产物之间的微小质量差（可能是每个原子一两个中子/质子），目前能够达到全部湮灭的方法只有把当中一半的燃料转化为反物质。
• 湮灭产生的巨大能量必须全部收集，不能产生耗散。
• 收集的能量必须全部转化成地球的机械能，这意味着不能有喷出的工质，因为这相当于有能量转化成工质机械能，同时整个体系必须不能发热，也发动机必须全部常温，不然会有热量作为热能耗散。

这当中做到任何一点都是科技的极大进步，而但凡物质不能全部湮灭（而是采取核聚变），能量难以完全收集与转化，则能量的利用率只会千万倍地往下降，消耗的燃料也必将回归上一小节当中的恐怖水平。

结语

经过上文的分析，我们发现其实流浪地球的世界观本身还是有一点不科学的，尤其是地球发动机，简直是Bug一般的存在。但是这并不是说流浪地球就不是一部好片子，因为几乎每一部科幻片都有能让人杠一杠的地方。因此写这篇文章只是为了分享一下自己想法，也帮助科普一下一些轨道力学的知识。当然在推式子的过程中也和同学们复习了一下指对数函数的知识，真是极好的。

安利

什么？你问我自己怎么会想到这些？

因为本人是Kerbal Space Program玩家啊#(滑稽)

在这里给大家安利一下这款“游戏”，游戏本身非常硬核，需要玩家自己设计并且发射航天器进入轨道，自己规划轨道计算燃料，自己操控交会对接，模拟了近乎真实的物理系统，内有一个类似太阳系的星系。上述的计算，霍曼转移，燃料量的计算都是游戏当中的必修课是的现在玩游戏都要用计算器。个人认为比一些纯粹娱乐向的手游之类的更能从中学到知识，但同时因为这款游戏玩起来消耗大量脑细胞与时间，一般每周末玩两个小时就够了。

（我自己搭的Jool飞掠器的截图，Jool对应太阳系的木星）

从动手设计到发射成功现实当中花了一个半小时，游戏当中是两年半，失败五次。本质新手劝退向游戏

当然，虽然这款游戏Steam上要七十多块，但是破解版烂大街。欢迎有兴趣的同学来玩。

数据来源

大部分公式，数据均来源于英文维基百科，少部分来自NASA官网。

公式全部基于万有引力定律与能量守恒定律，绝大部分的公式可以在高中知识下进行推导。

简介

后缀数组是字符串当中常见的一种数据结构，它将一个字符串的所有后缀进行排序，得到以下三个数组：

1. ：排序后排名为的后缀的起始位置
2. ：起始位置为的后缀的排名
3. ：起始位置的后缀和起始位置的后缀的最长公共前缀长度

利用这三个数组我们可以有力地解决字符串当中的大部分问题。

构造

显然如果我们朴素地去构造长度为的字符串的后缀树组的话，时间复杂度为。这对于长字符串来说是完全无法接受的。然而利用倍增的思想我们可以在之内快速构造后缀树组，虽然像DC3，SAIS之类的算法也存在，但是倍增的编码难度更低，也足以通过大部分的题目（毕竟可能很多时候瓶颈不在于构造）。

的构造代码如下：

void build_sa() {
    static int a[N], b[N], d[N], tmp[N], buc[N];
    copy(s + 1, s + 1 + n, d + 1);
    sort(d + 1, d + 1 + n);
    int *end = unique(d + 1, d + 1 + n);
    for (int i = 1; i <= n; i++) a[i] = lower_bound(d + 1, end, s[i]) - d;
    for (int i = 1; i <= n; i++) buc[a[i]]++;
    for (int i = 1; i <= n; i++) buc[i] += buc[i - 1];
    for (int i = 1; i <= n; i++) rk[i] = buc[a[i] - 1] + 1;
    for (int t = 1; t <= n; t *= 2) {
        copy(rk + 1, rk + 1 + n, a + 1);
        for (int i = 1; i <= n; i++) b[i] = i + t > n ? 0 : rk[i + t];
        
        fill(buc, buc + n + 1, 0);
        for (int i = 1; i <= n; i++) buc[b[i]]++;
        for (int i = 1; i <= n; i++) buc[i] += buc[i - 1];
        for (int i = n; i; i--) tmp[buc[b[i]]--] = i;

        fill(buc, buc + n + 1, 0);
        for (int i = 1; i <= n; i++) buc[a[i]]++;
        for (int i = 1; i <= n; i++) buc[i] += buc[i - 1];
        for (int i = n; i; i--) sa[buc[a[tmp[i]]]--] = tmp[i];

        bool flag = true;
        for (int i = 1, last = 0; i <= n; i++) {
            int p = sa[i];
            if (!last) rk[p] = 1;
            else if (a[p] == a[last] && b[p] == b[last])
                flag = false, rk[p] = rk[last];
            else rk[p] = rk[last] + 1;
            last = p;
        }
        if (flag) break;
    } 
}

其思想如下：如果我们先把所有后缀按照它们前个字符进行排序，那基于排完序的结果我们就可以使用双关键字排序对所有后缀按照前个字符进行排序（原因是后缀的后缀还是后缀）,当然在实践过程中我们常使用桶排来加速这一个过程。

代码的最后一段是在互异的时候直接跳出，这是一个小优化。

计算

除了朴素比较之外，数组的计算可以在时间内完成。

观察：设，则我们有：

证明：

1. ：Trivial
2. ：设表示数组当中排在以为起始的后缀前的后缀的起始位置，以开始的后缀与以开始的后缀的LCP为。考虑统一截掉开头的一个字符，得到起始位置为与的两个后缀，其LCP为，在里面肯定排在前面，又因为后缀数组的性质，与的LCP一定不会小于和的LCP，因此不等式成立。

基于以上的观察，计算的代码如下：

for (int i = 1, k = 0; i <= n; i++) {
    if (rk[i] == 1) k = 0;
    else {
        if (k > 0) k--;
        int j = sa[rk[i] - 1];
        while (i + k <= n && j + k <= n && s[i + k] == s[j + k]) k++;
    }
    ht[rk[i]] = k;
}

NTT的思想与FFT类似，但是避免了浮点数运算带来的精度问题。

原根

定义：设且，使得的最小称为模的阶，记作。

定理：显然，对于所有，都有。

同时根据欧拉定理我们有

因此有。

定义：若则称为模意义下的原根。

存在原根当且仅当，其中为奇素数，。

原根的意义：如果存在原根，那恰能表示所有个与互素的数，且构成的简化剩余系。

原根的计算方法：原根只有个，因此只要枚举即可。

NTT

NTT中使用单位原根代替单位复数根进行计算，单位原根定义为：

其中为模素数的原根。

显然有与单位负数根类似的性质：

又因为以及一些二次剩余的知识，只可能是，又因为原根的性质 同时，观察的定义，易证（折半引理），原根具备所有FFT中利用的单位复数根的性质，因此也可以使用单位原根进行FFT，得到的算法就是NTT（Number Theoretic Transform），NTT适用于整数序列的变换，且精度与数值稳定性较FFT更优。

代码：

#include <cstdio>
#include <algorithm>

using namespace std;
const int MOD = 998244353;
const int N = 270000;
const int G = 3, INVG = 332748118;
typedef long long ll;
int n, m; 
ll a[N], b[N];

ll qpow(ll x, ll y) {
    ll ret = 1;
    for (; y; y >>= 1) {
        if (y & 1) ret = ret * x % MOD;
        x = x * x % MOD;
    }
    return ret;
}

void ntt(int len, ll *x, int dir = 1) {
    static ll X[N];
    for (int i = 0, k = 0; i < len; i++) {
        X[i] = x[k];
        int y = len >> 1;
        while (y & k) k ^= y, y >>= 1;
        k |= y;
    }
    for (int s = 2; s <= len; s <<= 1) {
        ll w_ = qpow(G, (MOD - 1) / s);
        if (dir == -1) w_ = qpow(w_, MOD - 2);
        for (int l = 0; l < len; l += s) {
            int mid = l + (s >> 1);
            ll w = 1;
            for (int i = 0; i < s >> 1; i++) {
                ll t = w * X[mid + i] % MOD;
                X[mid + i] = (X[l + i] - t) % MOD;
                X[l + i] = (X[l + i] + t) % MOD;
                w = w * w_ % MOD;
            }
        }
    }
    copy(X, X + len, x);
}

int main() {
    scanf("%d%d", &n, &m);
    for (int i = 0; i <= n; i++) scanf("%lld", &a[i]);
    for (int i = 0; i <= m; i++) scanf("%lld", &b[i]);
    int sz = 1; while (sz < n + m + 1) sz <<= 1;
    ntt(sz, a); ntt(sz, b);
    for (int i = 0; i < sz; i++) printf("%lld %lld\n", a[i], b[i]);
    for (int i = 0; i < sz; i++) a[i] = a[i] * b[i] % MOD;
    ntt(sz, a, -1); ll inv = qpow(sz, MOD - 2);
    for (int i = 0; i < n + m + 1; i++)
        printf("%lld ", (a[i] * inv % MOD + MOD) % MOD);
    puts("");
    return 0;
} 

常用的NTT素数

NTT所用的素数应该满足什么条件呢？观察的定义显然我们期望对于一定范围内的（比如说）都最好有。即的形式。

常见的NTT素数为：

1. 
2. 
3. 
4. 
5. 

这些素数的原根都包含。

任意模数NTT

有些时候题目要求任意模数的NTT变换，这个时候我们可以选取三个上述列表中模数分别进行NTT，然后通过CRT进行合并。具体上来说，对于一个数，若有：

我们考虑两两合并，先求解，不难得到（逆元在模意义下）。

令，再求解，此时有（逆元在模意义下）。

在计算的最后再对于指定的模数取模即可。

代码：

#include <cstdio>
#include <algorithm>

using namespace std;
typedef long long ll;
const ll G = 3;
const ll MOD1 = 998244353, MOD2 = 1004535809, MOD3 = 469762049;
const int N = 270000;
int n, m, P;

int qpow(int x, int y, int mod) {
    int ret = 1;
    for (; y; y >>= 1) {
        if (y & 1) ret = (ll)ret * x % mod;
        x = (ll)x * x % mod;
    }
    return ret;
}

struct mint {
    int m1, m2, m3;
    mint() {}
    mint(int x) : m1(x % MOD1), m2(x % MOD2), m3(x % MOD3) {}
    mint(int m1, int m2, int m3) : m1(m1), m2(m2), m3(m3) {}
    mint reduce() { 
        return mint(m1 + (m1 >> 31 & MOD1), // m1 >> 31 是符号位
                    m2 + (m2 >> 31 & MOD2), 
                    m3 + (m3 >> 31 & MOD3));
    }
    mint operator +(const mint &x) const {
        return mint(m1 + x.m1 - MOD1, 
                    m2 + x.m2 - MOD2, 
                    m3 + x.m3 - MOD3).reduce();
    }
    mint operator -(const mint &x) const {
        return mint(m1 - x.m1, m2 - x.m2, m3 - x.m3).reduce();
    }
    mint operator *(const mint &x) const {
        return mint((ll)m1 * x.m1 % MOD1, 
                    (ll)m2 * x.m2 % MOD2, 
                    (ll)m3 * x.m3 % MOD3);
    }
    mint inv() {
        return mint(qpow(m1, MOD1 - 2, MOD1), 
                    qpow(m2, MOD2 - 2, MOD2), 
                    qpow(m3, MOD3 - 2, MOD3));
    }
    int get() {
        const int INV_1 = qpow(MOD1, MOD2 - 2, MOD2);
        const int INV_2 = qpow((ll)MOD1 * MOD2 % MOD3, MOD3 - 2, MOD3);
        ll x = (ll)(m2 - m1 + MOD2) % MOD2 * INV_1 % MOD2 * MOD1 + m1;
        return ((ll)(m3 - x % MOD3 + MOD3) % MOD3 * INV_2 % MOD3 
                * ((ll)MOD1 * MOD2 % P) % P + x) % P;
    }
    static mint unit(int s) {
        return mint(qpow(G, (MOD1 - 1) / s, MOD1), 
                    qpow(G, (MOD2 - 1) / s, MOD2), 
                    qpow(G, (MOD3 - 1) / s, MOD3));
    }
} a[N], b[N];

void ntt(int len, mint *x, int dir = 1) {
    static mint X[N];
    for (int i = 0, k = 0; i < len; i++) {
        X[i] = x[k];
        int y = len >> 1;
        while (y & k) k ^= y, y >>= 1;
        k |= y;
    }

    for (int s = 2; s <= len; s <<= 1) {
        mint w_ = mint::unit(s);
        if (dir == -1) w_ = w_.inv();
        for (int l = 0; l < len; l += s) {
            int mid = l + (s >> 1);
            mint w = 1;
            for (int i = 0; i < s >> 1; i++) {
                mint t = w * X[mid + i];
                X[mid + i] = X[l + i] - t;
                X[l + i] = X[l + i] + t;
                w = w * w_;
            }
        }
    }
    if (dir == -1) {
        mint inv = mint(len, len, len).inv();
        for (int i = 0; i < len; i++) X[i] = X[i] * inv;
    }
    copy(X, X + len, x);
}

int main() {
    scanf("%d%d%d", &n, &m, &P);
    for (int i = 0; i <= n; i++) { int tmp; scanf("%d", &tmp); a[i] = tmp; }
    for (int i = 0; i <= m; i++) { int tmp; scanf("%d", &tmp); b[i] = tmp; }
    int sz = 1; while (sz < n + m + 1) sz <<= 1;
    ntt(sz, a); ntt(sz, b);
    for (int i = 0; i < sz; i++) a[i] = a[i] * b[i];
    ntt(sz, a, -1);
    for (int i = 0; i < n + m + 1; i++) printf("%d ", a[i].get());
    puts("");
    return 0;
}

问题

求满足：

的的值。

ExCRT

当互素的时候显然可以使用中国剩余定理（CRT）解决，时间复杂度为。

考虑不互素的情况，假设已经解决了前个方程，得到的解为，。

那么我们考虑解决

易证得到的解不会再和前个解产生冲突。

稍微移项一下：

用ExGCD解出即可。

这套方法称为ExCRT（拓展中国剩余定理），复杂度依旧是。全面优于CRT。

代码：

#include <cstdio>

using namespace std;
typedef long long ll;
const int N = 100010;
int n;
ll a[N], b[N];

ll exgcd(ll a, ll b, ll &x, ll &y) {
    if (b == 0) return x = 1, y = 0, a;
    ll d = exgcd(b, a % b, y, x);
    return y -= a / b * x, d;
}

ll smul(ll x, ll y, ll mod) {
    if (y < 0) return -qmul(x, -y, mod);
    ll ret = 0;
    for (; y; y >>= 1) {
        if (y & 1) ret = (ret + x) % mod;
        x = (x + x) % mod;
    }
    return ret;
}

ll excrt() {
    ll lcm = a[1], t = b[1], x, y;
    for (int i = 2; i <= n; i++) {
        ll d = exgcd(lcm, a[i], x, y);
        x = smul(x, (b[i] - t) / d, a[i]);
        t += x * lcm;
        lcm = lcm / d * a[i];
        t = (t % lcm + lcm) % lcm;
    }
    return t;
}

int main() {
    scanf("%d", &n);
    for (int i = 1; i <= n; i++) scanf("%lld%lld", &a[i], &b[i]);
    printf("%lld\n", excrt());
    return 0;
}

上次更新是在十月份？我到底摸了多久

卷积

对于两个序列，定义它们的卷积为：

单位根

定义复平面单位圆上的次单位根

易证其有如下的性质：

1. 
2. 
3. 

DFT

对于一长度为的序列，定义其离散傅里叶变换（DFT）为：

则有其离散傅里叶逆变换（IDFT）为： （标准定义里面和应该互换，但是不影响）

暴力证明：考虑的计算式：

设左边矩阵为，接下来从充分性上证明： 显然当时有，而时，由等比数列求和公式： 因此，得证。

上面的DFT有如下意义：

1. 将作为一组正交基的系数，计算和向量。（在此意义下逆变换就是将向量进行分解，因此逆变换的意义显然，无需上述证明）
2. 将作为以个多项式的系数，计算该多项式在点上的取值。

DFT是线性的，即和的DFT等于DFT的和，且常数倍的DFT是DFT的常数倍。

同时DFT的另一个重要性质被称为卷积性，即

证明：把式子写出来： 显然左边式子展开对应的任意一项都唯一地出现在右式内层求和的展开式中，因此两个式子总是等价的。

DFT将卷积的操作转化为的操作，因此可以用于加速卷积。

FFT

因为DFT对应点相乘这一操作本身是的，因此我们只要找到快速的DFT/IDFT算法，就可以避免地计算卷积。FFT就是这样一种基于分治的快速DFT算法，时间复杂度为。

首先假设我们需要计算一个序列的DFT（假设是的整数次幂），我们按照下标的奇偶性将这个序列分为两部分，其中且。那我们有：

看起来好像没有问题？我们最后一步推导其实是有一点小问题的，因为当依据我们的定义只有项，因此上面的推导只能让我们计算的值。对于的情况要如何处理呢？

首先，注意到：

因此我们稍微变一下： （其实本质上就是DFT的序列其实是可以周期拓展的）

这使得我们得以计算时DFT的值。因此，最为朴素的FFT实现就是基于递归不断将进行分治，并通过上式对于分治结果进行合并。时间复杂度自然就是。

然而事实一定没有那么简单。递归写法的FFT常数很大，很大程度上是因为合并本身不是in-place的，需要额外申请空间。更快的in-place FFT实现，Cooley–Tukey算法，是基于这样一个思路的：

假设每次分出来的向左走，向右走，那递归版的FFT就会画出一棵很漂亮的二叉树，其中叶子节点对应的下标和序列本身的下标顺序是有很大差别的。我们可以预先快速计算出叶子节点的顺序，然后直接从叶子开始就地一层一层地合并上去。

那叶子节点的顺序遵循怎样的规律呢？其实我们依据奇偶分组本质上是在依据二进制的末位分组，在这样的思想下，我们脑补一下发现这就是叶子节点的下标把原来序列的下标的二进制翻转一下（低位变高位，高位变低位）（即到叶子节点变成）（可以结合画图理解）。这样的话我们维护一个从高位到低位的二进制加法器就行了。

代码如下：

#include <cstdio>
#include <complex>
#include <cmath>

using namespace std;
typedef complex<double> cplx;
const double PI = 3.1415926535897932384626;

void fft(int len, cplx *x, cplx *X, int dir = 1) {
    for (int i = 0, k = 0; i < len; i++) {
        X[i] = x[k];
        int y = len >> 1;
        while (y & k) k ^= y, y >>= 1; // 维护倒序加法器
        k |= y; 
    }
    for (int s = 2; s <= len; s <<= 1) {
        cplx w_(cos(2 * PI / s), dir * sin(2 * PI / s));
        for (int l = 0; l < len; l += s) {
            int mid = l + (s >> 1);
            cplx w = 1;
            for (int i = 0; i < s >> 1; i++) { // 逐层向上合并（一次蝴蝶操作）
                cplx t = w * X[mid + i];
                X[mid + i] = X[l + i] - t;
                X[l + i] = X[l + i] + t;
                w *= w_; 
            }
        }
    }
}

电路图

绿色箭头为PWM开段电流方向，红色为闭段电流方向。

模型

符号定义

符号	含义
	PWM频率
	PWM周期
	PWM高电位时间
	PWM低电位时间
	电流（随时间变化）
	PWM周期开始时的电流
	PWM峰值电流
	PWM周期结束时的电流
	电池电压
	反电动势
	电感电压
	二极管压降
	电感
	马达内阻
	系统内阻

部分参数可以自行测定，测定方法见文末。

开段分析

首先，由基尔霍夫定律和欧姆定律，可得：

对上式进行整理，并结合电感电压公式： 我们把它看做一个关于的线性微分方程，初始条件为解得： 我们知道电流在高电压段逐步变大，因此，即： 出于简化式子的需要，令 注意到为时的电流，因此也称作稳态电流。

整理得到：

为了之后平均电流的计算需要我们顺便计算一波的积分：

闭段分析

再次由基尔霍夫定律及欧姆定律，有：

同样整理： 以初始条件求解这个微分方程解得： 类似地： 我们也令： 可见也是时的电流，也是稳态电流。

整理得到：

为了后文计算平均电流的需要，我们也计算积分：

注意到以上式子都依赖于和，前者就是占空比乘以PWM总周期长度，而后者则有些棘手，因为并不总是成立。

接下来的分析有两种情况：

1. 连续电流：PWM周期中的电流是连续的，即。
2. 不连续电流：电流在PWM闭段下降至被二极管截断，即。

两种情况需要分类讨论。

连续电流

由，将展开：

解出：

不连续电流

这个时候我们知道，由的计算式我们可以反解出：

从实现上来说，我们先假定（即假设没有二极管），计算出，检查它的正负性，来判断电流是否被截断。

平均电流计算

最后，我们结合计算出来的$ t_, t_$，以及开闭段计算出来的两个积分，我们终于可以计算PWM过程中的平均电流：

这也是我们模型计算的最终结果。

在实际的实现当中需要考虑更多的细节，比如说转速的方向，命令的方向等等……并不总是如以上的公式那么简单，是一项艰苦的体力劳动，具体可以参考In The Zone下的model.c。

逆向计算

很多时候我们不仅需要计算马达一定情况下的电流，更需要计算马达在某一状态下为了达到某一指定电流（扭矩）所需要的指令大小。例如：

1. 计算空载情况下马达到达某一速度所需要的指令值（精确的马达线性化）
2. 把PID的输出看做目标扭矩使用模型再反解出指令。

针对这个需求，因为马达曲线直观上总是单调的，因此理论上你总可以根据数学反推推出一个公式，但是由于模型当中牵涉到了分类讨论（电流的连续性），因此较为复杂。利用马达曲线的单调性，二分不失为一个更好的做法。由于指令的范围不大，我们至多只需要二分七八次即可。

注意：模型的计算在Cortex上还是较为耗时间的，实际当中根据当前电压计算一次完整的马达曲线需要4秒左右，优化常数，提高模型计算的效率一直是一个大问题，这也是为什么在实际应用中也许tanh线性化等计算更快的线性化的方法会更好的原因。

结果

这是真实的速度-指令曲线（空载）：

这是使用模型计算的速度-指令曲线：

两者几乎一致，说明建模正确。

参数测定

马达内阻

使用伏安法，我们首先需要测量出马达在127命令下堵转电流，127命令是为了尽可能消除电感的影响（即PWM全开），堵转是为了避免电能转化为机械能消耗掉，则根据开段分析当中的

此时，可反解出

反电动势常数

我们知道，其中是转速， 为反电动势常数，为了测定其值，我们测量出马达在127命令无负载运转电流和此时的转速（单位不论，但是需要保持统一），由

此时，解得 那么

简介

成套方法（repertoire method）可以用来求解一大类的递归式的封闭形式，由于和式也可以化为特殊的递归式，因而该方法对于处理和式也有很好的效果，具体而言，对于如下类似形式的递归式：

我们总是可以大胆地猜想：或许有着如下的封闭形式： 成套方法的精髓就是，试图用几组特殊的 （可以让我们相对容易地解出等的值），来推测等之间的关系，由于在上面的式子当中等都是线性组合的，因而最后我们会得到一组线性方程，求解这一组线性方程我们可以得出从而解出递归式的特殊形式。这是一个特殊→一般→特殊→一般→特殊的方法，十分巧妙。

例子：求解

这是一个和式，但是不打紧，我们可以把它处理成递归式。

令，那么有等价定义

可以看到从和式到递归式的转化是十分自然的。

这个递归式自然是满足最前面的形式的，其中，我们先着手求解一般情况，选取几组“特殊的”：

首先令，代入最开始的定义：

显然有，再把他们带入我们认为的“通用形式”： 很好，我们现在得到一组的线性关系了，还需要三组！

令

令： 最后，令： 最后，把这四组关系结合起来： 非常简单的线性方程组，解得： 这个时候，回到我们要解的递归式，其中,因此 求解到此结束，可以看到这个公式和我们的常识一致。虽然以上过程看似繁琐，但是注意我们已经求解了一大类的问题，因为是不会因而变化的！，现在，对于任何的二次多项式的和式或者递归式我们都可以飞快地计算出其封闭形式。更有用的是，不仅仅是二次三项式，如果我们以后要处理三次多项式（意味着我们要加入和），之前的的解也不会变化，因而只要继续带入即可！ （这是显然的）

注意：为了找出线性关系，我们也可以先找出等的一组特解，然后基于这个解计算，这和代入特殊然后反解是等价的！

思考：为什么可以这样做 ？

为什么几个看似毫不相关的可以指导我们正在处理的特殊情形，注意到我们的步骤：

1. 将特殊带入递归式的定义，解出
2. 将1的解带入假设的封闭形式，求出的关系

注意到，我们在第一步解出了等的值，说明此时我们的特殊本来就是以上的递归式定义可以表出的，而我们假设的通用形式是对任何一组参数通用的，因此求解出来的线性关系也是通用的。

思考：真的“万能”和“机械化”吗？

成套方法从刚刚的示例可以看到似乎是非常机械化的，但是其实不然，《具体数学》中有言，成套方法对于线性的递归式最为有效，因此前是不能有系数的（不然直觉上通用形式里就会出现棘手的指数项），但是这不妨碍我们处理和式，其次，对于一开始的“通用形式”的构造和特殊的构造其实不是那么简单的工作，对于多项式，我们可以很简单地构造多项式形式的“通用形式”，但是如果是：

最上面的通用形式就不管用了！我们需要自己寻找适合的通用定义和特殊函数值进行求解（这也是《具体数学》当中的习题），虽然找到之后的过程倒是非常地机械化没错。

思考：作为算法实现的复杂度？

我们仅在这里思考以上方法在多项式情况下的特例。

成套方法分为三步走：

1. 构造线性方程
2. 解线性方程
3. 回代

第一步，由于是多项式的情况，观察之前的特例并结合二项式定理可以得到， 如果我们带入的情况，对于对应的多项式系数（类似），其系数为（在时，在时系数自然就是0了）：

因此，配合杨辉三角递推组合数，这一部分的时间复杂度可以做到。

第二步，考虑到我们解的线性方程当中含有多项式，因此需要采用系数向量的表示方法，基本算数运算代价变为，但是注意到线性方程组列出时已经是简化阶梯形式，因此不需要高斯消元，只需要代入消元，时间复杂度是。

第三步的复杂度自然是。

因此，基于成套方法的算法的综合复杂度为。

题意

给定以及，要求最大化：

其中。

解法

我们可以二分答案，这样这个问题就转变为：能否求出一组解，使得？

对不等式进行变形：

我们计算的最大值，即构造使当且仅当，并判断最大值的正负性，在此基础上不断二分答案即可。时间复杂度为。

其实这个并没有什么卵用，基本没有人会卡。

对于求单个数的逆元我们有的算法，因此，如果我们要求的所有数的逆元的话，朴素的做法是可以做到的，但是其实我们有线性递推的算法，假设我们要求在模意义下的逆元，且的逆元已经求出：

设，其中，我们有：

两边同时乘上： 由于我们之前已经求得了的逆元，的逆元依上式可以直接递推，代码也很短：

inv[1] = 1;
for (int i = 2; i <= n; i++)
    inv[i] = (-(P / i) * inv[P % i] % P + P) % P;

但是真的没有卵用啊，动不动就是或者这种，就算你能线性递推时间和空间也要炸糊的。

简介

目标：对于一个位二进制表示的子集x，我们希望枚举其子集。

朴素枚举其子集的代码，复杂度：

for (int s = 0; s <= x; s++)
    if (s | x == x)
        // s 表示的二进制是x的子集

快速枚举其子集的代码，复杂度，其中是x表示的集合：

for (int s = x; s > 0; s = (s - 1) & x)
    // s表示的二进制是x的子集
// 以上算法不！会！处！理！空！集！，额外处理！

正确性证明：

代码中的& x保证s一定是s的子集，而(s - 1)保证枚举出来的s递减至0，详细证明不在此赘述脑补即可.

拓展

在很多的状压DP当中，我们不仅要枚举一个集合的子集，而是要枚举的所有二进制数表示的集合的子集，如果使用朴素地暴力做法，那么复杂度显然是，可以过的数据。

那么如果我们使用刚刚介绍的优化方法，复杂度会变成多少呢？

分析：总体复杂度为：

逆用二项式定理： 因此，算法的复杂度为（而且是最优的），可以通过的数据。