DumpOperation has unexpected global phase for double-controlled-R1 gate

[fedimser]

Describe the bug

Std.Diagnostics.DumpOperation, when applied to R1 gate with 2 controls, produces incorrect matrix.

To Reproduce

First, run this code:

operation Op(qs: Qubit[]) : Unit is Adj {
    Controlled Z([qs[0], qs[1]], (qs[2]));
}
operation Main() : Unit {
    Std.Diagnostics.DumpOperation(3, Op);
}

You will get this correct output:

 1.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖
 0.0000+0.0000𝑖,  1.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖
 0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  1.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖
 0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  1.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖
 0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  1.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖
 0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  1.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖
 0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  1.0000+0.0000𝑖,  0.0000+0.0000𝑖
 0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖, −1.0000+0.0000𝑖

Then, run this code:

operation Op(qs: Qubit[]) : Unit is Adj {
    Controlled R1([qs[0], qs[1]], (Std.Math.PI(), qs[2]));
}
operation Main() : Unit {
    Std.Diagnostics.DumpOperation(3, Op);
}

You will get:

 0.9239−0.3827𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖
 0.0000+0.0000𝑖,  0.9239−0.3827𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖
 0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.9239−0.3827𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖
 0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.9239−0.3827𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖
 0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.9239−0.3827𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖
 0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.9239−0.3827𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖
 0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.9239−0.3827𝑖,  0.0000+0.0000𝑖
 0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖,  0.0000+0.0000𝑖, −0.9239+0.3827𝑖

This is the correct matrix, multiplied by the global phase 0.9239−0.3827𝑖.

The second snippet should return the same matrix as the first, because by definition, R1(π)=diag(1, e^iπ)=diag(1,-1)=Z.

Expected behavior

Matrices from both snippets should match.

Screenshots

N/A

System information

qsharp 1.12.1

Additional context

• Documentation for DumpOperation does say that global phase is expected when there are other qubits allocated. However, in this example there are no other qubits besides the three qubits on which this operation acts. If this is intended behaviour, please update the documentation.
• I discovered this bug when migrating tests for my quantum_decomp library to use modern Q# in tests (and use DumpUnitary instead of calling DumpMachine many times). It worked for 2x2 and 4x4 matrices, but failed for 8x8 matrices. Temporary fix is to compare matrices up to global phase.