整数分解： 反转乘法

• 下载 RevMult_20230930 - 18.1 KB
• 下载 xPrompt - 741.6 KB

背景

本文是系列文章的一部分。每篇文章都重点介绍了整数因式分解的一个方面。

• 整数因式分解：试除法算法[^] 重点介绍了试除法算法的变体及其各自的效率。
• 整数因式分解：可怕的素数列表[^] 重点介绍了一种通过压缩处理大量素数的方法。
• 整数因式分解：优化小因子检查[^] 重点介绍了一种比除法更快地检查小因子方法。
• 整数因式分解：逆转乘法[^] 重点介绍了另一种优化试除法算法的方法。

xPrompt 下载和 prg 文件

xPrompt.zip 包含 Windows 版 xHarbour 脚本引擎的二进制文件。这是运行 prg 代码所需的应用程序。

用法

• 将 prg 文件复制到与 xPrompt 相同的目录中
• 启动 xPrompt
• 键入“do FileName.prg”来运行代码

xHarbour 是免费软件，可以从 xHarbour.org[^] 下载。

xHarbour 是 xBase 语言家族的一部分：xBase - 维基百科[^]。

ub 文件

.ub 文件是 UBASIC 程序文件。UBASIC 是一个 DOS 可执行文件，我在运行旧 DOS 硬盘副本的 VirtualBox 中使用它。

• UBASIC - 维基百科[^]

我为什么要使用 Basic 解释器？

因为 UBASIC 整数原生支持多达 2600 位数字，无需库，无需任何东西。

引言

当开始研究整数因式分解时，尝试所有可能的因子是第一个想法，该算法称为试除法。

问题在于，每次除法如果没有得到因子，所做的工作都会丢失。这让我很恼火，于是我寻找了另一种方法。

情况分析

对试除法算法的研究得出了一些观察结果

• 基本上，该算法是关于通过除法测试所有可能的因子。
• 我们不关心除法的结果。
• 我们不关心余数，只关心它是否为零。
• 复合整数是至少两个整数相乘的结果。

对于除法，唯一的优化是尽可能多地移除非素数。

我很恼火，因为整数除法需要大量工作来检查一个可能的因子，而且我不喜欢在检查因子后所有工作都丢失了。所以我寻找方法来改进这部分。

然后，在研究乘法方法时产生了想法。计算的某些部分可以在因子之间共享。

重排序与征服

关键在于利用乘法的一个特性：乘积的任何最低有效部分都由两个因数（被乘数和乘数）的相同最低有效部分完全定义。这意味着对于共享相同最低有效部分的任何因数，可以共享该最低有效部分的计算。只需在加法和乘法所组成的移位层面进行操作。

为了利用这一点，必须重新排序因子以将具有相同左侧部分的因子分组。然后，我们正在探索一棵树。当因子被重新排序时，有时我们很幸运，有时则不然。尝试分解一个素数也是最坏的情况。

方程重写

需要重写方程以展示其工作原理。

A compound integer can be written as:
PQ = P * Q
In order to do the reverse multiplication, I rewrite it as:
PQ = PQh * Scale + Pl * Ql
where
PQh is ( PQ - Pl * Ql ) / Scale
Scale is a power of 2 (power of 10 for demonstration)
Pl is the least significant digits of P
Ql is the least significant digits of Q 

被乘数 乘数 ---------- 部分积 部分积 ---------- 完全积 待分解整数 剩余部分

可以看出，每个黄色尝试中都重用了绿色部分。请注意，此动画使用十进制进行演示，这使事情复杂化。

动画来自下载中的 RMDemo10.prg。

来自 SPLIT_00.UB 的未优化 UBASIC 代码，SPLIT_05.UB 包含优化版本；但代码更长。两个版本都是递归的。

 1000   '
 1010   '	Spliter
 1020   fnSpliter(Produit)
 1030   ' PQ est le nombre … factoriser
 1040   '
 1050   ' vérifier les facteurs de la base
 1060   if Produit@2=0 then return(2)
 1070   if Produit@5=0 then return(5)
 1080   '
 1090   return(fnSpliter_calc(Produit,0,0,1))
 1100   '
 1110   '	Lancement du calcul
 1120   fnSpliter_calc(Pq,P,Q,Coef)
 1130   local Coefn,X,Y,Pn,Qn,PQn,Diviseur
 1140   Coefn=Coef*10
 1150   X=0
 1160      Y=0
 1170      while (Pq-(P*X+Q*Y+X*Y*Coef))@10<>0
 1180         Y+=1:if Y>9 then goto 1300
 1190      wend
 1200      inc Tour:gosub *Disp(Pq,P,Q,Coef,X,Y)
 1210      Pn=P+Y*Coef
 1220      Qn=Q+X*Coef
 1230      PQn=(Pq-(P*X+Qn*Y))\10
 1240      if PQn>0 then
 1250        :Diviseur=fnSpliter_calc(PQn,Pn,Qn,Coefn)
 1260        :if Diviseur<>1 then return(Diviseur) endif
 1270      :elseif PQn=0 then
 1280        :if Pn<>1 and Qn<>1 then return(Pn) endif
 1290      endif
 1300   X+=1:if X<=9 then goto 1160
 1310   return(1) 

基数 2

使用基数 2 简化了事情，算法是一种二叉树探索。

主要算法

待分解的数字是 PQ。

• 检查 PQ 是否为 1：在这种情况下，1 是一个因子。
• 检查 PQ 是否为偶数：在这种情况下，2 是一个因子。
• 算法维持方程为：PQ = PQh * Scale + Pl * Ql
• 从根节点开始，设置 PQh = (PQ-1)/2, Pl = 1, Ql = 1 和 Scale = 2。
• 对于每个节点
  • 当 PQh 为 0 且 Pl <> 1 且 Ql <> 1 时，找到因子。
  • 如果 PQh 小于 Pl 和 Ql 中较小的一个，则节点是死胡同。
  • 如果 PQh 是偶数，则生成 2 个新子节点，新位将是 (0,0) 和 (1,1)。Scale*2。PQh 相应更新。
  • 如果 PQh 是奇数，则生成 2 个新子节点，新位将是 (0,1) 和 (1,0)。Scale*2。PQh 相应更新。
• 当探索结束且未找到因子时，数字为素数。

不幸的是，代码相当长，小片段没有太大意义。

因式分解示例

下载中的程序 RMTree2.prg 提供了因式分解的日志，并生成一个 XML 文件，用于以树状显示探索过程。

Integer to factorize: 53 = 0b110101
Status   Action  Parent  NodeId  Equation
Root     Process      0       0  0b110101 * 0b1 + 0b0 * 0b0
New 11   Process      0       1  0b11010 * 0b10 + 0b1 * 0b1
New 11   Push         1       2  0b1011 * 0b100 + 0b11 * 0b11
New 00   Process      1       3  0b1101 * 0b100 + 0b01 * 0b01
New 01   Symetry      3       4  0b110 * 0b1000 + 0b001 * 0b101
New 10   Process      3       5  0b110 * 0b1000 + 0b101 * 0b001
Tail P   Reject       5       6  0b11 * 0b10000 + 0b0101 * 0b0001
Tail Q   Process      5       7  0b11 * 0b10000 + 0b0101 * 0b0001
Tail Q   Process      7       8  0b1 * 0b10000 + 0b10101 * 0b0001
Tail Q   Reject       8       9  0b0 * 0b10000 + 0b100101 * 0b0001
New 01   Symetry      2      10  0b100 * 0b1000 + 0b011 * 0b111
New 10   Process      2      11  0b100 * 0b1000 + 0b111 * 0b011
Tail Q   Reject      11      12  0b10 * 0b10000 + 0b0111 * 0b0011
Factorization 53 = 1 * 53

Integer to factorize: 251 = 0b11111011
Status   Action  Parent  NodeId  Equation
Root     Process      0       0  0b11111011 * 0b1 + 0b0 * 0b0
New 11   Process      0       1  0b1111101 * 0b10 + 0b1 * 0b1
New 01   Symetry      1       2  0b111110 * 0b100 + 0b01 * 0b11
New 10   Process      1       3  0b111110 * 0b100 + 0b11 * 0b01
New 11   Push         3       4  0b11011 * 0b1000 + 0b111 * 0b101
New 00   Process      3       5  0b11111 * 0b1000 + 0b011 * 0b001
New 01   Push         5       6  0b1111 * 0b10000 + 0b1011 * 0b0001
New 10   Process      5       7  0b1110 * 0b10000 + 0b1001 * 0b0011
Tail P   Reject       7       8  0b111 * 0b100000 + 0b01001 * 0b00011
Tail Q   Process      7       9  0b111 * 0b100000 + 0b01001 * 0b00011
Tail Q   Reject       9      10  0b10 * 0b100000 + 0b101001 * 0b00011
Tail P   Reject       6      11  0b10 * 0b100000 + 0b01011 * 0b10001
Tail Q   Process      6      12  0b111 * 0b100000 + 0b11011 * 0b00001
Tail Q   Process     12      13  0b11 * 0b100000 + 0b111011 * 0b00001
Tail Q   Process     13      14  0b1 * 0b100000 + 0b1011011 * 0b00001
Tail Q   Reject      14      15  0b0 * 0b100000 + 0b1111011 * 0b00001
New 01   Push         4      16  0b1011 * 0b10000 + 0b1111 * 0b0101
New 10   Process      4      17  0b1010 * 0b10000 + 0b1101 * 0b0111
Tail Q   Reject      17      18  0b101 * 0b100000 + 0b01101 * 0b00111
Tail Q   Reject      16      19  0b11 * 0b100000 + 0b11111 * 0b00101
Factorization 251 = 1 * 251

Graphs

为了获得这些图表，我使用了 yWorks 的免费软件 yEd，网址是 yEd - Graph Editor[^]。

我喜欢它的一个巧妙小功能是自动布局，这样我在生成过程中就无需操心了。

请注意，生成的文件经过过度简化，仅包含必要的信息，以便在使用 yEd 美化后获得树状图。

使生成的图表可读

• 在 xPrompt 目录中打开生成的文件 MyFile.graphml
• 菜单 > 工具 > 将节点适配到标签 > 确定
• 菜单 > 布局 > 树形 > 有向 > 有向选项卡 > 方向 = 从上到下 > 确定

基准测试

这就是乘法逆转如何随着数字演变。可以看出，工作量遵循平方根曲线。

对于定时基准测试，我需要将算法移植到 C++ 甚至汇编语言，以利用每一次优化。

Number	Nodes	-		>>		&		|		Stack	TD_SR
11		6		4		9		3		3		0		2
31		8		7		12		3		6		1		4
53		13		10		22		6		8		1		6
71		13		8		23		8		7		1		7
97		21		21		35		8		16		2		8
101		17		17		28		6		13		1		9
1009	53		55		98		26		47		7		30
2003	54		59		98		21		51		11		43
3001	86		73		158		39		61		13		53
4001	90		90		167		43		74		12		62
5003	86		84		160		35		71		19		69
6007	98		89		181		41		75		19		76
7001	123		120		219		42		98		22		82
8009	132		132		238		49		108		22		88
9001	135		133		243		52		107		21		93
10007	123		105		231		58		 90		22		99
20011	174		156		323		73		130		34		140
30011	214		210		391		85		172		39		172
40009	268		259		495		117		213		42		199
49999	264		226		507		135		198		49		222
1000003	1193	990		2325	680		854		198		999
9000011	3548	3359	6614	1499	2845	719		2999

列为

• 节点是二叉树中的节点。
• - 是减法。
• >> 是移位。
• & 是按位 AND。
• | 是按位 OR。
• 堆栈是递归中的堆栈使用情况。
• TD_SR 是试除法（平方根版本）中的除法次数。

关注点

从“试除法”算法和一些修改开始，我最终得到了一种新算法。

该算法只使用非常快的处理器操作，如加/减和移位，这些操作与多精度（大整数）配合得非常好。由于该算法探索二叉树，多线程效率很高。如果堆栈使用成为瓶颈，可以将代码重写为非递归方式。

这是一项原创工作。我未能找到有关此算法的任何参考文献，如果您知道此类工作，请随时发表评论。

链接

• 整数因式分解：整数因式分解 - 维基百科[^]
• 试除法：试除法 - 维基百科[^]
• xHarbour: xHarbour.org[^]
• xBase 家族语言：xBase - 维基百科[^]
• yWorks - 图形专家[^]
• UBASIC - 维基百科[^]

历史

• 20世纪90年代：算法发现，首次实验
• 2021年1月16日：第一版
• 2021年1月31日：修正和改进
• 2023年9月30日：图像更改为 svg