Function getRandomUnitary

tests/utils/random.cpp:348–376

qmatrix getRandomUnitary(int numQb) {
    DEMAND( numQb >= 1 );

    // create Z ~ random complex matrix (distribution not too important)
    size_t dim = getPow2(numQb);
    qmatrix matrZ = getRandomMatrix(dim);
    qmatrix matrZT = getTranspose(matrZ);

    // create Z = Q R (via QR decomposition) ...
    qmatrix matrQT = getOrthonormalisedRows(matrZ);
    qmatrix matrQ = getTranspose(matrQT);
    qmatrix matrR = getZeroMatrix(dim);

    // ... where R_rc = (columm c of Z) . (column r of Q) = (row c of ZT) . (row r of QT) = <r|c>
    for (size_t r=0; r<dim; r++)
        for (size_t c=r; c<dim; c++)
            matrR[r][c] = getInnerProduct(matrQT[r], matrZT[c]);

    // create D = normalised diagonal of R
    qmatrix matrD = getZeroMatrix(dim);
    for (size_t i=0; i<dim; i++)
        matrD[i][i] = matrR[i][i] / std::abs(matrR[i][i]);

    // create U = Q D
    qmatrix matrU = matrQ * matrD;

    DEMAND( isApproxUnitary(matrU) );
    return matrU;
}


qmatrix getRandomDiagonalUnitary(int numQb) {