# Prime factorization algorithm by using the multiprocessing module of Python

I may introduce a Python code of prime factorization which is from my personal project that I'm working on.

import timeit
import gmpy2
import math

A = gmpy2.mpz(11111111111)
# Input a positive integer to factor.
# Put it inside the parentheses as "A = gmpy2.mpz(....)"

def factor(A):
B = gmpy2.mpz(math.sqrt(2 * A + 0.25) - 1)
D = gmpy2.mpz(A - (B * B + B) / 2)
while (D > 0):
B += 1
D = gmpy2.mpz(A - (B * B + B) / 2)
n = gmpy2.mpz(0)
while (D != 0 and B <= A):
if (D > 0):
B += 1
D -= B
else:
n += 1
D += n
if B > A:
return output(0, 0, 0)
else:
if (B - n) % 2 == 0:
E = gmpy2.mpz((B - n) / 2)
F = gmpy2.mpz(B + n + 1)
else:
E = gmpy2.mpz(B - n)
F = gmpy2.mpz((B + n + 1) / 2)
return output(A, E, F)

def output(A, E, F):
if A == 0:
print(f"Initial Value Error: Reset the B value")
else:
if A % E != 0 or A % F != 0 or A != E * F:
print(f"[Error Occurred]  {A}  !=  {E}  *  {F}")
else:
print(f"{A}  =  {E}  *  {F}")
return 0

if __name__ == '__main__':
if A >= 2:
timer_start = timeit.default_timer()
factor(A)
timer_stop = timeit.default_timer()
running_time = round(timer_stop - timer_start, 6)
print("running time: ", running_time, " seconds")
else:
print("undefined for A<2")


If you are interested in the theoretical background of this code, please consider visiting my other question in CS stackexchange.

Since you can observe the procedure of finite sum of integers at the following part

while (D != 0 and B <= A):
if (D > 0):
B += 1
D -= B
else:
n += 1
D += n


as somewhat of a natural step, I've made a guess that wouldn't the computation time become faster by trying the following modification.

import timeit
import gmpy2
import math
import multiprocessing

A = gmpy2.mpz(math.pow(2,29)-1)
# Input a positive integer to factor.
# Put it inside the parentheses as "A = gmpy2.mpz(....)"

k = 6
# Input the number of cores to be used as "k=..."
# [!!Caveat!!] Don't use too much number of cores which exceeds the
# availability of CPU resources of your personal computer environment.
# It can damage the CPU !!!

def factor(n):
if A < 2:
print("undefined for A<2")
else:
B = gmpy2.mpz(math.sqrt(2 * A + 0.25) - 1)
D = gmpy2.mpz(A - (B * B + B) / 2)
while (D > 0):
B += 1
D = gmpy2.mpz(A - (B * B + B) / 2)
D += gmpy2.mpz(n * (n + 1) / 2)
while (D != 0 and B <= A):
if (D > 0):
B += 1
D -= B
else:
D += k * n + (k * (k + 1)) / 2
n += k
if B > A:
return output(0, 0, 0)
else:
if (B - n) % 2 == 0:
E = gmpy2.mpz((B - n) / 2)
F = gmpy2.mpz(B + n + 1)
else:
E = gmpy2.mpz(B - n)
F = gmpy2.mpz((B + n + 1) / 2)
return output(A, E, F)

def output(A, E, F):
timer_stop = timeit.default_timer()
running_time = round(timer_stop - timer_start, 6)
if A == 0:
print(f"Initial Value Error: Reset the B or k value \n")
else:
if A % E != 0 or A % F != 0 or A != E * F:
print(f"[Error Occurred]  {A}  !=  {E}  *  {F} \n")
else:
print(f"[running time: {running_time} seconds]  {A}  =  {E}  *  {F} \n")
return 0

if __name__ == '__main__':
n = []
timer_start = timeit.default_timer()
with multiprocessing.Pool(processes=k) as pool:
for x in range(0, k):
y = gmpy2.mpz(x)
n.append(y)
results = pool.map(factor, n)


!!!!Caveat!!! Don't use the code above unless you are acquainted with the CPU resources in your computer.

The code aims to add integers by k, instead of 1, for which k denotes the number of cores being used. I suspect that it is natural to assume that executing some addition of integers until meeting its target value, which is 0 in our code, shall be done faster if the added size of integers becomes bigger. For example, it is quite obvious that adding integers as 4+4+4... is faster to reach 100 than adding integers as 1+1+1...

However, my numerical experiments have shown that the latter code of the multiprocessing module is nowhere faster than the previous code of single-core. So I'm here to ask for help or any kind of advice from professionals and seniors who are experienced in dealing with parallel computation of algorithms.

I may sort my questions as follows.

1. Let's say the code of the multiprocessing module prints out any kind of result that correctly demonstrates the value of "A = a * b". In that case, even if the other cores are still running, we assume that the procedure of the algorithm is done.

In this regard, why is the heuristic computation time of the multiprocessing module version slower than the computation time of the single-core version?

2. As I've argued above, the speed of reaching the target value is assumed to be faster if the size of the constant being added becomes bigger. For example, 2+2+2... reaches 1000 faster than 1+1+1... and 4+4+4... reaches 1000 faster than 2+2+2... and 8+8+8... reaches 1000 faster than 4+4+4...

My personal computer CPU doesn't allow more than 8 cores to be concurrently operated. So I've couldn't check more than that. Still, trying my best, I've observed that the running time becomes slower when the number of used cores become near 8. Meanwhile, I've observed that the running time of the multiprocessing module version code reached its fastest peak when the number of used cores are 4 or 6. (However, none of them were faster than the running time of the single-core version code)

Would there be any plausible explanation for what I've observed? That is, why the computation time meets its peak at 4 cores or 6 cores instead of 8 cores?

3. Would the running time of the multiprocessing version code become faster eventually if the number of used cores becomes as much as one desires?

For extra questions, any tips or advice to increase the speed of the computation would be grateful. I would like to welcome any kind of opinions that could let me think differently about approaching this matter. Also, please let me know if any errors are found with my code. Before trying gmpy2 module, I encountered some overflow errors when doing computations with the original Python. I've not encountered any errors yet with my current code. So please let me know of any kind of error or counterexample if you discover it.

Thank you for reading my question.

Edit

The question has been updated. If we can port the Python code into C or C++ codes as what the answerer has suggested to the question, then we get the result what we desire.

Still, I would like to know other kinds of solutions aside from the referenced answer if there is any suitable and easy-going method to embody.

The original questions attached on this page are also still valid. If there's one, many people would learn a lot about the interaction between hardware and software by the answer, especially for those of the explanation on question 2 and 3.

If somebody has read all the referenced sources, I sincerely appreciate it. Thank you.

• If the order is not essential, could you try imap_unordered instead of map, requiring ordered results can drastically slow down parallel code. Maybe sorting the results after could still be faster. Jan 1 at 18:27
• @Caridorc: Thank you very much! You are literally almost the first and the only who gave some interest to this project! Really appreciate it. Now I am trying various methods to improve it and I also will test your method too. I will update it if there is some improvement.
– MYUN
Jan 1 at 18:34
• Also try reading here: stackoverflow.com/questions/23816546/… to choose the optimal number of workers automatically, to improve both performance on your system and portabiltiy. Jan 1 at 18:44
• @Caridorc: I thought I've read through all the relevant Stackoverflow sources that I needed,, but not yours! Thank you very much xD
– MYUN
Jan 1 at 18:46
• @Caridorc: I've tried your suggestion and in my numerical experiments, the imap_unordered method seemed to be promising when the size of argument got bigger. However, for a low size of arguements, I've couldn't observe a drastic change. So I'm yet leaving the code as it was. But still I'm working on the method which is introduced at the attacked link!
– MYUN
Jan 2 at 0:45

I may have arrived at some development.

Consider there are "30000 * 30000" number of cores. For each of those cores we distribute s=0..29999, r=0..29999, where n = s (mod 30000) and B = r (mod 30000).

Then, we put a code as follows before we start the while loop.

n = gmpy2.mpz(s)
D += (n * (n + 1)) / 2
B += r
D -= r * B - (r * (r-1)) / 2


For our original code

import timeit
import gmpy2
import math

A = gmpy2.mpz(12421772708041)
# Input a positive integer to factor.
# Put it inside the parentheses as "A = gmpy2.mpz(....)"

def factor(A):
B = gmpy2.mpz(math.sqrt(math.pow(2,1) * A + 0.25) - 1)
D = gmpy2.mpz(A - (B * B + B) / 2)
while (D > 0):
B += 1
D = gmpy2.mpz(A - (B * B + B) / 2)
print(f"interim B: {B}")
n = gmpy2.mpz(0)
while (D != 0 and B <= A):
if (D > 0):
B += 1
D -= B
else:
n += 1
D += n
if B > A:
return output(0, 0, 0)
else:
print(f"[B={B}, n={n}]")
if (B - n) % 2 == 0:
E = gmpy2.mpz((B - n) / 2)
F = gmpy2.mpz(B + n + 1)
else:
E = gmpy2.mpz(B - n)
F = gmpy2.mpz((B + n + 1) / 2)
return output(A, E, F)

def output(A, E, F):
if A == 0:
print(f"Initial Value Error: Reset the B value")
else:
if A % E != 0 or A % F != 0 or A != E * F:
print(f"[Error Occurred]  {A}  !=  {E}  *  {F}")
else:
print(f"{A}  =  {E}  *  {F}")
return 0

if __name__ == '__main__':
if A >= 2:
timer_start = timeit.default_timer()
factor(A)
timer_stop = timeit.default_timer()
running_time = round(timer_stop - timer_start, 6)
print("running time: ", running_time, " seconds")
else:
print("undefined for A<2")


it took around 1.4 second.

If we put s=9840 and r=5028, then

import timeit
import gmpy2
import math

A = gmpy2.mpz(12421772708041)

# Input a positive integer to factor.
# Put it inside the parentheses as "A = gmpy2.mpz(....)"

def factor(A):
B = gmpy2.mpz(math.sqrt(math.pow(2, 1) * A + 0.25) - 1)
D = gmpy2.mpz(A - (B * B + B) / 2)
while (D > 0):
B += 1
D = gmpy2.mpz(A - (B * B + B) / 2)
print(f"interim B: {B}")
n = gmpy2.mpz(9840)  # 8529840 = 9840 (mod 30000)
D += (n * (n + 1)) / 2
B += 5028  # 9879358-4984330 = 5028 (mod 30000)
D -= 5028 * B - (5028 * 5027) / 2
while (D != 0 and B <= A):
if (D > 0):
B += 30000
D -= 30000 * B - 449985000
else:
D += 30000 * n + 450015000
n += 30000
if B > A:
return output(0, 0, 0)
else:
print(f"[B={B}, n={n}]")
if (B - n) % 2 == 0:
E = gmpy2.mpz((B - n) / 2)
F = gmpy2.mpz(B + n + 1)
else:
E = gmpy2.mpz(B - n)
F = gmpy2.mpz((B + n + 1) / 2)
return output(A, E, F)

def output(A, E, F):
if A == 0:
print(f"Initial Value Error: Reset the B value")
else:
if A % E != 0 or A % F != 0 or A != E * F:
print(f"[Error Occurred]  {A}  !=  {E}  *  {F}")
else:
print(f"{A}  =  {E}  *  {F}")
return 0

if __name__ == '__main__':
if A >= 2:
timer_start = timeit.default_timer()
factor(A)
timer_stop = timeit.default_timer()
running_time = round(timer_stop - timer_start, 6)
print("running time: ", running_time, " seconds")
else:
print("undefined for A<2")


it took less than 1 millisecond.

If we can make automatic copies of numerous similar codes, doing the same task, as what I have referred to at this question

then we are done.

As more cores we operate concurrently, the faster we get. At the answer from the question which I referred to, he says with confidence that if we port the code into C or C++ then we can do the job what we desire.

I'm not an expertise at handling all those C++ stuffs. If somebody could embody what the answer has suggested, please consider posting your answer on this page. I'm sure it will be a valuable source of study for numerous people.

• Shy of $10^9$ cores? Will be a month or two. Jan 3 at 7:54
• Looking back on this topic after several months, I now see that I've underestimated the difficulty of managing the exponential number of cores. I'm pretty sure that managing such a large number of cores is an unrealistic task at the moment.
– MYUN
Nov 2 at 16:46