#13703: special matrices
----------------------------------+-----------------------------------------
Reporter: jason | Owner: jason, was
Type: enhancement | Status: needs_work
Priority: minor | Milestone: sage-5.9
Component: linear algebra | Resolution:
Keywords: | Work issues: documentation, doctests,
checking
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 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 )
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")
}}
The matrixform patch below first implements a decorator to take care of an
optional first argument of the ring, as well as an optional keyword
argument of 'sparse'. Then it implements most of these using that
decorator, so the boilerplate code is abstracted away.
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:11>
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.
