- C++代码的空间复杂度优化
- 近世代数基本知识+理论推导
- 除掉因子2、5
- 欧拉函数的计算
- 求ord10t
- 求公因子
- 计算用的CalCycleDiv类
- 测试的main函数
- 结语
- 附录：子模块的UnitTest
在文章【C/C++无聊练手（一）】从「计算1/N的循环节」到「计算M/N的循环节」的理论推导&代码实现中，我们轻而易举地计算了千万量级输入的分数循环节。
诚然，上面的代码的确能用，时间复杂度是 O ( N ) mathcal{O}left(Nright) O(N) 。然而空间复杂度近似于 5 N 5N 5N ，常数 5 5 5 十分大。这使得我们在计算 1000000007 这种十亿级别的输入时，内存会被直接吃完。
让我们反思一下，空间复杂度其实分为 N + 4 N N+4N N+4N ，那 4 N 4N 4N 是因为我们用了 int 类型的 N N N 长度数组来计算循环节长度。如果我们通过质因数分解+欧拉定理的方式简化计算过程，那么 4 N 4N 4N 这部分空间复杂度将会被优化为 O ( log N ) mathcal{O}left( log Nright) O(logN) 。
因为这一部分需要用到质因数分解，所以我单独为质因数分解写了一串代码，放到上一篇文章【C/C++无聊练手（二）】将质因数分解封装为类+快速幂求幂模函数+代码正确性单元测试（本文作为下文的基础类）中
因为这部分涉及到了一些近世代数的理论知识，所以我没有在第一篇文章进行介绍，而是在本篇文章与上一篇文章中进行。
近世代数基本知识+理论推导定义1：设 a a a 、 n n n 为互素的整数， a ≠ 0 aneq 0 a=0 且 N ≠ 1 Nneq1 N=1，使得 a x ≡ 1 ( m o d n ) a^xequiv 1left(mathrm{mod} nright) ax≡1(mod n) 的最小正整数 x x x 称为 a a a 模 n n n 的阶，记为 o r d n a mathrm{ord}_na ordna
我们考虑最简分式
M
/
N
M/N
M/N ，并设
0
<
M
<
N
0 我们设
N
=
2
s
1
5
s
2
t
N=2^{s_1}5^{s_2}t
N=2s15s2t ，其中正整数
t
t
t 的因子中不含
2
2
2 或
5
5
5 。那么有限小数部分
p
=
max
{
s
1
,
s
2
}
p=maxleft{s_1, s_2right}
p=max{s1, s2} ，循环节长度为
o
r
d
t
10
mathrm{ord}_{t}10
ordt10 ，也就是说，
(
p
,
p
+
o
r
d
t
10
]
left(p, p+mathrm{ord}_{t}10right]
(p, p+ordt10] 为最小循环节部分。 显然，
t
t
t 与
p
p
p 都很好求，现在问题在于如何求
o
r
d
t
10
mathrm{ord}_{t}10
ordt10 。根据欧拉公式
a
ϕ
(
n
)
≡
1
(
m
o
d
n
)
a^{phileft(nright)}equiv 1left(mathrm{mod} nright)
aϕ(n)≡1(mod n) ，得到
o
r
d
10
t
mathrm{ord}_{10}t
ord10t 是
ϕ
(
t
)
phileft(tright)
ϕ(t) 的因子，即
o
r
d
t
10
∣
ϕ
(
t
)
mathrm{ord}_{t}10|phileft(tright)
ordt10∣ϕ(t) ，此时
o
r
d
t
10
mathrm{ord}_{t}10
ordt10 的搜索范围成功缩小到了
ϕ
(
t
)
phileft(tright)
ϕ(t) 的全体因子。 举个例子，假设
ϕ
(
t
)
=
2
a
3
b
5
c
phileft(tright)=2^a3^b5^c
ϕ(t)=2a3b5c ，它的全体因子个数为
(
a
+
1
)
(
b
+
1
)
(
c
+
1
)
left(a+1right)left(b+1right)left(c+1right)
(a+1)(b+1)(c+1) ，分别为
(
2
0
,
2
1
,
⋯
2
a
)
×
(
3
0
,
3
1
,
⋯
3
b
)
×
(
5
0
,
5
1
,
⋯
5
c
)
left(2^0,2^1,cdots 2^aright)timesleft(3^0,3^1,cdots 3^bright)timesleft(5^0,5^1,cdots 5^cright)
(20,21,⋯2a)×(30,31,⋯3b)×(50,51,⋯5c) ，其中符号
×
times
× 代表直和。 具体而言，
ϕ
(
63
)
=
36
=
2
2
3
2
phileft(63right)=36=2^23^2
ϕ(63)=36=2232 ，此时恰好相当于
a
=
2
,
b
=
2
,
c
=
0
a=2,b=2,c=0
a=2,b=2,c=0 ，相当于
ϕ
(
t
)
=
36
phileft(tright)=36
ϕ(t)=36 ，那么 其中
φ
(
63
)
varphileft(63right)
φ(63) 共有
9
9
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
2^03^0,2^03^1,2^03^2,2^13^0,2^13^1,2^13^2,2^23^0,2^23^1,2^23^2
2030,2031,2032,2130,2131,2132,2230,2231,2232 。对于每一个因子
s
=
2
a
3
b
s=2^a3^b
s=2a3b ，此时，我们只需要计算一下
1
0
s
(
m
o
d
63
)
10^{s}left(mathrm{mod} 63right)
10s(mod 63) ，看谁算出来是
1
1
1 。 比如这里，最小的一组是
a
=
1
,
b
=
1
a=1,b=1
a=1,b=1 ，即
1
0
2
1
3
1
≡
1
(
m
o
d
63
)
10^{2^13^1}equiv 1left(mathrm{mod} 63right)
102131≡1(mod 63) ，那么循环节长度就是
2
1
3
1
=
6
2^13^1=6
2131=6 。 事实上，这里不需要搜索所有的因子，也就是说复杂度不是
O
(
(
a
+
1
)
⋅
(
b
+
1
)
⋅
(
c
+
1
)
)
mathcal{O}left(left(a+1right)cdotleft(b+1right)cdotleft(c+1right)right)
O((a+1)⋅(b+1)⋅(c+1)) ，而是
O
(
(
a
+
1
)
+
(
b
+
1
)
+
(
c
+
1
)
)
mathcal{O}left(left(a+1right)+left(b+1right)+left(c+1right)right)
O((a+1)+(b+1)+(c+1)) ，具体算法将在后面对应部分讲到。 等下，好像欧拉函数
ϕ
(
t
)
phileft(tright)
ϕ(t) 还不知道怎么算呢。别急，这里其实和质因数分解是一体的。 假设有质因数分解
t
=
p
1
s
1
p
2
s
2
⋯
p
n
s
n
t=p_1^{s_1}p_2^{s_2}cdots p_n^{s_n}
t=p1s1p2s2⋯pnsn ，那么
ϕ
(
t
)
=
t
×
(
p
1
−
1
p
1
)
(
p
2
−
1
p
2
)
⋯
(
p
n
−
1
p
n
)
phileft(tright)=ttimes 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)
ϕ(t)=t×(p1p1−1)(p2p2−1)⋯(pnpn−1) 。 大家应该看出来了，这里面最复杂的内容在于质因数分解类，因此，我将这部分内容单独放到了一篇文章中进行——【C/C++无聊练手（二）】将质因数分解封装为类+快速幂求幂模函数+代码正确性单元测试（本文作为下文的基础类），其中顺便实现了快速幂模函数。 简要分析一下复杂度—— 然后总结一下循环节长度的计算步骤： 回忆一下，我们要算的数是
1
/
N
1/N
1/N ，但实际上是先进行了分解
N
=
2
s
1
5
s
2
t
N=2^{s_1}5^{s_2}t
N=2s15s2t ，得到有限小数部分
p
=
max
{
s
1
,
s
2
}
p=maxleft{s_1, s_2right}
p=max{s1, s2} ，真正计算的是
t
t
t 。所以我们先编写一个函数，输入整数
N
N
N ，返回 Factor 类的
t
t
t ，以及小数部分
p
p
p 。 首先，根据我这篇文章【C/C++无聊练手（二）】将质因数分解封装为类+快速幂求幂模函数+代码正确性单元测试（本文作为下文的基础类）编写好 Factor.hpp 和 Factor.cpp ，作为质因数分解的基础。 根据欧拉函数定义式，假设有质因数分解
t
=
p
1
s
1
p
2
s
2
⋯
p
n
s
n
t=p_1^{s_1}p_2^{s_2}cdots p_n^{s_n}
t=p1s1p2s2⋯pnsn ，那么
ϕ
(
t
)
=
t
×
(
p
1
−
1
p
1
)
(
p
2
−
1
p
2
)
⋯
(
p
n
−
1
p
n
)
phileft(tright)=ttimes 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)
ϕ(t)=t×(p1p1−1)(p2p2−1)⋯(pnpn−1) ，我们编写如下函数—— 回忆前面所说的—— 将
ϕ
(
t
)
phileft(tright)
ϕ(t) 质因数分解为
ϕ
(
t
)
=
p
1
s
1
p
2
s
2
⋯
p
n
s
n
phileft(tright)=p_1^{s_1}p_2^{s_2}cdots p_n^{s_n}
ϕ(t)=p1s1p2s2⋯pnsn ，令
s
=
p
1
a
1
p
2
a
2
⋯
p
n
a
n
s=p_1^{a_1}p_2^{a_2}cdots p_n^{a_n}
s=p1a1p2a2⋯pnan ，其中
a
1
⩽
s
1
,
⋯
,
a
n
⩽
s
n
a_1leqslant s_1, cdots ,a_nleqslant s_n
a1⩽s1,⋯,an⩽sn 。找出满足
1
0
s
≡
1
(
m
o
d
t
)
10^{s}equiv 1left(mathrm{mod} tright)
10s≡1(mod t) 最小的一组
(
a
1
,
⋯
,
a
n
)
left(a_1,cdots ,a_nright)
(a1,⋯,an) ，此时的
s
s
s 记为所求的
o
r
d
10
t
mathrm{ord}_{10}t
ord10t 通过递归的方法，我们很容易找出这一组最小的
(
a
1
,
a
2
,
⋯
a
n
)
left(a_1,a_2,cdots a_nright)
(a1,a2,⋯an) 。直接看代码的注释和代码本身吧，懒得讲了—— 我们不希望
M
/
N
M/N
M/N 中的
M
M
M 和
N
N
N 有公因子。求公因子的办法可以用辗转相除法。 和文章【C/C++无聊练手（一）】从「计算1/N的循环节」到「计算M/N的循环节」的理论推导&代码实现很像，不过是按照近世代数理论优化后的实现，所以内容比较模块化。 使用 main 函数对结果进行测试，发现历时 5.485s 后， 1000000007 的循环节被计算出，长度为 1000000006 位（因为太长了，没有执行显示函数） 事实上，绝大多数时间都被浪费在那10亿次循环上了，实际上，我们还能进一步对代码进行优化。比如一次循环里不再是 *=10 ，而是 *=1000000000 ，这样刚好不会溢出，还能提升 9 倍的效率。还可以结合快速幂引入多线程的思想，充分利用多核的性能，这里就不再写代码演示了（反正也没人看，写这个本来就是我自娱自乐练手的）。本期练手系列文章到此结束，撒花
(
2
0
,
⋯
,
2
a
)
×
(
3
0
,
⋯
,
3
b
)
=
(
2
0
,
2
1
,
2
2
)
×
(
3
0
,
3
1
,
3
2
)
=
(
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
)
(1.1)
begin{aligned} &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^03^0,2^03^1,2^03^2,2^13^0,2^13^1,2^13^2,2^23^0,2^23^1,2^23^2 right)\ end{aligned} tag{1.1}
==(20,⋯,2a)×(30,⋯,3b)(20,21,22)×(30,31,32)(2030,2031,2032,2130,2131,2132,2230,2231,2232)(1.1)
除掉因子2、5
// main.cpp(part 1)
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;
}
欧拉函数的计算
// main.cpp(part 2)
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
// main.cpp(part 3)
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);
}
}
求公因子
// main.cpp(part 4)
inline uint32 GCD(uint32 a,uint32 b) {
return b>0 ? GCD(b, a % b):a;
}
计算用的CalCycleDiv类
// 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.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;
}
结语
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 =
}
