```# Here is a simplified program for Shor's Factoring Algorithm.
# Consider it useful for instructional purposes only.
# It is very inefficient, I know. The Fourier transform is brute
# force, using neither the FFT nor any optimizations. The "quantum
# entanglement" is of course carried out classically.<p>

#! /usr/bin/ruby

require 'complex'
include Math

# This is Shor's factoring algorithm as explained in the class on Quantum Mechanics and
# Quantum Computation. It does not use the QFT, but simply the Fourier transform matrix
# without any optimizations. Then the program looks for the greatest common divisor among
# the indices of all the non-zero elements in the output. This works fine, because the order
# of the transform is an integer multiple of the order of the function f(x) = a^x (mod N),
# the principle disadvantage being that the transform is slow.

# 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

# Build the input array as it would look after measuring f(x) = a^x (mod N)
# and entanglement.
# Input inp -- initially a blank array; filled with the entangled state on return
# Input rand_ind -- a random index with 0 <= rand_ind < r
# Input r -- the order of a (mod N)
# Input m -- the length of inp
# Return -- none
def build_sft_input(inp, rand_ind, r, m)
m.times { |i| inp[i] = 1.0 if i % r == rand_ind }
end

# The slow Fourier transform
# Input inp -- transform input
# Input outp-- initially a blank array; contains the transform on output
# Return -- none
def sft(inp, outp)
m, sum = inp.length, nil
omg = exp(2 * Complex::I * PI / m)
puts "M = #{m}"
m.times do |i|
sum = Complex(0.0, 0.0)
m.times do |j|
sum += omg**(i * j) * inp[j]
end
outp << sum
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

# 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

# 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

# Collect the (non-zero) qubit indices after the Fourier transform
# Input outp -- the Fourier transform
# Input output_indices -- the qubits in the transform
# Return -- none
def build_output_indices_array(outp, output_indices)
outp.each_index { |i| output_indices << i if outp[i].polar > Epsilon }
end

Uncertainty, Epsilon, inp, outp, output_indices, m_over_r = 0.01, 0.0001, [], [], [], nil
have_p, have_q = false, false

begin
p, q = nil, nil
loop do
if !have_p
STDOUT.write "Enter first prime (<= 47): "
p = gets().to_i()
if p < 2 or p > 47
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 (<= 47): "
q = gets().to_i()
if q < 2 or q > 47
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
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
loop do
STDOUT.write "Enter M / r: "
m_over_r = gets().to_i()
break if m_over_r > 0
puts "#{m_over_r} is invalid"
end
m = m_over_r * r
inp = Array.new(m, 0.0)
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_sft_input(inp, rand_ind, r, m)
sft(inp, outp)
#               outp.each { |el| printf "|(%8.4f, %9.4f)| = %8.4f\n", el.real, el.image, el.polar }
build_output_indices_array(outp, output_indices)
p output_indices
puts "The are #{output_indices.length} elements in the Fourier transform."
gcd_count = nil
loop do
STDOUT.write "Enter the number of those elements you would like to use for the gcd calculation: "
gcd_count = gets().to_i()
break if gcd_count >= 1 and gcd_count <= output_indices.length
puts "#{gcd_count} is invalid"
end
running_gcd = 0
gcd_count.times do
an_index = rand(output_indices.length)
# The 0 entry need not be used in the gcd calculations.
next if an_index == 0
running_gcd = gcd(running_gcd, output_indices[an_index])
output_indices.delete_at(an_index)
end
r = m / running_gcd
puts "The gcd result is #{running_gcd}, resulting in 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
break
end
end
break
end