From adb4bc82d22373575f3fb96966fd0703781b3cf6 Mon Sep 17 00:00:00 2001
From: Michael Orlitzky <michael@orlitzky.com>
Date: Thu, 2 Feb 2023 09:19:58 -0500
Subject: [PATCH] mjo/random.py: new module for random_unitary_matrix().

---
 mjo/all.py    |  1 +
 mjo/random.py | 70 +++++++++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 71 insertions(+)
 create mode 100644 mjo/random.py

diff --git a/mjo/all.py b/mjo/all.py
index 78550cc..2a2cbc9 100644
--- a/mjo/all.py
+++ b/mjo/all.py
@@ -12,3 +12,4 @@ from mjo.hurwitz import Octonions
 from mjo.orthogonal_polynomials import *
 from mjo.polynomial import *
 from mjo.symbol_sequence import *
+from mjo.random import *
diff --git a/mjo/random.py b/mjo/random.py
new file mode 100644
index 0000000..f08c201
--- /dev/null
+++ b/mjo/random.py
@@ -0,0 +1,70 @@
+r"""
+Generate random things.
+"""
+
+from sage.matrix.constructor import matrix
+from sage.rings.real_lazy import ComplexLazyField
+from sage.rings.all import ZZ
+
+def random_unitary_matrix(F,n):
+    r"""
+    Generate a random unitary matrix of size ``n`` with entries
+    in the field ``F``.
+
+    INPUT:
+
+    - ``F`` -- a field; specifically, a subfield of the complex
+      numbers having characteristic zero.
+
+    - ``n`` -- integer; the size of the random matrix you want.
+
+    OUTPUT:
+
+    A random ``n``-by-``n`` unitary matrix with entries in ``F``. A
+    ``ValueError`` is raised if ``F`` is not an appropriate field.
+
+    REFERENCES:
+
+    - Hans Liebeck and Anthony Osborne. The Generation of All Rational
+      Orthogonal Matrices. The American Mathematical Monthly, Vol. 98,
+      No. 2. (February, 1991), pp. 131-133.
+
+    SETUP::
+
+        sage: from mjo.random import random_unitary_matrix
+
+    TESTS::
+
+        sage: n = ZZ.random_element(10)
+        sage: U = random_unitary_matrix(QQ,n)
+        sage: U.is_unitary()
+        True
+        sage: U.base_ring() is QQ
+        True
+
+    ::
+
+        sage: n = ZZ.random_element(10)
+        sage: K = QuadraticField(-1,'i')
+        sage: U = random_unitary_matrix(K,n)
+        sage: U.is_unitary()
+        True
+        sage: U.base_ring() is K
+        True
+    """
+    if not F.is_field():
+        raise ValueError("F must be a field")
+    elif not F.is_subring(ComplexLazyField()):
+        raise ValueError("F must be a subfield of the complex numbers")
+    elif not F.characteristic().is_zero():
+        raise ValueError("F must have characteristic zero")
+
+    I = matrix.identity(F,n)
+    A = matrix.random(F,n)
+    S = A - A.conjugate_transpose()
+    U = (S-I).inverse()*(S+I)
+    D = matrix.identity(F,n)
+    for i in range(n):
+        if ZZ.random_element(2).is_zero():
+            D[i,i] *= F(-1)
+    return D*U
-- 
2.43.2