function ym=shor(M,a)
% This function performs the quantum Fourier transform of Shor's algorithm
%  by simulating the changes that occur in the quantum state as the
%  algorithm is executed
% ym=shor(M,a)
% Inputs
%   M = The number we are factoring
%   a = Arbitrary number such that 0 < a < M
% Output
%   ym = The probabilities of measuring each outcome of the QFT
m=ceil(log2(M*M)); % M^2 <= 2^m < 2*M^2
L=ceil(log2(M));
em=2^m;
eL=2^L;
emL=em*eL; % Total number of quantum states
% State space ranges from |0,0> to |2^m-1,2^L-1>
y=zeros(1,emL);
yt=zeros(1,emL);
% Create superposition of states |x,a^x mod M> for all x where 0 <= x < 2^m
% Since a^0 = 1, the first state in the superposition is |0,1>
axM = 1;
y(1+axM)=1/sqrt(em); % |0,1> is 2nd element; |0,0> is 1st (MATLAB uses 1-based indices)
for x=1:em-1
    axM=mod(axM*a, M);
    y(1+axM+x*eL)=1/sqrt(em);
end
% Perform quantum Fourier transform on the first m qubits
% The last L qubits are unchanged
% Therefore, the actual state transition matrix is the tensor product
%  of the QFT matrix and a (2^L)x(2^L) identity matrix
% Because the full matrix is very large, we use a for loop to do the
%  matrix multiplication instead of building the full matrix
for u=1:eL
    yt(u:eL:emL)=fft(y(u:eL:emL))/sqrt(em);
end
% Plot and return the probabilities of measuring each possible state of the
% first m qubits; the peaks are used to determine the period of a^x mod M
ym=zeros(1,em);
for x=0:em-1
    ym(1+x)=norm(yt(1+eL*x:eL*(x+1)))^2;
end
plot(0:em-1,ym);
return;
