share/quantum/quantum_computing.mac

   1 /*
   2    Quantum Computing Simulator for Maxima
   3    Copyright (C) 2021 Jaime E. Villate
   4
   5    Time-stamp: "2023-08-03 08:43:43 villate"
   6
   7    This program is free software; you can redistribute it and/or
   8    modify it under the terms of the GNU General Public License
   9    as published by the Free Software Foundation; either version 2
  10    of the License, or (at your option) any later version.
  11
  12    This program is distributed in the hope that it will be useful,
  13    but WITHOUT ANY WARRANTY; without even the implied warranty of
  14    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15    GNU General Public License for more details.
  16
  17    You should have received a copy of the GNU General Public License
  18    along with this program; if not, write to the Free Software
  19    Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
  20    MA  02110-1301, USA
  21 */
  22
  23 /*
  24    linsert
  25    Inserts an expression or list "e" into a list "list" at position "pos".
  26    The list could be empty and "pos" must be an integer between 1 and the
  27    length of "list" plus 1.
  28 */
  29 linsert(e, list, pos) := block(
  30   if not(listp(e)) then e: [e],
  31   append(rest(list, pos-length(list)-1), e, rest(list, pos-1)))$
  32
  33 /*
  34    lreplace
  35    If e is a list of length n, the elements in the positions pos, pos+1,
  36    ..., pos+n-1 of the list "list" are replaced by e, or the first elements
  37    of e if the end of "list" is reached.
  38    If e is an expression, the element in position "pos" of "list" is
  39    replaced by that expression.
  40    "pos" must be an integer between 1 and the length of "list".
  41 */
  42 lreplace(e, list, pos) := block([r: copylist(list)],
  43   if not(listp(e)) then e: [e],
  44   for i:0 thru min(length(e)-1, length(list)-pos) do r[pos+i]: e[i+1],
  45   r)$
  46
  47 /* binlist
  48    Should be called with one or two arguments, binlist(k) or binlist(k,n),
  49    where k and n are natural numbers. With only one argument, a list of
  50    binary digits is returned, which are the representation of k in binary.
  51    With two arguments, zeros are included at the beginning in order
  52    to return a list with exactly n binary digits. Notice that for the result
  53    to represent a possible state of m qubits, n should be equal to 2^m and
  54    k should be between 0 and 2^m-1.
  55 */
  56 binlist([n]):= block([m:n[1], q, l:[], k],
  57   do([m,q]: divide(m, 2), l: cons(q,l), if m=0
  58     then (
  59       if length(n)>1
  60       then (
  61         k: length(l),
  62         for i:1 thru n[2]-k do l:cons(0,l)),
  63       return(l))))$
  64
  65 /* binlist2dec
  66    Given a list with n binary digits, it returns the decimal number it
  67    represents.
  68 */
  69 binlist2dec(l) := block([m:l[1]],
  70   for i:2 thru length(l) do m:m*2+l[i],
  71   m)$
  72
  73 /* tprod
  74    Returns the tensor product of n matrices or n lists.
  75 */
  76 tprod([ops]) := block(
  77   if listp(ops[1])
  78   then list_matrix_entries(apply(tprod,map(matrix,ops)))
  79   else lreduce(kronecker_product, ops))$
  80
  81 /* normalize
  82    Returns the normalized version of a quantum state given as a list q.
  83 */
  84 normalize(q) := block([p:0],
  85   for b in q do p: p+cabs(b)^2,
  86   for i:1 thru length(q) do q[i]: rootscontract(radcan(q[i]/sqrt(p))),
  87   q)$
  88
  89 /* qubits
  90    Returns a list with 2^n elements, representing the state of n qubits.
  91    The argument can be a single natural number n, or n binary digits.
  92    If given only one natural number, greater than 1, it returns the ground
  93    state for a system of n qubits.
  94    If given n binary digits (n can be 1), it returns the state of a system
  95    of n qubits with those binary values.
  96 */
  97 qubits([b]) := block([l],
  98   if (length(b) = 1) and (first(b) > 1)
  99   then cons(1, makelist(0,2^first(b)-1))
 100   else (
 101     l: makelist(0,2^length(b)),
 102     l[binlist2dec(b)+1]: 1,
 103     l))$
 104
 105 /* gate
 106    Applies a matrix U, acting on mu qubits on the list q corresponding
 107    to a state of mq qubits (mq >= mu).
 108    When U acts on several qubits (mu > 1), and no more arguments are
 109    given, the matrix will act on qubits 1, 2, ..., mu. A third
 110    argument can be given, a integer from 1 through mq+1-mu, in which
 111    case the matrix will act on qubits i through i+mu-1.
 112    When U acts on only one qubit (mu = 1), and no more arguments are
 113    given, the matrix will act in all of the qubits. If more arguments
 114    are given, they should be different integer numbers between 1 and
 115    nq, and the matrix will act in the qubits corresponding to those
 116    numbers.
 117    U can be a symbol representing one of the predefined matrices.
 118    It modifies the list q and returns that new list.
 119 */
 120 gate(U,q,[bits]) := block([nq: length(q), mq, nu, mu, ln, lk, j, l, aux],
 121   if symbolp(U) then U: qmatrix[U],
 122   nu: length(U[1]),
 123   mq: round(float(log(nq)/log(2))),
 124   mu: round(float(log(nu)/log(2))),
 125   if length(bits) = 0
 126   then (
 127     if mu > 1
 128     then bits:[1]
 129     else bits: makelist(j, j, 1, mq)),
 130   for i in bits do (
 131     aux: makelist(0,nq),
 132     for n:1 thru nq do (
 133       ln: binlist(n-1, mq),
 134       j: binlist2dec(rest(rest(ln, i-1), mu+i-mq-1)) + 1,
 135       for k:1 thru nu do (
 136         lk: binlist(k-1, mu),
 137         l: binlist2dec(lreplace(lk, ln, i)) + 1,
 138         aux[n]: aux[n] + U[j][k]*q[l])),
 139     for n:1 thru nq do q[n]: aux[n]),
 140   q)$
 141
 142 /* gate_matrix
 143    Creates the matrix corresponding to the action of a single qbit matrix U
 144    acting on one of several of the qubits in a state of n qubits.
 145    The first argument must be a 2 by 2 matrix or a symbol representing one
 146    of the predefined matrices.
 147    There must be a second argument, a natural number, corresponding to the
 148    total number of qbits.
 149    If no more arguments are given, the matrix will act in all of the qubits.
 150    If there are more arguments, they should be one or more different
 151    natural numbers, between 1 and n, and the matrix will act in the qubits
 152    corresponding to those numbers.
 153    It returns a 2^n by 2^n matrix.
 154 */
 155 gate_matrix(U,n,[bits]) := block([lU],
 156   if symbolp(U) then U: qmatrix[U],
 157   if length(bits) = 0
 158   then apply(tprod, makelist(U, n))
 159   else (
 160     lU: makelist(ident(2), n),
 161     for i in bits do lU[i]: U,
 162     apply(tprod, lU)))$
 163
 164 /* controlled
 165    Applies a matrix U, acting on mu qubits, on qubits i through i+mu-1
 166    on the state q of a system of mq qubits (mq > mu), when the value of
 167    the c'th qubit in q equals 1.
 168    i should be an intenger between 1 and mq+1-mu.
 169    c should be an integer within 1 and mq, excluding the qubits to be
 170    modified (i through i+mu-1)
 171    U can be a symbol representing one of the predefined matrices.
 172    It modifies the list q and returns that new list.
 173 */
 174 controlled(U,q,c,i) := block([nq: length(q), mq, nu, mu, ln, lk, j, l, aux],
 175   if symbolp(U) then U: qmatrix[U],
 176   aux: makelist(0,nq),
 177   nu: length(U[1]),
 178   mq: round(float(log(nq)/log(2))),
 179   mu: round(float(log(nu)/log(2))),
 180   for n:1 thru nq do (
 181     ln: binlist(n-1, mq),
 182     if (ln[c] = 1)
 183     then (
 184       j: binlist2dec(rest(rest(ln, i-1), mu+i-mq-1)) + 1,
 185       for k:1 thru nu do (
 186         lk: binlist(k-1, mu),
 187         l: binlist2dec(lreplace(lk, ln, i)) + 1,
 188         aux[n]: aux[n] + U[j][k]*q[l]))
 189     else (
 190       aux[n]: q[n])),
 191   for n:1 thru nq do q[n]: aux[n],
 192   q)$
 193
 194 /* CNOT
 195    Changes the value of the j'th qubit, in a state q of m qubits, when
 196    the value of the i'th qubit equals 1.
 197    It modifies the list q and returns its modified value.
 198 */
 199
 200 CNOT(q,i,j) := block([lst: [], n: length(q), d, m, k],
 201   m: round(float(log(length(q))/log(2))),
 202   d: 2^m/2^j,
 203   if (m > 2)
 204   then (
 205     for b:0 thru n/4-1 do (
 206      lst: binlist(b, m-2),
 207       if (i < j)
 208       then lst: linsert(0, linsert(1,lst,i), j)
 209       else lst: linsert(1, linsert(0,lst,j), i),
 210       k: binlist2dec(lst)+1,
 211       [q[k],q[k+d]]: [q[k+d],q[k]]))
 212   else (
 213     [q[4-i],q[4-i+d]]: [q[4-i+d],q[4-i]]),
 214   q)$
 215
 216 /* qswap
 217    Interchanges the states of bits i and j in the qubits state q.
 218    It modifies the list q and returns its modified value.
 219 */
 220
 221 qswap(q, i, j) := (CNOT(q, i, j), CNOT(q, j, i), CNOT(q, i, j))$
 222
 223 /* toffoli
 224    Changes the value of the k'th qubit, in a state q of m qubits, when
 225    the value of the i'th snf j'th qubitd equal 1.
 226    It modifies the list q and returns its modified value.
 227 */
 228
 229 toffoli(q,i,j,k) := block([lst: [], n: length(q), c, d, m],
 230   m: round(float(log(length(q))/log(2))),
 231   d: 2^m/2^k,
 232   if (m > 3)
 233   then (
 234     for b:0 thru n/8-1 do (
 235      lst: binlist(b, m-3),
 236       if (i < j)
 237       then (
 238         if (i < k)
 239         then (
 240           if (j < k)
 241           then lst: linsert(0, linsert(1, linsert(1,lst,i), j), k)
 242           else lst: linsert(1, linsert(0, linsert(1,lst,i), k), j))
 243         else lst: linsert(1, linsert(0, linsert(1,lst,j), k), i))
 244       else if (i < k)
 245       then lst: linsert(0, linsert(1, linsert(1,lst,j), i), k)
 246       else (
 247         if (j < k)
 248         then lst: linsert(1, linsert(1, linsert(0,lst,k), i), j)
 249         else lst: linsert(1, linsert(1, linsert(0,lst,k), j), i)),
 250       c: binlist2dec(lst)+1,
 251       [q[c],q[c+d]]: [q[c+d],q[c]]))
 252   else (
 253     [q[8-d],q[8]]: [q[8],q[8-d]]),
 254   q)$
 255
 256 /* qmeasure
 257    Measures the value of one or more qubits in a system of n qubits with
 258    state q.
 259    It requires 1 or more arguments. The first argument must be the state q.
 260    If no more arguments are given, all the n qubits will be measured.
 261    If more arguments are given, they should be different natural numbers
 262    between 1 and n, and the state of those qubits will be measured.
 263    It returns a list with the values of the qubits measured (either 0 or 1),
 264    in the same order that the were requested or in ascending order if no
 265    second argument was given.
 266    It modifies the list q, reflecting the collapse of the quantum state
 267    after the measurement.
 268 */
 269 qmeasure(q,[bits]) := block([r, p, l:[], bl, n, m, v, result: [],
 270   test:true],
 271   n: length(q),
 272   m: round(float(log(n)/log(2))),
 273   if length(bits) = 0 then bits: makelist(j, j, 1, m),
 274   for i in bits do (
 275     r: random(1.0),
 276     p: 0,
 277     for j:1 thru n do (
 278       if not(test) then return(),
 279       if constantp(q[j])
 280       then (
 281         bl: binlist(j-1,m),
 282         if bl[i] = 0 then p: p+cabs(q[j])^2)
 283       else test:false),
 284     if test
 285     then (
 286       if r<p then v:0 else (v:1, p: 1-p),
 287       for j: 1 thru n do (
 288         bl:binlist(j-1,m),
 289         if bl[i]=v
 290         then (
 291           if atom(p)
 292           then q[j]: q[j]/sqrt(p)
 293           else q[j]: rootscontract(radcan(q[j]/sqrt(p))))
 294         else q[j]: 0),
 295       result: endcons(v, result))),
 296   result)$
 297
 298 /* qdisplay
 299    Represents the state q of a system of n qubits as a linear combination
 300    of the computational states with n binary digits.
 301    It returns an expression including strings and symbols.
 302 */
 303 qdisplay(q):=block([n, m, s, g:1, v, sg: ["-","",""], test:true],
 304   n:length(q),
 305   m:round(float(log(n)/log(2))),
 306   for e in q do
 307   /* Find the GCD only where the non-zero components are expressions
 308   **** Currently disabled ****
 309   if ((imagpart(e) = 0) and atom(e) and (sign(e) # 'zero)) then test: false,
 310   if test
 311   then g: 1/lreduce('gcd, q)
 312   else g: 1, */
 313   if (g = 1) and (length(q) > 1) then s: "" else s: "(",
 314   for i:1 thru n do (
 315     if (g = 1)
 316     then v: q[i]
 317     else v: q[i]*g,
 318     if ((imagpart(v) # 0) or (sign(v) # 'zero))
 319     then (
 320       if (numberp(v) or ((imagpart(v) = 0) and atom(v) and (sign(v) # 'pnz)))
 321       then (
 322         if ((imagpart(v^2-1) = 0) and (sign(v^2-1) = 'zero))
 323         then s:concat(s, sg[round((v+3)/2)], printf(false,"|~v,'0b>",m,i-1))
 324         elseif (sign(v) = 'neg)
 325         then s:concat(s, sg[3], printf(false,"~a|~v,'0b>",v,m,i-1))
 326         else s:concat(s, sg[2], printf(false,"~a|~v,'0b>",v,m,i-1)))
 327       elseif atom(v)
 328       then s: concat(s, sg[2], printf(false,"~a|~v,'0b>",v,m,i-1))
 329       else s: concat(s, sg[2], printf(false,"(~a)|~v,'0b>",v,m,i-1)),
 330       sg: [" -"," +"," "])),
 331   if (g # 1)
 332   then (
 333     if atom(g)
 334     then s: concat(s, printf(false,")/~a",g))
 335     else s: concat(s, printf(false,")/(~a)",g))),
 336   s)$
 337
 338 /*
 339    Predefined matrices.
 340 */
 341 qmatrix[I]: matrix([1,0],[0,1])$
 342 qmatrix[X]: matrix([0,1],[1,0])$
 343 qmatrix[Y]: matrix([0,-%i],[%i,0])$
 344 qmatrix[Z]: matrix([1,0],[0,-1])$
 345 qmatrix[H]: matrix([1,1],[1,-1])/sqrt(2)$
 346 qmatrix[S]: matrix([1,0],[0,%i])$
 347 qmatrix[T]: matrix([1,0],[0,exp(%i*%pi/4)])$
 348 Rx(ang) :=  matrix([cos(ang/2),-%i*sin(ang/2)],[-%i*sin(ang/2),cos(ang/2)])$
 349 Ry(ang) :=  matrix([cos(ang/2),-sin(ang/2)],[sin(ang/2),cos(ang/2)])$
 350 Rz(ang) :=  matrix([cos(ang/2)-%i*sin(ang/2),0],[0,cos(ang/2)+%i*sin(ang/2)])$