# 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")