[算法笔记] 数学基础

完整代码在这里。

1. 求sqrt(2)的值

基于零点定理。

考虑函数 (f(x)=x^2) ，给定 (l=1, r=2)，则：

l = mid if f(mid)<2
r = mid if f(mid)>2

重复二分直到 l 和 r 足够靠近（即小于设置的精度 EPS）

代码：

#include <iostream>
#include <math.h>
using namespace std;
const double eps = 1e-8;
double square(double x)
{
   return x * x;
}
int main()
{
   double l = 1, r = 1.5;
   double mid;
   while ((r - l) >= eps)
   {
      mid = (l + r) / 2;
      if (square(mid) > 2)
      {
         r = mid;
      }
      else
      {
         l = mid;
      }
   }
   cout << mid << endl;
   cout << sqrt(2) << endl;
}

2. 快速幂

给定 (a, b, m) ，求出 (a^b\%m) 的值（模 (m) 是为了防止溢出）。

比较暴力的做法就是循环 b 次乘以 a 。

伪代码：

for i=1 to b
{
	a = (a*a)%m
}

快速幂解法：

• 递归快速幂：当 (b) 是偶数，那么先计算 (k = a^{b/2})，显然对于 (k) 可以采取递归处理，那么最后只需要计算 (k^2)。如果 (b) 是奇数，那么 (b-1) 是偶数，即计算 (a*f(a,b-1,m))。算法复杂度为 (O(logb)) 。

• 迭代快速幂：将 (b) 看作是二进制数，例如 (b=13=(1101)_2) ，显然有 (a^{13} = a^8×a^4×a^1) ，由于二进制数中每个 bit 代表的十进制都可以由前一个 bit 的十进制计算出来（例如，(a^8 = a^4 × a^4)）。迭代快速幂主要思想是遍历 (b) 的每一个 bit 位，同时记住当前 bit 的十进制，累积之。算法复杂度 (O(32)) 。

需要特别注意的地方：

1. 输入 b=1
2. 输入 a>m

代码：

/*
    Input:a,b,m
    Output: a^(b)%m
 */
#include <iostream>
using namespace std;
uint64_t binaryPowerRec(uint64_t a, uint64_t b, uint64_t m)
{
    if (b == 0)
        return 1;
    if (b == 1)
        return a;
    else if (b & 1)
    {
        //return (binaryPowerRec(a, b / 2, m) * binaryPowerRec(a, b / 2, m) * a) % m;
        return (binaryPowerRec(a, b - 1, m) * a) % m;
    }
    else
    {
        int k = binaryPowerRec(a, b / 2, m);
        return (k * k) % m;
    }
}
uint64_t binaryPowerItor(uint64_t a, uint64_t b, uint64_t m)
{
    if (b == 0)
        return 1;
    uint64_t k = 1;
    while (b)
    {
        if (b & 1)
        {
            k = (k * a) % m;
        }
        a = (a * a) % m;
        b >>= 1;
    }
    return k;
}
int main()
{
    uint64_t a, b, m;
    a = 31, b = 31, m = 32;
    a = a % m; //if a>m
    cout << binaryPowerRec(a, b, m) << endl;
    cout << binaryPowerItor(a, b, m) << endl;
}

3. GCD 和 LCM

GCD 是最大公约数， LCM 是最小公倍数。

GCD 基于欧几里得算法：

[gcd(a,b) = gcd(b,a\%b) ]

证明也比较简单:

GCD证明

设 (a=kb+r)，其中 (k) 是 (a) 除以 (b) 的商，(r) 是 (a) 除以 (b) 的余数。
显然有 (r = a - kb) 成立。
设 (d) 是 (a) 和 (b) 的一个公约数， 则 (d) 也是 (r = a - kb = a \% b) 的一个约数。
因此：(d) 是 (a,b) 的一个公约数，也是 (b,a\%b) 的一个公约数。
由于 (d) 的任意性，所以 { (a,b) 的公约数} (subseteq) { (b, a\%b) 的公约数}。
同理可证：{ (b, a\%b) 的公约数} (subseteq) { (a,b) 的公约数}
因此：{ (a,b) 的公约数} (=) { (b, a\%b) 的公约数}
显然：(max({ a,b 的公约数}) = max({b,a\%b的公约数}))
即：(gcd(a,b) = gcd(b, a\%b))

代码实现：

#include <iostream>
using namespace std;

//最大公约数， require a>b
int gcd(int a, int b)
{
    return (b == 0) ? a : gcd(b, a % b);
}

//最小公倍数
int lcm(int a, int b)
{
    return a / gcd(a, b) * b;
}
int main()
{
    cout << gcd(10, 3) << endl;
    cout << lcm(10, 5) << endl;
    cout << lcm(lcm(3, 9), 12) << endl;
    cout << lcm(-9, -3) << endl;
}

注意 LCM ，如果是 a*b/gcd(a,b) ，a*b 可能会溢出。

4. 分数四则运算

实现分数的四则运算，《算法笔记》一书比较基础（确实不该看这部分内容，太简单了）。就当复习一下C++基础知识吧。

/*
    1. 当分数为负，约定分子为负
    2. 当分数为0， up=1,down=0
    3. 分子分母不可约
 */
#include <iostream>
#include <algorithm>
using namespace std;
//最大公约数
int gcd(int a, int b)
{
    return (b == 0) ? a : gcd(b, a % b);
}

//最小公倍数
int lcm(int a, int b)
{
    return a / gcd(a, b) * b;
}
class Fraction
{
public:
    int up, down;
    Fraction(int u = 0, int d = 1)
    {
        up = u, down = d;
        reduction();
    }
    void reduction()
    {
        if (down < 0)
        {
            up = -up, down = -down;
        }
        if (up == 0)
        {
            down = 1;
            return;
        }
        int a = abs(up), b = abs(down);
        if (a < b)
            swap(a, b);
        int k = gcd(a, b);
        up /= k;
        down /= k;
    }
    void print()
    {
        cout << up << '/' << down << endl;
    }
    void operator=(const Fraction &f)
    {
        up = f.up, down = f.down;
    }
    Fraction operator+(const Fraction &f) const
    {
        Fraction c(up * f.down + down * f.up, down * f.down);
        return c;
    }
    Fraction operator-(const Fraction &f) const
    {
        Fraction c(up * f.down - down * f.up, down * f.down);
        return c;
    }
    Fraction operator*(const Fraction &f) const
    {
        Fraction c(up * f.up, down * f.down);
        return c;
    }
    Fraction operator/(const Fraction &f) const
    {
        Fraction c(up * f.down, down * f.up);
        return c;
    }
};

int main()
{
    Fraction a(1, 2), b(1, 4);
    Fraction c;
    c = a + b, c.print();
    c = a - b, c.print();
    c = b - a, c.print();
    c = a * b, c.print();
    c = a / b, c.print();
}

5. 素数表

求解素数的经典算法就不说了。

在这里的要求是输出 (0 - n) 范围内的素数，普通算法复杂度为 (O(nsqrt{n})) 。

还有另一类经典算法叫「筛选法」，(Eratosthenes) 筛选法和欧拉筛选法。前者复杂度 (O(nln(n))) ，后者复杂度为 (O(n)) 。

(Eratosthenes) 筛选法举例：

2 3 4 5 6 7 8 9 10 (2 是素数，那么去除 2 的倍数)
2 3 5 7 9 (3 是素数，那么去除 3 的倍数)
2 3 5 7 (5 是素数，以此类推)

若果一个数没有被前面的数「筛选掉」（例如上面的 5 ），那么这个数一定是素数。

算法复杂度：

[O(sum_{i=2}^{n} frac{n}{i}) = O(nln(n)) ]

代码：

#include <iostream>
#include <vector>
#define N 100
using namespace std;
vector<int> table;
bool nonprime[N] = {false};
void findPrime()
{
    for (int i = 2; i < N; i++)
    {
        if (!nonprime[i])
        {
            table.push_back(i);
            for (int j = i + i; j < N; j += i)
                nonprime[j] = true;
        }
    }
}
int main()
{
    findPrime();
    for (int x : table)
        cout << x << " ";
}

6. 质因子分解

将一个正整数 (n) 写成若干个质数的乘积的形式。

对于任意正整数 (n) ，如果它存在 1 和本身之外的因子，那么一定是在 (sqrt{n}) 的左右成对出现。

在这里，就可推出二级结论：一个整数 (n) ，如果存在 ([2,n]) 内的因子，那么最有只有一个质因子大于 (sqrt{n})。

规定输入与输输出如下：

Input: n = 12
Output: 12 = 2^2 * 3

在这里，使用结构体 struct Factor 存储一个「质因子项」，其包括质因子和它出现的次数。

struct Factor
{
    int x;
    int count;
    Factor(int k = 0, int c = 0)
    {
        x = k, count = c;
    }
};

完整代码：

vector<Factor> decomposeFactor(int n)
{
    vector<Factor> result;
    int k = sqrt(n);
    for (int i = 2; i <= k; i++)
    {
        if (n % i == 0)
        {
            Factor f(i, 0);
            while (n % i == 0)
            {
                f.count++;
                n /= i;
            }
            result.push_back(f);
        }
    }
    if (n != 1)
    {
        result.push_back(Factor(n, 1));
    }
    return result;
}
int main()
{
    int n = 97532468;
    cin >> n;

    if (n == 1 || n == 0)
    {
        cout << n << "=" << n << endl;
        return 0;
    }

    auto v = decomposeFactor(n);
    int len = v.size();
    cout << n << "=";
    cout << v.front().x;
    if (v.front().count != 1)
    {
        cout << "^" << v.front().count;
    }
    for (int i = 1; i < len; i++)
    {
        cout << "*" << v[i].x;
        if (v[i].count != 1)
        {
            cout << '^' << v[i].count;
        }
    }
    cout << endl;
    return 0;
}

7. 大整数运算

实现大整数 (A quad op quad B) 的运算，比如「9999999999999999999 * 9」这样的等式，比较经典的模拟型的题目。

数据结构：

class BigInteger
{
    vector<int> bits;//存放各个位的数字
    //实现API
    BigInteger operator+(BigInteger &b);
    BigInteger operator-(BigInteger &b);
    BigInteger operator*(int k);
    BigInteger operator/(int k);
}

对于加法和减法，有 2 种处理办法：

• 对于A和B，先「对准长度」，再从低位往高位处理（本次采取此法）
• 对于A和B，先对较短的一部分相加，再处理剩余部分

加减法要特别注意的地方就是：处理进位和借位。

加法代码：

BigInteger operator+(BigInteger &b)
{
    // 方法0：开始对准长度
    // 方法1：先算短的，再加剩余的
    BigInteger c;
    int maxlen = max(bits.size(), b.bits.size());
    b.bits.resize(maxlen), bits.resize(maxlen);
    int t = 0, flag = 0;
    for (int i = 0; i < maxlen; i++)
    {
        t = bits[i] + b.bits[i] + flag;
        c.bits.push_back(t % 10);
        flag = t / 10;
    }
    c.bits.push_back(flag);
    format(), b.format(), c.format();
    return c;
}

减法代码：

BigInteger operator-(BigInteger &b)
{
    //require: a>=b
    BigInteger c;
    int maxlen = max(bits.size(), b.bits.size());
    bits.resize(maxlen), b.bits.resize(maxlen);
    int t = 0, flag = 0;
    for (int i = 0; i < maxlen; i++)
    {
        t = bits[i] - b.bits[i] - flag;
        if (t >= 0)
            c.bits.push_back(t), flag = 0;
        else
            c.bits.push_back(t + 10), flag = 1;
    }
    format(), b.format(), c.format();
    return c;
}

对于乘除法，我们不实现 BigInteger op BigInteger ，而是先实现 BigInteger op int。

对于乘法：

   1 4 7
×    3 5
---------
   2 4 5
 1 4 0
 3 5
---------
 5 1 4 5

过程：

1. (7 × 35 = 245)，当前位为 (5)，进位 (24)
2. (4×35+24=164)，当前位为 (4)，进位为 (16)
3. (1×35+16=51) ，当前位为 (1)，进位为 (5)
4. 当前位为进位 (5)

乘法代码：

BigInteger operator*(int k)
{
    BigInteger c;
    auto &v = bits;
    int len = v.size();
    int t = 0, flag = 0;
    for (int i = 0; i < len; i++)
    {
        t = v[i] * k + flag;
        c.bits.push_back(t % 10);
        flag = (t / 10);
    }
    while (flag)
    {
        c.bits.push_back(flag % 10);
        flag = flag / 10;
    }
    return c;
}

现在考虑如何实现 BigInteger1 * BigInteger2 :

BigInteger c("0");
for k in BigInteger2.bits: # high bits to low bits
	c = (c*10) + (BigInteger1 * k)

其中，BigInteger = BigInteger * int 和 BigInteger = BigInteger + BigInteger 均已实现。

完整代码:

BigInteger operator*(BigInteger b)
{
    BigInteger c("0");
    BigInteger t;
    int len = b.bits.size();
    for (int i = len - 1; i >= 0; i--)
    {
        t = ((*this) * b.bits[i]);
        c = (c * 10) + t;
    }
    return c;
}

下面考虑 BigInteger = BigInteger / int 的实现。

除法的过程：

  _0_1_7_6_
7/ 1 2 3 4
     7
-----------
     5 3
     4 9
-----------
       4 4
       4 2
-----------
         2

算法主要思想：

• 1 和 7 比较，不够除，该位商为0，余数 (r=1)
• 当前被除数为 (r=10r+2=12) ，商为1，余数 (r=5)
• 当前被除数为 (r=10r+3=53) ，商为7，余数 (r=4)
• 当前被除数为 (r=10r+4=44) ，商为6，余数 (r=2)

除法代码：

BigInteger operator/(int k)
{
    BigInteger c;
    int len = bits.size();
    int r = 0;
    auto &v = bits;
    c.bits.resize(len);
    for (int i = len - 1; i >= 0; i--)
    {
        r = r * 10 + v[i];
        if (r < k)
            c.bits[i] = 0;
        else
        {
            c.bits[i] = r / k;
            r = r % k;
        }
    }
    c.format();
    return c;
}

至于如何实现 BigInteger = BigInteger1 / BigInteger2 ？暂时没想到好的办法，初步想法是把 BigInteger2 分解为 2 个不超过 int的因子，即：BigInteger2 = int1 * int2 ，然后令：

BigInteger c = BigInteger1
c = c / int1
c = c / int2
return c

显然，当除数为一个超过 int 范围的质数的时候，此法无能为力，因此不是一个好的算法。如果有好的想法，欢迎留言。

8. n!​ 的质因子个数问题

现考虑一个关于 (n!) 的问题：求 (n!) 中有多少个值为 (p) 的质因子。

比如：(6!) 中有 4 个质因子 2 。

对于这个问题，比较直观的做法是计算从 [1,n] 的每个数字各有多少个质因子 (p) ，然后累计结果。时间复杂度为 (O(nlogn)) 。

代码如下：

int calFactor(int n, int p)
{
    int ans = 0;
    for (int i = 2; i <= n; i++)
    {
        int t = i;
        while (t % p == 0)
        {
            ans++;
            t /= p;
        }
    }
    return ans;
}

现在给出一种 (O(logn)) 的解法。

考虑 (10!) 中质因子 2 的个数。显然：

• (10!) 中有因子 (2^1) 的个数为 (10/2=5) .
• (10!) 中有因子 (2^2) 的个数为 (10/4=2) .
• (10!) 中有因子 (2^3) 的个数为 (10/8=1) .

所以 (10!) 中质因子 2 的个数为 5+2+1=8。

实际上，该过程可以推广为：(n!) 中有 $(frac{n}{p} + frac{n}{p^2} + frac{n}{p^3} + ...) $ 个质因子 (p) 。其中除法均为向下取整。

(该推论可以从“分数约分”的角度来证明)。

代码实现：

int calFactor2(int n, int p)
{
    int ans = 0;
    while (n)
    {
        ans += (n / p);
        n /= p;
    }
    return ans;
}

利用该算法，可以快速计算出 (n!) 末尾有多少个零：末尾 0 的个数等于 (n!) 中因子 10 的个数，而这又等于质因子 5 的个数（为什么不是质因子 2 的个数呢？因此 4, 6, 8 中也含有 2 ）。亦即计算：calFactor2(n,5)。

9. 组合数求解

众所周知，组合数具有下列性质：

[C_n^{m} = frac{n!}{m!(n-m)!} \ C_n^{m} = C_n^{n-m} \ C_n^{0} = C_n^{n} = 1 ]

当然我们可以按照式子 1 去暴力求解，时间复杂度为 (O(n)) ，但是显然这种方式需要存储三个阶乘数值，容易溢出。

但其实，组合数还有下面这一条性质：

[C_n^{m} = C_{n-1}^{m} + C_{n-1}^{m-1} ]

顺便一提该公式的含义：(C_n^{m}) 表示从 (n) 个里取 (m) 个，现在如果第 (n) 个必取出，那么有 (C_{n-1}^{m-1}) 的取法；如果第 (n) 个必然不取出，那么就有 (C_{n-1}^{m}) 种取法。

在计算机编程中，我们是十分喜欢递推公式的，因为有递推公式就意味着可以通过循环迭代进行优化。

对于上式的第一项，递归边界是 n==m；对于第二项，递归的边界是 m==0。

递归形式实现：

//方法1：递推公式-递归形式
uint64_t CombineRec(int n, int m)
{
    return (m == 0 || m == n) ? 1 : CombineRec(n - 1, m) + CombineRec(n - 1, m - 1);
}

显然，(C_{n-1}^{m-1}) 和 (C_{n-1}^{m}) 必然包含了重复的计算。（类似于斐波那契数列的过程）

通过二维数组记录计算过程，递归的第二种形式：

//方法1：递推公式-消除重复计算
const int N = 100;
uint64_t table[N][N] = {{0}};
uint64_t combineRec2(int n, int m)
{
    if (n == m || m == 0)
        return 1;
    if (table[n][m])
        return table[n][m];
    return (table[n][m] = combineRec2(n - 1, m) + combineRec2(n - 1, m - 1));
}

计算所有结果并存放到二维数组：

//方法1：递推公式-迭代形式
uint64_t caltable(int n, int m)
{
    for (int i = 0; i < N; i++)
        table[i][0] = table[i][i] = 1;
    for (int i = 2; i < N; i++)
    {
        for (int j = 0; j <= (i / 2); j++)
        {
            table[i][j] = table[i][i - j] = table[i - 1][j] + table[i - 1][j - 1];
        }
    }
    return table[n][m];
}

此外，我们可以组合数定义式有以下规律：

[egin{equation} egin{aligned} C_n^{m} &= frac{n!}{m!(n-m)!} \ &= frac{(n-m+1)(n-m+2)...(n)}{m!} \ &= frac{(n-m+1)(n-m+2)...(n-m+m)}{1*2*3...*m} \ &= frac{frac{frac{frac{n-m+1}{1}*(n-m+2)}{2}*...}{...}*(n-m+m)}{m} end{aligned} end{equation} ]

显然可以发现 k*(n-m+i)/i 这个规律。

//方法2：定义式变形
uint64_t combineNumber(int n, int m)
{
    uint64_t k = 1;
    for (int i = 1; i <= m; i++)
    {
        k = k * (n - m + i) / i;
    }
    return k;
}

10. 组合数模 p 求解

基于上一小节的组合数求解，很容易改写出以下算法：

int table[MAX][MAX] = {{0}};
const int p = 13;
int C1(int n, int m)
{
    if (m == 0 || m >= n)
        return 1;
    if (table[n][m])
        return table[n][m];
    return (table[n][m] = (C1(n - 1, m) + C1(n - 1, m - 1)) % p);
}
int C2(int n, int m)
{
    for (int i = 0; i < n; i++)
        table[i][0] = table[i][i] = 1;
    for (int i = 2; i <= n; i++)
    {
        for (int j = 0; j <= (i / 2); j++)
        {
            table[i][j] = table[i][i - j] = (table[i - 1][j] + table[i - 1][j - 1]) % p;
        }
    }
    return table[n][m];
}

上述算法对 (p) 的素数性没有要求，且允许 (p) 在 int 的范围内即可。

前面写了这么多，总不能到这里就结束了，对吧。下面接着看基于质因数分解的方法。

现在将组合数 (C_n^{m}) 进行质因子分解，设分解结果为：

[C_n^{m} = p_1^{c_1} * p_2^{c_2} * p_3^{c_3} * ... * p_k^{c_k} ]

那么：

[egin{equation} egin{aligned} C_n^{m}\%p &= (p_1^{c_1} * p_2^{c_2} * p_3^{c_3} * ... * p_k^{c_k}) \% p \ &= (p_1^{c_1} \% p)*(p_2^{c_2} \% p)*(p_3^{c_3} \% p)*...*(p_k^{c_k} \% p) \ end{aligned} end{equation} ]

(p_i^{c_i} \% p) 可通过快速幂算法来实现，时间复杂度是 (O(1)) 。

那么怎么将 (C_n^{m}) 进行质因子分解呢？考虑到 (C_n^{m} = frac{n!}{m!(n-m)!}) ，公式中存在三个阶乘，这就是上面提到的 n!的质因子个数问题 。

遍历不超过 (n) 的所有质数 (p_i) ，然后计算 (n!, m!, (n-m)!) 中分别含质因子 (p_i) 的个数 (x_i, y_i, z_i)，就可以得到 (C_n^{m}) 中质因子 (p_i) 的个数为 (x_i-y_i-z_i) 。

先来回忆一下需要用到的辅助函数。

• (n!) 的质因子 (p) 的个数，复杂度 (O(log_p{n}))

int countFactorial(int n, int p)
{
    int k = 0;
    while (n)
    {
        k += n / p;
        n /= p;
    }
    return k;
}

• 快速幂，复杂度 (O(1))

uint64_t power(int a, int b, int mod)
{
    if (b == 0)
        return 1;
    uint64_t n = 1;
    while (b)
    {
        if (b & 1)
            n = (n * a) % mod;
        b >>= 1;
        a = (a * a) % mod;
    }
    return n;
}

• 计算 (n) 以内的素数表，复杂度 (O(nln(n)))

vector<int> getPrimeTable(int n)
{
    bool nonprime[MAX] = {0};
    vector<int> v;
    for (int i = 2; i <= n; i++)
    {
        if (!nonprime[i])
        {
            v.push_back(i);
            for (int j = i + i; j <= n; j += i)
                nonprime[j] = true;
        }
    }
    return v;
}

• 组合数模 (p)

uint64_t C3(int n, int m)
{

    vector<int> primes = getPrimeTable(n);
    int x, y, z;
    uint64_t ans = 1;
    for (int p : primes)
    {
        x = countFactorial(n, p);
        y = countFactorial(m, p);
        z = countFactorial(n - m, p);
        ans = ans * power(p, x - y - z, MODP) % MODP;
    }
    return ans;
}

假设素数表里面有 (k(k≤n)) 个素数。

现在来计算一下总的时间复杂度：(O(nln(n)) + O(klog_p{n}) = O(nlogn))

11. Lucas 定理

Lucas 定理其实还是用于求解组合数模 (p) 的问题。

如果 (p) 是素数，将 (m,n) 表示为 (p) 进制：

[m = m_kp^k + ... + m_0 \ n = n_kp^k + ... + n_0 ]

Lucas 定理：

[C_n^{m} = C_{n_k}^{m_k} * ... * C_{n_0}^{m_0} ]

例如，对于 (C_8^3 \% 5) 来说，(n=8, m=3, p=5)，那么：

[m = (03)_5 \ n = (13)_5 ]

所以：

[C_8^3\%5 = (C_3^3 * c_1^0) quad MOD quad 5 = 1 ]

代码实现：

uint64_t lucas(int n, int m)
{
    if (m == 0)
        return 1;
    return (C3(n % MODP, m % MODP) * lucas(n / MODP, m / MODP)) % MODP;
}

对于 C3(n % MODP, m % MODP) 这一项，显然是对应于公式的 (C_{n_0}^{m_0})，其次是递归项 lucas(n / MODP, m / MODP) 表示的是 (p) 进制表示中，从低位到高位，可利用 2 进制这一特殊情况去代入理解。

总结

此文章写了有 2 天了，为了避免文章过长，将下面内容放在下一次的博文（提醒自己不要忘 = =）：

• 扩展欧几里得算法
• 基于扩展欧几里得算法的组合数模 (p) 问题的求解