Phase estimation

 Since we have the inverse QFFT, we can move to the phase estimation.

In the example, we are going to estimate the eigenvalues of the following matrix:

.

This matrix (or operator) can be written as a tensor product of 

and

.

.

This matrix is obviously unitary and can be used as the U operator in the diagram.
This matrix must be used with a control qubit, and I created a new "operation" matrixOperation in the following code.
The results agree with the example, and it seems that no need to repeat the measurement many times. The following table shows the number of measured bit strings out of 1,000 trials.

BTW, MS-Q# has QFT function which is "Adjoint"able (invertible), and the results are completely identical to mine. And also, I realize that the R1 gate should also take Adjoint, and I should use it instead of multiplying -1.0 to theta.

namespace phaseEstimation {

    open Microsoft.Quantum.Canon;

    open Microsoft.Quantum.Intrinsic;

    open Microsoft.Quantum.Convert;

    open Microsoft.Quantum.Arrays;

    open Microsoft.Quantum.Measurement;

    open Microsoft.Quantum.Math;

    open Microsoft.Quantum.Arithmetic as Arithmetic;

    

    @EntryPoint()

    operation SayHello() : Unit {

        let nQubit = 3; // Number of control qubits

        let maxNum = PowI(2,nQubit);

        mutable counts = new Int[maxNum];

        mutable theta  = new Double[maxNum];

        mutable revSeq = new Int[nQubit];

        for iQubit in 0..(nQubit-1){

            set theta w/= iQubit <- calcTheta(IntAsDouble(iQubit+1));

        }

        for iQubit in 0..(nQubit-1){

            set revSeq w/= iQubit <- nQubit-1-iQubit;

        }

        for trial in 1..1000 {

        use  (resultBits,vec) = (Qubit[nQubit],Qubit[2]) {

            ApplyToEach(Reset,resultBits); // Zero-clear the control qubits

            ApplyToEach(H,resultBits);

            //X(vec[0]);

            //X(vec[1]);

            for iQubit in 0..(nQubit-1){

                (Controlled matrixOperation)([resultBits[iQubit]],(iQubit,vec));

            }

            // Start Inverse QFFT

            let numLoopSwap = Floor(IntAsDouble(nQubit)*0.5);

            for iQubit in revSeq{

                for jQubit in 0..(nQubit-1)-iQubit-1{

                    let kQubit = revSeq[(nQubit-1)-jQubit];

                    let controlQubit = nQubit-kQubit-1;

                    let rotationIndex= nQubit-(kQubit+iQubit)-1;

                    (Controlled R1)([resultBits[controlQubit]],(-1.0*theta[rotationIndex],resultBits[iQubit]));

                }

                H(resultBits[iQubit]);

            }

            //for iQubit in 0..(numLoopSwap-1) {

            //    SWAP(resultBits[iQubit],resultBits[nQubit-1-iQubit]);

            //}              

            //let register = Arithmetic.BigEndian(resultBits);

            //(Adjoint QFT)(register);         

            let res = MultiM(resultBits);

            let bit = ResultArrayAsInt(res);

            set counts w/= bit <- counts[bit]+ 1;

            ResetAll(resultBits);

            ResetAll(vec);

        }

        }

        for j in 0..(maxNum-1)  {

           Message(IntAsString(counts[j]));

        }

    }

    function calcTheta(l : Double) : Double {

        return 2.0*PI()/(PowD(2.0,l));

    }

    operation matrixOperation(iQubit : Int, q : Qubit[]): Unit is Ctl{

       //Message(IntAsString(iQubit)+" "+IntAsString(PowI(2,iQubit)));

        for i in 0..(PowI(2,iQubit)-1){

            S(q[0]);

            T(q[1]);

        }

        return ();

    }

}