Ticket #13703: trac-13703-matrixform.patch

File trac-13703-matrixform.patch, 5.2 KB (added by jason, 7 years ago)

take 2, apply this patch *or* the other patch, but not both.

  • sage/matrix/constructor.py

    # HG changeset patch
    # User Jason Grout <jason.grout@drake.edu>
    # Date 1363386055 18000
    # Node ID 1dd0b3418dbf9af6be25cc4982842521c1cd2653
    # 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  
    7676    else:
    7777        return lambda func: matrix_method(func, name=name)
    7878
     79from sage.misc.decorators import decorator_defaults
     80@decorator_defaults
     81def matrix_form(f, defaults = None):
     82    if defaults is None:
     83        defaults = {}
     84    defaults.update(sparse=False)
     85    from functools import wraps
     86    @matrix_method
     87    @wraps(f)
     88    def r(*args, **kwds):
     89        args = list(args)
     90        ring = None
     91        opts = {}
     92        if rings.is_Ring(args[0]):
     93            ring = args.pop(0)
     94        for k in defaults:
     95            opts[k] = kwds.pop(k, defaults[k])
     96        rows,cols,entries = f(*args, **kwds)
     97        if ring is not None:
     98            return matrix(ring, rows, cols, entries, **opts)
     99        else:
     100            return matrix(rows, cols, entries, **opts)
     101    r._entries = f
     102    return r
     103
    79104def _matrix_constructor(*args, **kwds):
    80105    """
    81106    Create a matrix.
     
    13591384
    13601385
    13611386@matrix_method
    1362 def diagonal_matrix(arg0=None, arg1=None, arg2=None, sparse=True):
     1387def diagonal_matrix(arg0=None, arg1=None, arg2=None, sparse=True, offset=0):
    13631388    r"""
    13641389    Return a square matrix with specified diagonal entries, and zeros elsewhere.
    13651390   
     
    15431568   
    15441569    # Reconcile matrix size and number of entries
    15451570    try:
    1546         nentries = len(entries)
     1571        nentries = len(entries) + abs(offset)
    15471572    except TypeError:
    15481573        raise TypeError('unable to determine number of entries for diagonal matrix construction')
    15491574    # 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))
     1575    if nrows is not None and nentries > nrows:
     1576
     1577        raise ValueError('number of diagonal matrix entries (%s) offset at diagonal (%s) exceeds the requested matrix size (%s)' % (len(entries), offset, nrows))
    15521578    if nrows is None:
    15531579        nrows = nentries
    15541580   
     
    15681594    # Create a "diagonal" dictionary for matrix constructor
    15691595    # If nentries < nrows, diagonal is effectively padded with zeros at end
    15701596    w = {}
    1571     for i in range(len(v)):
    1572         w[(i,i)] = v[i]
    1573    
     1597    if offset>=0:
     1598        for i in range(len(v)):
     1599            w[(i,i+offset)] = v[i]
     1600    else:
     1601        offset = -offset
     1602        for i in range(len(v)):
     1603            w[(i+offset,i)] = v[i]
     1604
    15741605    # Ship ring, matrix size, dictionary to matrix constructor
    15751606    if ring is None:
    15761607        return matrix(nrows, nrows, w, sparse=sparse)
     
    40764107    entries.update({(j,j):aa, (j,i):bb, (i,j):-bb, (i,i):aa})
    40774108    return matrix(entries, nrows=dim, ring=ring)
    40784109
    4079 
    4080 
     4110@matrix_form
     4111def hilbert_matrix(n):
     4112    return n, n, lambda i,j: 1/(i+j+1)
     4113
     4114@matrix_form
     4115def vandermonde_matrix(v):
     4116    return len(v), len(v), lambda i,j: v[i]^j
     4117
     4118@matrix_form
     4119def toeplitz_matrix(c,r):
     4120    return len(c), len(r), lambda i,j: c[i-j] if i>=j else r[j-i]
     4121
     4122@matrix_form
     4123def hankel_matrix(c,r):
     4124    entries=c+r[1:]
     4125    return len(c), len(r), lambda i,j: entries[i+j]
     4126
     4127@matrix_form
     4128def circulant_matrix(E):
     4129    return toeplitz_matrix._entries(E[0:1]+E[-1:0:-1], E)
     4130
     4131@matrix_form
     4132def skew_circulant_matrix(E):
     4133    return hankel_matrix._entries(E, E[-1:]+E[:-1])
     4134
     4135#Hadamard matrices:
     4136def legendre_symbol(x):
     4137    """
     4138    Extend the built in legendre_symbol function to handle prime power
     4139    fields.  Assume x is an element of a finite field as well
     4140    """
     4141    if x==0:
     4142        return 0
     4143    elif x.is_square():
     4144        return 1
     4145    else:
     4146        return -1
     4147
     4148@matrix_method
     4149def jacobsthal_matrix(p,n):
     4150    """
     4151    See http://en.wikipedia.org/wiki/Paley_construction for a way to
     4152    use jacobsthal matrices to construct hadamard matrices
     4153    """
     4154    if n == 1:
     4155        elts = GF(p).list()
     4156    else:
     4157        elts = GF(p^n,'a').list()
     4158    return matrix(len(elts), lambda i,j: legendre_symbol(elts[i]-elts[j]))
     4159
     4160@matrix_method
     4161def paley_matrix(p,n):
     4162    """See http://en.wikipedia.org/wiki/Paley_construction"""
     4163    mod = p^n%4
     4164    if mod == 3:
     4165        # Paley Type 1 construction
     4166        ones = vector([1]*p^n)
     4167        QplusI = jacobsthal(p,n)
     4168        # Q+=I efficiently
     4169        for i in range(p^n):
     4170            QplusI[i,i]=-1
     4171        return block_matrix(2,[
     4172        [1,ones.row()],
     4173        [ones.column(), QplusI]])
     4174    elif mod == 1:
     4175        # Paley Type 2 construction
     4176        ones = vector([1]*p^n)
     4177        QplusI = jacobsthal(p,n)
     4178        QminusI = copy(QplusI)
     4179        for i in range(p^n):
     4180            QplusI[i,i]=1
     4181            QminusI[i,i]=-1
     4182        SplusI = block_matrix(2,[[1,ones.row()],[ones.column(), QplusI]])
     4183        SminusI = block_matrix(2,[[-1,ones.row()], [ones.column(), QminusI]])
     4184        return block_matrix(2,[[SplusI,SminusI],[SminusI,-SplusI]])       
     4185    else:
     4186        raise ValueError("p^n must be congruent to 1 or 3 mod 4")