#13703: special matrices
----------------------------------+-----------------------------------------
Reporter: jason | Owner: jason, was
Type: enhancement | Status: new
Priority: minor | Milestone: sage-5.9
Component: linear algebra | Resolution:
Keywords: | Work issues:
Report Upstream: N/A | Reviewers:
Authors: | Merged in:
Dependencies: | Stopgaps:
----------------------------------+-----------------------------------------
Description changed by jason:
Old description:
> It would be great to have a matrices namespace to put special matrix
> commands, sort of like the graphs.* namespace or groups.* namespace.
> Here are some starter definitions:
>
> {{{
> def hilbert(R,n): return matrix(R, n, lambda i,j: 1/(i+j+1))
> def vandermonde(R, v): return matrix(R, len(v), lambda i,j: v[i]^j)
> def toeplitz(R,c,r): return matrix(R, len(c), len(r), lambda i,j: c[i-j]
> if i>=j else r[j-i])
> def hankel(R,c,r): entries=c+r[1:]; return matrix(R, len(c), len(r),
> lambda i,j: entries[i+j])
> def circulant(R,E): return hankel(R, E, E[-1:]+E[:-1])
>
> #Hadamard matrices:
> def legendre_symbol(x):
> """Extend the built in legendre_symbol function to handle prime power
> fields. Assume x is an element of a finite field as well"""
> if x==0:
> return 0
> elif x.is_square():
> return 1
> else:
> return -1
>
> def jacobsthal(p,n):
> """See http://en.wikipedia.org/wiki/Paley_construction for a way to
> use jacobsthal matrices to construct hadamard matrices"""
> if n == 1:
> elts = GF(p).list()
> else:
> elts = GF(p^n,'a').list()
> return matrix(len(elts), lambda i,j:
> legendre_symbol(elts[i]-elts[j]))
> def paley_matrix(p,n):
> """See http://en.wikipedia.org/wiki/Paley_construction""";
> mod = p^n%4
> if mod == 3:
> # Paley Type 1 construction
> ones = vector([1]*p^n)
> QplusI = jacobsthal(p,n)
> # Q+=I efficiently
> for i in range(p^n):
> QplusI[i,i]=-1
> return block_matrix(2,[
> [1,ones.row()],
> [ones.column(), QplusI]])
> elif mod == 1:
> # Paley Type 2 construction
> ones = vector([1]*p^n)
> QplusI = jacobsthal(p,n)
> QminusI = copy(QplusI)
> for i in range(p^n):
> QplusI[i,i]=1
> QminusI[i,i]=-1
> SplusI = block_matrix(2,[[1,ones.row()],[ones.column(), QplusI]])
> SminusI = block_matrix(2,[[-1,ones.row()], [ones.column(),
> QminusI]])
> return block_matrix(2,[[SplusI,SminusI],[SminusI,-SplusI]])
> else:
> raise ValueError("p^n must be congruent to 1 or 3 mod 4")
>
> }}
>
> Additionally, we could use scipy to create more matrices (or do it
> ourselves): http://docs.scipy.org/doc/scipy/reference/linalg.html
> #special-matrices
New description:
It would be great to have a matrices namespace to put special matrix
commands, sort of like the graphs.* namespace or groups.* namespace. Here
are some starter definitions:
{{{
def hilbert(R,n): return matrix(R, n, lambda i,j: 1/(i+j+1))
def vandermonde(R, v): return matrix(R, len(v), lambda i,j: v[i]^j)
def toeplitz(R,c,r): return matrix(R, len(c), len(r), lambda i,j: c[i-j]
if i>=j else r[j-i])
def hankel(R,c,r): entries=c+r[1:]; return matrix(R, len(c), len(r),
lambda i,j: entries[i+j])
def circulant(R, E): return toeplitz(R, E[0:1]+E[-1:0:-1], E)
def skew_circulant(R,E): return hankel(R, E, E[-1:]+E[:-1])
#Hadamard matrices:
def legendre_symbol(x):
"""Extend the built in legendre_symbol function to handle prime power
fields. Assume x is an element of a finite field as well"""
if x==0:
return 0
elif x.is_square():
return 1
else:
return -1
def jacobsthal(p,n):
"""See http://en.wikipedia.org/wiki/Paley_construction for a way to
use jacobsthal matrices to construct hadamard matrices"""
if n == 1:
elts = GF(p).list()
else:
elts = GF(p^n,'a').list()
return matrix(len(elts), lambda i,j: legendre_symbol(elts[i]-elts[j]))
def paley_matrix(p,n):
"""See http://en.wikipedia.org/wiki/Paley_construction""";
mod = p^n%4
if mod == 3:
# Paley Type 1 construction
ones = vector([1]*p^n)
QplusI = jacobsthal(p,n)
# Q+=I efficiently
for i in range(p^n):
QplusI[i,i]=-1
return block_matrix(2,[
[1,ones.row()],
[ones.column(), QplusI]])
elif mod == 1:
# Paley Type 2 construction
ones = vector([1]*p^n)
QplusI = jacobsthal(p,n)
QminusI = copy(QplusI)
for i in range(p^n):
QplusI[i,i]=1
QminusI[i,i]=-1
SplusI = block_matrix(2,[[1,ones.row()],[ones.column(), QplusI]])
SminusI = block_matrix(2,[[-1,ones.row()], [ones.column(),
QminusI]])
return block_matrix(2,[[SplusI,SminusI],[SminusI,-SplusI]])
else:
raise ValueError("p^n must be congruent to 1 or 3 mod 4")
}}
Additionally, we could use scipy to create more matrices (or do it
ourselves): http://docs.scipy.org/doc/scipy/reference/linalg.html#special-
matrices
(thanks to pascal on sage-support for correcting the circulant code above:
https://groups.google.com/d/msg/sage-support/RnKjQ9n2YB0/vfCEvIV_HZUJ )
--
--
Ticket URL: <http://trac.sagemath.org/sage_trac/ticket/13703#comment:7>
Sage <http://www.sagemath.org>
Sage: Creating a Viable Open Source Alternative to Magma, Maple, Mathematica,
and MATLAB
--
You received this message because you are subscribed to the Google Groups
"sage-trac" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
To post to this group, send email to [email protected].
Visit this group at http://groups.google.com/group/sage-trac?hl=en.
For more options, visit https://groups.google.com/groups/opt_out.
