To use QDD, add the following header to your program :

#include "qdd.h"

# --- Location of QDD/BDD package installation
QDDLOC = ..

# --- Location of QDD/BDD library
LIBDIR = $(QDDLOC)/lib

# --- Location of QDD/BDD package header
INCDIR = $(QDDLOC)/include

# --- Linker flags for use with QDD (order and location *is* important)
QDDFLAGS = -lqdd -lbdd -lm -lg++

# --- Compiler flags
CFLAGS = -g -O3 -Wall

# --- Your C++ compiler
CPP = g++

all: shornuf

shornuf: shornuf.o
	$(CPP) shornuf.o -L$(LIBDIR) $(QDDFLAGS) -o shornuf

shornuf.o: shornuf.cc
	$(CPP) $(CFLAGS) -c shornuf.cc -I$(INCDIR) -o shornuf.o

void muxaddbit (int a0, int a1, QuantumBit & sel, QuantumBit & b, QuantumBit & ci_sum, QuantumBit & co, QuantumBit & e, int dir) {

// Note that the arguments here are object references.
// We don't allow QuantumBits/Registers to be copied at all.
// This helps guarantee quantum consitency in the model.

// We have to cheat to implement reversible operations :(

  if (dir > 0) {

// Allocate a single quantum bit M with an initial value of 0.

    QuantumBit M(0);


// The C++ "^=" operator has been overloaded to be the quantum XOR
// operator.  The C++ infix "&" allows accumulation of conditions for
// use in the quantum XOR.  In the following line, if the classical
// control signal a1 equals 1, then we will conditionally invert the
// quantum bit M in the quantum state where the quantum bits
// e and sel are true.


    if (a1) M ^= e & sel;

// Conditionally invert M if a0=1 when e is true and sel is false.

    if (a0) M ^= e & (! sel);

// .. etc 

    ci_sum ^= e & M;
    ci_sum ^= e & b;

    co  ^= e & M & (! ci_sum);
    co  ^= e & (! M) & b & (! ci_sum);
    co  ^= e & M & b & ci_sum ;

// Undo the operations performed on M

    if (a0) M ^= e & (! sel);

    if (a1) M ^= e & sel;

// By the time we get here, M should be in a pure state ..
// in fact, it should be in its original (false) state.
// (which I could/should have tested for, I guess :)

    try {
      M.unMixed();
    }
    catch (QuantumException e) {

// This is really a violation of quantum computing rules .. You cannot
// "tell" if a bit is in a pure state.  Nonetheless, such tests are 
// considered necessary during software development .. at least 
// it helped me :)


      cerr << "muxaddbit() : Mixed Bits\n" ;
      throw e;
    }

// At this point, C++ helps us do quantum garbage collection (automagically)
// Since M is in a pure state, we can reuse it later

  }

  else {

// The reverse operation.

    QuantumBit M(0);

    if (a1) M ^= e & sel;

    if (a0) M ^= e & (! sel);

    co  ^= e & M & b & ci_sum ;
    co  ^= e & (! M) & b & (! ci_sum);
    co  ^= e & M & (! ci_sum);

    ci_sum ^= e & b;
    ci_sum ^= e & M;

    if (a0) M ^= e & (! sel);

    if (a1) M ^= e & sel;

    try {
      M.unMixed();
    }
    catch (QuantumException e) {
      cerr << "muxaddbit(-1) : Mixed Bits\n" ;
      throw e;
    }

  }

}

Integer shornuf(Integer Number) {

// We use g++ BIGNUM Integers to make life easier when factoring really big numbers.

  Integer Factor;
  Integer f1,f2;
  Integer A;
  Integer Cv;
  Integer Xv;

  int size = ilog2(Number);

  cout << "Bits : " << size << "\n";

// Here we allocate a quantum register with size bits.

  QuantumRegister F(size);
  Integer Fv;

  QuantumRegister X(size);
  Integer Pv;

  int attempt = 1;
  int found = 0;
  while ((attempt < 8) && (! found)) {

    cout << "Attempt # " << attempt << "\n" ;

    do {
      
      A =  (rand(size) % Number) | 3 ;
    
    } while (gcd(Number,A) != 1) ;

    cout << "Random A : " << A << "\n";

// Here we set the value of the X register to all zeros.

    X.Set(0);

// We then perform an X.Mix() to place the register in a superposition of states.

    X.Mix();
            
    F.Set(0);

// Perform Quantum Exponentiation (Wow!)

    expn(A,Number,X,F, 1); // F = a^X mod n

// Perform a classical sampling of F.  This collapses the entire quantum
// state so that it is consistant with the value measured in F.

    Fv = F.Sample();

// Measure the period of X.  Unfortunately, a bit of a quantum cheat.

    Pv = X.Period();

    f1 = gcd((powmod(A,Pv/2,Number) + 1) % Number , Number) ;
    f2 = gcd((powmod(A,Pv/2,Number) + Number - 1) % Number , Number) ;

    Xv = X.Sample();

    Cv = powmod(A,Xv,Number);

    assert(Cv == Fv, "Mismatch");

    f1 = gcd(Number,f1);
    f2 = gcd(Number,f2);

    Factor = (f1 > f2 ? f1 : f2);

    found = (Factor > 1) && (Factor < Number);

    attempt++;

  }

  return(Factor);

}