We are currently solving the following issue within our repository. Here is the issue text:
--- BEGIN ISSUE ---
TwoQubitWeylDecomposition: failed to diagonalize M2 -- lack of normalization?
While writing [random test for 3586](https://github.com/Qiskit/qiskit-terra/pull/3586/files#diff-4dd440b6693c3194f4667ac2cdb99094), I noticed that a set this range wrongly (see https://github.com/Qiskit/qiskit-terra/pull/3586/commits/6e2b2709d8ffe702df8bce9a04089d8a9fb55245). However it seem that the previous range discover some crashes in `weyl_coordinates`.

### Steps to reproduce the problem

Currently, `two_qubit_decompose` is non-determinisic (see https://github.com/Qiskit/qiskit-terra/pull/3585), making reproduction complicated. Here are the situations I found:

```python
def test_cx_equivalence_1cx_random(rnd):
    qr = QuantumRegister(2, name='q')
    qc = QuantumCircuit(qr)

    qc.u3(rnd[0], rnd[1], rnd[2], qr[0])
    qc.u3(rnd[3], rnd[4], rnd[5], qr[1])

    qc.cx(qr[1], qr[0])

    qc.u3(rnd[6], rnd[7], rnd[8], qr[0])
    qc.u3(rnd[9], rnd[10], rnd[11], qr[1])

    sim = UnitarySimulatorPy()
    U = execute(qc, sim).result().get_unitary()
    two_qubit_cnot_decompose.num_basis_gates(U)


for _ in range(100):
    test_cx_equivalence_1cx_random(rnd=(0.0, 0.0, 0.0, 20895.999830149118, 4097.624991416554,
                                        0.00112247467014881, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0))
```

```python
def test_cx_equivalence_2cx_random(rnd):
    qr = QuantumRegister(2, name='q')
    qc = QuantumCircuit(qr)

    qc.u3(rnd[0], rnd[1], rnd[2], qr[0])
    qc.u3(rnd[3], rnd[4], rnd[5], qr[1])

    qc.cx(qr[1], qr[0])

    qc.u3(rnd[6], rnd[7], rnd[8], qr[0])
    qc.u3(rnd[9], rnd[10], rnd[11], qr[1])

    qc.cx(qr[0], qr[1])

    qc.u3(rnd[12], rnd[13], rnd[14], qr[0])
    qc.u3(rnd[15], rnd[16], rnd[17], qr[1])

    sim = UnitarySimulatorPy()
    U = execute(qc, sim).result().get_unitary()
    two_qubit_cnot_decompose.num_basis_gates(U)


for _ in range(100):
    test_cx_equivalence_2cx_random(rnd=(0.042351722707424384, 73735.43772411175,
                                        0.01757621764727491, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
                                        0.004490852354911424, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
                                        0.0))
```
```python
def test_cx_equivalence_3cx_random(rnd):
    qr = QuantumRegister(2, name='q')
    qc = QuantumCircuit(qr)

    qc.u3(rnd[0], rnd[1], rnd[2], qr[0])
    qc.u3(rnd[3], rnd[4], rnd[5], qr[1])

    qc.cx(qr[1], qr[0])

    qc.u3(rnd[6], rnd[7], rnd[8], qr[0])
    qc.u3(rnd[9], rnd[10], rnd[11], qr[1])

    qc.cx(qr[0], qr[1])

    qc.u3(rnd[12], rnd[13], rnd[14], qr[0])
    qc.u3(rnd[15], rnd[16], rnd[17], qr[1])

    qc.cx(qr[1], qr[0])

    qc.u3(rnd[18], rnd[19], rnd[20], qr[0])
    qc.u3(rnd[21], rnd[22], rnd[23], qr[1])

    sim = UnitarySimulatorPy()
    U = execute(qc, sim).result().get_unitary()
    two_qubit_cnot_decompose.num_basis_gates(U)

for _ in range(100):
    test_cx_equivalence_3cx_random(rnd=(0.17845153804438854, 8128.06054498254, 12287.99999713898,
                                        0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
                                        1983.8820662647765, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0))
```

All this cases trigger (eventually) `qiskit.exceptions.QiskitError: 'TwoQubitWeylDecomposition: failed to diagonalize M2'`.

Because they do not raise when the range is shorter (see https://github.com/Qiskit/qiskit-terra/pull/3586/commits/6e2b2709d8ffe702df8bce9a04089d8a9fb55245) I assume is lack of normalization.

--- END ISSUE ---

Below are some code segments, each from a relevant file. One or more of these files may contain bugs.
--- BEGIN FILES ---
Path: `qiskit/quantum_info/synthesis/weyl.py`
Content:
```
1 # -*- coding: utf-8 -*-
2 
3 # This code is part of Qiskit.
4 #
5 # (C) Copyright IBM 2017, 2019.
6 #
7 # This code is licensed under the Apache License, Version 2.0. You may
8 # obtain a copy of this license in the LICENSE.txt file in the root directory
9 # of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
10 #
11 # Any modifications or derivative works of this code must retain this
12 # copyright notice, and modified files need to carry a notice indicating
13 # that they have been altered from the originals.
14 # pylint: disable=invalid-name
15 
16 """Routines that compute  and use the Weyl chamber coordinates.
17 """
18 
19 import numpy as np
20 import scipy.linalg as la
21 from qiskit.exceptions import QiskitError
22 
23 _B = (1.0/np.sqrt(2)) * np.array([[1, 1j, 0, 0],
24                                   [0, 0, 1j, 1],
25                                   [0, 0, 1j, -1],
26                                   [1, -1j, 0, 0]], dtype=complex)
27 _Bd = _B.T.conj()
28 
29 
30 def weyl_coordinates(U):
31     """Computes the Weyl coordinates for
32     a given two-qubit unitary matrix.
33 
34     Args:
35         U (ndarray): Input two-qubit unitary.
36 
37     Returns:
38         ndarray: Array of Weyl coordinates.
39 
40     Raises:
41         QiskitError: Computed coordinates not in Weyl chamber.
42     """
43     pi2 = np.pi/2
44     pi4 = np.pi/4
45 
46     U = U / la.det(U)**(0.25)
47     Up = _Bd.dot(U).dot(_B)
48     M2 = Up.T.dot(Up)
49 
50     # M2 is a symmetric complex matrix. We need to decompose it as M2 = P D P^T where
51     # P ∈ SO(4), D is diagonal with unit-magnitude elements.
52     # D, P = la.eig(M2)  # this can fail for certain kinds of degeneracy
53     for _ in range(3):  # FIXME: this randomized algorithm is horrendous
54         M2real = np.random.randn()*M2.real + np.random.randn()*M2.imag
55         _, P = la.eigh(M2real)
56         D = P.T.dot(M2).dot(P).diagonal()
57         if np.allclose(P.dot(np.diag(D)).dot(P.T), M2, rtol=1.0e-13, atol=1.0e-13):
58             break
59     else:
60         raise QiskitError("TwoQubitWeylDecomposition: failed to diagonalize M2")
61 
62     d = -np.angle(D)/2
63     d[3] = -d[0]-d[1]-d[2]
64     cs = np.mod((d[:3]+d[3])/2, 2*np.pi)
65 
66     # Reorder the eigenvalues to get in the Weyl chamber
67     cstemp = np.mod(cs, pi2)
68     np.minimum(cstemp, pi2-cstemp, cstemp)
69     order = np.argsort(cstemp)[[1, 2, 0]]
70     cs = cs[order]
71     d[:3] = d[order]
72 
73     # Flip into Weyl chamber
74     if cs[0] > pi2:
75         cs[0] -= 3*pi2
76     if cs[1] > pi2:
77         cs[1] -= 3*pi2
78     conjs = 0
79     if cs[0] > pi4:
80         cs[0] = pi2-cs[0]
81         conjs += 1
82     if cs[1] > pi4:
83         cs[1] = pi2-cs[1]
84         conjs += 1
85     if cs[2] > pi2:
86         cs[2] -= 3*pi2
87     if conjs == 1:
88         cs[2] = pi2-cs[2]
89     if cs[2] > pi4:
90         cs[2] -= pi2
91 
92     return cs[[1, 0, 2]]
93
```
--- END FILES ---

Please first localize the bug based on the issue statement, and then generate a patch according to the `git diff` format fenced by three backticks.

Here is an example:

```diff
diff --git a/examples/server_async.py b/examples/server_async.py
--- a/examples/server_async.py
+++ b/examples/server_async.py
@@ -313,4 +313,4 @@
 
 
 if __name__ == "__main__":
-    asyncio.run(run_async_server("."), debug=True)
+    asyncio.run(run_async_server(), debug=True)
diff --git a/examples/server_sync.py b/examples/server_sync.py
--- a/examples/server_sync.py
+++ b/examples/server_sync.py
@@ -313,5 +313,5 @@
 
 
 if __name__ == "__main__":
-    server = run_sync_server(".")
+    server = run_sync_server()
     server.shutdown()

```

diff --git a/qiskit/quantum_info/synthesis/weyl.py b/qiskit/quantum_info/synthesis/weyl.py
--- a/qiskit/quantum_info/synthesis/weyl.py
+++ b/qiskit/quantum_info/synthesis/weyl.py
@@ -54,10 +54,13 @@
         M2real = np.random.randn()*M2.real + np.random.randn()*M2.imag
         _, P = la.eigh(M2real)
         D = P.T.dot(M2).dot(P).diagonal()
-        if np.allclose(P.dot(np.diag(D)).dot(P.T), M2, rtol=1.0e-13, atol=1.0e-13):
+        if np.allclose(P.dot(np.diag(D)).dot(P.T), M2, rtol=1.0e-10, atol=1.0e-10):
             break
     else:
-        raise QiskitError("TwoQubitWeylDecomposition: failed to diagonalize M2")
+        raise QiskitError("TwoQubitWeylDecomposition: failed to diagonalize M2. "
+                          "Please submit this output to "
+                          "https://github.com/Qiskit/qiskit-terra/issues/4159 "
+                          "Input %s" % U.tolist())
 
     d = -np.angle(D)/2
     d[3] = -d[0]-d[1]-d[2]

{"golden_diff": "diff --git a/qiskit/quantum_info/synthesis/weyl.py b/qiskit/quantum_info/synthesis/weyl.py\n--- a/qiskit/quantum_info/synthesis/weyl.py\n+++ b/qiskit/quantum_info/synthesis/weyl.py\n@@ -54,10 +54,13 @@\n         M2real = np.random.randn()*M2.real + np.random.randn()*M2.imag\n         _, P = la.eigh(M2real)\n         D = P.T.dot(M2).dot(P).diagonal()\n-        if np.allclose(P.dot(np.diag(D)).dot(P.T), M2, rtol=1.0e-13, atol=1.0e-13):\n+        if np.allclose(P.dot(np.diag(D)).dot(P.T), M2, rtol=1.0e-10, atol=1.0e-10):\n             break\n     else:\n-        raise QiskitError(\"TwoQubitWeylDecomposition: failed to diagonalize M2\")\n+        raise QiskitError(\"TwoQubitWeylDecomposition: failed to diagonalize M2. \"\n+                          \"Please submit this output to \"\n+                          \"https://github.com/Qiskit/qiskit-terra/issues/4159 \"\n+                          \"Input %s\" % U.tolist())\n \n     d = -np.angle(D)/2\n     d[3] = -d[0]-d[1]-d[2]\n", "issue": "TwoQubitWeylDecomposition: failed to diagonalize M2 -- lack of normalization?\nWhile writing [random test for 3586](https://github.com/Qiskit/qiskit-terra/pull/3586/files#diff-4dd440b6693c3194f4667ac2cdb99094), I noticed that a set this range wrongly (see https://github.com/Qiskit/qiskit-terra/pull/3586/commits/6e2b2709d8ffe702df8bce9a04089d8a9fb55245). However it seem that the previous range discover some crashes in `weyl_coordinates`.\r\n\r\n### Steps to reproduce the problem\r\n\r\nCurrently, `two_qubit_decompose` is non-determinisic (see https://github.com/Qiskit/qiskit-terra/pull/3585), making reproduction complicated. Here are the situations I found:\r\n\r\n```python\r\ndef test_cx_equivalence_1cx_random(rnd):\r\n    qr = QuantumRegister(2, name='q')\r\n    qc = QuantumCircuit(qr)\r\n\r\n    qc.u3(rnd[0], rnd[1], rnd[2], qr[0])\r\n    qc.u3(rnd[3], rnd[4], rnd[5], qr[1])\r\n\r\n    qc.cx(qr[1], qr[0])\r\n\r\n    qc.u3(rnd[6], rnd[7], rnd[8], qr[0])\r\n    qc.u3(rnd[9], rnd[10], rnd[11], qr[1])\r\n\r\n    sim = UnitarySimulatorPy()\r\n    U = execute(qc, sim).result().get_unitary()\r\n    two_qubit_cnot_decompose.num_basis_gates(U)\r\n\r\n\r\nfor _ in range(100):\r\n    test_cx_equivalence_1cx_random(rnd=(0.0, 0.0, 0.0, 20895.999830149118, 4097.624991416554,\r\n                                        0.00112247467014881, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0))\r\n```\r\n\r\n```python\r\ndef test_cx_equivalence_2cx_random(rnd):\r\n    qr = QuantumRegister(2, name='q')\r\n    qc = QuantumCircuit(qr)\r\n\r\n    qc.u3(rnd[0], rnd[1], rnd[2], qr[0])\r\n    qc.u3(rnd[3], rnd[4], rnd[5], qr[1])\r\n\r\n    qc.cx(qr[1], qr[0])\r\n\r\n    qc.u3(rnd[6], rnd[7], rnd[8], qr[0])\r\n    qc.u3(rnd[9], rnd[10], rnd[11], qr[1])\r\n\r\n    qc.cx(qr[0], qr[1])\r\n\r\n    qc.u3(rnd[12], rnd[13], rnd[14], qr[0])\r\n    qc.u3(rnd[15], rnd[16], rnd[17], qr[1])\r\n\r\n    sim = UnitarySimulatorPy()\r\n    U = execute(qc, sim).result().get_unitary()\r\n    two_qubit_cnot_decompose.num_basis_gates(U)\r\n\r\n\r\nfor _ in range(100):\r\n    test_cx_equivalence_2cx_random(rnd=(0.042351722707424384, 73735.43772411175,\r\n                                        0.01757621764727491, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,\r\n                                        0.004490852354911424, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,\r\n                                        0.0))\r\n```\r\n```python\r\ndef test_cx_equivalence_3cx_random(rnd):\r\n    qr = QuantumRegister(2, name='q')\r\n    qc = QuantumCircuit(qr)\r\n\r\n    qc.u3(rnd[0], rnd[1], rnd[2], qr[0])\r\n    qc.u3(rnd[3], rnd[4], rnd[5], qr[1])\r\n\r\n    qc.cx(qr[1], qr[0])\r\n\r\n    qc.u3(rnd[6], rnd[7], rnd[8], qr[0])\r\n    qc.u3(rnd[9], rnd[10], rnd[11], qr[1])\r\n\r\n    qc.cx(qr[0], qr[1])\r\n\r\n    qc.u3(rnd[12], rnd[13], rnd[14], qr[0])\r\n    qc.u3(rnd[15], rnd[16], rnd[17], qr[1])\r\n\r\n    qc.cx(qr[1], qr[0])\r\n\r\n    qc.u3(rnd[18], rnd[19], rnd[20], qr[0])\r\n    qc.u3(rnd[21], rnd[22], rnd[23], qr[1])\r\n\r\n    sim = UnitarySimulatorPy()\r\n    U = execute(qc, sim).result().get_unitary()\r\n    two_qubit_cnot_decompose.num_basis_gates(U)\r\n\r\nfor _ in range(100):\r\n    test_cx_equivalence_3cx_random(rnd=(0.17845153804438854, 8128.06054498254, 12287.99999713898,\r\n                                        0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,\r\n                                        1983.8820662647765, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0))\r\n```\r\n\r\nAll this cases trigger (eventually) `qiskit.exceptions.QiskitError: 'TwoQubitWeylDecomposition: failed to diagonalize M2'`.\r\n\r\nBecause they do not raise when the range is shorter (see https://github.com/Qiskit/qiskit-terra/pull/3586/commits/6e2b2709d8ffe702df8bce9a04089d8a9fb55245) I assume is lack of normalization.\n", "before_files": [{"content": "# -*- coding: utf-8 -*-\n\n# This code is part of Qiskit.\n#\n# (C) Copyright IBM 2017, 2019.\n#\n# This code is licensed under the Apache License, Version 2.0. You may\n# obtain a copy of this license in the LICENSE.txt file in the root directory\n# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.\n#\n# Any modifications or derivative works of this code must retain this\n# copyright notice, and modified files need to carry a notice indicating\n# that they have been altered from the originals.\n# pylint: disable=invalid-name\n\n\"\"\"Routines that compute  and use the Weyl chamber coordinates.\n\"\"\"\n\nimport numpy as np\nimport scipy.linalg as la\nfrom qiskit.exceptions import QiskitError\n\n_B = (1.0/np.sqrt(2)) * np.array([[1, 1j, 0, 0],\n                                  [0, 0, 1j, 1],\n                                  [0, 0, 1j, -1],\n                                  [1, -1j, 0, 0]], dtype=complex)\n_Bd = _B.T.conj()\n\n\ndef weyl_coordinates(U):\n    \"\"\"Computes the Weyl coordinates for\n    a given two-qubit unitary matrix.\n\n    Args:\n        U (ndarray): Input two-qubit unitary.\n\n    Returns:\n        ndarray: Array of Weyl coordinates.\n\n    Raises:\n        QiskitError: Computed coordinates not in Weyl chamber.\n    \"\"\"\n    pi2 = np.pi/2\n    pi4 = np.pi/4\n\n    U = U / la.det(U)**(0.25)\n    Up = _Bd.dot(U).dot(_B)\n    M2 = Up.T.dot(Up)\n\n    # M2 is a symmetric complex matrix. We need to decompose it as M2 = P D P^T where\n    # P \u2208 SO(4), D is diagonal with unit-magnitude elements.\n    # D, P = la.eig(M2)  # this can fail for certain kinds of degeneracy\n    for _ in range(3):  # FIXME: this randomized algorithm is horrendous\n        M2real = np.random.randn()*M2.real + np.random.randn()*M2.imag\n        _, P = la.eigh(M2real)\n        D = P.T.dot(M2).dot(P).diagonal()\n        if np.allclose(P.dot(np.diag(D)).dot(P.T), M2, rtol=1.0e-13, atol=1.0e-13):\n            break\n    else:\n        raise QiskitError(\"TwoQubitWeylDecomposition: failed to diagonalize M2\")\n\n    d = -np.angle(D)/2\n    d[3] = -d[0]-d[1]-d[2]\n    cs = np.mod((d[:3]+d[3])/2, 2*np.pi)\n\n    # Reorder the eigenvalues to get in the Weyl chamber\n    cstemp = np.mod(cs, pi2)\n    np.minimum(cstemp, pi2-cstemp, cstemp)\n    order = np.argsort(cstemp)[[1, 2, 0]]\n    cs = cs[order]\n    d[:3] = d[order]\n\n    # Flip into Weyl chamber\n    if cs[0] > pi2:\n        cs[0] -= 3*pi2\n    if cs[1] > pi2:\n        cs[1] -= 3*pi2\n    conjs = 0\n    if cs[0] > pi4:\n        cs[0] = pi2-cs[0]\n        conjs += 1\n    if cs[1] > pi4:\n        cs[1] = pi2-cs[1]\n        conjs += 1\n    if cs[2] > pi2:\n        cs[2] -= 3*pi2\n    if conjs == 1:\n        cs[2] = pi2-cs[2]\n    if cs[2] > pi4:\n        cs[2] -= pi2\n\n    return cs[[1, 0, 2]]\n", "path": "qiskit/quantum_info/synthesis/weyl.py"}], "after_files": [{"content": "# -*- coding: utf-8 -*-\n\n# This code is part of Qiskit.\n#\n# (C) Copyright IBM 2017, 2019.\n#\n# This code is licensed under the Apache License, Version 2.0. You may\n# obtain a copy of this license in the LICENSE.txt file in the root directory\n# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.\n#\n# Any modifications or derivative works of this code must retain this\n# copyright notice, and modified files need to carry a notice indicating\n# that they have been altered from the originals.\n# pylint: disable=invalid-name\n\n\"\"\"Routines that compute  and use the Weyl chamber coordinates.\n\"\"\"\n\nimport numpy as np\nimport scipy.linalg as la\nfrom qiskit.exceptions import QiskitError\n\n_B = (1.0/np.sqrt(2)) * np.array([[1, 1j, 0, 0],\n                                  [0, 0, 1j, 1],\n                                  [0, 0, 1j, -1],\n                                  [1, -1j, 0, 0]], dtype=complex)\n_Bd = _B.T.conj()\n\n\ndef weyl_coordinates(U):\n    \"\"\"Computes the Weyl coordinates for\n    a given two-qubit unitary matrix.\n\n    Args:\n        U (ndarray): Input two-qubit unitary.\n\n    Returns:\n        ndarray: Array of Weyl coordinates.\n\n    Raises:\n        QiskitError: Computed coordinates not in Weyl chamber.\n    \"\"\"\n    pi2 = np.pi/2\n    pi4 = np.pi/4\n\n    U = U / la.det(U)**(0.25)\n    Up = _Bd.dot(U).dot(_B)\n    M2 = Up.T.dot(Up)\n\n    # M2 is a symmetric complex matrix. We need to decompose it as M2 = P D P^T where\n    # P \u2208 SO(4), D is diagonal with unit-magnitude elements.\n    # D, P = la.eig(M2)  # this can fail for certain kinds of degeneracy\n    for _ in range(3):  # FIXME: this randomized algorithm is horrendous\n        M2real = np.random.randn()*M2.real + np.random.randn()*M2.imag\n        _, P = la.eigh(M2real)\n        D = P.T.dot(M2).dot(P).diagonal()\n        if np.allclose(P.dot(np.diag(D)).dot(P.T), M2, rtol=1.0e-10, atol=1.0e-10):\n            break\n    else:\n        raise QiskitError(\"TwoQubitWeylDecomposition: failed to diagonalize M2. \"\n                          \"Please submit this output to \"\n                          \"https://github.com/Qiskit/qiskit-terra/issues/4159 \"\n                          \"Input %s\" % U.tolist())\n\n    d = -np.angle(D)/2\n    d[3] = -d[0]-d[1]-d[2]\n    cs = np.mod((d[:3]+d[3])/2, 2*np.pi)\n\n    # Reorder the eigenvalues to get in the Weyl chamber\n    cstemp = np.mod(cs, pi2)\n    np.minimum(cstemp, pi2-cstemp, cstemp)\n    order = np.argsort(cstemp)[[1, 2, 0]]\n    cs = cs[order]\n    d[:3] = d[order]\n\n    # Flip into Weyl chamber\n    if cs[0] > pi2:\n        cs[0] -= 3*pi2\n    if cs[1] > pi2:\n        cs[1] -= 3*pi2\n    conjs = 0\n    if cs[0] > pi4:\n        cs[0] = pi2-cs[0]\n        conjs += 1\n    if cs[1] > pi4:\n        cs[1] = pi2-cs[1]\n        conjs += 1\n    if cs[2] > pi2:\n        cs[2] -= 3*pi2\n    if conjs == 1:\n        cs[2] = pi2-cs[2]\n    if cs[2] > pi4:\n        cs[2] -= pi2\n\n    return cs[[1, 0, 2]]\n", "path": "qiskit/quantum_info/synthesis/weyl.py"}]}