3
* This file implements the shor algorithm.
6
* 0.1 initial version, implementation of the classical version of
7
* the factoring algorithm, and some matrix operations using
12
/* A short description due to Matthew Hayward :
15
Steps to Shor's Algorithm
17
Shor's algorithm for factoring a given integer n can be broken into
20
Step 1) Determine if the number n is a prime, a even number, or an
21
integer power of a prime number. If it is
22
we will not use Shor's algorithm. There are efficient classical methods
23
for determining if a integer n belongs
24
to one of the above groups, and providing factors for it if it is. This
25
step would be performed on a classical computer.
27
Step 2) Pick a integer q that is the power of 2 such that "n^2 <= q < 2*n^2".
28
This step would be done on a classical computer.
30
Step 3) Pick a random integer x that is coprime to n. When two numbers
31
are coprime it means that their greatest common divisor is 1. There
32
are efficient classical methods for picking such an x. This step would be
33
done on a classical computer.
35
Step 4) Create a quantum register and partition it into two sets,
36
register 1 and register 2. Thus the state of our quantum computer can
37
be given by: |reg1, reg2>. Register 1 must have enough qubits to represent
38
integers as large as q - 1. Register 2 must have enough qubits to represent
39
integers as large as n - 1. The calculations for how many qubits are
40
needed would be done on a classical computer.
42
Step 5) Load register 1 with an equally weighted superposition of all
43
integers from 0 to q - 1. Load register 2 with all zeros. This
44
operation would be performed by our quantum computer. The total state
45
of the quantum memory register at this point is:
47
"1/Sqrt(q) Sum(a,0,q-1) |a,0>"
49
Step 6) Now apply the transformation "x^a mod n" to for each number stored in
50
register 1 and store the result in register 2. Due to quantum parallelism
51
this will take only one step, as the quantum computer will
52
only calculate "x^|a> mod n", where |a> is the superposition of states created in
53
step 5. This step is performed on the quantum computer. The state of the
54
quantum memory register at this point is:
56
"1/Sqrt(q) Sum(a,0,q-1) |a,x^a mod n>"
58
Step 7) Measure the second register, and observe some value k. This has
59
the side effect of collapsing register one into a equal superposition of
60
each value a between 0 and q-1 such that
64
This operation is performed by the quantum computer. The state of the
65
quantum memory register after this step is:
67
"1/Sqrt(||A||) Sum(a' elements of A) |a',k>"
69
Where A is the set of a's such that "x^a mod n = k", and ||A|| is the
70
number of elements in that set.
72
Step 8) Now compute the discrete Fourier transform on register one. The
73
discrete Fourier transform when applied to a state |a>changes it in the
76
"|a> = 1/Sqrt(q) Sum(c=0,q-1) |c> *Exp(2*Pi*I*a*c/q)"
78
This step is performed by the quantum computer in one step through quantum
79
parallelism. After the discrete Fourier transform our register is in
82
"1/Sqrt(||A||) Sum(a' element of A) 1/Sqrt(q) Sum(c=0,q-1) |c,k> *Exp(2*Pi*I*a'*c/q)"
84
Step 9) Measure the state of register one, call this value m, this integer
85
m has a very high probability of being a multiple of q/r, where r is the
86
desired period. This step is performed by the quantum computer.
88
Step 10) Take the value m, and on a classical computer do some post
89
processing which calculates r based on knowledge of m and q. There
90
are many ways to do this post processing. This post processing is done on a
93
Step 11) Once you have attained r, a factor of n can be determined by
94
taking "Gcd(x^(r/2)+1,n)" and "Gcd(x^(r/2)-1,n)"
95
. If you have found a factor of n, then stop, if not
96
go to step 4. This final step is done on a classical computer.
98
Step 11 contains a provision for what to do if Shor's algorithm failed
99
to produce factors of n. There are a few reasons why Shor's algorithm
100
can fail, for example the discrete Fourier transform could be measured
101
to be 0 in step 9, making the post processing in step 10 impossible.
102
The algorithm will sometimes find factors of 1 and n, which is not
103
useful either. For these reasons step 11 must be able to jump back
104
to step four to start over. (Williams, Clearwater)
113
/* ClassicPeriod : this is the computation that the Shor
114
algorithm is supposed to perform more efficiently.
115
It tries to find the period of n^(p+1) mod m. That is, it
116
will return the p for which n mod m = n^p mod m
117
Example: ClassicPeriod(4,15) -> 2
122
Local(i,first,period);
127
While(Mod(i,m)!=first)
135
/* ClassicFactors does the same as Factors, but using a period
136
Example: ClassicFactors(4,15,ClassicPeriod(4,15)) -> {5,3}
138
ClassicFactors(n,m,period):=
142
{Gcd(res+1,m),Gcd(res-1,m)};
147
/* Rot: tensor defining a rotation in a plane spanned by two
148
* arbitrary base vectors. Rot(n,m,i,j,a) gives the matrix element
149
* (i,j) of a rotation matrix that rotates the space in the plane (n,m)
153
5 # Rot(ii_IsInteger,_ii,jj_IsInteger,_jj,_angle) <-- 1;
154
6 # Rot(ii_IsInteger,_ii,jj_IsInteger,kk_IsInteger,_angle) <-- 0;
155
10 # Rot(_ii,_jj,_ii,_ii,_angle) <-- Cos(angle);
156
10 # Rot(_ii,_jj,_ii,_jj,_angle) <-- -Sin(angle);
157
10 # Rot(_ii,_jj,_jj,_ii,_angle) <-- Sin(angle);
158
10 # Rot(_ii,_jj,_jj,_jj,_angle) <-- Cos(angle);
159
20 # Rot(ii_IsInteger,jj_IsInteger,kk_IsInteger,_kk,_angle)_
160
(kk!= ii And kk != jj) <-- 1;
162
30 # Rot(ii_IsInteger,jj_IsInteger,kk_IsInteger,ll_IsInteger,_angle)_
163
(kk != ii And kk != jj) <-- 0;
164
30 # Rot(ii_IsInteger,jj_IsInteger,kk_IsInteger,ll_IsInteger,_angle)_
165
(ll != ii And ll != jj) <-- 0;
168
/* Two examples of use of Rot:
169
* - f1(i,a) returns the i-th vector element
170
* after rotating a vector X in an aight-dimensional space in the
173
* - f2(i,a1,a2) does the same, rotating over (5,3) and then (2,5)
175
f1(i,a):=Eval(TExplicitSum(8)TSum({j})Rot(2,5,i,j,a)*X(j));
176
f2(i,a1,a2):=Eval(TExplicitSum(16)TSum({j,k})Rot(2,5,i,j,a1)*Rot(5,3,j,k,a2)*X(k));
180
/* Mir: mirror a vector. Mir(i,j,k) is the matrix that converts a vector
181
* into one with the i-th component changed from X(i) to -X(i).
184
5 # Mir(_ii,_ii,_ii) <-- -1;
185
6 # Mir(_ii,_jj,_jj) <-- 1;
186
7 # Mir(_ii,jj_IsInteger,kk_IsInteger) <-- 0;
188
/* Example, mirror the i-th component, and show the j-th component
189
f3(i,j):=Eval(TExplicitSum(8)TSum({k})Mir(i,j,k)*X(k));
192
/* InsertBit inserts a bit into a number num, at position index. the
193
* first bit has index 0.
196
InsertBit(_index,_bit,_num) <--
198
Local(low,high,lowmask);
199
lowmask := (1<<index)-1;
200
low := BitAnd(num,lowmask);
202
BitOr(BitOr(low,high),bit*(1<<index));
205
/* Had is an implementation of the Hadamard gate. It rotates the bit n
206
* into a superposition of states (0,1) by rotating around angle Pi/4.
208
Had(_n,_Nbits,_ii,_jj) <--
210
Local(list,result,i1,i2,kk,in,kkk);
212
list:= (0 .. (2^(Nbits-1) - 1));
213
kkk:=MakeVector(kk,2^(Nbits-1)-1);
214
kk:=Concat({ii},kkk,{jj});
218
i1:=InsertBit(n,0,item)+1;
219
i2:=InsertBit(n,1,item)+1;
220
result:=result*Rot(i1,i2,kk[in],kk[in+1],Pi/4);
223
TExplicitSum(2^Nbits)(TSum(kkk)result);
226
/* Example: apply one Hadamard gate to a vector space on qubit 1.
227
f4(i):=Eval(TExplicitSum(8)TSum({j})Had(1,3,i,j)*X(j));
230
/* Applying the Hadamard gate twice: a quantum not gate on qubit 1.
231
f5(i):=Eval(TExplicitSum(4)TSum({j,k})Had(1,2,i,j)*Had(1,2,j,k)*X(k));
234
/* This example prepares the state in its ground state:
235
* XX(1)=1 and XX(n!=1)=0.
236
* Applying the Hadamard gate to all the bits. This prepares the system in
237
* a state where each combination is equally probable. Try it out with f6(n),
238
* which should return 1/2 for any integer (1 .. 4)
241
f6(i):=Eval(TExplicitSum(4)TSum({j,k})Had(0,2,i,j)*Had(1,2,j,k)*XX(k));
247
/* CN: Controlled Not gate. If the control bit is set, the target
250
CN(_target,_control,_Nbits,_ii,_jj) <--
252
Local(list,result,i1,i2,kk,in,kkk,cntl);
254
list:= (0 .. (2^(Nbits-1) - 1));
255
kkk:=MakeVector(kk,2^(Nbits-2)-1);
256
kk:=Concat({ii},kkk,{jj});
262
i1:=InsertBit(target,0,item)+1;
263
i2:=InsertBit(target,1,item)+1;
264
If (BitAnd(cntl,i1-1)!= 0,
266
result:=result*Rot(i1,i2,kk[in],kk[in+1],Pi/2);
271
TExplicitSum(2^Nbits)(TSum(kkk)result);
274
/* Example using CN: 3-qubit space, with bit 0 the control bit
275
* and bit 1 th target bit. This f7(1)=X(1), f7(2)=-X(4), f7(3)=X(3)
277
f7(i):=Eval(TExplicitSum(8)TSum({j})CN(1,0,3,i,j)*X(j));
284
10 # ShorFactor(_M,_x)_(Mod(M , 2) = 0) <--
286
Echo({"Number is even. Factors found : ",M," = 2*",M>>1});
289
20 # ShorFactor(M_IsPrime,_x) <--
291
Echo({"Numer is prime. It cannot be factored."});
295
25 # ShorFactor(M_IsPrimePower,_x) <--
297
Echo({M," is a prime power (not handled yet)."});
301
30 # ShorFactor(_M,_x)_
304
factor := MathGcd(M,x);
305
factor != 1 And factor != M;
309
factor := MathGcd(M,x);
310
Echo({"Factor found since GCD(",M,",",x,")=",factor});
313
40 # ShorFactor(_M,_x) <--
315
Local(q,inbits,outbits);
316
q:= (1<<CountBits(M^2));
317
inbits:=CountBits(q-1);
318
outbits:=CountBits(M-1);
321
shorinitstate(i):=Eval(InState(inbits,i));
323
shormodexpstate(i):=Eval(ModExpState(inbits,i,x,M));
324
nr:= (1<<(inbits+outbits));
327
If(shormodexpstate(i) != 0, ShowModExp(inbits,i,x,M));
331
/* TODO old way of calculating */
332
ClassicFactors(x,M,ClassicPeriod(x,M));
336
10 # InState(_inbits,i_IsInteger)_((i>>inbits) = 0) <-- 1/Sqrt(1<<inbits);
337
20 # InState(_inbits,i_IsInteger) <-- 0;
340
ModExpState(_inbits,i_IsInteger,_x,_M) <--
344
in := BitAnd(i,(1<<inbits)-1);
345
out:= (i-in)>>inbits;
346
/* Echo({out,aMod(x^in,M)}); */
347
If(out = Mod(x^in,M), 1/Sqrt(1<<inbits), 0);
350
ShowModExp(_inbits,i_IsInteger,_x,_M) <--
354
in := BitAnd(i,(1<<inbits)-1);
355
out:= (i-in)>>inbits;
added
removed
inprogress/shor2
