#13703: special matrices
----------------------------------+-----------------------------------------
       Reporter:  jason           |         Owner:  jason, was
           Type:  enhancement     |        Status:  new       
       Priority:  minor           |     Milestone:  sage-5.5  
      Component:  linear algebra  |    Resolution:            
       Keywords:                  |   Work issues:            
Report Upstream:  N/A             |     Reviewers:            
        Authors:                  |     Merged in:            
   Dependencies:                  |      Stopgaps:            
----------------------------------+-----------------------------------------

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])
> def jacobsthal(p,n):
>     """See http://en.wikipedia.org/wiki/Paley_construction for a way to
> use jacobsthal matrices to construct hadamard matrices"""
>     elts = GF(p^n).list()
>     return matrix(len(elts), lambda i,j:
> legendre_symbol(elts[i]-elts[j],p))
> }}}
>
> 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 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

--

Comment (by jason):

 More code, this time for Hadamard matrices.  Thanks to wikipedia and a
 paper by Alex Kramer summarizing these ideas.  I can't seem to find a
 source for the Paley type 2 matrix construction, but it seems to check out
 for small values of p^n, and it looks like a reasonable modification of
 the wikipedia method.

-- 
Ticket URL: <http://trac.sagemath.org/sage_trac/ticket/13703#comment:5>
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 post to this group, send email to [email protected].
To unsubscribe from this group, send email to 
[email protected].
For more options, visit this group at 
http://groups.google.com/group/sage-trac?hl=en.

Reply via email to