```# This is a complete implementation of the Shor factoring algorithm.

#! /usr/bin/ruby

# Complex
require 'complex'
# log
include Math

# This version of the Shor factoring algorithm is complete. It uses the QFT with size
# a power of 2 and does the continued fraction expansion of the index of one of the non-trivial
# elements of the transform followed by a search through the expansion convergent denominators for
# the largest one less than N. It continues with as many indices as necessary until r, the order of
# f(x) = a^x (mod N), is found, or terminates if it exhausts the output indices without having
# found a value for r, i.e., a solution.

# Euclid's algorithm
# Input x_in -- one of the numbers to find the greatest common denominator of
# Input y_in -- one of the numbers to find the greatest common denominator of
# Return -- gcd(x_in, y_in)
def gcd(x_in, y_in)
x, y = x_in.abs, y_in.abs
if (x > 0)  or (y > 0)
if y == 0
y = x
elsif x > 0
x, y = y, x if x < y
r = nil
x, y = y, r while (r = x % y) > 0
end
y
end
end

# Fast and efficient algorithm for computing b^e (mod m)
# Input b -- the base, or number to be exponentiated
# Input e -- the exponent to raise the base, b, to
# Input m -- the modulus of the exponentiation problem, b^e (mod m)
# Return -- b^e (mod m)
def pow_mod(b, e, m)
result = 1
while e > 0
result = (result * b) % m if e & 1 == 1
e >>= 1
b = (b * b) % m
end
result
end

# Miller's test for primality with respect to a given base
# Input n -- the number to test for primality
# Input base -- the base for the test
# Return -- true if n is prime to base, false otherwise
def millers_test(n, base)
j, rv = 0, false
d = n_minus_1 = n - 1
while d % 2 != 1
d /= 2
j += 1
end
base_to_power = pow_mod(base, d, n)
if base_to_power == 1 or base_to_power == n_minus_1
rv = true
else
j.times do
base_to_power = pow_mod(base_to_power, 2, n)
if base_to_power == n_minus_1
rv = true
break
end
end
end

rv
end

# Rabin-Miller probabilistic test for primality
# Input n -- the number to test for primality
# Input uncertainty -- the level of uncertainty of the prime test
# Return -- true if believed to be prime, false otherwise
def rabin_miller(n, uncertainty)
rv = true
if n != 2
count, base = (log(uncertainty) / log(0.25) + 1).to_int, 2
count.times do
if !millers_test(n, base)
rv = false
break
end
base += 1
end
end

rv
end

# Given a fraction's numerator and denominator, find its continued fraction coefficients
# Input num -- the numerator of the fraction
# Input den -- the denominator of the fraction
# Return -- the array of continued fraction coefficients
def fraction_to_cf(num, den)
a = []
a <<  num / den

until num == 0
num, den = den, num
i = num / den
a << i
num -= i * den
end
a
end

def reverse_and_swap(n, count, s, swapped)
r = 0
mask = 1 << (count - 1)
0.upto(count - 1) do |i|
r += mask if (n[i] == 1)
end
if r != n
if !swapped.has_key?(n)
swapped[r] = 't'
s[2 * n + 1], s[2 * r + 1] = s[2 * r + 1], s[2 * n + 1]
s[2 * n + 2], s[2 * r + 2] = s[2 * r + 2], s[2 * n + 2]
end
end
end

def realft(samples, i_sign)
n = samples.length - 1
c1 = 0.5
theta = 2 * PI  / n
if i_sign == 1
c2 = -0.5
jns_four1(samples, 1)
else
c2 = 0.5
theta = -theta
end
hold = sin(0.5 * theta)
wpr, wpi = -2.0 * hold * hold, sin(theta)
wr, wi = 1.0 + wpr, wpi
np3 = n + 3
2.upto(n / 4) do |i|
i1 = i + i - 1
i2 = i1 + 1
i3 = np3 - i2
i4 = i3 + 1
h1r, h1i =  c1 * (samples[i1] + samples[i3]), c1 * (samples[i2] - samples[i4])
h2r, h2i = -c2 * (samples[i2] + samples[i4]), c2 * (samples[i1] - samples[i3])
samples[i1] =  h1r + wr * h2r - wi * h2i
samples[i2] =  h1i + wr * h2i + wi * h2r
samples[i3] =  h1r - wr * h2r + wi * h2i
samples[i4] = -h1i + wr * h2i + wi * h2r
hold = wr
wr += wr * wpr -   wi * wpi
wi += wi * wpr + hold * wpi
end
h1r = samples
if i_sign == 1
samples = h1r + samples
samples = h1r - samples
else
samples = c1 * (h1r + samples)
samples = c1 * (h1r - samples)
jns_four1(samples, -1)
end
end

def jns_four1(samples, i_sign)
len = len2 = samples.length
bits = 0
while len2 > 2
len2 >>= 1
bits += 1
end
swapped = Hash.new
0.upto(len / 2 - 1) do |n|
reverse_and_swap(n, bits, samples, swapped)
end
swapped = nil
n = len - 1
m_max = 2
while n > m_max
i_step = 2 * m_max
theta = 2 * PI  * i_sign / m_max
hold = sin(0.5 * theta)
wpr, wpi = -2.0 * hold * hold, sin(theta)
wr, wi = 1.0, 0.0
1.step(m_max, 2) do |m|
m.step(n, i_step) do |i|
j = i + m_max
tempr = wr * samples[j]     - wi * samples[j + 1]
tempi = wr * samples[j + 1] + wi * samples[j]
samples[j]     = samples[i]     - tempr
samples[j + 1] = samples[i + 1] - tempi
samples[i]     += tempr
samples[i + 1] += tempi
end
hold = wr
wr += wr * wpr - wi   * wpi
wi += wi * wpr + hold * wpi
end
m_max = i_step
end
end

# Routine period_finder fakes "quantum entanglement" in O(r) time,
# where r is the repeat period.
# Input a -- the base for f(x) = a^x (mod N)
# Input n -- N in f(x) = a^x (mod N)
# Return -- r, the order of a (mod N), where f(x + r) = f(x) (mod N); also a^r = 1 (mod N)
def period_finder_1(a, n)
x, f = 0, 1
loop do
x += 1
f = (f * a) % n
break if f == 1
end
x
end

# Find the smallest power of 2 >= n * n
# Input n -- the number to be factored
# Return -- n * n <= 2^q < 2 * n * n
def find_q(n)
y = 2 * log(n) / log(2.0)
2**(y.to_i() + 1)
end

# The FFT used, realft, only outputs the first half of the transform plus the Nyquist value
# Input inp -- the output of the FFT, which is only half the transform, and
#              with the Nyquist value in the second position
# Return -- the output array as Complexes
def build_outp(inp)
rv = []
nyq = inp
inp = 0.0

0.step(inp.length - 1, 2) { |i| rv << Complex(inp[i], inp[i + 1]) }
rv << Complex(nyq, 0.0)
k = rv.length - 2
k.downto(1) { |i| rv << rv[i].conjugate }
rv
end

# Collect the (non-zero) qubit indices after the Fourier transform
# Input outp -- the Fourier transform
# Input max_factor -- the factor for the threshold of non-zero interference
# Return output_indices -- the qubits in the transform
def build_output_indices_array(outp, max_factor)
rv = []
almost_mx = max_factor * outp.max.polar
# Skip the first element, which is always 0.
1.upto(outp.length - 1) { |i| rv << i if outp[i].polar > almost_mx }
rv
end

# From the continued fraction coefficients recursively generate the denominators of the convergents
# Input a -- array of continued fraction coefficients
# Input n -- the number to factor
# Return -- the minimum value of the period, r, from the continued fraction convergent denominators
def find_r_min(a, n)
k, k_back_1, k_back_2 = nil, 0, 1
done, rv = false, nil

a.each do |el|
k = el * k_back_1 + k_back_2
if k > n
rv = k_back_1
done = true
break
end
k_back_2, k_back_1 = k_back_1, k
end
rv = k if !done
rv
end

# Try integer multiples of r_min for r until a^r = 1 (mod n)
# Input r_min -- the smallest possible value of the period, r
# Input a -- the base
# Input n -- the modulus
# Return -- an array, the first element of which is true only if an r < n could be found
#           and the second is the order of a mod n, i.e., a^r = 1 (mod n), if found
def find_r(r_min, a, n)
rv, found = r_min, true
until pow_mod(a, rv, n) == 1
rv += r_min
if rv >= n
found = false
break
end
end
return found, rv
end

Uncertainty, max_factor, inp, outp= 0.01, 0.9, nil, nil
have_p, have_q, all_done = false, false, false

begin
p, q = nil, nil
loop do
if !have_p
STDOUT.write "Enter first prime (< 50): "
p = gets().to_i()
if p < 2 or p >= 50
puts "#{p} is out of range"
redo
elsif (p != 2 and p != 3 and p != 5 and p != 7) and !rabin_miller(p, Uncertainty)
puts "#{p} is not believed to be prime"
redo
else
have_p = true
end
end
if !have_q
loop do
STDOUT.write "Enter second prime (< 50): "
q = gets().to_i()
if q < 2 or q >= 50
puts "#{q} is out of range"
redo
elsif rabin_miller(q, Uncertainty)
have_q = true
break
elsif q == 2 or q == 3 or q == 5 or q == 7
have_q = true
break
end
puts "#{q} is not believed to be prime"
end
end
n = p * q
a = 0
loop do
STDOUT.write "Enter a base number < #{n}: "
a = gets().to_i()
break if a > 1 and a < n
puts "#{a} is out of range"
end

g = gcd(a, n)
if g > 1
puts "N has the non trivial factor #{g}"
break
else
r = period_finder_1(a, n)
puts "Order of #{a} (mod #{n}) = #{r} (using fake \"quantum entanglement\")"
if r % 2 == 1
puts 'Index is odd: need to try another base.'
redo
elsif pow_mod(a, r / 2, n) == n - 1
puts 'Trivial square root: need to try another base'
redo
else
q = find_q(n)
puts "q = #{q}"
puts "There are #{r} possible values for the measurement of f(x)"
rand_ind = rand(r)
puts "The program randomly picked f(#{rand_ind}) = #{a}^#{rand_ind} (mod #{n}) = #{pow_mod(a, rand_ind, n)}"
# Build the input array as it would look after measuring f(x) = a^x (mod N) and entanglement.
inp = Array.new(q) { |i| i % r == rand_ind ? 1.0: 0.0 }
inp.unshift(-1.0)
realft(inp, 1)
inp.shift
outp = build_outp(inp)
output_indices = build_output_indices_array(outp, max_factor)
puts "The are #{output_indices.length} non-zero elements in the Fourier transform."
while !output_indices.empty?
rand_ind = rand(output_indices.length)
c = output_indices.delete_at(rand_ind)
puts "The program randomly picked #{c}."
cf_array =  fraction_to_cf(c, q)
r_min = find_r_min(cf_array, n)
puts "From the continued fraction expansion the minimum r is #{r_min}."
found, r = find_r(r_min, a, n)
if found
puts "After testing integer multiples of #{r_min}, we find r = #{r}."
#                       Ruby's power modulus is either not working or not implemented!
#                       sr = a**(r / 2) % n
sr = pow_mod(a, r / 2, n)
printf "%d is a non-trivial square root of 1 (mod %d).\n", sr, n
p, q = gcd(sr + 1, n), gcd(sr - 1, n)
printf "The factors of %d are %d and %d.\n", n, p, q
all_done = true
break
else
puts 'Could not find the period from the chosen random index; trying another one.'
end
end
"Exhausted all elements in the QFT. Run the program again, please." if !all_done
end
break
end
break
end