【C/C++无聊练手（二）】将质因数分解封装为类+快速幂求幂模函数+代码正确性单元测试（本文作为下文的基础类）

前言

在上篇文章【C/C++无聊练手（一）】从「计算1/N的循环节」到「计算M/N的循环节」的理论推导&代码实现中，我们对于 $1/N$ 循环节的问题实现了用 C/C++ 的计算。但在文末，笔者留下一个未解决的问题——

对于这个代码，其实我们还不够满意。因为只能计算千万级别的输入显然不够coooooool！我们希望挑战一个高难度任务，让我们能计算 1000000007 这种十亿级别的输入。不过我们会遇到一大困难——结果的长度有 1000000006 位，而我们代码的空间复杂度近似于 $5N$ ，常数 $5$ 十分大。这使得我们在计算 1000000007 这种十亿级别的输入时，内存会被直接吃完。
下下篇文章中，我们将通过近世代数的理论，将代码的空间复杂度优化到 $N$ （事实上，可以优化到 $N/2$ ，不过这样写起来麻烦些）。为了实现这样的优化，我们需要编写一个质因数分解的类，作为中间的辅助类，而这将放到下篇文章中进行

对于具体的优化方法，笔者暂且按下不表（其实就是根据欧拉定理进行优化），这部分内容将在后续的第三篇文章中提到。但笔者可以透露的是：本方法基于大量质因数分解。因此，笔者将质因数分解的代码封装为了类，放到本文单独介绍。

因为这部分代码较长，我也为每部分代码写了一个单元测试，以保证代码的正确性。

简短的理论介绍

根据算术基本定理，每个 $⩾2$ 的正整数都有唯一的质因数分解形式。例如 $9828=2^2\times 3^3\times 7^1\times 13^1$ 。因此，我们可以将 $2^2,3^3,7^1,13^1$ 每个都看作元组 $\left(2,2\right),\left(3,3\right),\left(7,1\right),\left(13,1\right)$ ，用这 $4$ 个元组来表示原正整数 $9828$ ，这一思想就是质因数分解。

根据我们问题需要，我们只需要对 $\left[1,2^{32}\right)$ 的整数进行质因数分解即可。我们知道：质数的个数是无限的。那么对于 $\left[1,2^{32}\right)$ 的整数，我们需要用无限的元组来存储它们的质因数分解吗？并不是！我们只需要 $9$ 个元组就够了。

其中的原理是——

$$\begin{gathered} 2\times 3\times 5\times 7\times 11\times 13\times 17\times 23\times 29<2^{32} \\ 2\times 3\times 5\times 7\times 11\times 13\times 17\times 23\times 29\times 31>2^{32}\text{\hspace{1em}(1.1)} \end{gathered}$$

这里就不细讲了。

Factor.hpp头文件

将整个快速幂求幂模函数+质数的幂类+因式分解类打包成一个 Factor.hpp 文件，使得整个工程的逻辑更加清晰。

// Factor.hpp
#ifndef FACTOR_HPP
#define FACTOR_HPP

#include <iostream>
#include <cmath>                    // 用sqrt函数

typedef unsigned char      uint8;   // 使用uint8表示<256的正整数，用以节省空间
typedef unsigned long      uint32;  // 无符号32位整数
typedef unsigned long long uint64;  // 无符号64位整数

#define MAX8        (uint8) 255             // 2^8-1
#define MAX32       (uint32)4294967295      // 2^32-1
#define UNIT_TEST   0                       // 0: 不测试；1: 测试1；2: 测试1~2；3: 测试1~3；4: 测试1~4

/**
 * @brief 计算整数乘方，不做结果越界检查，请小心
 * @param base  底数，范围为uint32
 * @param pow   乘方，要求pow>0，因为pow<32，故用uint8
 *              因为pow通常很小，这里不采用快速幂
 */
uint32 MyPow (uint32 base, uint8 pow);

/**
 * @brief 计算(base^pow) % mod
 * @param base  底数（在本任务中，实际约定base=10）
 *              注意：要求该参数实际数据范围在uint32内！
 *              此处设置为uint64是为了防止溢出
 *              因为uint32 * uint32 < uint64
 * @param pow   乘方，一般是一个很大的数，所以要用快速幂
 * @param mod   模，一般是一个很大的数
 *              所以要求base和ans都用uint64储存，防止溢出
 */
uint32 MyPowMod (uint64 base, uint32 pow, uint32 mod);


/**
 * @brief PrimePow用于表示一个素数的m_prime的m_pow次幂。
 * 
 *        例如PrimePow(3, 5) == 3^5 == 243，m_pow不应当为0。
 *        本类主要作为质因数分解结果储存容器，不主动执行计算。
 */
class PrimePow {
private:
    uint32  m_prime;    // 素数
    uint8   m_pow;      // 该素数的幂

public:
    // 构造函数
    PrimePow() {};                  // 空构造，别删！
    PrimePow(uint32 prime, uint8 pow) : m_prime (prime)               // 赋值构造
                                      , m_pow   (pow) { }

    // 显示方法
    uint32 ToInt    ()      {return MyPow(m_prime, m_pow);}
    uint32 GetPrime ()      {return m_prime;}
    uint8  GetPow   ()      {return m_pow;}
    void   DispAll  ()      {std::cout << "(" << m_prime << "^" << (int)m_pow << ")";}

    // 操作接口（因为该接口由上层Factor控制，溢出检查由上层执行，此处不做溢出检查）
    /**
     * @brief 尝试乘以PrimePow类元素。
     *        乘成功返回1，自己prime小返回0，自己prime大返回2
     */
    uint8  TryMultiple (PrimePow pp);
    /**
     * @brief 尝试除以PrimePow类元素。
     *        除成功返回1，不匹配返回0，除完m_pow == 0返回2
     */
    uint8  TryDiv      (PrimePow pp);
};

/**
 * @brief Factor类是质因数分解类。
 * 
 *        构造时只是将N传入，不参与计算，要求N不等于0，注意N非负！
 *        调用Factorization()后执行质因数分解，结果储存在m_factor中
 */
class Factor {
private:        // 数据成员
    uint32      m_N;            // 整数本身，记为N
    PrimePow    m_factor[9];    // 存因式分解后的结果，因为2*3*5*7*11*13*17*23*29*31>2^32，所以长度为9就够了
    uint8       m_len;          // m_factor前m_len个元素才是有效的，初始m_len == 255代表未执行因式分解

protected:      // 这些函数用于TryMultiple和TryDiv调用，不应当直接调用，因为Insert和Del都不改变N
    
    bool       InsertPP      (PrimePow pp, uint8 pos);  // 尝试插入PrimePow因子到位置pos，超长返回false
    bool       DelPP         (uint8 pos);               // 尝试删除PrimePow因子，找不到返回false

public:
    // 构造对象（不进行因式分解）
    Factor()        {m_len = 255;}  // 空构造，别删！
    Factor(uint32 N);               // 只构造对象，不因式分解，N不能为0

    // 显示方法，参数获取接口
    void       DispFactor    ();    // 显示素因子形式，未执行Factorization()返回自身
    void       DispCalResult ();    // 显示N == 素因子形式
    uint32     GetN          () {return m_N;}
    uint8      GetLen        () {return m_len;}
    PrimePow   GetPP         (uint8  k) {return m_factor[k];}
    void       SetN          (uint32 N);

    // 操作接口
    void       Factorization ();                // 执行因式分解，结果储存到m_factor和m_len中
    bool       TryMultiple   (PrimePow pp);     // 尝试乘以PrimePow类元素，溢出返回false
    bool       TryDiv        (PrimePow pp);     // 尝试除以PrimePow类元素，不匹配返回false
    bool       TryMultiple   (Factor   ft);     // 尝试乘以Factor类元素，溢出返回false
    bool       TryDiv        (Factor   ft);     // 尝试除以Factor类元素，不匹配返回false
};

#if   UNIT_TEST >  0
void UnitTest();
#endif

#endif

准备工作1：求幂、求幂模函数

如果我们将一个数因式分解为了 $125=5^3$ ，我们得到了元组 $\left(5,3\right)$ ，此时我们想根据元组计算这个数本身是什么，就应当使用 pow 函数求幂。虽然 C++ 的 cmath 库里给了 pow 函数，然而输出输出参数都是浮点数，并且考虑这里的幂次通常很小（ $<32$ ），没必要使用快速幂，所以我手动实现了一个 pow 函数。

再考虑到我们在求方程的阶的时候经常用到 $10^t\left(\mathrm{m}\mathrm{o}\mathrm{d}N\right)$ ，其中 $t$ 和 $N$ 都可能很大，于是我们采用快速幂的方法来实现求幂模函数。

具体代码如下，具体各项参数类型的设置理由，我已经在 Factor.hpp 的注释里说得很详细了——

// Factor.cpp(part 1)
#include "Factor.hpp"

/****** MyPow|MyPowMod函数 ******/

uint32 MyPow (uint32 base, uint8 pow) {
    uint32 ans = base;

    for (uint8 i = 1; i < pow; ++i) {
        ans *= base;
    }
    return ans;
}

uint32 MyPowMod (uint64 base, uint32 pow, uint32 mod) {
    uint64 ans    = 1;

    while (pow > 0) {
        if (pow & 1) {      // if pow % 2 == 1
            ans = (ans * base) % mod;
        }
        base = (base * base) % mod;
        pow /= 2;
    }

    return (uint32)ans;
}

准备工作2：PrimePow类

我们定义一个名为 PrimePow 的类来储存质因数，其中有 m_prime 和 m_pow 两个成员，分别表示底数和幂次。该类配备有各种显示函数，还有 TryMultiple 和 TryDiv 两个函数接口，以供上层的 Factor 类调用。

具体类的定义见 Factor.hpp ，大部分函数已经内联实现，只有两个函数在 Factor.cpp 里实现——

// Factor.cpp(part 2)
#include "Factor.hpp"

/****** PrimePow 类 ******/

uint8 PrimePow::TryMultiple (PrimePow pp) {
    if (m_prime == pp.GetPrime()) {
        m_pow += pp.GetPow();
        return 1;
    }
    else if (m_prime < pp.GetPrime()) {     // 底数不匹配则返回0或者2
        return 0;
    }
    else {
        return 2;
    }
}

uint8 PrimePow::TryDiv (PrimePow pp) {
    // 底数不匹配或被除数的pow比除数的pow小则返回false
    if (m_prime != pp.GetPrime() || m_pow < pp.GetPow()) {
        return 0;
    }
    else {  // 注意：此处返回2表示m_pow == 0，此时上层类应当主动删掉该节点
        m_pow -= pp.GetPow();
        return (m_pow == 0 ? 2 : 1);
    }
}

准备工作3：Factor类

我们定义一个名为 Factor 的类来储存质因数，其中有 m_N 和 PrimePow[9] 、 m_len 三个成员， m_N 就是记录下输入参数 N ， m_len 记录 PrimePow[9] 中有多少个是有效的。

我们引入了 InsertPP 和 DelPP 两个保护类成员，它们不应当被类外调用，因为它们只对 PrimePow[9] 、 m_len 执行操作，不会对 m_N 执行操作，所以只应当被上层的 TryMultiple 和 TryDiv 调用。

其余公有显示接口不必多说，主要提一下操作接口。我们构造一个 Factor 类的对象时，并不会执行质因数分解，只有在执行了 Factorization() 才会执行质因数分解。同时，如果只想把这个库用于质因数分解使用，可以调用 SetN 接口，从而不必反复构造该类的对象。

整个类的接口的调用层次是 Factor::TryDiv(Factor) --> Factor::TryDiv(PrimePow) --> PrimePow::TryDiv(PrimePow) 。

// Factor.cpp(part 3)
#include "Factor.hpp"

/****** Factor 类 ******/

Factor::Factor (uint32 N) {
    // 实际上，负数也不应当被放入Factor类，因为只是我自己用这个库，所以就没管这个事情了
    try {
        if (N == 0)                         // 异常处理
            throw -1;
    }
    catch (int errNum) {
        std::cout << "0不应当作为Factor类的参数，错误！" << std::endl;
    }

    m_N   = N;
    m_len = 255;                // m_len == 255代表未进行因式分解
}

void Factor::SetN (uint32 N) {
    // 实际上，负数也不应当被放入Factor类，因为只是我自己用这个库，所以就没管这个事情了
    try {
        if (N == 0)                         // 异常处理
            throw -1;
    }
    catch (int errNum) {
        std::cout << "0不应当作为Factor类的参数，错误！" << std::endl;
    }

    m_N   = N;
    m_len = 255;                // m_len == 255代表未进行因式分解
}

void Factor::DispFactor () {
    if (m_len == 255) {       // 如果还没执行Factorization()，则返回自身
        std::cout << m_N << "（尚未执行质因数分解）" << std::endl;
        return;
    }
    for (int i = 0; i < m_len; i++) {
        m_factor[i].DispAll();
    }
    std::cout << std::endl;
}

void Factor::DispCalResult () {
    std::cout << m_N << " = ";
    DispFactor();
}


void Factor::Factorization () {
    uint32  N       = m_N;
    uint32  sqrtN   = sqrt(N);
    uint8   counter = 0;

    m_len = 0;                          // 初始化m_len

    while (N % 2 == 0) {                // 单独算2
        N /= 2;
        ++counter;
    }
    if (counter != 0) {
        m_factor[m_len] = PrimePow(2, counter);     // k^counter
        sqrtN = sqrt(N);
        ++m_len;
    }

    counter = 0;
    for (uint32 k = 3;k <= sqrtN;k += 2) {  // 只判断奇数
        if (N % k == 0) {               // 如果k是因子
            while (N % k == 0) {
                N /= k;
                ++counter;
            }
            m_factor[m_len] = PrimePow(k, counter);     // k^counter
            sqrtN = sqrt(N);
            ++m_len;
            counter = 0;
        }
    }
    if (N != 1) {                       // 如果N不能整除sqrt(N)及以下的数，那么N是素数
        m_factor[m_len] = PrimePow(N, 1);
        ++m_len;
    }
}


bool Factor::InsertPP (PrimePow pp, uint8 pos) {
    if (m_len == MAX8) {
        std::cout << "未执行Factorization()不能对" << m_N << "进行InsertPP()" << std::endl;
        return false;
    }
    if (m_len >= 9 || pos > m_len) {
        std::cout << "InsertPP()执行错误，m_len=" << m_len << ", pos=" << pos << std::endl;
        return false;
    }
    for (int k = m_len - 1;k >= pos;--k) {      // 这里要一定int k！不然pos=0就死循环了
        m_factor[k + 1] = m_factor[k];  // 右移
    }
    ++m_len;
    m_factor[pos] = pp;
    return true;
}


bool Factor::DelPP (uint8 pos) {
    if (m_len == MAX8) {
        std::cout << "未执行Factorization()不能对" << m_N << "进行DelPP()" << std::endl;
        return false;
    }
    if (m_len == 0 || pos > m_len) {
        std::cout << "DelPP()执行错误，m_len=" << m_len << ", pos=" << pos << std::endl;
        return false;
    }
    for (uint8 k = pos + 1;k < m_len;++k) {
        m_factor[k - 1] = m_factor[k];  // 左移
    }
    --m_len;
    return true;
}


bool Factor::TryMultiple (PrimePow pp) {
    uint8 multState;
    
    if ((uint64)m_N * (uint64)pp.ToInt() > MAX32) { // 溢出
        std::cout << "TryMultiple(pp)乘法溢出，N=" << m_N << ", pp.toInt=" << pp.ToInt() << std::endl;
        return false;
    }
    
    m_N *= pp.ToInt();
    for (uint8 k = 0;k < m_len;++k) {
        multState = m_factor[k].TryMultiple(pp);
        if (multState == 1) {       // 乘法成功
            return true;
        }
        else if (multState == 2) {  // 自己的prime > pp的prime，应当插入
            return InsertPP(pp, k);
        }
    }
    return InsertPP(pp, m_len);     // 所有因子都比pp小，插入到最后
}


bool Factor::TryDiv (PrimePow pp) {
    uint8 multState;

    for (uint8 k = 0;k < m_len;++k) {
        multState = m_factor[k].TryDiv(pp);
        if (multState == 1) {       // 除法成功
            m_N /= pp.ToInt();
            return true;
        }
        else if (multState == 2) {  // 除完m_pow为0，删掉该因子
            m_N /= pp.ToInt();
            return DelPP(k);
        }
    }

    std::cout << "TryDiv(pp)除法执行失败，N=" << m_N << ", pp.toInt=" << pp.ToInt() << std::endl;
    return false;
}

bool Factor::TryMultiple (Factor ft) {
    if ((uint64)m_N * (uint64)ft.GetN() > MAX32) { // 溢出
        std::cout << "TryMultiple(ft)乘法溢出，N=" << m_N << ", ft.m_N=" << ft.GetN() << std::endl;
        return false;
    }
    for (uint8 k = 0;k < ft.GetLen();++k) {
        // 尝试乘法，如果第k次乘法失败，那么把之前的都除回去
        if (!TryMultiple(ft.GetPP(k))) {
            for (uint8 j = 0;j < k;++j) {
                TryDiv(ft.GetPP(j));
            }
            return false;
        }
    }
    return true;
}

bool Factor::TryDiv (Factor ft) {
    if ((uint64)m_N % (uint64)ft.GetN() != 0) { // 不匹配
        std::cout << "TryDiv(ft)除法溢出，N=" << m_N << ", ft.m_N=" << ft.GetN() << std::endl;
        return false;
    }
    for (uint8 k = 0;k < ft.GetLen();++k) {
        // 尝试除法，如果第k次除法失败，那么把之前的都乘回去
        if (!TryDiv(ft.GetPP(k))) { 
            for (uint8 j = 0;j < k;++j) {
                TryMultiple(ft.GetPP(j));
            }
            return false;
        }
    }
    return true;
}

正片：合并到一起进行单元测试——质因数分解类

上面写了这么长的代码，如何保证它们的正确性呢？我们在 Factor.cpp 代码加入如下单元测试代码，并将 Factor.hpp 的 #define UNIT_TEST 0 改为 #define UNIT_TEST 4 ，准备执行单元测试——

// Factor.cpp(part 4)
#if   UNIT_TEST >  0
void UnitTest() {
#if   UNIT_TEST >= 1
    std::cout << "unit test1: uint32 MyPow (uint32 base, uint8 pow);" <<std::endl;
    std::cout << MyPow(2, 31) << std::endl;         // expected: 2147483648
    std::cout << MyPow(3, 20) << std::endl;         // expected: 3486784401
    std::cout << MyPow(7, 10) << std::endl;         // expected: 282475249
    std::cout << MyPow(5, 1)  << std::endl;         // expected: 5
#endif
#if UNIT_TEST >= 2
    std::cout << "unit test2: uint32 MyPowMod (uint64 base, uint32 pow, uint32 mod);" <<std::endl;
    std::cout << MyPowMod(MAX32-2, 114514, MAX32)   << std::endl;   // expected: 262144
    std::cout << MyPowMod(10,      114514, 1919810) << std::endl;   // expected: 1199790
    std::cout << MyPowMod(114514,  1926,   817)     << std::endl;   // expected: 723
    std::cout << MyPowMod(13,      8,      1024)    << std::endl;   // expected: 33
    std::cout << MyPowMod(2,       31,     MAX32)   << std::endl;   // expected: 2147483648
    std::cout << MyPowMod(3,       20,     MAX32)   << std::endl;   // expected: 3486784401
    std::cout << MyPowMod(7,       10,     MAX32)   << std::endl;   // expected: 282475249
    std::cout << MyPowMod(5,       1,      MAX32)   << std::endl;   // expected: 5
#endif
#if UNIT_TEST >= 3
    std::cout << "unit test3: PrimePow" << std::endl;
    PrimePow    pp1(5, 2);
    PrimePow    pp2(5, 3);
    PrimePow    pp3(3, 2);
    pp1.DispAll();                                        // expected: (5^2)
    std::cout << (int)pp1.TryMultiple(pp2) << std::endl;  // expected: 1
    pp1.DispAll();                                        // expected: (5^5)
    std::cout << (int)pp1.TryDiv(pp2)      << std::endl;  // expected: 1
    pp1.DispAll();                                        // expected: (5^2)
    pp3.DispAll();                                        // expected: (3^2)
    std::cout << (int)pp1.TryMultiple(pp3) << std::endl;  // expected: 2
    pp1.DispAll();                                        // expected: (5^2)
    std::cout << (int)pp1.TryDiv(pp2)      << std::endl;  // expected: 0
    pp1.DispAll();                                        // expected: (5^2)
    // 注意，实际上应当避免出现下面的情况，执行完后应当立即调用上层类方法删除pp1
    std::cout << (int)pp1.TryDiv(pp1)      << std::endl;  // expected: 2
    pp1.DispAll();                                        // expected: (5^0)
    std::cout << std::endl;
#endif
#if UNIT_TEST >= 4
    std::cout << "unit test4: Factor" << std::endl;
    Factor      factor1(2*2*3*5*7);
    Factor      factor2(3486784401);
    factor1.DispFactor();                             // expected: 420（尚未执行质因数分解）
    factor1.DispCalResult();                          // expected: 420 = 420（尚未执行质因数分解）
    factor2.Factorization();
    factor2.DispCalResult();                          // expected: 3486784401 = (3^20)
    std::cout << (int)factor1.GetLen() << std::endl;  // expected: 255
    factor1.Factorization();
    factor1.DispFactor();                             // expected: (2^2)(3^1)(5^1)(7^1)
    factor1.DispCalResult();                          // expected: 420 == (2^2)(3^1)(5^1)(7^1)
    factor1.TryMultiple(pp2);
    factor1.DispCalResult();                          // expected: 52500 = (2^2)(3^1)(5^4)(7^1)
    factor2.TryMultiple(pp3);                         // expected: TryMultiple(pp)乘法溢出，N=3486784401, pp.toInt=9
    factor2.DispCalResult();                          // expected: 3486784401 = (3^20)
    PrimePow    pp4(5, 2);
    factor1.TryDiv(pp4);
    factor1.DispCalResult();                          // expected: 2100 = (2^2)(3^1)(5^2)(7^1)
    factor1.TryDiv(pp2);                              // expected: TryDiv(pp)除法执行失败，N=2100, pp.toInt=125
    factor1.DispCalResult();                          // expected: 2100 = (2^2)(3^1)(5^2)(7^1)
    factor1.TryDiv(pp4);
    factor1.DispCalResult();                          // expected: 84 = (2^2)(3^1)(7^1)
    factor2.TryMultiple(factor1);                     // expected: TryMultiple(ft)乘法溢出，N=3486784401, ft.m_N=84
    factor2.DispCalResult();                          // expected: 3486784401 = (3^20)
    Factor      factor3(2*3*11);                      
    factor1.TryDiv(factor3);                          // expected: TryDiv(ft)除法溢出，N=84, ft.m_N=66
    factor1.DispCalResult();                          // expected: 84 = (2^2)(3^1)(7^1)
#endif

    // 用于回显终端，输入任意键退出
    getchar();
}
#endif

然后新建一个文件，比如就叫 main.cpp 吧，执行如下代码（注意：要把在 Factor.hpp 代码倒数第三行加入上面的单元测试代码，并将 #define UNIT_TEST 0 改为 #define UNIT_TEST 4 ！）进行单元测试。

// main.cpp
#include "Factor.h"

int main() {
    std::cout << "Hello world!" << std::endl;
    
    UnitTest();

    return 0;
}

输出结果——

Hello world!
unit test1: uint32 MyPow (uint32 base, uint8 pow);
2147483648
3486784401
282475249
5
unit test2: uint32 MyPowMod (uint64 base, uint32 pow, uint32 mod);
262144
1199790
723
33
2147483648
3486784401
282475249
5
unit test3: PrimePow
(5^2)1
(5^5)1
(5^2)(3^2)2
(5^2)0
(5^2)2
(5^0)
unit test4: Factor
420（尚未执行质因数分解）
420 = 420（尚未执行质因数分解）
3486784401 = (3^20)
255
(2^2)(3^1)(5^1)(7^1)
420 = (2^2)(3^1)(5^1)(7^1)
52500 = (2^2)(3^1)(5^4)(7^1)
TryMultiple(pp)乘法溢出，N=3486784401, pp.toInt=9
3486784401 = (3^20)
2100 = (2^2)(3^1)(5^2)(7^1)
TryDiv(pp)除法执行失败，N=2100, pp.toInt=125
2100 = (2^2)(3^1)(5^2)(7^1)
84 = (2^2)(3^1)(7^1)
TryMultiple(ft)乘法溢出，N=3486784401, ft.m_N=84
3486784401 = (3^20)
TryDiv(pp)除法执行失败，N=84, pp.toInt=0
84 = (2^2)(3^1)(7^1)

单元测试正确，收工