Multiplication modulo $2^n$

For multiplying a state $|j\rangle$ with a constant $a$ we can follow the scheme to multiply two numbers in binary representation. As an example, multiply $5$ with $3$. $5$ has binary representation $0101$ and $3$ has $0011$. So, the procedure for $5\cdot 3$ is

  00011 * 1 (1*3)
+ 00110 * 0 (0*6)
+ 01100 * 1 (1*12)
+ 11000 * 0 (0*24)
= 01111 (15)

Now, if we have a controlled add operation, we can use the qubits of the first register (in this case representing $5$) as control bits and the other register as the constant to add. The single terms in the sum can be efficiently pre-computed classically as follows

summands <- function(x, n, N) {
  b <- as.integer(intToBits(x))
  ret <- c()
  for(i in c(1:N)) {
    s <- 0
    for(j in c(1:N)) {
      s <- s+as.integer(b[j])*2^(i+j-2)
    }
    ret[i] <- s %% n
  }
  return(ret)
}

Example

x <- 3
summands(3, 2^3, 3)

[1] 3 6 4

Now we need a controlled add operation. Here we can build on our add operation using the qft for which we need a controlled phase shift operation first

cRtheta <- function(bits, theta=0.) {
  cqgate(bits=bits, gate=methods::new("sqgate", bit=as.integer(bits[2]),
                                      M=array(as.complex(c(1, 0, 0, exp(1i*theta))),
                                              dim=c(2,2)), type="Rt"))
}

cadd <- function(c, bits, x, y) {
  n <- length(bits)
  z <- cqft(c=c, x=x, bits=bits)
  for(i in c(1:n)) {
    z  <- cRtheta(bits=c(c, bits[i]), theta = 2*pi*y/2^(n-i+1)) * z
  }
  z <- cqft(c=c, x=z, inverse=TRUE, bits=bits)
  return(invisible(z))
}

Let’s check whether it works as expected (keep in mind that the register is 3 qubit wide plus 1 control qubit, so addition is modulo $2^3=8$):

basis <- c()
for(i in c(0:(2^4-1))) basis[i+1] <- paste0("|", i %/% 2, ">|", i %% 2, ">")
x <- H(1)*qstate(4, basis=basis)
c <- 1
bits <- c(2:4)
z <- cadd(c=c, bits=bits, x=x, y=5)
z

   ( 0.7071068 )    * |0>|0> 
 + ( 0.7071068 )    * |5>|1> 

z <- cadd(c=c, bits=bits, x=z, y=2)
z

   ( 0.7071068 )    * |0>|0> 
 + ( 0.7071068 )    * |7>|1> 

z <- cadd(c=c, bits=bits, x=z, y=8)
z

   ( 0.7071068 )    * |0>|0> 
 + ( 0.7071068 )    * |7>|1> 

Equipped with this functionality, we can finally perform a binary multiplication. Note that we need two registers, the first one to store the initial value and the second one to store the final result of the multiplication

mult <- function(reg1, reg2, x, y, swap=TRUE) {
  stopifnot(length(reg1) == length(reg2))
  n <- length(reg2)
  s <- summands(y, 2^n, n)
  for(i in c(1:n)) {
    x <- cadd(c=reg1[i], bits=reg2, x=x, y=s[i])
  }
  if(swap) {
    for(i in c(1:n)) {
      x <- SWAP(c(reg1[i], reg2[i])) * x
    }
  }
  return(invisible(x))
}

With this we can perform a reversible multiplication, which is why we we introduced the SWAP operations at the end. They interchange the two registers. The result is the following

basis <- c()
for(i in c(0:(2^3-1))) {
  for(j in c(0:(2^3-1))) {
    basis[i*2^3+j + 1] <- paste0("|", i, ">|", j, ">")
  }
}
x <- X(2)*qstate(6, basis=basis)
x

   ( 1 )    * |0>|2> 

reg1 <- c(1:3)
reg2 <- c(4:6)
z <- mult(reg1, reg2, x=x, y=3)
z <- X(5) * z
z

   ( 1 )    * |0>|6> 

z <- mult(reg1, reg2, x=z, y=3)
z

   ( 1 )    * |6>|2> 

Let’s be a bit more precise here for the multiplication of, say $a$ and $b$ in two registers, both $n$ qubits wide. Starting with both registers in state $|0\rangle$, we can first bring the result register to state $|1\rangle$ using a NOT operation, i.e. $|0\rangle|1\rangle$. Now, we multiply by $a$, which leaves us with state $$|a \mod 2^n \rangle|1\rangle$$ which we can apply the NOT gate again to reset the result register to state $|0\rangle$ and then we swap the two registers to arrive at $$|0\rangle|a \mod 2^n\rangle\,.$$ Now we multiply by $b$ getting us to state $$|a\times b \mod 2^n\rangle|a \mod 2^n\rangle\,.$$ Multiply the result register with the inverse of $a \mod n$ and apply the NOT gate getting us to $$|0\rangle|a\times b \mod 2^n\rangle\,.$$ The inverse modulo $2^n$ can be computed efficiently in a classical way by the extended Euclidean algorithm

eEa <- function(a, b) {
  if(a == 0) return(c(b, 0, 1))
  res <- eEa(b %% a, a)
  return(c(res[1], res[3] - (b %/% a) * res[2], res[2]))
}

moduloinverse <- function(a, n) {
  res <- eEa(a=a, b=n)
  if(res[1] != 1) stop("inverse does not exist!")
  return(res[2] %% n)
}

If $a$ and $2^n$ are not coprime, the inverse does not exist. However, for the application we have in mind this is not an issue.

Working modulo $N$

So far we have assumed that we work modulo $2^n$, where $n$ was dictated by the number of qubits. However, this is not the realistic case. We have to write an adder modulo $N<2^n$, i.e. $|x\rangle\to|x + y\mod N\rangle$. We can implement this by subtracting $N < 2^n$ whenever needed. We will follow the convention that if $x\geq N$ the operation will be $|x\rangle\to|x\rangle$. Moreover, we will assume $x,y\geq 0$.

To find out, when this subtraction is needed, is a bit tricky. We want to add $y$ to $x$. To decide beforehand, whether or not we have to subtract $N$, we have to check whether $x < N-y$. If not, we have to subtract $N$. If we subtract $N-y$ from this state $|x\rangle$, the most significant qubit indicates whether there occurred an overflow. Using a CNOT gate we can store this info in one ancilla bit $c_1$. Then we add $N-y$ again to retain the original state. Such an operation can be implemented as follows in a controlled manner (control bit $c$, bits the bits in state $|x\rangle$ where x is stored, $a$ an ancilla bit and $y$ the value to compare with.

cis.less <- function(c, bits, x, c1, a, y) {
  ## add ancilla bit as most significant bit to bits
  b <- c(bits, a)
  n <- length(b)
  ## cadd works modulo 2^n
  z <- cadd(c=c, bits=b, x=x, y=2^n-y)
  ## 'copy' overflow bit
  z <- CNOT(c(a, c1)) * z
  ## add back, resetting ancilla a to |0>
  z <- cadd(c=c, bits=b, x=z, y=y)
  return(z)
}

This routine will set the qubit $|c_1\rangle$ to 1 if $|x_\mathrm{bits}\rangle$ is smaller than $y$ and leave it at zero otherwise. It uses $|a\rangle$ as ancilla bit and $|c\rangle$ as control bit. Here an example

basis <- c()
for(i in c(0:(2^6-1))) {
  basis[i + 1] <-
    paste0("|", i %/% 8 ,">|a=",
    (i %/% 4) %% 2, ">|c1=", (i%/%2) %% 2,
    ">|c=", i%%2, ">")
}

x <- H(1)*qstate(6, basis=basis)
z <- cadd(c=1, bits=c(4,5,6), x=x, y=5)
z

   ( 0.7071068 )    * |0>|a=0>|c1=0>|c=0> 
 + ( 0.7071068 )    * |5>|a=0>|c1=0>|c=1> 

## 5 < 7 -> c1 = 1
v <- cis.less(c=1, bits=c(4,5,6), x=z, c1=2, a=3, y=7)
v

   ( 0.7071068 )    * |0>|a=0>|c1=0>|c=0> 
 + ( 0.7071068 )    * |5>|a=0>|c1=1>|c=1> 

## 5 > 3 -> c1 = 0
w <- cis.less(c=1, bits=c(4,5,6), x=z, c1=2, a=3, y=3)
w

   ( 0.7071068 )    * |0>|a=0>|c1=0>|c=0> 
 + ( 0.7071068 )    * |5>|a=0>|c1=0>|c=1> 

## 5 < 9 -> c1 = 1
w <- cis.less(c=1, bits=c(4,5,6), x=z, c1=2, a=3, y=9)
w

   ( 0.7071068 )    * |0>|a=0>|c1=0>|c=0> 
 + ( 0.7071068 )    * |5>|a=0>|c1=1>|c=1> 

Now, recall that if $x\geq N$ we want the operation to leave the state unchanged. So, we need two cis.less operations, one to check whether $x < N$, which we store in ancilla qubit $c_1$ and another one to check whether $x < N-y$ stored in $c_2$. Note that the combination $c_1=0, c_2=1$ is not possible. The implementation looks as follows:

caddmodN <- function(c, bits, c1, c2, a, x, y, N) {
  stopifnot(length(a) == 1 && length(c1) == 1 &&
            length(c2) == 1 &&
            length(unique(c(c1, c2, a))) == 3)
  y <- y %% N
  ## set c1=1 if x < N
  z <- cis.less(c=c, bits=bits, x=x, c1=c1, a=a, y=N)
  ## set c2=1 if x < N - y
  z <- cis.less(c=c, bits=bits, x=z, c1=c2, a=a, y=N-y)
  
  ## if c1 and not c2, x = x + y - N
  z <- X(c2) *( CCNOT(c(c1, c2, a)) * (X(c2) * z))
  z <- cadd(c=a, bits=bits, x=z, y=y - N)
  z <- X(c2) * (CCNOT(c(c1, c2, a)) * (X(c2) * z))
  
  ## if c1 and c2 add x = x + y
  z <- CCNOT(c(c1, c2, a)) * z
  z <- cadd(c=a, bits=bits, x=z, y=y)
  z <- CCNOT(c(c1, c2, a)) * z
  
  ## reset c1,2
  z <- cis.less(c=c, bits=bits, x=z, c1=c2, a=a, y=y)
  z <- CNOT(c(c1, c2)) * z
  z <- cis.less(c=c, bits=bits, x=z, c1=c1, a=a, y=N)
  return(invisible(z))
}

For the reset part: in the first step we flip $c_2$ if $x + y \mod N < y$. This can only be true if $x>N-y$. If $c_2$ was 1, it’s zero now and the other way around. The next CNOT gate flips $c_2$ to $|0\rangle$. The last step resets $c_1$ to $|0\rangle$.

You can also see from the routine above that, if $x\geq N$ then caddmodN leaves the state unchanged.

Example

basis <- c()
for(i in c(0:(2^7-1))) {
  basis[i + 1] <-
    paste0("|", i %/% 16 , ">|a=", (i %/% 8) %% 2,
           ">|c2=", (i %/% 4) %% 2,
           ">|c1=", (i%/%2) %% 2, ">|c=", i%%2, ">")
}

x <- X(1)*qstate(7, basis=basis)
x

   ( 1 )    * |0>|a=0>|c2=0>|c1=0>|c=1> 

bits <- c(5,6,7)
c <- 1
c1 <- 2
c2 <- 3
a <- 4
N <- 5
z <- caddmodN(c=c, bits=bits, c1=c1, c2=c2, a=a, x=x, y=3, N=N) # 0 + 3 mod 5
z

   ( 1 )    * |3>|a=0>|c2=0>|c1=0>|c=1> 

z <- caddmodN(c=c, bits=bits, c1=c1, c2=c2, a=a, x=z, y=1, N=N) # 3 + 1 mod 5
z

   ( 1 )    * |4>|a=0>|c2=0>|c1=0>|c=1> 

z <- caddmodN(c=c, bits=bits, c1=c1, c2=c2, a=a, x=z, y=6, N=N) # 4 + 6 mod 5
z

   ( 1 )    * |0>|a=0>|c2=0>|c1=0>|c=1> 

Controlled Multiplier modulo $N$

Now, like the mult function above a version performing $|x\rangle\to|x+y\mod N\rangle$. Here we also include the un-computation of the second register. reg1 is the result register.

cmultmodN <- function(c, reg1, reg2, ancillas, x, y, N) {
  stopifnot(length(reg1) == length(reg2))
  ## need 4 ancilla registers
  stopifnot(length(ancillas) == 4 &&
            length(unique(ancillas)) == 4)
  n <- length(reg2)
  ## precompute terms in the sum
  s <- summands(y, N, n)
  ## start with |x>|0>
  for(i in c(1:n)) {
    x <- CCNOT(c(c, reg1[i], ancillas[4])) * x
    x <- caddmodN(c=ancillas[4], bits=reg2,
                  c1=ancillas[1], c2=ancillas[2],
                  a=ancillas[3],
                  x=x, y=s[i], N=N)
    x <- CCNOT(c(c, reg1[i], ancillas[4])) * x
  }
  ## now |x>|xy mod N>
  for(i in c(1:n)) {
    x <- CSWAP(c(c, reg1[i], reg2[i])) * x
  }
  ## now |xy mod N>|x>
  ## -y_inv mod N
  yinv <- N - moduloinverse(a=y, n=N)
  s <- summands(yinv, N, n)
  for(i in c(1:n)) {
    x <- CCNOT(c(c, reg1[i], ancillas[4])) * x
    x <- caddmodN(c=ancillas[4], bits=reg2,
                  c1=ancillas[1], c2=ancillas[2],
                  a=ancillas[3],
                  x=x, y=s[i], N=N)
    x <- CCNOT(c(c, reg1[i], ancillas[4])) * x
  }
  ## finally |xy mod N>|0>
  return(invisible(x))
}

For the un-computation of the second register, we start in the state (only the two registers reg1 and reg2) $$|x\cdot y \mod N\rangle|x\rangle\,.$$ Now we determine $yy_\mathrm{inv} = 1\mod N$ the modular inverse of $y$ (which is only possible, if $y$ and $N$ are co-prime). If we set $y' = xy\mod N$ it follows $y_\mathrm{inv}y' = y_\mathrm{inv}yx\mod N=x\mod N$. Thus, $x=y_\mathrm{inv}y'\mod N$. So, if we perform the following trafo $$|y'\rangle|x\rangle\to|y'\rangle|x-y_\mathrm{inv}y'\mod N\rangle = |y'\rangle|0\rangle\,.$$ we obtain $|x\cdot y \mod N\rangle|0\rangle$.

Example with two 3 qubit registers, which is starting to become slow, because in total we need 11 qubits

basis <- c()
for(i in c(0:(2^11-1))) {
  basis[i + 1] <-
    paste0("|reg1=", i %/% (32*2^3) , ">|reg2=", (i %/% 32) %% 2^3 ,
           "|anc=", (i %/% 16) %% 2,
           (i %/% 8) %% 2, (i %/% 4) %% 2,
           (i%/%2) %% 2, ">|c=", i%%2, ">")
}
x <- CNOT(c(1,10)) * (H(1)*qstate(11, basis=basis))
x

   ( 0.7071068 )    * |reg1=0>|reg2=0|anc=0000>|c=0> 
 + ( 0.7071068 )    * |reg1=2>|reg2=0|anc=0000>|c=1> 

c <- 1
ancillas <- c(2:5)
reg2 <- c(6:8)
reg1 <- c(9:11)
N <- 5
z <- cmultmodN(c=c, reg1=reg1, reg2=reg2,
               ancillas=ancillas, x=x, y=3, N=N)
z

   ( 0.7071068 )    * |reg1=0>|reg2=0|anc=0000>|c=0> 
 + ( 0.7071068 )    * |reg1=1>|reg2=0|anc=0000>|c=1> 

z <- cmultmodN(c=c, reg1=reg1, reg2=reg2,
               ancillas=ancillas, x=z, y=3, N=N)
z

   ( 0.7071068 )    * |reg1=0>|reg2=0|anc=0000>|c=0> 
 + ( 0.7071068 )    * |reg1=3>|reg2=0|anc=0000>|c=1> 

z <- cmultmodN(c=c, reg1=reg1, reg2=reg2,
               ancillas=ancillas, x=z, y=3, N=N)
z

   ( 0.7071068 )    * |reg1=0>|reg2=0|anc=0000>|c=0> 
 + ( 0.7071068 )    * |reg1=4>|reg2=0|anc=0000>|c=1> 

z <- cmultmodN(c=c, reg1=reg1, reg2=reg2,
               ancillas=ancillas, x=z, y=3, N=N)
z

   ( 0.7071068 )    * |reg1=0>|reg2=0|anc=0000>|c=0> 
 + ( 0.7071068 )    * |reg1=2>|reg2=0|anc=0000>|c=1> 

This seems to work up to this point.

cexpomodN <- function(c, reg1, reg2, ancillas, x, y, a, N) {
  stopifnot(length(reg1) == length(reg2))
  ## need 4 ancilla registers
  stopifnot(length(ancillas) == 4 &&
            length(unique(ancillas)) == 4)
  for(i in c(1:a)) {
    x <- cmultmodN(c=c, reg1=reg1, reg2=reg2,
                   ancillas=ancillas, x=x, y=y, N=N)
  }
  return(invisible(x))
}

x <- CNOT(c(1,10)) * (H(1)*qstate(11, basis=basis))
x

   ( 0.7071068 )    * |reg1=0>|reg2=0|anc=0000>|c=0> 
 + ( 0.7071068 )    * |reg1=2>|reg2=0|anc=0000>|c=1> 

x <- cexpomodN(c, reg1, reg2, ancillas, x, y=3, a=4, N)
x

   ( 0.7071068 )    * |reg1=0>|reg2=0|anc=0000>|c=0> 
 + ( 0.7071068 )    * |reg1=2>|reg2=0|anc=0000>|c=1> 

cexpomodN2 <- function(c, reg1, reg2, ancillas, x, y, a, N) {
  stopifnot(length(reg1) == length(reg2))
  ## need 4 ancilla registers
  stopifnot(length(ancillas) == 4 &&
            length(unique(ancillas)) == 4)
  ab <- as.integer(intToBits(a))
  n <- max(which(ab == 1))
  y2 <- y %% N
  for(i in c(1:n)) {
    if(ab[i] == 1) {
      x <- cmultmodN(c=c, reg1=reg1, reg2=reg2,
                     ancillas=ancillas, x=x, y=y2, N=N)
    }
    y2 <- ((y2%%N) * (y2%%N)) %% N # y2=y^(2^i) mod N
  }
  return(invisible(x))
}

x <- CNOT(c(1,10)) * (H(1)*qstate(11, basis=basis))
x

   ( 0.7071068 )    * |reg1=0>|reg2=0|anc=0000>|c=0> 
 + ( 0.7071068 )    * |reg1=2>|reg2=0|anc=0000>|c=1> 

x <- cexpomodN2(c, reg1, reg2, ancillas, x, y=3, a=4, N)
x

   ( 0.7071068 )    * |reg1=0>|reg2=0|anc=0000>|c=0> 
 + ( 0.7071068 )    * |reg1=2>|reg2=0|anc=0000>|c=1>