【C/C++无聊练手（三）】计算1/1000000007的循环节（长度为1000000006）

C++代码的空间复杂度优化

在文章【C/C++无聊练手（一）】从「计算1/N的循环节」到「计算M/N的循环节」的理论推导&代码实现中，我们轻而易举地计算了千万量级输入的分数循环节。

诚然，上面的代码的确能用，时间复杂度是 $\mathcal{O}\left(N\right)$ 。然而空间复杂度近似于 $5N$ ，常数 $5$ 十分大。这使得我们在计算 1000000007 这种十亿级别的输入时，内存会被直接吃完。

让我们反思一下，空间复杂度其实分为 $N+4N$ ，那 $4N$ 是因为我们用了 int 类型的 $N$ 长度数组来计算循环节长度。如果我们通过质因数分解+欧拉定理的方式简化计算过程，那么 $4N$ 这部分空间复杂度将会被优化为 $\mathcal{O}\left(\log N\right)$ 。

因为这一部分需要用到质因数分解，所以我单独为质因数分解写了一串代码，放到上一篇文章【C/C++无聊练手（二）】将质因数分解封装为类+快速幂求幂模函数+代码正确性单元测试（本文作为下文的基础类）中

因为这部分涉及到了一些近世代数的理论知识，所以我没有在第一篇文章进行介绍，而是在本篇文章与上一篇文章中进行。

近世代数基本知识+理论推导

定义1：设 $a$ 、 $n$ 为互素的整数， $a\ne 0$ 且 $N\ne 1$，使得 $a^x\equiv 1\left(\mathrm{m}\mathrm{o}\mathrm{d}n\right)$ 的最小正整数 $x$ 称为 $a$ 模 $n$ 的阶，记为 ${\mathrm{o}\mathrm{r}\mathrm{d}}_{n}a$

我们考虑最简分式 $M/N$ ，并设 $0<M<N$ ，那么它的循环节长度和 $1/N$ 的循环节长度相同。

我们设 $N=2^{s_1}{5}^{s_2}t$ ，其中正整数 $t$ 的因子中不含 $2$ 或 $5$ 。那么有限小数部分 p = max ⁡ { s 1 ,   s 2 } p=\max\left\{s_1,\ s_2\right\} p=max{s1​, s2​} ，循环节长度为 ${\mathrm{o}\mathrm{r}\mathrm{d}}_{t}10$ ，也就是说， $\left(p,p+{\mathrm{o}\mathrm{r}\mathrm{d}}_{t}10\right]$ 为最小循环节部分。

显然， $t$ 与 $p$ 都很好求，现在问题在于如何求 ${\mathrm{o}\mathrm{r}\mathrm{d}}_{t}10$ 。根据欧拉公式 $a^{\varphi \left(n\right)}\equiv 1\left(\mathrm{m}\mathrm{o}\mathrm{d}n\right)$ ，得到 ${\mathrm{o}\mathrm{r}\mathrm{d}}_{r}10$ 是 $\varphi \left(t\right)$ 的因子，即 ${\mathrm{o}\mathrm{r}\mathrm{d}}_{t}10|\varphi \left(t\right)$ ，此时 ${\mathrm{o}\mathrm{r}\mathrm{d}}_{t}10$ 的搜索范围成功缩小到了 $\varphi \left(t\right)$ 的全体因子。

举个例子，假设 $\varphi \left(t\right)=2^a{3}^b{5}^c$ ，它的全体因子个数为 $\left(a+1\right)\left(b+1\right)\left(c+1\right)$ ，分别为 $\left(2^0,2^1,\cdots 2^a\right)\times \left(3^0,3^1,\cdots 3^b\right)\times \left(5^0,5^1,\cdots 5^c\right)$ ，其中符号 $\times$ 代表直和。

具体而言， $\varphi \left(63\right)=36=2^2{3}^2$ ，此时恰好相当于 $a=2,b=2,c=0$ ，相当于 $\varphi \left(t\right)=36$ ，那么
$$\begin{array}{rl}& \left(2^0,\cdots ,2^a\right)\times \left(3^0,\cdots ,3^b\right)\\ =& \left(2^0,2^1,2^2\right)\times \left(3^0,3^1,3^2\right)\\ =& \left(2^0{3}^0,2^0{3}^1,2^0{3}^2,2^1{3}^0,2^1{3}^1,2^1{3}^2,2^2{3}^0,2^2{3}^1,2^2{3}^2\right)\end{array}\text{\hspace{1em}(1.1)}$$

其中 $\phi \left(63\right)$ 共有 $9$ 个因子，分别是 $2^0{3}^0,2^0{3}^1,2^0{3}^2,2^1{3}^0,2^1{3}^1,2^1{3}^2,2^2{3}^0,2^2{3}^1,2^2{3}^2$ 。对于每一个因子 $s=2^a{3}^b$ ，此时，我们只需要计算一下 $10^s\left(\mathrm{m}\mathrm{o}\mathrm{d}63\right)$ ，看谁算出来是 $1$ 。

比如这里，最小的一组是 $a=1,b=1$ ，即 $10^{2^1{3}^1}\equiv 1\left(\mathrm{m}\mathrm{o}\mathrm{d}63\right)$ ，那么循环节长度就是 $2^1{3}^1=6$ 。

事实上，这里不需要搜索所有的因子，也就是说复杂度不是 $\mathcal{O}\left(\left(a+1\right)\cdot \left(b+1\right)\cdot \left(c+1\right)\right)$ ，而是 $\mathcal{O}\left(\left(a+1\right)+\left(b+1\right)+\left(c+1\right)\right)$ ，具体算法将在后面对应部分讲到。

等下，好像欧拉函数 $\varphi \left(t\right)$ 还不知道怎么算呢。别急，这里其实和质因数分解是一体的。

假设有质因数分解 $t=p_1^{s_1}{p}_2^{s_2}\cdots p_n^{s_n}$ ，那么 $\varphi \left(t\right)=t\times \left(\frac{p_1-1}{p_1}\right)\left(\frac{p_2-1}{p_2}\right)\cdots \left(\frac{p_n-1}{p_n}\right)$ 。

大家应该看出来了，这里面最复杂的内容在于质因数分解类，因此，我将这部分内容单独放到了一篇文章中进行——【C/C++无聊练手（二）】将质因数分解封装为类+快速幂求幂模函数+代码正确性单元测试（本文作为下文的基础类），其中顺便实现了快速幂模函数。

简要分析一下复杂度——

• 质因数分解我就不打算用Pollard-Rho算法这么麻烦的东西了，反正 $N$ 不大，直接暴力试除法已经够用了。时间复杂度 $\mathcal{O}\left(\sqrt{N}\right)$ ，空间复杂度 $\mathcal{O}\left(\log N\right)$ 。
• 按照因子搜索复杂度，时间复杂度 $\mathcal{O}\left(\log N\right)$ ，空间复杂度 $\mathcal{O}\left(\log N\right)$ 。
• 纯粹计算 $10^x\left(\mathrm{m}\mathrm{o}\mathrm{d}n\right)$ 也可以用快速幂实现，我们可以为每个因子提前算好，时间复杂度 $\mathcal{O}\left(\log N\right)$ ，空间复杂度 $\mathcal{O}\left(1\right)$ 。

然后总结一下循环节长度的计算步骤：

• 对于一个数 $N=2^{s_1}{5}^{s_2}t$ ，其有限小数部分 p = max ⁡ { s 1 ,   s 2 } p=\max\left\{s_1,\ s_2\right\} p=max{s1​, s2​}
• 对于 $t$ ，计算 $\varphi \left(t\right)$
• 将 $\varphi \left(t\right)$ 质因数分解为 $\varphi \left(t\right)=p_1^{s_1}{p}_2^{s_2}\cdots p_n^{s_n}$ ，令 $s=p_1^{a_1}{p}_2^{a_2}\cdots p_n^{a_n}$ ，其中 $a_1⩽s_1,\cdots ,a_n⩽s_n$ 。找出满足 $10^s\equiv 1\left(\mathrm{m}\mathrm{o}\mathrm{d}t\right)$ 最小的一组 $\left(a_1,\cdots ,a_n\right)$ ，此时的 $s$ 记为所求的 ${\mathrm{o}\mathrm{r}\mathrm{d}}_{t}10$

除掉因子2、5

回忆一下，我们要算的数是 $1/N$ ，但实际上是先进行了分解 $N=2^{s_1}{5}^{s_2}t$ ，得到有限小数部分 p = max ⁡ { s 1 ,   s 2 } p=\max\left\{s_1,\ s_2\right\} p=max{s1​, s2​} ，真正计算的是 $t$ 。所以我们先编写一个函数，输入整数 $N$ ，返回 Factor 类的 $t$ ，以及小数部分 $p$ 。

// main.cpp(part 1)
/**
 * @brief 设N=(2^s1)(5^s2)t
 *        除掉所有因子2^s1和5^s2，返回引用p=max{s1, s2}
 *        返回Factor t
 */
Factor Mod2and5 (uint32 N, uint32& p) {
    Factor      tmpT(N);
    tmpT.Factorization();

    uint32      counter2 = 0;
    uint32      counter5 = 0;
    uint8       tmpLen   = tmpT.GetLen();
    uint32      tmpPrime;

    if (tmpLen == 0){
        p = 0;
        return tmpT;
    }
    for (uint8 k = 0;k < tmpLen;++k) {
        tmpPrime = tmpT.GetPP(k).GetPrime();
        if (tmpPrime == 2) {
            counter2 = tmpT.GetPP(k).GetPow();
            tmpT.TryDiv(tmpT.GetPP(k));
        }
        if (tmpPrime == 5) {
            counter5 = tmpT.GetPP(k).GetPow();
            tmpT.TryDiv(tmpT.GetPP(k));
        }
    }

    p = counter2 > counter5 ? counter2 : counter5;
    return tmpT;
}

欧拉函数的计算

首先，根据我这篇文章【C/C++无聊练手（二）】将质因数分解封装为类+快速幂求幂模函数+代码正确性单元测试（本文作为下文的基础类）编写好 Factor.hpp 和 Factor.cpp ，作为质因数分解的基础。

根据欧拉函数定义式，假设有质因数分解 $t=p_1^{s_1}{p}_2^{s_2}\cdots p_n^{s_n}$ ，那么 $\varphi \left(t\right)=t\times \left(\frac{p_1-1}{p_1}\right)\left(\frac{p_2-1}{p_2}\right)\cdots \left(\frac{p_n-1}{p_n}\right)$ ，我们编写如下函数——

// main.cpp(part 2)
/**
 * @brief 计算整数N的EulerPhi(N)，返回一个Factor类型
 *        设N=(p1^{s1})*...(pn^{sn})
 *        EulerPhi(N)=(p1^{s1-1})*...(pn^{sn-1})(p1-1)...(pn-1)
 */
Factor EulerPhi (uint32 N) {
    Factor   tmpN(N);
    Factor   tmpMult;       // 空Factor，每次使用时记得手动SetN并Factorization
    PrimePow pp[9];         // 存(p1^1)到(pn^1)

    tmpN.Factorization();
    uint8    tmpLen = tmpN.GetLen();

    // 把每个因子取出来一个，放到pp
    for (uint8 k = 0;k < tmpLen;++k) {
        pp[k] = PrimePow(tmpN.GetFactor(k).GetPrime(), 1);
    }
    // 除掉每个pp
    for (uint8 k = 0;k < tmpLen;++k) {
        tmpN.TryDiv(pp[k]);
    }
    // 乘以每个(pp[k]-1)
    for (uint8 k = 0;k < tmpLen;++k) {
        tmpMult.SetN(pp[k].ToInt() - 1);
        tmpMult.Factorization();
        tmpN.TryMultiple(tmpMult);
    }

    return tmpN;
}

求ord10t

回忆前面所说的——

将 $\varphi \left(t\right)$ 质因数分解为 $\varphi \left(t\right)=p_1^{s_1}{p}_2^{s_2}\cdots p_n^{s_n}$ ，令 $s=p_1^{a_1}{p}_2^{a_2}\cdots p_n^{a_n}$ ，其中 $a_1⩽s_1,\cdots ,a_n⩽s_n$ 。找出满足 $10^s\equiv 1\left(\mathrm{m}\mathrm{o}\mathrm{d}t\right)$ 最小的一组 $\left(a_1,\cdots ,a_n\right)$ ，此时的 $s$ 记为所求的 ${\mathrm{o}\mathrm{r}\mathrm{d}}_{t}10$

通过递归的方法，我们很容易找出这一组最小的 $\left(a_1,a_2,\cdots a_n\right)$ 。直接看代码的注释和代码本身吧，懒得讲了——

// main.cpp(part 3)
/**
 * @brief 计算ord10t
 *        设EulerPhi(t)=(p1^s1)(p2^s2)...(pn^sn)=s
 *        减少s1，计算10^s(mod t)，直到再少后10^s(mod t)不为1为止（记此时为a1）
 *        减少s2到a2，同样的操作，做完就s3，一直到sn
 *        返回(p1^a1)(p2^a2)...(pn^an)
 */
Factor CalOrd10T (Factor ft, uint8 now, uint32 mod) {
    if (now == ft.GetLen()) {   // 算完了
        return ft;
    }
    PrimePow tmpPP(ft.GetPP(now).GetPrime(), 1);    // 准备好打算除以的因子

    ft.TryDiv(tmpPP);                               // ft /= tmpPP
    if (MyPowMod(10, ft.GetN(), mod) == 1) {    // 如果余数为1，那么继续除以自己
        return CalOrd10T(ft, now, mod);
    }
    else {      // 如果余数为不为1，那么就不应该除，赶快乘回去，然后尝试除下一个now
        ft.TryMultiple(tmpPP);
        ++ now;
        return CalOrd10T(ft, now, mod);
    }
}

求公因子

我们不希望 $M/N$ 中的 $M$ 和 $N$ 有公因子。求公因子的办法可以用辗转相除法。

// main.cpp(part 4)
/**
 * @brief 辗转相除法求公因子，输入可以有0
 */
inline uint32 GCD(uint32 a,uint32 b) {
    return b>0 ? GCD(b, a % b):a;
}

计算用的CalCycleDiv类

和文章【C/C++无聊练手（一）】从「计算1/N的循环节」到「计算M/N的循环节」的理论推导&代码实现很像，不过是按照近世代数理论优化后的实现，所以内容比较模块化。

// main.cpp(part 5)
class CalCycleDiv{
public:
    CalCycleDiv(uint32 N, uint32 M=1);                        // 构造对象，并执行计算

    void displayAsCyclicDecimal ();                     // 显示该分数的循环小数形式
    void displayAllCyclicDecimal();                     // 显示分子+分母+该分数的循环小数形式
    uint32  getM         (){return m_M;}                // 返回m_M
    uint32  getN         (){return m_N;}                // 返回m_N
    uint32  getP         (){return m_loopP;}            // 返回m_loopP
    uint32  getQ         (){return m_loopQ;}            // 返回m_loopQ
    uint32  getLoopLength(){return m_loopQ - m_loopP;}  // 返回(m_loopQ - m_loopP)，即循环节长度
    uint8   GetDecimal   (uint32 N){return m_decimal[N];}

private:
    uint32  m_N;             // N
    uint32  m_M;             // M
    uint8  *m_decimal;       // 小数部分
    uint32  m_loopP;         // 循环节范围(m_loopP, m_loopQ]
    uint32  m_loopQ;
};

CalCycleDiv::CalCycleDiv (uint32 N, uint32 M)
    : m_N (N)               // m_N = N; m_M = M;
    , m_M (M) {

    try {
	    if (N == 0)                         // 异常处理
            throw -1;
    }
	catch (int errNum) {
		std::cout << "CalCycleDiv的分母N为0，错误！" << std::endl;
	}

    DWORD time_start = GetTickCount();       // 计时
    
    M = M % N;                              // 只需要计算小数部分
    uint32 tmpGCD = (M == 0 ? N : GCD(N, M));
    M /= tmpGCD;                            // 除掉公因数
    N /= tmpGCD;
    
    // 正式开始计算
    Factor tmpT  = Mod2and5(N, m_loopP);    // N == (2^s1)(5^s2)(t), m_loopP计算完后会被赋值为max{s1, s2}！

    Factor tmpEP = EulerPhi(tmpT.GetN());   // tmpEP = EulerPhi(t)

    m_loopQ      = m_loopP + CalOrd10T(tmpEP, 0, tmpT.GetN()).GetN();

    m_decimal    = new uint8[m_loopQ];      // 小数部分
    
    uint64 tmpM = M;
    uint64 tmpN = N;

    for (uint32 k = 0; k < m_loopQ; ++k) {
        tmpM        *= 10;
        m_decimal[k] = tmpM / tmpN;
        tmpM        %= tmpN;
    }
    
    std::cout << m_M << "/" << m_N << "循环节计算时长：" << GetTickCount() - time_start << "ms" << std::endl;
}

// 显示该分数的循环小数形式
void CalCycleDiv::displayAsCyclicDecimal () {
    std::cout << m_M/m_N << ".";              // 整数部分
    
    for (uint32 i = 0;i < m_loopP;++i)          // 非循环小数部分
        std::cout << (int) m_decimal[i];

    if (m_loopP != m_loopQ) {                 // 循环小数部分
        std::cout << "[";
        for (uint32 i = m_loopP ;i < m_loopQ;++i)
            std::cout << (int) m_decimal[i];
        std::cout << "]";
    }
    std::cout << std::endl;
}

// 显示该分数的循环小数形式
void CalCycleDiv::displayAllCyclicDecimal () {
    std::cout << m_M << " / " << m_N << " == ";
    displayAsCyclicDecimal();
}

测试的main函数

使用 main 函数对结果进行测试，发现历时 5.485s 后， 1000000007 的循环节被计算出，长度为 1000000006 位（因为太长了，没有执行显示函数）

// main.cpp(part 6)
int main() {
                                                        // expected: 1/7循环节计算时长：0ms
    CalCycleDiv(7,  1  ).displayAllCyclicDecimal();     // expected: 1 / 7 == 0.[142857]
                                                        // expected: 29/14循环节计算时长：0ms
    CalCycleDiv(14, 29 ).displayAllCyclicDecimal();     // expected: 29 / 14 == 2.0[714285]
                                                        // expected: 0/14循环节计算时长：0ms
    CalCycleDiv(14, 0  ).displayAllCyclicDecimal();     // expected: 0 / 14 == 0.[0]
                                                        // expected: 0/1循环节计算时长：0ms
    CalCycleDiv(1,  0  ).displayAllCyclicDecimal();     // expected: 0 / 1 == 0.[0]
                                                        // expected: 11/1循环节计算时长：0ms
    CalCycleDiv(1,  11 ).displayAllCyclicDecimal();     // expected: 11 / 1 == 11.[0]
                                                        // expected: 100/5循环节计算时长：0ms
    CalCycleDiv(5,  100).displayAllCyclicDecimal();     // expected: 100 / 5 == 20.[0]
    CalCycleDiv ccd(1000000007, 1);                     // expected: 1/1000000007循环节计算时长：6985ms
    std::cout << ccd.getLoopLength() << std::endl;      // expected: 1000000006
    for (int k = 0; k < 100; ++k) {
        std::cout << (int)ccd.GetDecimal(k);
    }   // expected: 0000000009999999930000000489999996570000024009999831930001176489991764570057648009596463932824752470
    std::cout << std::endl;
    

    getchar();
    return 0;
}

结语

事实上，绝大多数时间都被浪费在那10亿次循环上了，实际上，我们还能进一步对代码进行优化。比如一次循环里不再是 *=10 ，而是 *=1000000000 ，这样刚好不会溢出，还能提升 9 倍的效率。还可以结合快速幂引入多线程的思想，充分利用多核的性能，这里就不再写代码演示了（反正也没人看，写这个本来就是我自娱自乐练手的）。本期练手系列文章到此结束，撒花

附录：子模块的UnitTest

void UnitTest2() {
    uint32 p1  = 0;
    Factor ft1 = Mod2and5(125, p1);
    ft1.DispCalResult();                        // expected: 1 =
    std::cout << p1 << std::endl;               // expected: 3
    ft1 = Mod2and5(250, p1);
    ft1.DispCalResult();                        // expected: 1 =
    std::cout << p1 << std::endl;               // expected: 3
    ft1 = Mod2and5(65535, p1);
    ft1.DispCalResult();                        // expected: 13107 = (3^1)(17^1)(257^1)
    std::cout << p1 << std::endl;               // expected: 1
    ft1 = Mod2and5(20, p1);
    ft1.DispCalResult();                        // expected: 1 =
    std::cout << p1 << std::endl;               // expected: 2
    ft1 = Mod2and5(50, p1);
    ft1.DispCalResult();                        // expected: 1 =
    std::cout << p1 << std::endl;               // expected: 2
    ft1 = Mod2and5(4294967295, p1);
    ft1.DispCalResult();                        // expected: 858993459 = (3^1)(17^1)(257^1)(65537^1)
    std::cout << p1 << std::endl;               // expected: 1
    ft1 = Mod2and5(127, p1);
    ft1.DispCalResult();                        // expected: 127 = (127^1)
    std::cout << p1 << std::endl;               // expected: 0

    EulerPhi(2).DispCalResult();                // expected: 1 =
    EulerPhi(3).DispCalResult();                // expected: 2 = (2^1)
    EulerPhi(4).DispCalResult();                // expected: 2 = (2^1)
    EulerPhi(2*2*2*3).DispCalResult();          // expected: 8 = (2^3)
    EulerPhi(3*3*3).DispCalResult();            // expected: 18 = (2^1)(3^2)
    EulerPhi(2*3*7).DispCalResult();            // expected: 12 = (2^2)(3^1)
    EulerPhi(2*2*3*7).DispCalResult();          // expected: 24 = (2^3)(3^1)
    EulerPhi(65535).DispCalResult();            // expected: 32768 = (2^15)
    EulerPhi(65537).DispCalResult();            // expected: 65536 = (2^16)
    EulerPhi(1000000007).DispCalResult();       // expected: 1000000006 = (2^1)(500000003^1)
    
    Factor ft2(63);
    ft2.Factorization();
    ft2.DispCalResult();                        // expected: 63 = (3^2)(7^1)
    Factor ft3(EulerPhi(ft2.GetN()));
    ft3.Factorization();
    ft3.DispCalResult();                        // expected: 36 = (2^2)(3^2)
    Factor ft4 = CalOrd10T(ft3, 0, ft2.GetN());
    ft4.DispCalResult();                        // expected: 6 = (2^1)(3^1)
    Factor ft5(1000000007);
    ft5.Factorization();
    ft5.DispCalResult();                        // expected: 1000000007 = (1000000007^1)
    Factor ft6(EulerPhi(ft5.GetN()));
    ft6.Factorization();
    ft6.DispCalResult();                        // expected: 1000000006 = (2^1)(500000003^1)
    Factor ft7 = CalOrd10T(ft6, 0, ft5.GetN());
    ft7.DispCalResult();                        // expected: 1000000006 = (2^1)(500000003^1)
    Factor ft8(9);
    ft8.Factorization();
    ft8.DispCalResult();                        // expected: 9 = (3^2)
    Factor ft9(EulerPhi(ft8.GetN()));
    ft9.Factorization();
    ft9.DispCalResult();                        // expected: 6 = (2^1)(3^1)
    Factor ft10 = CalOrd10T(ft9, 0, ft8.GetN());
    ft10.DispCalResult();                       // expected: 1 =
    Factor ft11(1);
    ft11.Factorization();
    ft11.DispCalResult();                       // expected: 1 =
    Factor ft12(EulerPhi(ft11.GetN()));
    ft12.Factorization();
    ft12.DispCalResult();                       // expected: 1 =
    Factor ft13 = CalOrd10T(ft12, 0, ft11.GetN());
    ft13.DispCalResult();                       // expected: 1 =
}