CS 596 Quantum Computing
The Feynman Computer Simulator
RSA Algorithm Example
Shor's Algorithm Simulation
%main program
%n=#of qbits
n=4
I=sqrt(-1);
BV=Basis(n);
%input the Hamiltonian
%A=[0 1;1 0]^(1/2);
A=1/sqrt(2)*[1 1;1 -1];
%H1=kron(eye(2^(n-1)),A);
H1=kron(kron(eye(4),A),A);
%Feynman's Hamiltonian
c=[0 0;1 0]; a=[0 1;0 0]; E=eye(2);
c1=kron(kron(kron(E,c),E),E);
a0=kron(kron(kron(a,E),E),E);
c2=kron(kron(kron(E,E),c),E);
a1=kron(kron(kron(E,a),E),E);
%H=c1*a0*H1+c2*a1*H1;
H=c1*a0*H1;
H=H+H';
%calculate the Unitary matrix
deltat=1;
U=expm(-I*deltat*H);%U=U^(1/2);
%Initial Vector
p=zeros(2^n,1);
p(9)=1;
%Evolve the system
for
i=1:5Up=U*p;
T=[real(Up),imag(Up), BV, Up.*conj(Up)];
fprintf(1,
' time = %3.0f\n',i)fprintf(1,
'%10.4f + %7.4f I ket[%1.0f %1.0f %1.0f %1.0f] %6.4f\n',T.') %read the cursor Up=readq(1,Up);Up=readq(2,Up);
%Up=readq(3,Up); %Up=readq(4,Up); T=[real(Up),imag(Up), BV, Up.*conj(Up)];fprintf(1,
' Read the Cursor time = %3.0f\n',i)fprintf(1,
'%10.4f + %7.4f I ket[%1.0f %1.0f %1.0f %1.0f] %6.4f\n',T.')p=Up;
end
function
BV=Basis(n)%generate the basis vectors
m=2^n;
BV=zeros(m,n); P=[0;1];
BV(:,1)=kron(P,ones(2^(n-1),1));
for
i=2:(n-1)BV(:,i)=kron(kron(ones(2^(i-1),1),P),ones(2^(n-i),1));
end
BV(:,n)=kron(ones(2^(n-1),1),P);
function
W=readq(i,Up)A=[1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0
1 1 1 1 0 0 0 0 1 1 1 1 0 0 0 0
1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0
1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0];
M=[A(i,:);ones(1,16)-A(i,:)];
pr=M*(Up.*conj(Up));
if
rand > pr(1)k=2;
else k=1;end
W=(1/sqrt(pr(k)))*M(k,:)'.*Up;
clear
%choose two large primes and calculate product
p=811
q=613
n=p*q
pause
%choose an integer d coprime to (p-1)(q-1)
d=3211
pmqm=(p-1)*(q-1)
gcd(d,pmqm)
pause
%Compute e the modular inverse
fp=factor(p-1)
fq=factor(q-1)
pause
phi=fix(pmqm*(1-1/2)*(1-1/3)*(1-1/5)*(1-1/17))
pause
e=modp(3211,phi-1,pmqm)
pause
%Broadcast e and n
% M is the integer message
M=123456
pause
%Encrypt the message
E=modp(M,e,n)
pause
%Send and then receiver decodes
Mdc=modp(E,d,n)
function
z=modp(s,p,a)%Calculates s^p mod a for p < 2^20 ~ 10^6
c=zeros(1,20); c(1)=s;
for
i=2:20c(i)=mod(c(i-1)*c(i-1),a);
end
b=zeros(1,20);
for
i=19:-1:0r=p-2^i;
if r>=0b(i+1)=1; p=r;
endend
z=1;
for
i=1:20 if b(i)==1z=mod(z*c(i),a);
endend
% Shor's Algorithm on a conventional computer
%N = number to be factored
N=35
%choose a number relatively prime to N
x=23
gcd(x,N)
%choose a number ~ N^2
q=1024
a=1:q; fna=zeros(1,q);
fna(1)=x;
for
i=2:(q)fna(i)=mod(x*fna(i-1),N);
end
plot(a(1:25),fna(1:25))
grid on
pause
ffna=abs(ifft(fna));
plot(a,ffna)
axis([0,q,0,1])
ffna(1:500)
f=input(
'Enter Frequency ')[Q,NU,D]=cfrac(f,1024);
Q,[NU;D]
[Q,NU,D]=cfrac(f+1,1024);
Q,[NU;D]
[Q,NU,D]=cfrac(f+2,1024);
Q,[NU;D]
[Q,NU,D]=cfrac(f+3,1024);
Q,[NU;D]
r=input(
'Enter Period ')gcd(x^(r/2)-1,N)
gcd(x^(r/2)+1,N)
Continued Fraction Expansion Function
function
[Q,N,D] = cfrac(p,q)%Calculates the Continued Fraction Expansion of p/q assuming that p<q
r=1; i=0;
while
r>0i=i+1;
Q(i)=floor(q/p);
r=q-Q(i)*p;
q=p; p=r;
end
n=length(Q);
for
i=n:-1:2a=Q(i); b=1;
for j=i-1:-1:1c=a;
a=Q(j)*a+b;
b=c;
end N(i)=b; D(i)=a;end
N(1)=1; D(1)=Q(1);
end