| 79 | from sage.misc.decorators import decorator_defaults |
| 80 | @decorator_defaults |
| 81 | def matrix_form(f, defaults = None): |
| 82 | if defaults is None: |
| 83 | defaults = {} |
| 84 | defaults.update(sparse=False) |
| 85 | from functools import wraps |
| 86 | @matrix_method |
| 87 | @wraps(f) |
| 88 | def r(*args, **kwds): |
| 89 | args = list(args) |
| 90 | ring = None |
| 91 | opts = {} |
| 92 | if rings.is_Ring(args[0]): |
| 93 | ring = args.pop(0) |
| 94 | for k in defaults: |
| 95 | opts[k] = kwds.pop(k, defaults[k]) |
| 96 | rows,cols,entries = f(*args, **kwds) |
| 97 | if ring is not None: |
| 98 | return matrix(ring, rows, cols, entries, **opts) |
| 99 | else: |
| 100 | return matrix(rows, cols, entries, **opts) |
| 101 | r._entries = f |
| 102 | return r |
| 103 | |
4079 | | |
4080 | | |
| 4110 | @matrix_form |
| 4111 | def hilbert_matrix(n): |
| 4112 | return n, n, lambda i,j: 1/(i+j+1) |
| 4113 | |
| 4114 | @matrix_form |
| 4115 | def vandermonde_matrix(v): |
| 4116 | return len(v), len(v), lambda i,j: v[i]^j |
| 4117 | |
| 4118 | @matrix_form |
| 4119 | def toeplitz_matrix(c,r): |
| 4120 | return len(c), len(r), lambda i,j: c[i-j] if i>=j else r[j-i] |
| 4121 | |
| 4122 | @matrix_form |
| 4123 | def hankel_matrix(c,r): |
| 4124 | entries=c+r[1:] |
| 4125 | return len(c), len(r), lambda i,j: entries[i+j] |
| 4126 | |
| 4127 | @matrix_form |
| 4128 | def circulant_matrix(E): |
| 4129 | return toeplitz_matrix._entries(E[0:1]+E[-1:0:-1], E) |
| 4130 | |
| 4131 | @matrix_form |
| 4132 | def skew_circulant_matrix(E): |
| 4133 | return hankel_matrix._entries(E, E[-1:]+E[:-1]) |
| 4134 | |
| 4135 | #Hadamard matrices: |
| 4136 | def legendre_symbol(x): |
| 4137 | """ |
| 4138 | Extend the built in legendre_symbol function to handle prime power |
| 4139 | fields. Assume x is an element of a finite field as well |
| 4140 | """ |
| 4141 | if x==0: |
| 4142 | return 0 |
| 4143 | elif x.is_square(): |
| 4144 | return 1 |
| 4145 | else: |
| 4146 | return -1 |
| 4147 | |
| 4148 | @matrix_method |
| 4149 | def jacobsthal_matrix(p,n): |
| 4150 | """ |
| 4151 | See http://en.wikipedia.org/wiki/Paley_construction for a way to |
| 4152 | use jacobsthal matrices to construct hadamard matrices |
| 4153 | """ |
| 4154 | if n == 1: |
| 4155 | elts = GF(p).list() |
| 4156 | else: |
| 4157 | elts = GF(p^n,'a').list() |
| 4158 | return matrix(len(elts), lambda i,j: legendre_symbol(elts[i]-elts[j])) |
| 4159 | |
| 4160 | @matrix_method |
| 4161 | def paley_matrix(p,n): |
| 4162 | """See http://en.wikipedia.org/wiki/Paley_construction""" |
| 4163 | mod = p^n%4 |
| 4164 | if mod == 3: |
| 4165 | # Paley Type 1 construction |
| 4166 | ones = vector([1]*p^n) |
| 4167 | QplusI = jacobsthal(p,n) |
| 4168 | # Q+=I efficiently |
| 4169 | for i in range(p^n): |
| 4170 | QplusI[i,i]=-1 |
| 4171 | return block_matrix(2,[ |
| 4172 | [1,ones.row()], |
| 4173 | [ones.column(), QplusI]]) |
| 4174 | elif mod == 1: |
| 4175 | # Paley Type 2 construction |
| 4176 | ones = vector([1]*p^n) |
| 4177 | QplusI = jacobsthal(p,n) |
| 4178 | QminusI = copy(QplusI) |
| 4179 | for i in range(p^n): |
| 4180 | QplusI[i,i]=1 |
| 4181 | QminusI[i,i]=-1 |
| 4182 | SplusI = block_matrix(2,[[1,ones.row()],[ones.column(), QplusI]]) |
| 4183 | SminusI = block_matrix(2,[[-1,ones.row()], [ones.column(), QminusI]]) |
| 4184 | return block_matrix(2,[[SplusI,SminusI],[SminusI,-SplusI]]) |
| 4185 | else: |
| 4186 | raise ValueError("p^n must be congruent to 1 or 3 mod 4") |