数学

在正式开始学习数学之前，让我们先来玩一个数学游戏~

幻方

经典的数学游戏——幻方

幻方是一种神奇的N x N矩阵：它由数字1，2，3，……，N x N构成，且每行，每列以及两条对角线的数字之和都相同。

相信我们从小学的时候就开始接触幻方了，记得从那个时候起就有这种填幻方的数学题。不过以前都是一通瞎填，毫无章法。而最近因为各种机缘巧合再一次碰到了幻方（~~死去的记忆开始攻击我~~）那今天就让我们正式来讨论一下幻方的解法。

（笔者能力有限，在这里我们暂时只考虑奇数阶的幻方填法）

起源

幻方最早起源于中国。又称方阵或者厅平方，宋代数学家杨辉称之为纵横图。相传，大约两千年前，夏禹治水之时，黄河中跃出一匹神马，马背上驮着一幅图，人称《河图》；又洛水河中浮出一只神龟，龟背上也有一张象征吉祥的图案，被人们称为《洛书》。他们发现，图案每一列，每一行及对角线，加起来的数字和都是一样的，这就是我们所称的幻方。

幻方的解法

幻方的解法有很多种，一个幻方的可行解也不止一个。本文我们主要探讨几种我们常见的解法。

罗伯法

图解

罗伯法是最常见的一种奇数阶幻方解法。又称楼梯法。其原理的数学证明相对来讲比较复杂，这里我们因篇幅考虑就不对原理做详细阐述，只讲方法。

罗伯特法可以总结出如下口诀：

1填顶行正中央，依序斜填逐步上，上填出框转底行，右填出框左边放，重复便在下格填，出角重复一个样。

罗伯法图解如下，以3阶幻方为例：

1填顶行正中央——数字1填在首行的最中间的格子。

依序斜填逐步上——向右上角斜行，依次填入数字。

上填出框转底行——如果出了上框边界，就将数字填在该列的最下方。

右填出框左边放——如果出了右框边界，就将数字填在最左行。

重复便在下格填——如果右上格子已有数字，就在该格子的下方填数。

出角重复一个样——如果朝右上角出框界，和“重复”的情况同步。

最后8和9按照右出界、上出界的方法填写。

代码

罗伯法的描述如下：

首先将1写在第一行中间。

之后，按如下方式从小到大依次填写每个数K（K = 2，3，……，N x N）：

若（K - 1）在第一行但不在最后一列，则将K填在最后一行，（K - 1）所在列的右一列；
若（K - 1）在最后一列但不在第一行，则将K填在第一列，（K - 1）所在行的上一行；
若（K - 1）在第一行最后一列，则将K填在（K - 1）的正下方；
若（K - 1）既不在第一行，也不在最后一列，如果（K - 1）的右上方还未填数，则将K填在（K - 1）的右上方，否则将K填在（K - 1）的正下方。

直接上代码（洛谷P2615）

#include<iostream>
using namespace std;
int a[1005][1005];
int n, x, y;
int main() {
    cin >> n;
    a[1][n / 2 + 1] = 1;  //1写在第一行中间
    x = 1, y = n / 2 + 1;    //记录位置
    for (int i = 2; i <= n * n; ++i) {
        if (x == 1 && y != n) {
            a[n][y + 1] = i;
            x = n, y++;
        }
        else if (x != 1 && y == n) {
            a[x - 1][1] = i;
            x--, y = 1;
        }
        else if (x == 1 && y == n) {
            a[x + 1][y] = i;
            x++;
        }
        else if (x != 1 && y != n && a[x - 1][y + 1] == 0) {
            a[x - 1][y + 1] = i;
            x--, y++;
        }
        else if (x != 1 && y != n && a[x - 1][y + 1] != 0) {
            a[x + 1][y] = i;
            x++;
        }
    }
    for (int i = 1; i <= n; ++i) {
        for (int j = 1; j <= n; ++j) {
            cout << a[i][j] << ' ';
        }
        cout << endl;
    }
    return 0;
}

杨辉法

我国数学家杨辉也早就给出了三阶幻方的解法，其方法为：九子斜排，上下对易，左右相更，四维挺进。（~~文言文终究看起来还是高大上一点，老祖宗牛逼！~~）

如下图：9个数斜着排，上下的两个数1,9，左右的两个数3,7互相换一下，四个角上的2,4,6,8就移到那四个角上去，这样就填好了一个三阶幻方了。

巴舍法

巴舍法只不过是杨辉法的变形，在你完成杨辉法的“九子斜排”后（如上图中左侧图所示），将1向下3格填入，9向上3格填入，7向右3格填入，3向左3格填入即可。即在n x n方格外的数字向其对称方位移动n格。

幻方通项公式的推导

终于来到最后的通项公式推导啦，让我们以错位补角法为基础来推导幻方的通项公式吧！

之前有提过，罗伯法也称“楼梯法”，而所谓的“错位补角法”本质上也是罗伯法，而通项公式的推导就是从这个方法入手。

错位补角

先来看“错位补角”（来自百度百科）：

1、对于所有的奇数阶幻方，1 - n x n从小到大填入方格中，构造一张有序数表。以5阶为例：

2、横错位，将方格横向错位，每行错位数为 n-行数，即第一行横向移动n-1位，第二行横向移动n-2位…直到形成一个左低右高的楼梯。

3、横补角，以中间行为基准，将突出的数字补回本行所缺的方格内，4，5补到1的前，10补到6前，16补到20后，21，22补到25后。从而重新得到一个n*n方格。

4、竖错位，将方格纵向错位，每列错位数为 n-列数，即第一列横向移动n-1位，第二列横向移动n-2位…直到形成一个左低右高的楼梯。

5、竖补角，以中间列为基准，将突出的数字补回本列所缺的方格内，17，23补到4上，24补到5上，2补到21下，3，9补到22下。从而重新得到一个无序数表，及得到结果。

公式推导

我们的目的是利用错位补角法来推导出通项公式，那么我们首先来捋一下思路：

错位补角法先是让我们构造了一张有序数表（见错位补角法的第一步），我们规定行数从上往下数是0到n-1，列数从左往右数也是0到n-1。因为数表是有序的缘故，我们可以很容易地得到数表中元素和其所在行列数的关系，即第 i 行第 j 列的元素值是 i * n + j + 1 （1式）。例如有序数表中13的位置是第2行第2列（注意这里的行和列全都是从0开始数），那么 i = 2， j = 2，n = 5，而2 * 5 + 2 + 1 = 13，可以看到，其结果确实是和当前位置的元素值对应。

而后这个有序的数表经过错位补角变成了一张无序的数表（见错位补角法的最后一步）。假设我们现在要利用输出语句直接将最后的无序数表给输出出来，我们需要先知道这两个数表相同元素的位置关系。换句话来讲，当你需要输出幻方中第 i 行第 j 列的元素，你就需要知道这个元素在原来有序数表中是在第几行第几列，之后我们就可以利用（1式）直接计算出相对应的元素值。我们假设幻方中第 i 行第 j 列的元素在原来的有序数表中式第 i’ 行和第 j’ 列，那么我们需要找到这二者位置上的映射关系。

那么，如何找出这个映射关系呢？

我们利用（i，j）来表示在数表当中的第 i 行第 j 列。

我们知道，原有序数表当中（i’ ，j’）经过横错位、横补角、竖错位、竖补角的运算之后得到无序数表中对应的位置（i，j）。而我们想在知道（i，j）的情况下计算（i’，j’），显而易见，我们只需要对错位补角法进行逆运算就可以还原出（i，j）了。

我们先来进行第一步逆运算。（~~让我们开始进行华丽的逆转吧！~~）

观察错位补角法第五步和第三步的数表，可以发现，第五步到第三步，元素的列位置不发生变化，而行位置的变化从左到右依次是：向上2格，向上1格，不变，向下1格，向下2格。记向下为正，向下为负，则行变换的数值可以记为-2，-1，0，+1，+2。而对应的列数是0，1，2，3，4，故有行变换的数值 = 对应列数 - （n - 1）/2 ，第三步行数 = 第五步行数 + 行变换的数值。而考虑到计算出来的行数有可能是负数，但行数是大于0且小于等于n-1的，所以我们需要将其修正为对应位置的正数表示，这也很简单，只需要利用取模运算就好了，即第三步行数 = （第五步行数 + 对应列数 - （n - 1）/2 + n）% n，整理成第三步行数 = （第五步行数 + 对应列数 +（n + 1）/2 ）% n，所以，可以得到第一步逆运算为：（i，j）-> （（i + j + （n + 1）/2）% n，j ）。

接下来进行第二步逆运算。（~~逆转！故技重施！~~）

我们记（i + j + （n + 1）/2）% n = tmpi （2式），则上述的逆运算表示为（i，j）-> （tmpi，j ）。

观察错位补角法第三步和第一步的数表，可以发现，第三步到第一步，元素的行位置不发生变化，而列位置的变化从上到下依次是：向左2格，向左1格，不变，向右1格，向右2格。记向右为正，向左为负，则列变换的数值可以记为-2，-1，0，+1，+2。而对应的行数是0，1，2，3，4，故有列变换的数值 = 对应行数 - （n - 1）/2，第一步列数 = 第三步列数 + 列变换的数值。而考虑到计算出来的列数有可能是负数，但列数是大于0且小于等于n-1的，所以我们需要将其修正为对应位置的正数表示，故技重施，利用取模运算即可。即第一步的列数 = （第三步的列数 + 对应行数 - （n - 1）/2 + n）% n，整理成第一步列数 = （第三步列数 + 对应行数 +（n + 1）/2 ）% n，所以可以得到第二步的逆运算为：（tmpi，j）->（tmpi，（tmpi + j +（n + 1）/2）% n）= （i’，j’）。至此，我们要的 i’ 和 j’ 就计算出来了。

整理一下，结合（2式）可得：经过映射后的（i’，j’）= （ i + j + （n + 1）/2）% n，（（i + j + （n + 1）/2）% n + j +（n + 1）/2 ）% n ）。

但是这条公式实在是泰长辣！而我们这一章是关于数学的，数学嘛，讲究的就是精简，这条公式的后半段有着很明显的化简空间！

先来看公式，对于取模运算，有如下性质：（a + b + …… + z）% n = （a % n + b % n + …… + z % n）% n（3式）

而该公式的后半段（（i + j + （n + 1）/2）% n + j +（n + 1）/2 ）% n代入（3式）可以拆成（j % n + （（i + j + （n + 1）/2）% n）% n + （（n + 1）/2）% n）% n其中，（（i + j + （n + 1）/2）% n）% n和（i + j + （n + 1）/2）% n 没有区别，故我们直接利用后者代替前者，得（j % n + （i + j + （n + 1）/2）% n + （（n + 1）/2）% n）% n，逆用（3式）得原式 = （j + i + j + （n + 1）/2 + （n + 1）/2）% n，整理有：原式 = （2*j + i + n + 1）% n，最终，大功告成！

（i’，j’）= （i + j + （n + 1）/2）% n，（2*j + i + n + 1）% n ）。

最后，代入（1式），可得幻方中第 i 行第 j 列的元素值为 i’ * n + j’ + 1。

总结

记$a_{ij}$为幻方第 i 行第 j 列的元素，在罗伯法的解下，则有 $a_{ij}$ = （（i + j + （n + 1）/2）% n）* n + （ 2*j + i + n + 1）% n + 1

推导完毕，华丽收场。

高精度

高精度加法

int a[100], b[100], c[100];
string s1, s2;
cin >> s1 >> s2;    //输入两个数字
int len = max(s1.size(), s2.size());
for (int i = 0; i < s1.size(); ++i) {
    a[s1.size() - 1 - i] = s1[i] - '0';    //逆序输入
}
for (int i = 0; i < s2.size(); ++i) {
    b[s2.size() - 1 - i] = s2[i] - '0';    //逆序输入
}
//处理相加
for (int i = 0; i < len; ++i) {
    c[i] += a[i] + b[i];
    if (c[i] >= 10) {
        c[i + 1] += (c[i] / 10);
        c[i] = c[i] % 10;
    }
}
if (c[len] > 0) len++; //判断是否导致最终的位数增加
for (int i = len - 1; i >= 0; --i) {   //逆序输出
    cout << c[i];
}

高精度减法

int a[100], b[100], c[100];
string s1, s2;
cin >> s1 >> s2;    //输入两个数字
//判断相减之后是否为负数
if (s1.size() < s2.size() || s1.size() == s2.size() && s1 < s2) {
    swap(s1, s2);    //交换s1和s2，保证使用s1-s2
    cout << "-";
}
int len = max(s1.size(), s2.size());
for (int i = 0; i < s1.size(); ++i) {
    a[s1.size() - 1 - i] = s1[i] - '0';    //逆序输入
}
for (int i = 0; i < s2.size(); ++i) {
    b[s2.size() - 1 - i] = s2[i] - '0';    //逆序输入
}
//处理相减
for (int i = 0; i < len; ++i) {
    if (a[i] < b[i]) { //不够减借一位
        a[i + 1]--;
        a[i] += 10;
    }
    c[i] = a[i] - b[i];
}
//删除前导0
while (c[len - 1] == 0 && len > 1) {
    len--;
}
for (int i = len - 1; i >= 0; --i) {   //逆序输出
    cout << c[i];
}

高精度乘低精度

int a[100], b, c[100];
string s;
cin >> s;    //输入大数字
cin >> b;    //输入小数字
int len = s.size();
for (int i = 0; i < s.size(); ++i) {
    a[s.size() - 1 - i] = s1[i] - '0';    //逆序输入
}
//处理相乘
for (int i = 0; i < len; ++i) {
    c[i] += a[i] * b;
    if (c[i] >= 10) {
        c[i + 1] += c[i] / 10;
        c[i] = c[i] % 10;
    }
}
//判断最终数组长度
while (c[len] > 0) {
    if (c[len] >= 10) {
        c[len + 1] += c[len] / 10;
        c[len] = c[len] % 10;
    }
    len++;
}
//删除前导0
while (c[len - 1] == 0 && len > 1) {
    len--;
}
for (int i = len - 1; i >= 0; --i) {   //逆序输出
    cout << c[i];
}

高精度乘高精度

int a[100], b[100], c[100];
string s1, s2;
cin >> s1 >> s2;
for (int i = 0; i < s1.size(); ++i) {
    a[s1.size() - 1 - i] = s1[i] - '0';    //逆序输入
}
for (int i = 0; i < s2.size(); ++i) {
    b[s2.size() - 1 - i] = s2[i] - '0';    //逆序输入
}
//处理相乘
for (int i = 0; i < s1.size(); ++i) {
    for (int j = 0; j < s2.size(); ++j) {
        int k = i + j;
        c[k] += a[i] * b[j];
        if (c[k] >= 10) {
            c[k + 1] += c[k] / 10;
            c[k] %= 10;
        }
    }
}
//判断进位进到哪里，两个数相乘，位数最多是x+y位，所以从x+y+1那里判断
int len = s1.size() + s2.size() + 1;
if (c[len - 1] > 0) {
    len++;
}
//删除前导0
while (c[len - 1] == 0 && len > 1) {
    len--;
}
for (int i = len - 1; i >= 0; --i) {   //逆序输出
    cout << c[i];
}

低精度除法高精商

a除以b，要求输出小数点后n位

int a, b, n, c[100];
c[0] = a / b;   //整数部分
int t = a % b;
for (int i = 1; i <= n; ++i) {
    c[i] = t * 10 / b;
    t = t * 10 - c[i] * b;
}
//先输出小数点前的数字以及小数点
cout << c[0] << '.';
//然后再来输出小数点后面的数
for (int i = 1; i <= n; ++i) {
    cout << c[i];
}

高精度除法低精度商

//除数用s1存放，被除数用int b存放，余数用int t存放，商用s2存放
string s1, s2;
int b, t, x;
cin >> s1 >> b;
//高精度除法用正序转存s1
for (int i = 0; i < s1.size(); ++i) {
    a[i] = s1[i] - '0';
}
//商暂时存在数组c中，长度存在int x中
t = x = 0;
//代入计算的时候要注意余数t的参与
for (int i = 0; i < s1.size(); ++i) {
    c[i] = (t * 10 + a[i]) / b; //记得加上上一个数作除法之后留下的余数
    t = (t * 10 + a[i]) % b;
}
//处理前导0
while (c[x] == 0 && x < s1.size()) {
    x++;
}
//将商存到s2中
for (int i = x; i < s1.size(); ++i) {
    s2 += c[i] + '0';
}
//输出商和余数
cout << s2 << ' ' << t;

辗转相除法最大公约数

时间复杂度为O（logb）

定理

当a与b都为正整数且a>b时，记gcd(a,b)为a与b的最大公约数，则有gcd（a,b）=gcd（b, a mod b）

证明

a可以表示成 a = kb + r（a，b，k，r皆为正整数，且r不为0）

假设d是a，b的一个公约数，则有d|a，d|b，即a和b都可以被d整除。（x|y意为kx = y，k为正整数）

而r = a - kb，两边同时除以d，r/d = a/d - kb/d，由等式右边可知m = r/d为整数，因此d|r

因此d也是b，a mod b的公约数

故（a,b）与（b, a mod b）的公约数相等，则其最大公约数也相等，得证。

举例

假如需要求 1997 和 615 两个正整数的最大公约数,用欧几里得算法，是这样进行的：
1997 / 615 = 3 （余 152）
615 / 152 = 4（余7）
152 / 7 = 21（余5）
7 / 5 = 1 （余2）
5 / 2 = 2 （余1）
2 / 1 = 2 （余0）
至此，最大公约数为1。

代码

递归版

int gcd(int n, int m) {
    if (n < m) swap(n, m);    //保证n>m
    if (n % m == 0) return m;
    return gcd(m, n % m);
}

循环版

int gcd(int a, int b) {
    int c;
    if (a < b) swap(a, b);    //保证a>b
    c = a % b;
    while (c) {
        a = b;
        b = c;
        c = a % b;
    }
    reutrn b;
}

内置函数

C++可以使用内置函数__gcd（a,b）来求两数的最大公约数，使用时需包含头文件algorithm。

蔡勒公式

用于知道年月日求对应是周几

定理

$$
D=[\frac{c}{4}]-2c+y+[\frac{y}{4}]+[\frac{13(m+1)}{5}]+d-1
$$

$$
W=D\mod 7
$$

其中：

W是星期数
D是辅助计算数（意为当前日期到原点日期一共经过多少天）
c是年份前两位
y是年份后两位
m是月份。m的取值范围是3至14，某年的1、2月看作是上一年的13、14月
d是日数
[]是取整运算（向下取整）
mod 是求余运算

证明

证明见faterazer博客

举例

计算 1994 年 12 月 13 日是星期几。显然 c = 19，y = 94，m = 12，d = 13，代入公式：

$$
D=[\frac{19}{4}]-2*19 +94+[\frac{94}{4}]+[\frac{13(12+1)}{5}]+13-1=128
$$

$$
W=128\mod7 = 2
$$

补充

在计算机编程中，W 的计算结果有可能是负数。我们需要注意，数学中的求余运算和编程中的求余运算不是完全相同的，数学上余数不能是负数，而编程中余数可以是负数。因此，在计算机中 W 是负数时，我们需要进行修正。修正方法十分简单：让 W 加上一个足够大的 7 的整数倍，使其成为正数，得到的结果再对 7 取余即可。比如 -15，我可以让其加上 70，得到 55，再除以 7 余 6，通过余数可知这一天是星期六。
此公式只适用于格里高利历（也就是现在的公历）。

海伦公式

已知三角形三条边边长a，b，c，欲求三角形面积。

先求出半周长p：

$$
p=\frac{a+b+c}{2}
$$

则三角形面积S有：

$$
S=\sqrt{p(p-a)(p-b)(p-c)}
$$

素数筛法

素数概念

素数（Prime number），又称质数，指在大于1的自然数中，除了1和该数自身外，无法被其他的自然数整除的数。1不是素数。

朴素筛法

最朴素的筛法，一个一个试。

//朴素筛法
bool isprime(int x) {
    int sq = sqrt(x); //预处理
    for (int i = 2; i <= sq + 1 && i < x; ++i) {
        if (x % i == 0) return false;
    }
    return true;
}

埃氏（Eratosthenes）筛法

假设要筛2-n内的素数，则先将2的倍数从里面剔除，再将3的倍数从里面剔除，以此类推……（小学课本里面就已经记录了这种素数筛法，~~这下我连小学生都不如了。。~~）时间复杂度为O（nloglogn），已经非常接近线性了。时间复杂度相关证明见seh_sjlj的博客。

//埃氏筛法
const int maxn = 1e6 + 10;
bool mark[maxn];
int N, sum;  //从2筛到N
cin >> N;
for (int i = 2; i <= sqrt(N); ++i) {
    //从没有被标记的数开始剔除其倍数
    if (!mark[i]) {
        //从2*i到(i-1)*i的数会被2 ~ i-1筛掉，所以从i*i开始
        for (int j = i * i; j <= N; j += i) mark[j] = 1;
    }
}
for (int i = 2; i <= N; ++i) {
    if (!mark[i]) ++sum;
}
cout << sum;  //输出总数

欧拉（Euler）筛法

欧拉筛法是埃氏筛法的改进，埃氏筛法终究会出现一个数被多个数筛掉的情况。例如因为120 = 2^3 x 3 x 5，因为2，3，5是120的质因子，所以120会被2筛一次，被3筛一次，被5筛一次，共3次。

而欧拉筛法保证了每一个合数都被其最小质因子筛去，保证不会重复筛除。故遍历一次就好，时间复杂度为O（n）。（欧拉我神！）

算法步骤视频->传送门

mark[0] = mark[1] = 1;    //特判1和0
for (int i = 2; i <= N; ++i) {
    if (!mark[i]) {     //存放素数
        primes[++pp] = i;
        ++sum;
    }
    for (int j = 1; primes[j] * i <= N; ++j) { //保证数组不会越界
        mark[primes[j] * i] = 1;    //最小质因子筛合数
        if (i % primes[j] == 0) break;
    }
}
cout << sum;  //输出总数

快速幂算法

在计算a的b次方的时候，一种朴素的算法是利用for循环进行b次操作，这样子的时间复杂度是O（b），在实际操作中复杂度还是很高。而数学性质告诉我们，a的b次方可以转化为a*a的b/2次方，所以我们不妨增大a的值，减少b的值来达到缩减计算速度的效果。快速幂算法时间复杂度是O（logb）。

普通版

typedef unsigned long long ull;
ull a, b, p;
ull quick_pow(ull a, ull b) {   //b是指数
    ull ans = 1;
    while (b) {
        if (b % 2 == 0) {    //指数是偶数
            b /= 2;
            a = a * a;
        }
        else {  //指数是奇数
            --b;
            ans = ans * a;
        }
    }
    return ans;
}

利用位运算可以简化操作，其中，奇偶性可以通过&1来判断，除以2可以通过>>1来进行。

优化版

ull quick_pow(ull a, ull b) {   //b是指数
    ull ans = 1;
    while (b) {
        if (b & 1) ans = ans * a;//奇数
        b >>= 1;  //b除以2
        a = a * a;
    }
    return ans;
}

模运算

有些时候需要进行模运算操作，例如一道经典的题目就是让我们求a^b mod p的值，我们利用求余运算的运算性质 a*b mod p = （a mod p * b mod p） mod p 来进行计算即可。

ull quick_pow(ull a, ull b) {   //b是指数
    ull ans = 1;
    while (b) {
        if (b & 1) ans = ans * a % p;//奇数
        b >>= 1;  //b除以2
        a = a * a % p;
    }
    return ans % p;
}

负进制转换

一般情况下的10进制转2进制我们是采用短除法，而倘若10进制转-2进制的时候我们需要注意，因为我们需要除以-2的缘故，导致了余数可能出现负数。例如 -15 / -2 = 7 …… -1，进制转换要求我们余数不能为负数，所以当余数出现负数的时候，我们直接让商+1，即 -15 / -2 = 8 …… 1，这样就可以满足要求。当然了，这里空白太小，写不下相对应的数学证明。（~~实际上就是我找不到，找到了估计也看不懂~~）

直接上代码，洛谷P1017

#include<iostream>
using namespace std;
int ans[1 << 21];
int main() {
    int n, r, len = 0;
    cin >> n >> r;
    cout << n << "=";
    while (n) {  //10进制转换为r进制
        int tmp1 = n % r, tmp2 = n / r;
        if (tmp1 < 0) tmp2++; //余数小于0，商加1
        ans[len++] = n - r * tmp2;
        n = tmp2;
    }
    for (int i = len - 1; i >= 0; --i) {
        if (ans[i] >= 10) {
            cout << (char)('A' + ans[i] - 10);    //特殊处理16进制
            continue;
        }
        cout << ans[i];
    }
    cout << "(base" << r << ")";
    return 0;
}

组合数求解方法

先上公式

$$
C_n^m=1,(m=0或m=n)
$$

$$
C_n^m=C_{n-1}^m+C_{n-1}^{m-1}，(n>m>0)
$$

所以，我们可以使用动态规划来求解组合数，直接上代码。

//初始化
for (int i = 0; i <= n; ++i) {
    dp[i][i] = dp[i][0] = 1;
}
for (int i = 0; i <= n; ++i) {
    for (int j = 1; j < i; ++j) {
        dp[i][j] = dp[i - 1][j] + dp[i - 1][j - 1];
    }
}

观察到每一行的组合数都只需要用到上一行组合数的数值，所以可以进行状态压缩，注意倒序处理。

for (int i = 0; i <= n; ++i) {
    for (int j = i; j >= 0; --j) {    //倒序处理
        if (j == i || j == 0) dp[j] = 1;
        else dp[j] = dp[j] + dp[j - 1];
    }
}

卡特兰（Catalan）数

讲解视频

卡特兰数最初是用来解决凸包划分三角形的问题，从该问题延伸出来的很多问题也是非常经典的组合数学问题。

数列形式

其数列形式为：1，1，2，5，14，42，132，429……

定理

设$a_n$为卡特兰数的通项公式，则其定义式为：

$$
a_n=a_0a_{n-1}+a_1a_{n-2}+a_2a_{n-3}+…+a_{n-1}a_0=\sum_{k=0}^{n-1}a_ka_{n-1-k}（其中a_0=1，a_1=1）
$$

整理可得：

$$
C(n)=a_n=\frac{C_{2n}^{n}}{n+1}
$$

另外分别有两个变式分别为：

$$
C(n)=C_{2n}^{n}-C_{2n}^{n+1}
$$

$$
C(n)=C(n-1)\frac{4n-2}{n+1}
$$

这四条公式皆为卡特兰数的通项公式，不同的情境下有不同用法。

举例

卡特兰数所涉及到的问题都非常经典，常见的问题有：

凸包划分问题：一个凸n+2边形用共n-1条对角线分割成三角形，一共有多少种方法？
圆分割问题：圆上有2n个点，两个点能连成一条弦，在保证弦不交叉的前提下，共有多少种连法？
栈问题：1-n共n个数在栈中进行出栈入栈后一共可以有多少种排列方式？
括号匹配问题：n对括号，有多少种括号匹配的括号序列？（转化成栈）
二叉树问题：n+1个叶子节点能够构成多少种不同形状的满二叉树？
黑白球问题：2n个盒子，现有n个白球和n个黑球，每个盒子只放一个球，保证任意前m个盒子里黑球数量都大于白球，共有多少种放法？

以上问题的答案都是卡特兰数。详细证明略。

代码

一种是利用$C(n)=C_{2n}^n-C_{2n}^{n+1}$这条定义式来做，需要利用动态规划来求解组合数。

//初始化
dp[0] = 1;
for (int i = 0; i <= 2 * n; ++i) {
    for (int j = i; j >= 0; --j) {
        if (j == i || j == 0) dp[j] = 1;
        else dp[j] = dp[j] + dp[j - 1];
    }
}
cout << dp[n] - dp[n + 1];    //卡特兰数

该方法在求解小型数据的时候可以适用，在求解大数据题目的时候一般涉及到取模运算，取模运算下需要使用费马小定理进行处理，笔者目前的知识水平有限，暂时不考虑该方式下的取模问题。

第二种方法是利用定义式来做，代码也是很好理解。

int n;
cin >> n;
dp[0] = dp[1] = 1;
for (int i = 2; i <= n; ++i) {
    for (int j = 0; j <= i - 1; ++j) {
        dp[i] += dp[j] * dp[i - 1 - j];
        //    dp[i]%=mod; 如果题目需要取模则直接取模
    }
}
cout << dp[n];

两个for循环就可以搞定，并且该种方式可以适用于取模运算，但是数据量太大这种方法依旧会超时，具体优化需要用到逆元，笔者目前的知识水平有限，暂时不考虑该方式下的逆元优化问题。

斯特林（Stirling）公式

公式

斯特林公式是一条用来求取n的阶乘的近似值的公式，公式为：

$$
n! \thickapprox \sqrt{2\pi n}(\frac{n}{e})^n
$$

或者更加精确的：

$$
\lim\limits_{x\rightarrow+\infty} \frac{n!}{\sqrt{2\pi n}(\frac{n}{e})^n} = 1
$$

证明过程较复杂，这里不进行展示。

使用斯特林公式求n！的位数

要求n！的位数，在C++中，包含了头文件cmath之后，只需要ceil（log10（n！））即可。

代入上述的公式易得，对于位数 f 有：

$$
f = \frac{lg2\pi}{2} + \frac{lgn}{2} + nlgn - nlge
$$

在C++中，可以使用exp（1）来求e的值，使用M_PI来求$\pi$的值。

代码：

#include<iostream>
#include<cmath>
using namespace std;
const double e = exp(1);  //e的1次方
int main() {
    int n;  //输入n来表示n的阶乘
    while (cin >> n) {
        if (n == 1) {
            cout << "1" << '\n';    //特判1
            continue;
        }
        double a = log10(2.0 * M_PI) / 2.0;
        double b = log10(n) / 2.0;
        double c = n * 1.0 * log10(n);
        double d = n * log10(e);
        int ans = ceil(a + b + c - d);
        cout << ans << '\n';
    }
    return 0;
}

康托（Cantor）展开

概念

康托展开是一个全排列到一个自然数的双射，常用于构建哈希表时的空间压缩。康托展开的实质是计算当前排列在所有由小到大全排列中的顺序，因此是可逆的。

康托展开运算

已知一个由1-n的正整数任意组合而成的排列S，则该排列S在1-n的有序全排列中在S前面的排列个数X为：

$$
X=a_1(n-1)!+a_2(n-2)!+……+a_n0!
$$

其中，$a_i$表示原排列的在第 i 位以后比$a_i$小的数的个数$(1<=i<=n)$

举例

对于一个由1-5的正整数组合而成的排列S = {3，4，1，5，2}

首先观察第一个数，如果第一位数填1，则后续的4个数是2，3，4，5的一个排列，一共是4！种，因为是1开头，比3小，所以这4！个排列一定在S之前；同样的，如果第一位数填2，则后续的4个数是1，3，4，5的一个排列，也有4！种，因为是2开头，比3小，所以这4！个排列也在S之前。换句话来讲，$a_1$ = 2，因为此时在第一位数字3后面的4位数中有两个数比3小。

而对于第二个数4，按理讲第二位如果是1，2，3任意一个都可以满足排列在S之前，但是3已经在第一位固定住了，所以第二位只能填1或2，而后面就是剩下3个数的排列。换句话来讲，此时$a_2$ = 2，因为在第二位数字4后面有2个数比4小。

类比推理，可以得到最终X = 2x4！+ 2x3！+ 0x2！+ 1x1！+ 0x0！=61，故有61个排列在S前，所以S在由1-5组合而成的排列中位列第62名。

代码

#include<iostream>
using namespace std;
int n,num[10],fac[10];
void init_fac() {    
    fac[n] = 1;    //0!=1
    for (int i = 1; i <= n - 1; ++i) {
        fac[n-i] = fac[n-i+1] * i;    //如果需要取模，则在这里取模
    }
}
int cantor() {
    int rank = 0;
    for (int i = 1; i <= n; ++i) {
        int tmp = 0;
        for (int j = i + 1; j <= n; ++j) {
            if (num[j] < num[i]) ++tmp;    //统计比ai小的数的个数
        }
        rank += (tmp * fac[i]);    //如果需要取模，这里也要取模
    }
    return rank+1;    //有rank个排列在当前排列之前，所以当前排列位于第rank+1名
}
int main() {
    cin >> n;    //由1-n的正整数组合而成的排列
    init_fac();//初始化阶乘
    for (int i = 1; i <= n; ++i) {
        cin >> num[i];    //输入排列
    }
    cout << cantor();
    return 0;
}

康托展开逆运算

由上述已知：

$$
X=a_1(n-1)!+a_2(n-2)!+……+a_n0!
$$

康托展开逆运算需要在知道X的情况下推出相对应的排列。~~（康托展开，然后逆转！）~~

举例

要计算康托展开的逆运算，我们不妨先来看一个例子：

现在有一个数45613，我们知道，这个数可以写成以下形式：

$$
45613=410^4+510^3+610^2+110^1+3*10^0
$$

观察到$10^4$，$10^3$，……，$10^0$与$(n-1)!$，$(n-2)!$，……，$0!$有异曲同工之妙。

对于45613：

45613 ÷ $10^4$ = 4余5613，那么45613的第一位数就是4
5613 ÷ $10^3$ = 5余613，那么45613的第二位数就是5
……

这样子推下去，我们就可以把45613的各个位数给取出来。那么，按照这样子的方法，回到康托运算，如果我们知道了X，就可以类比推理出对应的排列了。

例如，如果知道排名是62位，那么X = 62 - 1 = 61。接下来：

61 ÷ 4！= 2余13，说明比首位小的数有2个，则首位是3（不要忘了$a_i$表示原排列的在第 i 位以后比$a_i$小的数的个数）
13 ÷ 3！=2余1，说明比第二位小的数有2个，则第二位是4
1 ÷ 2！=0余1，说明比第三位小的数有0个，则第三位是1
1 ÷ 1！=1余0，说明比第四位小的数有1个，则第四位是5
最后一位自然填2

故最终组合S = {3，4，1，5，2}

代码

#include<iostream>
#include<cstring>
using namespace std;
int n,m,num[10],fac[10];    //m为排列中数的个数
bool used[10];    //排列中有哪些数已经被用过
void init_fac() {
    fac[m] = 1;
    for (int i = 1; i <= m - 1; ++i) {
        fac[m - i] = fac[m - i + 1] * i;
    }
}
void decantor() {
    n--;    //还原排名
    int pointer = 1;
    for (int i = 1; i <= m; ++i) {
        int quotient = n / fac[i];    //商
        int remainder = n % fac[i];    //余数
        n %= fac[i];
        int tmp = 0;    //临时计数，记录比当前数小的数的个数
        for (int j = 1; j <= m; ++j) {    //填数
            if (!used[j]) ++tmp;
            if (tmp != quotient + 1) continue;
            used[j] = 1;    //使用这个数
            num[pointer++] = j;
            break;
        }
    }
}
void print() {
    for (int i = 1; i <= m; ++i) {
        cout << num[i] << ' ';
    }
    cout << '\n';
}
int main() {
    memset(used, 0, sizeof(used));
    cin >> n >> m;
    init_fac();    //初始化阶乘
    decantor();    //康托展开逆运算
    print();    //打印排列
    return 0;
}

全排列问题

给定一个数组 S = {1，2，3，……，n}，现欲打印其所有排列方式。

C++ 的 STL 中有一个函数叫做 next_permutation 可以求全排列。具体为：bool next_permutation(iterator start,iterator end)。对应的，还有一个求前排列的 prev_permutation：bool prev_permutation(iterator start,iterator end)。接下来本节主要探讨这两个函数的底层实现。

全排列推导

假设我们现在S的总数为3，即 S = {1，2，3}。易得其部分全排列有：123，132，213，231，312，321。我们定义非递减排列是最小排列，非递增排列是最大排列。（此处不使用递增和递减的原因是考虑到若一个数列当中有重复元素出现，则不能称为严格意义上的递增递减。）

我们观察132这个排列，132后面排上213是因为对于最后两位来讲，32已经是降序排列，是最大排列。对于后两位，已经没有其他比32更大的排列了，所以这个时候，我们需要增大排序位，从最后两位排列增大到三位的排列。那么，又该怎么排三位呢？观察到132变为了213，第一位从1变到了2，并且，13这个排列对于最后两位来讲，是一个最小排列。

对此，我们可以得到一个规律：对于数组S = {$a_1$，$a_2$，$a_3$，……，$a_{i-1}$，$a_i$，$a_{i+1}$，……，$a_n$}，如果从$a_i$到$a_n$已经是一个最大排列，则将$a_{i-1}$与 “在$a_i$到$a_n$中大于$a_{i-1}$的最小的数” 进行交换。并且，$a_i$到$a_n$要变为一个最小排列。

我们可以来思考如何实现这个规律：

如何定位到$a_i$呢？很简单，注意到$a_i$到$a_n$已经是一个最大排列，前面的数会违背最大排列的排列规则，所以，$a_{i-1}$就是第一个违背最大排列的数。（$a_i > a_{i-1}$）
如何找到$a_i$到$a_n$中大于$a_{i-1}$的最小的数呢？我们通过之前的规律了解到：此时的$a_i$到$a_n$已经是一个最大排列了，而最大排列是一个降序的排列。所以，$a_i$到$a_n$中大于$a_{i-1}$的最小的数就是从右往左扫描过去的第一个大于$a_{i-1}$的数。
那么我们又如何将$a_i$到$a_n$要变为一个最小排列呢？可以证明，当我们进行上述的交换操作之后，$a_i$到$a_n$依旧是降序排列，这个时候，我们只需要将$a_i$到$a_n$逆序，即可得到其最小排列。

了解到这里，就不难推出 next_permutation 和 prev_permutation 的底层原理了。