Main Program

word number;                    // number to be factored
word factor;                    // found factor
int width=duallog(number);      // length of N in bits
int nreg1=2*width,nreg2=width;  // width of registers
quBaseState qubase(nreg1+nreg2);// basestate
quWord reg1(nreg1,0,qubase);    // register 1
quWord reg2(nreg2,nreg1,qubase);// register 2
word x;                         // base value
word mreg1,mreg2;               // measurements of 1st and 2nd register
word pow;                       // pow^2==1 mod number
word a,b;                       // possible factors
word p,q;                       // fraction p/q for rational approximation
double qmax=1<<width;           // maximal period

while(1) {
  qubase.reset();                       // reseting state
  opFFT(nreg1)(reg1);                   // 1st Fourier transformation
  x=randcoprime(number);                // selecting random x
  opEXPN(nreg1,nreg2,x,number)(qubase); // modular exponentiation
  mreg2=reg2.measure().get-word();       // measure 2nd register
  opFFT(nreg1)(reg1);                   // 2nd Fourier transformation
  mreg1=reg1.select().getword();        // measure 1st register
  if(mreg1==0) continue;                // failed if measured zero

  // finding rational approximation for mreg1/rmax^2
  approx((double) mreg1/(qmax*qmax),(int) qmax,&p,&q);
  if(q%2) continue;                     // failed if q is odd.
  pow=powmod(x,q/2,number);             // pow = x^(q/2) mod number
  a=(pow+1)%number;                     // candidates with possible
  b=(pow+number-1)%number;              // common factors with number
  
  // testing for common factors with number                    
  if(a>1 && (factor=gcd(number,a))>1) break; 
  if(b>1 && (factor=gcd(number,b))>1) break;
};