Skip to content

Instantly share code, notes, and snippets.

@philzook58
Created March 23, 2020 03:05
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save philzook58/95b6b6ef97c4b2047d9833b7993fa1ef to your computer and use it in GitHub Desktop.
Save philzook58/95b6b6ef97c4b2047d9833b7993fa1ef to your computer and use it in GitHub Desktop.
syzygies
def matrix_to_eqs(m):
nrows, ncols = m.shape
gens = [sy.Dummy() for i in range(ncols)]
eqs = m @ sy.Matrix(gens)
return eqs, gens
def eqs_to_matrix(eqns, gens):
return sy.Matrix( [[ eq.coeff(g) for g in gens] for eq in eqns])
import sympy as sy
def spoly(f,g,*gens):
ltf = sy.LT(f,gens)
ltg = sy.LT(g,gens)
lcm = sy.lcm(ltf,ltg)
s = lcm / ltf * f - lcm / ltg * g
return s
#grobner tracking. Maintaining the relation of the grobner basis to the original
def extended_groebner(F, *gens):
n = len(F)
markers = [sy.Dummy() for i in range(n)]
Fext = [ f + a for a, f in zip(markers, F)]
gen_ext = list(gens) + markers
Gext = sy.groebner(Fext, gen_ext)
A = [[g.coeff(m) for m in markers ] for g in Gext]
G = [ sy.simplify(g - sum( [ m * aa for m,aa in zip(markers, a) ] )) for a,g in zip(A,Gext) ] #remove dummy parts
assert( sy.simplify(sy.Matrix(G) - sy.Matrix(A) * sy.Matrix(F)).is_zero )
# maybe assert buchberger criterion
return G, A
def reduce_basis(F,G,*gens):
B,rems = list(zip(*[sy.reduced(f,G, gens) for f in F]))
print(B)
print(rems)
assert( all([r == 0 for r in rems] )) # assuming G is a grobner basis
assert( sy.simplify(sy.Matrix(F) - sy.Matrix(B) * sy.Matrix(G)).is_zero )
return B
# generators for the syzygies of G. Schreyer's Theorem. Cox Little Oshea Theorem 3.2 chapter 5
def syzygy(G,*gens): # assuming G is groebner basis
n = len(G)
basis = []
I = sy.eye(n)
for i, gi in enumerate(G):
for j, gj in enumerate(G):
if i < j:
s = spoly(gi,gj,*gens)
print(s)
a,r = sy.reduced( s , G, gens )
assert(r == 0) # should be groebner basis
lti = sy.LT(gi,gens)
ltj = sy.LT(gj,gens)
lcm = sy.lcm(lti,ltj)
ei = I[:,i]
ej = I[:,j]
basis.append(lcm / lti * ei - lcm / ltj * ej - sy.Matrix(a))
assert( sy.simplify(sy.Matrix(G).T * basis[-1]).is_zero) # should all null out on G
return basis
x,y,z,s = sy.symbols("x y z s")
F = [x+y+z, x*y+y*z+z*x, x*y*z-1]
G, A = extended_groebner(F, x,y,z)
B = reduce_basis(F,G,x,y,z)
Gsyz = syzygy(G)
@isuruf
Copy link

isuruf commented Oct 27, 2020

@philzook58, I want to add some of this code to SymPy. Would you be willing to let SymPy use this code with the license https://github.com/sympy/sympy/blob/master/LICENSE#L1-L28 ?

@philzook58
Copy link
Author

philzook58 commented Oct 27, 2020

@isuruf Sure do whatever you want with it. I'd recommend giving it thorough testing though. I don't recall how confident I was on correctness. There is a related blog post I wrote https://www.philipzucker.com/computing-syzygy-modules-in-sympy/

@isuruf
Copy link

isuruf commented Oct 27, 2020

Thanks @philzook58. I'll definitely add some tests and compare with singular.

@isuruf
Copy link

isuruf commented Nov 2, 2020

@philzook58, it turns out that sympy already has the code to this.

import sympy as sp
x,y,z,s = sp.symbols("x y z s")
F = [x+y+z, x*y+y*z+z*x, x*y*z-1]
ring = sp.EX.old_poly_ring(x, y, z)
print(ring.free_module(1).submodule(*[[f] for f in F]).syzygy_module())

Thanks for pointing out the blog post

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment