#!/usr/bin/perl -w
die 'usage: ./shor.pl [number to factor]' unless @ARGV;
use strict;
use warnings;
use Quantum::Entanglement qw(:DEFAULT :complex :QFT);
$Quantum::Entanglement::destroy = 0;
my $num = $ARGV[0];
# do some early die'ing
die "$num is a multiple of two, here I am, brain the size..." unless $num %2;
die "$num is a non-integer, I only have whole numbers of fingers"
unless $num == int($num);
die "$num is less than 15" unless $num >= 15;
print "Performing initial classical steps:\n";
# work out q value
my $q_power = int(2* log($num) / log(2)) +1;
my $q = 2 ** $q_power;
# pick some x so that x is coprime to n.
my $x;
do {
$x = int(rand $num) + 1;
} until ($num % $x != 0 and $x > 2); #ok.. so this misses the point slightly
print "Using q:$q, x:$x\nStarting quantum steps\n";
# fill up a register with integers from 1..q
my $prob = 1/sqrt($q);
my $register1 = entangle(map {$prob, $_} (0..$q-1));
# apply transformation F = x**|a> mod n, store in register 2
# (need to do a p_func to avoid overflow while **)
sub power_mod {
my ($state, $x1, $num1) = @_;
my $rt = 1;
return 1 if $state == 0;
return 1 if $state == 1;
for (1..$state) {
$rt = ($rt * $x1) % $num1;
}
return $rt;
}
print "Performing F = x**|a> mod n\n";
my $register2 = p_func(\&power_mod, $register1, $x, $num);
# We now observe $register2, thus partially collapsing reg1
my $k = "$register2";
print "\$register2 collapsed to $k\n";
print "Finding period of F (this is where you wish for a QCD)\n";
# take a ft of the amplitudes of reg1, placing result in reg3
my $register3 = QFT($register1);
my $lqonr = "$register3"; # observe, this must be multiple of q/r
if ($lqonr == 0) {
print "Got period of '0', halting\n"; exit(0);
}
my $period = int($q / $lqonr + 0.5); # rounding
print "Period of F = x**|a> mod n is $period\n";
# now given the period, we need to work out the factor of n
# work out the two thingies:
if ($period % 2 != 0) {
print "$period is not an even number, doubling to";
$period *=2;
print " $period\n";
}
my $one = $x**($period/2) -1;
my $two = $x**($period/2) +1;
# one and two must have a gcd in common with n, which we now find...
print "$one * $two and $num might share a gcd (classical step)\n";
my ($max1, $max2) = (1,1);
for (2..$num) {
last if $_ > $num;
unless (($num % $_) || ($one % $_)) {
$max1 = $_;
}
unless (($num % $_) || ($two % $_)) {
$max2 = $_;
}
}
print "$max1, $max2 could be factors of $num\n";
__END__;
=head1 NAME
shor - A short demonstration of Quantum::Entanglement
=head1 SYNOPSIS
./shor.pl [number to factor (>14)]
=head1 DESCRIPTION
This program implements Shor's famous algorithm for factoring numbers. A
brief overview of the algorithm is given below.
=head2 The important maths
Given a number B which we are trying to factor, and some other number
which we have guessed, B, we can say that:
x**0 % n == 1 (as x**0 = 1, 1 % n =1)
There will also be some other number, B such that
x**r % n == 1
or, more specifically,
x**(kr) % n ==1
in other words, the function
F(a) = x**a % n
is periodic with period B.
Now, starting from
x**r = 1 % n
x**(2*r/2) = 1 % n
(x**(r/2))**2 - 1 = 0 % n
and, if r is an even number,
(x**(r/2) - 1)*(x**(r/2) + 1) = 0 mod n
or in nice short words, the term on the left is an integer multiple of B.
So long as x**(r/2) != +-1, at least one of the two brackets on the left
must share a factor with B.
Shor's alorithm provides a way to find the periodicity of the function F
and thus a way to calculate two numbers which share a factor with n, it
is then easy to use a classical computer to find the GCD and thus a
factor of B.
=head1 The steps of the algorithm
=head2 1. Remove early trivial cases
We have efficient classical methods for finding that 2 is a factor of 26,
so we do not need to use this method for this.
=head2 2. Pick an integer
Chose a number B so that C, this is done on a
classical computer. (This is the size we will use for our quantum register.)
=head2 3. Select at random a number coprime to n
Think of some number less than B so that B and B do not share
a common factor (if they do, we already know the answer...).
=head2 4. Fill a quantum register with integers from 0..q-1
This is where we create our first entangled variable, and is the first
non-classical step in this algorithm.
=head2 5. Calculate F, store in a second register
We now calculate C< F(a) = x**a % n> where a represents the superposition
of states in our first register, we store the result of this in our
second register.
=head2 6. Look at register2
We now look at the value of register two and get some value B, this forces
register1 into
a state which can only collapse into values satisfying the equation
x**a % n = k
The probability amplitudes for the remaining states are now all equal to zero,
note that we have not yet looked directly at register1.
=head2 7. Find period of register1
We now apply a fourier transform to the amplitudes of the states in
register1, storing the result as the probability amplitudes for a new
state with the values of register1. This causes there to be a high
probability that the register will collapse to a value which is some
multiple of C.
=head2 8. Observe register1
We now observe register1, and use the result to calculate a likely value
for B. From this we can easily calculate two numbers, one of which
will have a factor in common with n, by applying an efficient classical
algoirthm for finding the greatest common denominator, we will be able
to find a value which could be a factor of B.
=head1 Things to remember
This algorithm does not claim to produce a factor of our number the first
time that it is run, there are various conditions which will cause it
to halt mid-way, for instance, the FT step can give a result of 0 which
is clearly useless. The algorithm is better than any known classical one
because the expectation value of the time required to get a correct answer
is still O(n).
This also cannot factor a number which is prime (it being, as it were, prime)
and also cannot factor something which is a prime power (25, say).
=head1 COPYRIGHT
This code is copyright (c) Alex Gough (alex@rcon.org )2001. This is
free software, you may use, modify and redistribute it under the same
terms as Perl itself.
=head1 BUGS
This is slow, being run on classical computers, ah well.
=cut