# -*- mode: octave -*-
# no function file
1;

function retval = initqreg(qubits)
# initialize to |0>
  if (nargin != 1)
    usage("initqreg(qubits)");
  endif

  retval = zeros(2^qubits,1);
  retval(1) = 1;
endfunction

# this is the same as the builtin 'eye'
# but not as efficient; example function
function retval = identity(dim)
  if (nargin != 1)
    usage ("identity (dim)");
  endif

  retval = zeros(dim)
  for i = 1 : dim
    retval(i,i) = 1;
  endfor
endfunction

# returns an operator that does op on the specified pair of qubits
# op should be a two-qubit(4x4) unitary
function retval = twoqubitopon(numqubits,op,bita,bitb)
  if (nargin != 4)
    usage("oponbits(numqubits,op,bita,bitb)");
  endif

  t= numqubits;
  retval = swapop(t,t-2,bitb)*swapop(t,t-1,bita)*addlbitstoop(op,numqubits-2)*swapop(t,t-1,bita)*swapop(t,t-2,bitb);
endfunction

# returns an operator that does X
# on the bit-th bit of a line of
# qubits (index starting from 0!)
# all this hoodoo because I can't figure out
# how to do bitwise ops; really what I want
# to do is toggle the bit-th bit in the
# 2nd dimension
function retval = Xop(numqubits,bit)
  if (nargin != 2)
    usage ("Xop (numqubits,bit)");
  endif
  if (bit >= numqubits)
    error("bit must be less than numqubits");
  endif

  dim = 2^numqubits;
  retval = zeros(dim);
  v = 2^bit; vp = 2^(bit+1);
  for i = 0 : dim-1
    j = mod(i,vp);
    k = i - j;
    j += v;
    j = mod(j,vp);
    k += j;
    retval(i+1,k+1) = 1;
  endfor
endfunction

# returns an operator that does X on all bits for dim
function retval = Xops(numqubits)
  if (nargin != 1)
    usage ("Xops (numqubits)");
  endif

  dim = 2^numqubits;
  retval = zeros(dim);
  for i = 1 : dim
    retval(i,dim-i+1) = 1;
  endfor
endfunction

# do this just for one qubit for now
function retval = Xrot(theta)
  retval = [ cos(theta/2), -i*sin(theta/2);
	     -i*sin(theta/2), cos(theta/2) ];
endfunction

function retval = Yop(numqubits,bit)
  if (nargin != 2)
    usage ("Yop (numqubits,bit)");
  endif
  if (bit >= numqubits)
    error("bit must be less than numqubits");
  endif

  dim = 2^numqubits;
  retval = zeros(dim);
  v = 2^bit; vp = 2^(bit+1);
  for q = 0 : dim-1
    j = mod(q,vp);
    k = q - j;
    j += v;
    j = mod(j,vp);
    k += j;
    if (j >= v)
      retval(q+1,k+1) = -i;
    else
      retval(q+1,k+1) = i;
    endif
  endfor
endfunction

function retval = Yops(numqubits)
  if (nargin != 1)
    usage("Yops(numqubits)");
  endif

  dim = 2^numqubits;
  retval = eye(dim);
  for j = 0 : numqubits-1
    retval *= Yop(numqubits,j);
  endfor
endfunction

# do this just for one qubit for now
function retval = Yrot(theta)
  retval = [ cos(theta/2), -sin(theta/2);
	     sin(theta/2), cos(theta/2) ]
endfunction

# I'm pretty sure this isn't the most efficient way to do all of this,
# but it works.
function retval = Hadamard(numqubits,bit)
  if (nargin != 2)
    usage ("Hadamard (dim,bit)");
  endif
  if (bit >= numqubits)
    error("bit must be less than numqubits");
  endif

  dim = 2^numqubits;
  retval = zeros(dim);
  v = 2^bit; vp = 2^(bit+1);
# this is horrible code
  q = 0;
  while (q < dim)
    j = mod(q,vp);
    k = q - j;
    j += v;
    j = mod(j,vp);
    k += j;
    for r = 0 : v-1;
      retval(k+r+1,k+r+1) = 1;
      retval(q+r+1,k+r+1) = 1;
      retval(k+r+1,q+r+1) = 1;
      retval(q+r+1,q+r+1) = -1;
    endfor
    q += v;
  endwhile
  retval /= sqrt(2);
endfunction

function retval = Hadamards(numqubits)
  if (nargin != 1)
    usage("Hadamards(numqubits)");
  endif

  dim = 2^numqubits;
  retval = eye(dim);
  for j = 0 : numqubits-1
    retval *= Hadamard(numqubits,j);
  endfor
endfunction

# Z, phase, and pi/8 gates are all the same except for phase
function retval = eyeonebitmod(numqubits,bit,bitval)
  if (nargin != 3)
    usage ("eyeonebitmod (numqubits,bit,bitval)");
  endif
  if (bit >= numqubits)
    error("bit must be less than numqubits");
  endif

  dim = 2^numqubits;
  retval = zeros(dim);
  v = 2^bit; vp = 2^(bit+1);
  for q = 0 : dim-1
    j = mod(q,vp);
    k = mod(q,v);
    if (j == k)
      retval(q+1,q+1) = 1;
    else
      retval(q+1,q+1) = bitval;
    endif
  endfor
endfunction

function retval = Zop(numqubits,bit)
  if (nargin != 2)
    usage ("Zop (numqubits,bit)");
  endif
  if (bit >= numqubits)
    error("bit must be less than numqubits");
  endif

  retval = eyeonebitmod(numqubits,bit,-1);
endfunction

function retval = Zops(numqubits)
  if (nargin != 1)
    usage("Zops(numqubits)");
  endif

  dim = 2^numqubits;
  retval = eye(dim);
  for j = 0 : numqubits-1
    retval *= Zop(numqubits,j);
  endfor
endfunction

# do this just for one qubit for now
function retval = Zrot(theta)
  retval = [ e^(-i*theta/2), 0;
	     0, e^(i*theta/2) ]
endfunction

function retval = Phaseop(numqubits,bit)
  if (nargin != 2)
    usage ("Phaseop (numqubits,bit)");
  endif
  if (bit >= numqubits)
    error("bit must be less than numqubits");
  endif

  retval = eyeonebitmod(numqubits,bit,i);
endfunction

function retval = Phaseops(numqubits)
  if (nargin != 1)
    usage("Phaseops(numqubits)");
  endif

  dim = 2^numqubits;
  retval = eye(dim);
  for j = 0 : numqubits-1
    retval *= Phaseop(numqubits,j);
  endfor
endfunction

# Pi/8 gate
function retval = Top(numqubits,bit)
  if (nargin != 2)
    usage ("Top (numqubits,bit)");
  endif
  if (bit >= numqubits)
    error("bit must be less than numqubits");
  endif

  retval = eyeonebitmod(numqubits,bit,e^(i*pi/4));
endfunction

function retval = Tops(numqubits)
  if (nargin != 1)
    usage("Tops(numqubits)");
  endif

  dim = 2^numqubits;
  retval = eye(dim);
  for j = 0 : numqubits-1
    retval *= Top(numqubits,j);
  endfor
endfunction

# control-NOT for N control lines, N=numqubits-1
# low-order bit is the one toggled
# read "C to the N NOT"
function retval = ctonnot(numqubits)
  if (nargin != 1)
    usage("ctonnot(numqubits)");
  endif

  dim = 2^numqubits;
  retval = zeros(dim);
  for i = 1 : dim-2
    retval(i,i) = 1;
  endfor
  retval(dim-1,dim) = 1;
  retval(dim,dim-1) = 1;
endfunction

function retval = cnot()
  retval = ctonnot(2);
endfunction

function retval = ccnot()
  retval = ctonnot(3);
endfunction

function retval = swapop(numqubits,bita,bitb)
  if (nargin != 3)
    usage("swapop(numqubits,bita,bitb)");
  endif

  dim = 2^numqubits;
  twoa = 2^bita;
  twob = 2^bitb;
  retval = zeros(dim);
  for j = 0 : dim-1
    aval = mod(j,2^(bita+1)) - mod(j,twoa); # should mask to give only bit a
    bval = mod(j,2^(bitb+1)) - mod(j,twob); # should mask to give only bit b
    if ((aval & bval) | (!aval & !bval))
      retval(j+1,j+1) = 1;
    elseif (aval)
      retval(j+1,j-aval+twob+1) = 1;
    else
      retval(j+1,j-bval+twoa+1) = 1;
    endif
  endfor
endfunction

# for NMR
function retval = Jgate
  retval = sqrt(-i)*[1,0,0,0; 0,i,0,0; 0,0,i,0; 0,0,0,1]
endfunction

function retval = qftrot(d)
  retval = [ 1,0,0,0; 0,1,0,0; 0,0,1,0; 0,0,0, e^(i*pi/(2^d))]
endfunction

# j & k are zero-based
function retval = qftrotb(j,k,l)
  a = qftrot(k-j);
  b = addhbitstoop(a,l-2);
  c = swapop(l,k,1);
  d = swapop(l,j,0);
  retval = b*c*d
endfunction

function retval = addlbitstoop(op,lowbits)
  if (nargin != 2)
    usage("addlbitstoop(op,lowbits)");
  endif
  [ opsizer, opsizec ] = size(op);
  if (opsizer != opsizec)
    error("shouldn't op be square?");
  endif
  lbdim = 2^lowbits;
  retval = zeros(opsizer*lbdim);
  for j = 0 : opsizer-1
    for k = 0 : opsizec-1
      for l = 1 : lbdim
	retval(j*lbdim+l,k*lbdim+l) = op(j+1,k+1);
      endfor
    endfor
  endfor
endfunction

function retval = addhbitstoop(op,hibits)
  if (nargin != 2)
    usage("addhbitstoop(op,hibits)");
  endif
  [ opsizer, opsizec ] = size(op);
  if (opsizer != opsizec)
    error("shouldn't op be square?");
  endif
  hbdim = 2^hibits;
  retval = zeros(opsizer*hbdim);
  for j = 0 : opsizer-1
    for k = 0 : opsizec-1
      for l = 0 : hbdim-1
	retval(l*opsizer+j+1,l*opsizec+k+1) = op(j+1,k+1);
      endfor
    endfor
  endfor
endfunction

# adds one high control bit to an op
function retval = ctlop(op)
  [ opsizer, opsizec ] = size(op);
  if (opsizer != opsizec)
    error("shouldn't op be square?");
  endif
  for j = 1 : opsizer
    retval(j,j) = 1;
  endfor
  for j = 1 : opsizer
    for k = 1 : opsizec
      retval(j+opsizer,k+opsizec) = op(j,k);
    endfor
  endfor
endfunction

# adds one low control bit to an op
function retval = lctlop(op)
  [ opsizer, opsizec ] = size(op);
  if (opsizer != opsizec)
    error("shouldn't op be square?");
  endif
  for j = 1 : opsizer
    retval(2*j-1,2*j-1) = 1;
  endfor
  for j = 1 : opsizer
    for k = 1 : opsizec
      retval(2*j,2*k) = op(j,k);
    endfor
  endfor
endfunction


# Schonhage-Strassen multiplication, see Knuth .v. 2 sec. 4.3.3, pp. 308ish
# takes in two bignums (sort of) in arrays, in a basis already set up for this.
function retval = ssmult(u,v)
  if (nargin != 2)
    usage("ssmult(u,v)");
  endif
  [ usizer, usizec ] = size(u);
  if (usizec != 1)
    error("shouldn't u be a vector?");
  endif
  [ vsizer, vsizec ] = size(v);
  if (vsizec != 1)
    error("shouldn't v be a vector?");
  endif
  uhat = fft(u);
  vhat = fft(v);
  what = uhat .* vhat;
  w = ifft(what);
  retval = w;
endfunction

function retval = sseval(w,basis)
  if (nargin != 2)
    usage("sseval(w,basis)");
  endif
  [ wsizer, wsizec ] = size(w);
  if (wsizec != 1)
    error("shouldn't w be a vector?");
  endif
  retval = 0;
  for j = 1 : wsizer
    retval += w(j)*basis^(j-1);
  endfor
endfunction

function retval = fivegateccnot()
  srx = 1/2*[ 1+i, 1-i; 1-i, 1+i ];
  srxd = srx';		# dagger
#  srx = Xrot(-pi/2);
#  srxd = Xrot(pi/2);		# dagger
  csrx = ctlop (srx);
  csrxd = ctlop(srxd);
  csrx13 = ctlop(addhbitstoop (srx,1));
  csrx23 = addhbitstoop(ctlop(srx),1);
  csrxd23 = addhbitstoop(ctlop(srxd),1);
  cnot12 = ctlop (addlbitstoop (Xop(1,0),1));

  myfgcc = csrx13*cnot12*csrxd23*cnot12*csrx23;
  retval = myfgcc;
endfunction

# Vedral CARRY construct on non-neighbors
function retval = vcarry()
  srx = 1/2*[ 1+i, 1-i; 1-i, 1+i ];
  srxd = srx';		# dagger
#  srx = Xrot(-pi/2);
#  srxd = Xrot(pi/2);		# dagger
  csrx = ctlop (srx);
  csrxd = ctlop(srxd);
  csrx14 = ctlop(addhbitstoop (srx,2));
  csrx24 = addhbitstoop(ctlop(addhbitstoop (srx,1)),1);
  cnot13 = ctlop (addhbitstoop (addlbitstoop (Xop(1,0),1),1));
  csrxd34 = addhbitstoop (csrxd,2);
  cnot23 = addhbitstoop(ctlop (addlbitstoop (Xop(1,0),1),1),1);
  csrx34 = addhbitstoop (csrx,2);

  myv = csrx14*cnot13*csrxd34*cnot13*csrx24*cnot23*csrx34;
  retval = myv;
endfunction

# Vedral CARRY construct w/ neighbor-only gates
# 4 ctl-sqrt-X, 17 CNOT, latency of 19
# trailing 6 CNOTs to be optimized away
# leading 4 gates to be pipelined
# core latency 9
function retval = vcarryn()
  srx = 1/2*[ 1+i, 1-i; 1-i, 1+i ];
  srxd = srx';		# dagger
  csrx = ctlop (srx);
  csrxd = ctlop(srxd);
  csrxd34 = addhbitstoop (csrxd,2);
  csrx34 = addhbitstoop (csrx,2);
				# all of the CNOTs we need...
  cnot12 = addlbitstoop (ctlop (Xop(1,0),1),2);
  cnot21 = addlbitstoop (lctlop (Xop(1,0),1),2);
  cnot23 = addhbitstoop(addlbitstoop (ctlop (Xop(1,0),1),1),1);
  cnot32 = addhbitstoop(addlbitstoop (lctlop (Xop(1,0),1),1),1);

  myv = cnot12*cnot21*cnot12 *cnot32*cnot23*cnot32 *csrx34 \
      *cnot23*cnot32 *cnot12*cnot21*cnot12 *csrxd34 \
      *cnot32*cnot23 \
      *cnot12 *csrx34 *cnot32\
      *cnot23*cnot32 *csrx34;

  retval = myv;
endfunction

# prints the non-zero elements of a state vector
function qprint(a)
  printed = 0;
  str = "";
  [ rows, cols ] = size(a);
  for j = 1 : rows
    if (a(j) != 0)
      thisstr = idx2bitstr(j-1,log2(rows));
      if printed
	str = strcat(str, " + ");
      endif
      str = strcat(str, "(", com2str(a(j)), ")", thisstr);
      printed = 1;
    endif
  endfor
#  if (a(2) != 0)
#    if (printed)
#      str = strcat(str, " + (", com2str(a(2)), ")|1>");
#    else
#      str = strcat(str, "(", com2str(a(2)), ")|1>");
#    endif
#  endif
  disp (str)
endfunction

# prints the non-almost-zero elements of a state vector
# lets us supress uninteresting states either because of
# rounding errors or algorithm residue
function qprintwthreshold(a,threshold)
  printed = 0;
  str = "";
  [ rows, cols ] = size(a);
  for j = 1 : rows
    if (abs(a(j)) > threshold)
      thisstr = idx2bitstr(j-1,log2(rows));
      if printed
	str = strcat(str, " + ");
      endif
      str = strcat(str, "(", com2str(a(j)), ")", thisstr);
      printed = 1;
    endif
  endfor
#  if (a(2) != 0)
#    if (printed)
#      str = strcat(str, " + (", com2str(a(2)), ")|1>");
#    else
#      str = strcat(str, "(", com2str(a(2)), ")|1>");
#    endif
#  endif
  disp (str)
endfunction

function retval = idx2bitstr(j,len)
  j;
  printed = 0;
  retval = "|";
  for k = 0 : len-1
    k = len - k - 1;
    twok = 2^k;
    if (j >= twok)
      retval = strcat(retval, "1");
      j = j - twok;
      printed = 1;
    else
      retval = strcat(retval, "0");
    endif
  endfor
  retval = strcat(retval, ">");
endfunction

function retval = idx2bitstrsuppressleadzeros(j)
  j;
  maxqubits = 12;		# limit it to 12 bits; how do I make a global?
  printed = 0;
  retval = "|";
  for k = 0 : maxqubits
    k = maxqubits - k;
    twok = 2^k;
    if (j >= twok)
      retval = strcat(retval, "1");
      j = j - twok;
      printed = 1;
    else
      if printed | k == 0
	retval = strcat(retval, "0");
      endif
    endif
  endfor
  retval = strcat(retval, ">");
endfunction

function retval = noiseeq(a,b,thresh)
  if (nargin != 3)
    usage("noiseeq(a,b,thresh)");
  endif
  [ asizer, asizec ] = size(a);
  if (asizer != asizec)
    error("shouldn't a be square?");
  endif
  [ bsizer, bsizec ] = size(b);
  if (bsizer != bsizec)
    error("shouldn't b be square?");
  endif
  if (asizer != bsizer)
    error("stick to the same size a and b");
  endif
  for j = 1 : asizer
    for k = 1 : asizec
      if (abs(abs(a(j,k))-abs(b(j,k))) < thresh)
	retval(j,k) = 1;
      else
	retval(j,k) = 0;
      endif
    endfor
  endfor
endfunction

# a function to return the probability of measuring zero
function retval = zeroprob(sv,bitno)
  if (nargin != 2)
    usage("zeroprob(statevector,bitno)");
  endif
  [ svr, svc ] = size(sv);
  mask1 = 2^bitno;
  mask2 = 2^(bitno+1);
  retval = 0.0;
  for j = 0 : svr-1
    if (mod(j,mask1) == mod(j,mask2))
      # then we know the bit is zero
      retval += abs(sv(j+1))^2;
    endif
  endfor  
endfunction

# now a function to measure out a single qubit from a register...
# returns a reduced-dimension state vector
function retval = measure(sv,bitno)
  if (nargin != 2)
    usage("measure(sv,bitno)");
  endif
  zp = zeroprob(sv,bitno);
  szp = sqrt(zp);
  myrand = rand();
  if (myrand <= zp)
    printf("measured qubit %d: zero\n",bitno);
  else
    printf("measured qubit %d: one\n",bitno);
  endif
  [ svr, svc ] = size(sv);
  mask1 = 2^bitno;
  mask2 = 2^(bitno+1);
  retval = zeros(svr/2,1);
  for j = 0 : svr/2 - 1
    if (myrand <= zp)
      retval(j+1) = sv((j-mod(j,mask1))*2+mod(j,mask1)+1)/szp;
    else			# it's a one
      retval(j+1) = sv((j-mod(j,mask1))*2+mod(j,mask1)+mask1+1)/szp;
    endif
  endfor
endfunction

function retval = makerandomqubit()
  retval = initqreg(1);
  myrand = rand();
  retval = Xrot(myrand*2*pi)*retval;
  myrand = rand();
  retval = Zrot(myrand*2*pi)*retval;
endfunction

# take two unentangled registers and create a state vector
# for the combo; will still be mathematically separable,
# but this is where we start combining them
function retval = mergetwoqregisters(a,b)
  [ ar, ac ] = size(a);
  [ br, bc ] = size(b);
  retval = zeros(ar*br,1);
  for ari = 0 : ar-1
    for abi = 0 : br-1
      retval(ari*br+abi+1) = a(ari+1)*b(abi+1);
    endfor
  endfor
endfunction

function retval = makeeprpair()
  retval = initqreg(2);
  retval = Hadamard(2,1)*retval;
  retval = cnot()*retval;
endfunction

# reports the sum of the probability amplitudes; should
# always be 1.0, to within FP precision
function retval = totprob(a)
  retval = 0.0;
  [ ar, ac ] = size(a);
  for ari = 1 : ar
    retval = retval + abs(a(ari))^2;
  endfor
endfunction

function retval = blochvector(a)
  if (imag(a(1)) != 0.0)
    pla = a * abs(a(1))/a(1);	# phaseless a
  else
    pla = a;
  endif
  theta = acos(pla(1))*2;
  if (theta != 0.0)
    psi = log(a(2)/sin(theta/2))/i;
  else
    psi = 0.0;
  endif
  retval = [ cos(psi)*sin(theta) ; sin(psi)*sin(theta) ; cos(theta) ];
endfunction

function plot3blocharrowv(a)
  ar = real(a);
  u = [0 : pi/10 : 2*pi];
  v = [0 : pi/10 : 2*pi];
  plot3([ 0 ar(1) ], [ 0 ar(2) ], [ 0 ar(3) ], "k", 
	cos(u)' * cos(v), sin(u)' * cos(v), ones(21,1) * sin(v), "b", 
	(cos(u)' * cos(v)), (sin(u)' * cos(v))', (ones(21,1) * sin(v))', "b");
endfunction

function plot3blocharrowq(a)
  plot3blocharrowv(blochvector(a));
endfunction