Algorithms Used to Factor in Qwerty #

from qwerty import *

def qpe(precision, prep_eigvec, op,
        n_shots):
  @qpu[M,T](prep_eigvec, op)
  def kernel(
      prep_eigvec: qfunc[0,M],
      op: rev_qfunc[M][[...]]) \
      -> bit[T]:
    return '+'[T] + prep_eigvec() \
           | (std[T-1-j]+'1'+std[j]
              & op[[j]]
              for j in range(T)) \
           | fourier[T].measure \
             + discard[M]

  k_inst = kernel[[precision]]
  for meas in k_inst(shots=n_shots):
    yield meas.as_bin_frac()

from qpe import qpe
from qwerty import *
import sys

phi = float(sys.argv[1])
precision = int(sys.argv[2])

@qpu
def prep1() -> qubit:
    return '1'

@qpu[J](phi)
@reversible
def rot(phi: angle,
        q: qubit) -> qubit:
    return q | std >> \
      {'0', phase(2*pi*phi*2**J)*'1'}

print('Expected:', phi)
phi_got, = qpe(precision, prep1, rot,
               n_shots=1)
print('Actual:', float(phi_got))

import math
from qwerty import *
from qpe import qpe

def order_finding(epsilon, x, modN):
  L = math.ceil(math.log2(modN))
  t = 2*L + 1 + math.ceil(
        math.log2(2+1/(2*epsilon)))

  @qpu[M]
  def one() -> qubit[M]:
    return '0'[M-1] + '1'

  @classical[X,N,M,J]
  def xymodN(y: bit[M]) -> bit[M]:
    return X**2**J * y % N

  x_inv = pow(x, -1, modN)
  fwd = xymodN[[x,modN,L,...]]
  rev = xymodN[[x_inv,modN,L,...]]
  mult = fwd.inplace(rev)
  frac1, frac2 = qpe(t, one, mult, 2)

  def denom(frac):
    cf = cfrac.from_fraction(frac)
    for c in reversed(
        cf.convergents()):
      if c.denominator < modN:
        return c.denominator

  return math.lcm(denom(frac1),
                  denom(frac2))

import math
from qwerty import *
from order_finding import order_finding
import sys

def naive_classical(x, modN):
  for r in range(1, modN):
    if x**r % modN == 1:
      return r

err = 0.2
x = int(sys.argv[1])
modN = int(sys.argv[2])
if math.gcd(x, modN) != 1:
  raise ValueError('invalid x, modN')

print('Classical:',
      naive_classical(x, modN))
print('Quantum:',
      order_finding(err, x, modN))

import math
import random
from qwerty import *

from order_finding import \
  order_finding

def shors(epsilon, num):
  if num % 2 == 0:
    return 2

  x = random.randint(2, num-1)
  if (y := math.gcd(x, num)) > 1:
    print('Got lucky! Skipping '
          'quantum subroutine...')
    return y

  r = order_finding(epsilon, x, num)

  if r % 2 == 0 \
      and pow(x, r//2, num) != -1:
    if (gcd := math.gcd(x**(r//2)-1,
                        num)) > 1:
      return gcd
    if (gcd := math.gcd(x**(r//2)+1,
                        num)) > 1:
      return gcd

  raise Exception("Shor's failed")

from qwerty import *
from shors import shors

import sys

err = 0.2
num = int(sys.argv[1])
print('Nontrivial factor of', num,
      'is', shors(err, num))