Changes between Version 4 and Version 5 of Ticket #13703


Ignore:
Timestamp:
12/01/12 17:44:26 (7 years ago)
Author:
jason
Comment:

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 pn, and it looks like a reasonable modification of the wikipedia method.

Legend:

Unmodified
Added
Removed
Modified
  • Ticket #13703 – Description

    v4 v5  
    77def hankel(R,c,r): entries=c+r[1:]; return matrix(R, len(c), len(r), lambda i,j: entries[i+j])
    88def circulant(R,E): return hankel(R, E, E[-1:]+E[:-1])
     9
     10#Hadamard matrices:
     11def legendre_symbol(x):
     12    """Extend the built in legendre_symbol function to handle prime power fields.  Assume x is an element of a finite field as well"""
     13    if x==0:
     14        return 0
     15    elif x.is_square():
     16        return 1
     17    else:
     18        return -1
     19
    920def jacobsthal(p,n):
    1021    """See http://en.wikipedia.org/wiki/Paley_construction for a way to use jacobsthal matrices to construct hadamard matrices"""
    11     elts = GF(p^n).list()
    12     return matrix(len(elts), lambda i,j: legendre_symbol(elts[i]-elts[j],p))
    13 }}}
     22     if n == 1:
     23        elts = GF(p).list()
     24    else:
     25        elts = GF(p^n,'a').list()
     26    return matrix(len(elts), lambda i,j: legendre_symbol(elts[i]-elts[j]))
     27def paley_matrix(p,n):
     28    """See http://en.wikipedia.org/wiki/Paley_construction"""
     29    mod = p^n%4
     30    if mod == 3:
     31        # Paley Type 1 construction
     32        ones = vector([1]*p^n)
     33        QplusI = jacobsthal(p,n)
     34        # Q+=I efficiently
     35        for i in range(p^n):
     36            QplusI[i,i]=-1
     37        return block_matrix(2,[
     38        [1,ones.row()],
     39        [ones.column(), QplusI]])
     40    elif mod == 1:
     41        # Paley Type 2 construction
     42        ones = vector([1]*p^n)
     43        QplusI = jacobsthal(p,n)
     44        QminusI = copy(QplusI)
     45        for i in range(p^n):
     46            QplusI[i,i]=1
     47            QminusI[i,i]=-1
     48        SplusI = block_matrix(2,[[1,ones.row()],[ones.column(), QplusI]])
     49        SminusI = block_matrix(2,[[-1,ones.row()], [ones.column(), QminusI]])
     50        return block_matrix(2,[[SplusI,SminusI],[SminusI,-SplusI]])       
     51    else:
     52        raise ValueError("p^n must be congruent to 1 or 3 mod 4")
     53
     54}}
    1455
    1556Additionally, we could use scipy to create more matrices (or do it ourselves): http://docs.scipy.org/doc/scipy/reference/linalg.html#special-matrices