gusucode.com > qit_matlab_0.10.0工具箱源码程序 > qit/examples/shor_factorization.m

    function [p] = shor_factorization(N, cheat)
% SHOR_FACTORIZATION  Shor's factorization algorithm demo.
%  p = shor_factorization(N)       % simulate the full algorithm
%  p = shor_factorization(N, true) % cheat, avoid the quantum part
%
%  Simulates Shor's factorization algorithm, tries to factorize the integer N.
%
%  NOTE: This is a very computationally intensive quantum algorithm
%  to simulate classically, and probably will not run for any
%  nontrivial value of N (unless you choose to cheat, in which case
%  instead of simulating the quantum part we use a more efficient
%  classical algorithm for the order-finding).

%! P.W. Shor, "Algorithms For Quantum Computation: Discrete Logs and Factoring", Proc. 35th Symp. on the Foundations of Comp. Sci., 124 (1994).
%! M.A. Nielsen, I.L. Chuang, "Quantum Computation and Quantum Information" (2000), chapter 5.3.
% Ville Bergholm 2010


fprintf('\n\n=== Shor''s factorization algorithm ===\n\n')

global qit;

if (nargin < 2)
  cheat = false;
end
if (cheat)
  fprintf('(cheating)\n\n');
end

% number of bits needed to represent mod N arithmetic:
m = ceil(log2(N));
fprintf('Trying to factor N = %d (%d bits).\n', N, m);

% maximum allowed failure probability for the quantum order-finding part
epsilon = 0.25;

% number of index qubits required for the phase estimation
t = 2*m +1 +ceil(log2(2+1/(2*epsilon)));

fprintf('The quantum order-finding subroutine will need %d+%d qubits.\n\n', t, m);


% classical reduction of factoring to order-finding
while (true)
  a = 1+randi(N-2); % random integer, 1 < a < N
  fprintf('Random integer: a = %d\n', a)
  
  p = gcd(a, N);
  if (p ~= 1)
    % a and N have a nontrivial common factor p.
    % This becomes extremely unlikely as N grows.
    fprintf('Lucky guess, we found a common factor!\n')
    break
  end

  fprintf('Trying to find the period of f(x) = a^x mod N');

  if (cheat)
    % classical cheating shortcut
    r = find_order_cheat(a, N);
  else
    while (true)
      fprintf('.');
      % ==== quantum part of the algorithm ====
      [s1, r1] = find_order(a, N, t, m);
      [s2, r2] = find_order(a, N, t, m);
      % =============  ends here  =============

      if (gcd(s1, s2) == 1)
        % no common factors
        r = lcm(r1, r2);
        break;
      end
    end
  end
  fprintf('\n  =>  r = %d\n', r)

  % if r is odd, try again
  if (gcd(r, 2) == 2)
    % r is even

    x = mod_pow(a, r/2, N);
    if (mod(x, N) ~= N-1)
      % factor found
      p = gcd(x-1, N); % ok?
      if (p == 1)
        p = gcd(x+1, N); % no, try this
      end
      break
    else
      fprintf('a^(r/2) = -1 (mod N), try again...\n\n');
    end
  else
    fprintf('r is odd, try again...\n\n');
  end
end

fprintf('\nFactor found: %d\n', p)
end


function [s, r] = find_order(a, N, t, m)
% Quantum order-finding subroutine.
% Finds the period of the function f(x) = a^x mod N.

T = 2^t; % index register dimension
M = 2^m; % state register dimension

% applying f(x) is equivalent to the sequence of controlled modular multiplications in phase estimation
U = gate.mod_mul(a, M, N);
U = U.data; % FIXME phase_estimation cannot handle lmaps right now

% state register initialized in the state |1>
st = state(1, M);

% run the phase estimation algorithm
reg = phase_estimation(t, U, st, true); % use implicit measurement to save memory

% measure index register
[dummy, res] = measure(reg, 1);

num = res-1;

% another classical part
r = find_denominator(num, T, T+1);
s = num*r/T;
end


function d_1 = find_denominator(x, y, max_den)
% Finds the denominator q for p/q \approx x/y such that q < max_den
% using a continued fraction representation for x.
%
% We use floor and mod here, which could be both implemented using
% the classical Euclidean algorithm.

d_2 = 1;
d_1 = 0;

while (1)
  a = floor(x/y); % integer part == a_n
  temp = a*d_1 +d_2; % n:th convergent denumerator d_n = a_n*d_{n-1} +d_{n-2}

  if (temp >= max_den)
    return % stop
  end
  d_2 = d_1;
  d_1 = temp;
 
  temp = mod(x,y); % x - a*y; % subtract integer part

  if (temp == 0)
  %if (temp/y < 1 / (2*max_den^2))
    return
  end

  % invert the remainder (swap numerator and denominator)
  x = y;
  y = temp;
end
end


function r = find_order_cheat(a, N)
% Classical order-finding algorithm.

for r = 1:N
  if (mod_pow(a, r, N) == 1)
    return
  end
end
end


function res = mod_pow(a, x, N)
% Computes mod(a^x, N) using repeated squaring mod N.
% x must be a positive integer.

b = fliplr(dec2bin(x) - '0'); % exponent in little-endian binary
res = 1;
for k = 1:length(b)
  if (b(k))
    res = mod(res*a, N);
  end
  a = mod(a*a, N); % square it
end
end