Integer Factorization: Reversing the Multiplication

• Download RevMult_20230930 - 18.1 KB
• Download xPrompt - 741.6 KB

Background

This article is part of a set of articles. Each article highlights an aspect around Integer Factorization.

• Integer Factorization: Trial Division Algorithm[^] focuses on variants of Trial Division algorithm and their respective efficiency.
• Integer Factorization: Dreaded list of primes[^] focuses on a method to handle a large list of Primes by compression.
• Integer Factorization: Optimizing Small Factors Checking[^] focuses on a method to check small factors faster than division.
• Integer factorization: Reversing the Multiplication[^] focuses on another approach to optimize the Trial Division algorithm.

The xPrompt Download and prg Files

xPrompt.zip contains the binary of the xHarbour Scripting Engine for Windows. This is the app needed to run the prg code.

Usage:

• Copy prg file in the same directory than xPrompt
• Launch xPrompt
• Type "do FileName.prg" to run the code

xHarbour is freeware and can be dowloaded at xHarbour.org[^].

xHarbour is part of the xBase family languages: xBase - Wikipedia[^].

ub files

The .ub file are UBASIC programs files. UBASIC is a DOS executable, I use it in a VirtualBox which runs a copy of my old DOS HDD.

• UBASIC - Wikipedia[^]

Why do I use a basic interpreter?

Because UBASIC integers are up to 2600 digits, natively, no need of library, no nothing.

Introduction

When starting to play with Integer Factorization, trying all possible factors is the first idea, that algorithm is named Trial Division.

The problem is that with each division, if it doesn't lead to a factor, all the work done is lost. That annoyed me and I have searched another way.

Analyze of Situation

A study of Trial Division algorithm leads to a few observations:

• Basically, the algorithm is about testing all possible factors with a division.
• We don't care about the result of division.
• We don't care about the reminder, just if it is zero or not.
• A compound integer is the result of multiplication of at least 2 integers.

With divisions, the only optimization is about removing as many non primes as possible.

I was annoyed because an integer division was a lot of work to check a single possible factor, and I didn't like that all that work was lost after factor was checked. So I searched for ways to improve that part.

Then the idea came while studying the multiplication approach. Parts of the calculation can be mutualized across factors.

Reorder and Conquer

The twist is to take advantage of a feature of multiplication: Any least significant part of a product is fully defined by the same least significant part of both factors (multiplicand and multiplier). This mean that for any factor sharing the same least significant part, one can mutualize that least significant part calculation. One just has to work at the level of additions and shifts the multiplication is made of.

To take advantage of this, one has to reorder factors to group the factors with same left parts. And then, we are exploring a tree. As factors are reordered, sometimes we are lucky, and sometimes we are not. Trying to factorize a prime is also a worst case.

Equation Rewriting

The equation needs to be rewritten in order to show how it works.

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 

Multiplicand Multiplier ---------- Partial Product Partial Product ---------- Full Product Integer to factorize Remaining

As one can see, any green part is reused in each yellow try. Note that this animation is using base 10 for demonstration, it complicates things.

Animation comes from RMDemo10.prg in downloads.

Unoptimized UBASIC code from SPLIT_00.UB, SPLIT_05.UB contain an optimized version; but code is much longer. Both versions are recursive.

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)

Base 2

Using base 2 simplify things, and the algorithm is a binary tree exploration.

Main Algorithm

The number to factorize is PQ.

• Check if PQ is 1: In this case, 1 is a factor.
• Check if PQ is even: In this case, 2 is a factor.
• The algorithm maintains the equation is: PQ = PQh * Scale + Pl * Ql
• Start by root node with setting PQh = (PQ-1)/2, Pl = 1, Ql = 1 and Scale = 2.
• For each node:
  • Factor is found when PQh is 0 and Pl <> 1 and Ql <> 1.
  • If PQh is lower than smallest of Pl and Ql, node is a dead end.
  • If PQh is even, 2 new child nodes, new bits will be (0,0) and (1,1). Scale*2. PQh updated acordingly.
  • If PQh is odd, 2 new child nodes, new bits will be (0,1) and (1,0). Scale*2. PQh updated acordingly.
• Number is a prime when exploration ends with no factor found.

Unfortunately, the code is rather long and small snipsets are not of great interest.

Factorization Examples

Program RMTree2.prg from downloads gives a log of factorization, and generates an XML file used to show exploration as a tree.

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

To get the graphs, I use yEd, a freeware from yWorks at yEd - Graph Editor[^].

I like the neat little feature which is automatic placement, so that I don't have care of during generation.

Note that the files generated are over simplified and contain only necessary information to get the tree after beautifying with yEd.

Make the Generated Graph Readeable

• Open the generated file in directory of xPrompt MyFile.graphml
• Menu > Tools > Fit Node to Label > Ok
• Menu > Layout > Tree > Directed > Directed tab > Orientation = Top-Down > Ok

Benchmark

Here is how the Multiplication Reversion evolve with numbers. As one can see, the workload follows the Square Root curve.

For timed benchmark, I need to carry the algorithm to C++ or even assembler to take advantage of every optimization.

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

Columns are :

• Nodes are nodes in binary tree.
• - are subtractions.
• >> are Shifts.
• & are bitwise ANDs.
• | are bitwise ORs.
• Stack are stack usage in recursion.
• TD_SR are number of divisions in Trial Division (Square Root version).

Points of Interest

Starting with the 'Trial Division' algorithm and a twist, I ended with a new algorithm.

The algorithm uses only very fast processor operations like Add/Sub and shifts which scale very well with multi precision (Big Integers). Since the algorithm explores a binary tree, multi threading is efficient. If stack usage become a bottleneck, it is possible to rewrite code in a non recursive fashion.

This is an original work. I have not been able to find any reference about this algorithm, feel free to drop a comment if you know such work.

Links

• Integer factorization: Integer factorization - Wikipedia[^]
• Trial division: Trial division - Wikipedia[^]
• xHarbour: xHarbour.org[^]
• xBase family languages: xBase - Wikipedia[^]
• yWorks - The Diagramming Experts[^]
• UBASIC - Wikipedia[^]

History

• 1990s: Algorithm discovery, first experiments
• 16^th January, 2021: First version
• 31^st January, 2021: Corrections and improvement
• 30^th September, 2023: Changed Images to svg