function x=grover(n,p)
N=2^n; % The total number of items to search
assert(N>1);
assert(N<=1024);
assert(p>0);
assert(p<=N);
% Set up oracle
Uf=eye(N*2);
Uf(2*p-1,2*p-1)=0;
Uf(2*p,2*p)=0;
Uf(2*p,2*p-1)=1;
Uf(2*p-1,2*p)=1;
X=[0 1;1 0];
InX=kron(eye(N),X);
H=1;
for j=1:log2(N)+1
    H=kron(H,[1 1;1 -1]/sqrt(2));
end
W=repmat(2/N,N)-eye(N);
WxI=kron(W,eye(2));
% Set up a superposition of all states tensored with (|0> - |1>)/sqrt(2)
y=zeros(N*2,1);
y(1)=1;
y=InX*y;
y=H*y;
% Calculate optimal number of iterations
k=(pi/4)*sqrt(N)
for j=1:k
    % Apply to oracle to sign flip the desired item
    y=Uf*y;
    % Invert all states about the average
    y=WxI*y;
end
% Extract first n qubits (the oracle qubit never changes)
x=zeros(N,1);
for j=1:N
    x(j)=y(2*j-1)*sqrt(2);
end
bar(x);

