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])) |

| 27 | def 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 | }} |