modarith.qcl

set allow-redefines 1;

include "functions.qcl";
include "qufunct.qcl";

// conditional addition mod n for 1 integer and 1 qureg
// flag is set if a+b<n for invertability

qufunct addn(int a,int n,quconst b,quvoid flag,quvoid sum,quconst e) {
  qureg s=sum[0\#b-1];
  qureg f=sum[#b-1];
  qureg bb=b;                      // "abuse" sum and b as scratch
  lt(n-a,bb,f,s);                  // for the less-than operator
  CNot(flag,f & e);                // save result of comparison
  !lt(n-a,bb,f,s);                 // restore sum and b
  muxadd(2^#b+a-n,a,flag,b,sum,e); // add either a or a-n
}

// Conditional overwriting addition mod n: sum -> (a+sum) mod n

qufunct oaddn(int a,int n,qureg sum,quconst e) {
  qureg j[#sum];
  qureg f[1];
  
  addn(a,n,sum,f,j,e);             // junk -> a+b mod n
  Swap(sum,j);                     // swap junk and sum
  CNot(f,e);                       // toggle flag
  !addn(n-a,n,sum,f,j,e);          // uncompute b to zero
}

// Conditional Multiplication mod n of an integer a by the qureg b,
// prod <- ab mod n.

qufunct muln(int a,int n,quconst b,qureg prod,quconst e) {
  int i;
  
  for i=0 to #prod-1 {
    if bit(a,i) { CNot(prod[i],b[0] & e); }
  }
  for i=1 to #b-1 {
    oaddn(2^i*a mod n,n,prod,b[i] & e);
  }
}

// Conditional Overwriting multiplication mod n: b-> ab mod n

qufunct omuln(int a,int n,qureg b,quconst e) {
  qureg j[#b];

  if gcd(a,n)>1 { 
    exit "omuln: a and n have to be relativly prime"; 
  }
  muln(a,n,b,j,e);
  !muln(invmod(a,n),n,j,b,e);
  cxor(j,b,e);
  cxor(b,j,e);
}

// Modular exponentiation: b -> x^a mod n

qufunct expn(int a,int n,quconst b,quvoid ex) {
  int i;
  
  Not(ex[0]);                            // start with 1
  for i=0 to #b-1 {
    omuln(powmod(a,2^i,n),n,ex,b[i]);    // ex -> ex*a^2^i mod n
  }
}

set allow-redefines 0;