# Linear systems of graphs in Sage

Let $\Gamma$ be a graph. A divisor on $\Gamma$ is an element of the free group generated by the vertices $V$, ${\mathbb{Z}}[V]$.

We say that divisors $D$ and $D^\prime$ are linearly equivalent and write $D \sim D^\prime$ if $D^\prime-D$ is a principal divisor, i.e., if $D^\prime = D + \text{div}(f)$ for some function $f : V \rightarrow {\mathbb{Z}}$. Note that if $D$ and $D^\prime$ are linearly equivalent, they must have the same degree, since the degree of every principal divisor is 0. Divisors of degree 0 are linearly equivalent if and only if they determine the same element of the Jacobian. If $D$ is a divisor of degree 0, we denote by $[D]$ the element of the Jacobian determined by $D$. A divisor $D$ is said to be effective if $D(v) \geq 0$ for all vertices $v$. We write $D \geq 0$ to mean that $D$ is effective. The linear system associated to a divisor $D$ is the set

$|D| = \{ D^\prime \in \text{Div}(\Gamma ) : D^\prime \geq 0 \text{ and } D^\prime \sim D\},$

i.e., $|D|$ is the set of all effective divisors linearly equivalent to $D$. Note that if $D_1 \sim D_2$, then $|D_1| = |D_2|$. We note also that if $\text{deg}(D)<0$, then $|D|$ must be empty.

Sage can be used to compute the linear system of any divisor on a graph.

def linear_system(D, Gamma):
"""
Returns linear system attached to the divisor D.

EXAMPLES:
sage: Gamma2 = graphs.CubeGraph(2)
sage: Gamma1 = Gamma2.subgraph(vertices = ['00', '01'], edges = [('00', '01')])
sage: f = [['00', '01', '10', '11'], ['00', '01', '00', '01']]
sage: is_harmonic_graph_morphism(Gamma1, Gamma2, f)
True
sage: PhiV = matrix_of_graph_morphism_vertices(Gamma1, Gamma2, f); PhiV
[1 0 1 0]
[0 1 0 1]
sage: D = vector([1,0,0,1])
sage: PhiV*D
(1, 1)
sage: linear_system(PhiV*D, Gamma1)
[(2, 0), (1, 1), (0, 2)]
sage: linear_system(D, Gamma2)
[(0, 2, 0, 0), (0, 0, 2, 0), (1, 0, 0, 1)]
sage: [PhiV*x for x in linear_system(D, Gamma2)]
[(0, 2), (2, 0), (1, 1)]

"""
Q = Gamma.laplacian_matrix()
CS = Q.column_space()
N = len(D.list())
d = sum(D.list())
#print d
lin_sys = []
if d < 0:
return lin_sys
if (d == 0) and (D in CS):
lin_sys = [CS(0)]
return lin_sys
elif (d == 0):
return lin_sys
S = IntegerModRing(d+1)^N
V = QQ^N
for v in S:
v = V(v)
#print D-v,v,D
if D-v in CS:
lin_sys.append(v)
return lin_sys