# Give as input a RSA integer (N=pq) to find p or q
# of course this will simulate what a quantum computer would do in 1 cycle (Discrete Fourier Transform), however
# a normal computer takes ages, so do not input something bigger than 4 digits, or be patient.
# The explanation of this is in my blog: https://b3ck.blogspot.com/2019/03/quantum-factorization-intuition-and-how.html
# Example run:
# $ python3 fft.py 4141
# Discrete periods via FFT to use for factorization are: [4, 20, 20, 4, 4, 20, 20, 4, 20, 4, 20, 4, 20, 20, 4, 4, 4]
# Factor for 4141 is 101
# Eduardo Duarte
# toorandom [ at] g mail.com
import sys
#number you would like to factor
# Shor's Algorithm
N=int(sys.argv[1])
#A number that does not share primes with N (a prime number might work)
x=41
#Amount of times to append the sequence x^0, x^1, ..., x^N to a signal in order to make the period more visible
#when M is higher, the fourier is more explicit on the high frequency periods
M=8
flag=0
S=[]
#Calculate period by hand (slow, just to check)
for k in range(0,M*N+1):
S.append((x**k)%N)
## if((x**k)%N == 1 and k>0 and flag<1):
## print("Period: "+str(k%(M*N)))
## flag=1
#Calculating period and multiples of period via FFT
import numpy as np
S = np.array(S)
fourierS=np.abs(np.fft.fft(S))
freqdomain = np.fft.fftfreq(S.size,d=1)
l = list(zip(freqdomain,fourierS))
s = sorted(l,key=lambda t: t[1],reverse=True)
periods = []
#Obtaining the highest amplitude frequencies that reflect a discrete period
for k in range(0,100):
if(s[k][0]>0):
ifreq = 1/s[k][0]
if( abs(round(ifreq,0) - round(ifreq,1)) < 0.1):
if(round(ifreq)%2 == 0):
periods.append(int(round(ifreq)))
print("Discrete periods via FFT to use for factorization are: "+str(periods))
#Get the factor of N via periods
if len(periods) < 1:
print("No periodicity, maybe "+str(N)+" is prime")
exit(1)
from math import gcd
for k in sorted(periods,reverse=True):
a=gcd(N,x**int(round(k/2))+1)
b=gcd(N,x**int(round(k/2))-1)
f = max(a,b)
if(f>1):
print("Factor for "+str(N)+" is "+str(f) )
flag=1
break
if(flag):
exit(1)
else:
print("No period in discrete signal, maybe "+str(N)+" is prime")