mjo/random.py

   1 r"""
   2 Generate random things.
   3 """
   4
   5 from sage.matrix.constructor import matrix
   6 from sage.rings.real_lazy import ComplexLazyField
   7 from sage.rings.all import ZZ
   8
   9 def random_unitary_matrix(F,n):
  10     r"""
  11     Generate a random unitary matrix of size ``n`` with entries
  12     in the field ``F``.
  13
  14     INPUT:
  15
  16     - ``F`` -- a field; specifically, a subfield of the complex
  17       numbers having characteristic zero.
  18
  19     - ``n`` -- integer; the size of the random matrix you want.
  20
  21     OUTPUT:
  22
  23     A random ``n``-by-``n`` unitary matrix with entries in ``F``. A
  24     ``ValueError`` is raised if ``F`` is not an appropriate field.
  25
  26     REFERENCES:
  27
  28     - Hans Liebeck and Anthony Osborne. The Generation of All Rational
  29       Orthogonal Matrices. The American Mathematical Monthly, Vol. 98,
  30       No. 2. (February, 1991), pp. 131-133.
  31
  32     SETUP::
  33
  34         sage: from mjo.random import random_unitary_matrix
  35
  36     TESTS::
  37
  38         sage: n = ZZ.random_element(10)
  39         sage: U = random_unitary_matrix(QQ,n)
  40         sage: U.is_unitary()
  41         True
  42         sage: U.base_ring() is QQ
  43         True
  44
  45     ::
  46
  47         sage: n = ZZ.random_element(10)
  48         sage: K = QuadraticField(-1,'i')
  49         sage: U = random_unitary_matrix(K,n)
  50         sage: U.is_unitary()
  51         True
  52         sage: U.base_ring() is K
  53         True
  54     """
  55     if not F.is_field():
  56         raise ValueError("F must be a field")
  57     elif not F.is_subring(ComplexLazyField()):
  58         raise ValueError("F must be a subfield of the complex numbers")
  59     elif not F.characteristic().is_zero():
  60         raise ValueError("F must have characteristic zero")
  61
  62     I = matrix.identity(F,n)
  63     A = matrix.random(F,n)
  64     S = A - A.conjugate_transpose()
  65     U = (S-I).inverse()*(S+I)
  66     D = matrix.identity(F,n)
  67     for i in range(n):
  68         if ZZ.random_element(2).is_zero():
  69             D[i,i] *= F(-1)
  70     return D*U