Ticket #13703: trac-13703.patch

File trac-13703.patch, 4.4 KB (added by jason, 7 years ago)
  • sage/matrix/constructor.py

    # HG changeset patch
    # User Jason Grout <jason.grout@drake.edu>
    # Date 1363360163 18000
    # Node ID 8982bcf0e98fe80e02cae7cbac044328b9ad54be
    # Parent  f65732213e49b466c5288e32bb2f64204a241353
    Implement several more common matrix constructors, and give diagonal_matrix an offset argument.
    
    diff --git a/sage/matrix/constructor.py b/sage/matrix/constructor.py
    a b  
    13591359
    13601360
    13611361@matrix_method
    1362 def diagonal_matrix(arg0=None, arg1=None, arg2=None, sparse=True):
     1362def diagonal_matrix(arg0=None, arg1=None, arg2=None, sparse=True, offset=0):
    13631363    r"""
    13641364    Return a square matrix with specified diagonal entries, and zeros elsewhere.
    13651365   
     
    15431543   
    15441544    # Reconcile matrix size and number of entries
    15451545    try:
    1546         nentries = len(entries)
     1546        nentries = len(entries) + abs(offset)
    15471547    except TypeError:
    15481548        raise TypeError('unable to determine number of entries for diagonal matrix construction')
    15491549    # sometimes catches a negative size
    1550     if not nrows is None and nentries > nrows:
    1551         raise ValueError('number of diagonal matrix entries (%s) exceeds the requested matrix size (%s)' % (nentries, nrows))
     1550    if nrows is not None and nentries > nrows:
     1551
     1552        raise ValueError('number of diagonal matrix entries (%s) offset at diagonal (%s) exceeds the requested matrix size (%s)' % (len(entries), offset, nrows))
    15521553    if nrows is None:
    15531554        nrows = nentries
    15541555   
     
    15681569    # Create a "diagonal" dictionary for matrix constructor
    15691570    # If nentries < nrows, diagonal is effectively padded with zeros at end
    15701571    w = {}
    1571     for i in range(len(v)):
    1572         w[(i,i)] = v[i]
    1573    
     1572    if offset>=0:
     1573        for i in range(len(v)):
     1574            w[(i,i+offset)] = v[i]
     1575    else:
     1576        offset = -offset
     1577        for i in range(len(v)):
     1578            w[(i+offset,i)] = v[i]
     1579
    15741580    # Ship ring, matrix size, dictionary to matrix constructor
    15751581    if ring is None:
    15761582        return matrix(nrows, nrows, w, sparse=sparse)
     
    40764082    entries.update({(j,j):aa, (j,i):bb, (i,j):-bb, (i,i):aa})
    40774083    return matrix(entries, nrows=dim, ring=ring)
    40784084
    4079 
    4080 
     4085@matrix_method
     4086def hilbert_matrix(R,n):
     4087    return matrix(R, n, lambda i,j: 1/(i+j+1))
     4088
     4089@matrix_method
     4090def vandermonde_matrix(R, v):
     4091    return matrix(R, len(v), lambda i,j: v[i]^j)
     4092
     4093@matrix_method
     4094def toeplitz_matrix(R,c,r):
     4095    return matrix(R, len(c), len(r), lambda i,j: c[i-j] if i>=j else r[j-i])
     4096
     4097@matrix_method
     4098def hankel_matrix(R,c,r):
     4099    entries=c+r[1:]
     4100    return matrix(R, len(c), len(r), lambda i,j: entries[i+j])
     4101
     4102@matrix_method
     4103def circulant(R, E):
     4104    return toeplitz_matrix(R, E[0:1]+E[-1:0:-1], E)
     4105
     4106@matrix_method
     4107def skew_circulant_matrix(R,E):
     4108    return hankel_matrix(R, E, E[-1:]+E[:-1])
     4109
     4110#Hadamard matrices:
     4111def legendre_symbol(x):
     4112    """
     4113    Extend the built in legendre_symbol function to handle prime power
     4114    fields.  Assume x is an element of a finite field as well
     4115    """
     4116    if x==0:
     4117        return 0
     4118    elif x.is_square():
     4119        return 1
     4120    else:
     4121        return -1
     4122
     4123@matrix_method
     4124def jacobsthal_matrix(p,n):
     4125    """
     4126    See http://en.wikipedia.org/wiki/Paley_construction for a way to
     4127    use jacobsthal matrices to construct hadamard matrices
     4128    """
     4129    if n == 1:
     4130        elts = GF(p).list()
     4131    else:
     4132        elts = GF(p^n,'a').list()
     4133    return matrix(len(elts), lambda i,j: legendre_symbol(elts[i]-elts[j]))
     4134
     4135@matrix_method
     4136def paley_matrix(p,n):
     4137    """See http://en.wikipedia.org/wiki/Paley_construction"""
     4138    mod = p^n%4
     4139    if mod == 3:
     4140        # Paley Type 1 construction
     4141        ones = vector([1]*p^n)
     4142        QplusI = jacobsthal(p,n)
     4143        # Q+=I efficiently
     4144        for i in range(p^n):
     4145            QplusI[i,i]=-1
     4146        return block_matrix(2,[
     4147        [1,ones.row()],
     4148        [ones.column(), QplusI]])
     4149    elif mod == 1:
     4150        # Paley Type 2 construction
     4151        ones = vector([1]*p^n)
     4152        QplusI = jacobsthal(p,n)
     4153        QminusI = copy(QplusI)
     4154        for i in range(p^n):
     4155            QplusI[i,i]=1
     4156            QminusI[i,i]=-1
     4157        SplusI = block_matrix(2,[[1,ones.row()],[ones.column(), QplusI]])
     4158        SminusI = block_matrix(2,[[-1,ones.row()], [ones.column(), QminusI]])
     4159        return block_matrix(2,[[SplusI,SminusI],[SminusI,-SplusI]])       
     4160    else:
     4161        raise ValueError("p^n must be congruent to 1 or 3 mod 4")