1/*
2
3  GRAPHS - graph theory package for Maxima
4  Copyright (C) 2007-2011 Andrej Vodopivec <andrej.vodopivec@gmail.com>
5
6  This program is free software; you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation; either version 2 of the License, or
9  (at your option) any later version.
10
11  This program is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14  GNU General Public License for more details.
15
16  You should have received a copy of the GNU General Public License
17  along with this program; if not, write to the Free Software
18  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
19
20*/
21
22
23/*********************
24 *
25 * Chromatic polynomial
26 *
27 *********************/
28
29chromatic_polynomial(gr, x) :=
30if graph_size(gr)=0 then x^graph_order(gr)
31else  block(
32  [comp],
33  comp : connected_components(gr),
34  if length(comp)>1 then (
35    comp : map(lambda([u], chromatic_polynomial(induced_subgraph(u, gr), x)), comp),
36    expand(apply("*", comp)))
37  else if graph_size(gr)=graph_order(gr)-1 then x*(x-1)^graph_size(gr)
38  else if graph_size(gr)=graph_order(gr)*(graph_order(gr)-1)/2 then c_poly_complete(graph_order(gr),x)
39  else if min_degree(gr)[1]=2 and max_degree(gr)[1]=2 then c_poly_cycle[graph_order(gr)](x)
40  else block(
41    [g1, g2, u, v, p1, p2, e],
42    u : max_degree(gr)[2],
43    v : first(neighbors(u, gr)),
44    e : [min(u,v), max(u,v)],
45    g1 : copy_graph(gr),
46    g2 : copy_graph(gr),
47    remove_edge(e, g1),
48    contract_edge(e, g2),
49    p1 : chromatic_polynomial(g1, x),
50    p2 : chromatic_polynomial(g2, x),
51    p1-p2))$
52
53c_poly_cycle[1](x) := x$
54c_poly_cycle[3](x) := x*(x-1)*(x-2)$
55c_poly_cycle[n](x) := x*(x-1)^(n-1)-c_poly_cycle[n-1](x)$
56c_poly_complete(n,x) := apply("*", makelist(x-i, i, 0, n-1))$
57
58/*******************
59 *
60 *  Matching polynomial
61 *
62 *******************/
63
64matching_polynomial(gr, x) := (
65  if max_degree(gr)[1]<3 then
66    matching_polynomial_simple(gr, x)
67  else block(
68    [g1 : copy_graph(gr), g2 : copy_graph(gr), md, mv],
69    md : max_degree(g1),
70    mv : md[2],
71    md : md[1],
72    u : first(neighbors(mv, g1)),
73    remove_vertex(mv, g1),
74    remove_vertex(u, g1),
75    remove_edge([u, mv], g2),
76    matching_polynomial(g2, x) - matching_polynomial(g1, x)))$
77
78matching_polynomial_simple(gr, x) := block(
79  [conn, pol : 1, c, deg, u],
80  conn : connected_components(gr),
81  for c in conn do (
82    deg : apply(min,
83      args(map(lambda([u], vertex_degree(u, gr)), c))),
84    if deg=2 then pol : pol * cycle_poly(length(c), x)
85    else pol : pol * path_poly[length(c)](x)),
86  expand(pol))$
87
88cycle_poly(n, x) := path_poly[n](x) - path_poly[n-2](x)$
89path_poly[1](x) := x$
90path_poly[2](x) := x^2-1$
91path_poly[n](x) := x*path_poly[n-1](x) - path_poly[n-2](x)$
92
93/*******************
94 *
95 *  Tutte polynomial
96 *
97 *******************/
98
99tutte_polynomial(g, x, y) := block(
100  [non_bridge:false, tpzero:1, components: biconnected_components(g)],
101  /* Reduce to biconnected components */
102  if length(components)>1 then block(
103    [n_loops:0],
104    for v in vertices(g) do n_loops: n_loops+get_vertex_label(v, g, 0),
105    components: map(lambda([comp], induced_subgraph(comp, g)), components),
106    map(
107      lambda([gr],
108        for e in edges(gr) do
109        set_edge_weight(e, get_edge_weight(e, g), gr)),
110      components),
111    xreduce("*", map(lambda([gr], tutte_polynomial(gr, x, y)), components))*y^n_loops)
112  /* check for ``small'' graphs: */
113  /*   - point with loops */
114  else if graph_order(g)=1 then
115    y^get_vertex_label(first(vertices(g)), g, 0)
116  /*   - a multiedge with loops */
117  else if graph_order(g)=2 then block(
118    [e: first(edges(g))],
119    (x+xreduce("+", makelist(y^i, i, 1, get_edge_weight(e, g) - 1)))*
120      y^(get_vertex_label(e[1], g, 0) + get_vertex_label(e[2], g, 0)))
121  /* a cycle on n vertices */
122  else if first(max_degree(g))=2 and lmax(makelist(get_edge_weight(e, g), e, edges(g)))=1 then (
123    (y + xreduce("+", makelist(x^i, i, 1, graph_order(g)-1)))*
124      y^lsum(get_vertex_label(v, g, 0), v, vertices(g)))
125  /* The graph is biconnected - no edge is a bridge */
126  else (
127    /* choose the edge with one endpoint of minimum degree in the graph */
128    non_bridge: [second(min_degree(g))],
129    non_bridge: cons(first(neighbors(non_bridge[1], g)), non_bridge),
130    if non_bridge[1]>non_bridge[2] then non_bridge: reverse(non_bridge),
131    if non_bridge=false then block(
132      [tp:1],
133      tp: tp*x^graph_size(g),
134      for v in vertices(g) do
135        tp: tp*y^get_vertex_label(v, g, 0),
136      tp)
137    else block(
138      [g1: copy_graph(g), g2: copy_graph(g), mfactor:1, tp],
139      contract_edge(non_bridge, g2),
140      if get_edge_weight(non_bridge, g)=1 then (
141        remove_edge(non_bridge, g1))
142      else (
143        set_edge_weight(non_bridge, 1, g1),
144        mfactor: xreduce("+", makelist(y^i, i, 1, get_edge_weight(non_bridge, g)-1))),
145      for u in neighbors(non_bridge[2], g) do
146        if u#non_bridge[1] then
147          set_edge_weight([non_bridge[1], u],
148            get_edge_weight([non_bridge[1], u], g, 1, 0) +
149            get_edge_weight([non_bridge[2], u], g, 1, 0),
150            g2),
151      set_vertex_label(non_bridge[1],
152        get_vertex_label(non_bridge[1], g, 0) +
153        get_vertex_label(non_bridge[2], g, 0),
154        g2),
155      tp: tutte_polynomial(g1, x, y) + mfactor*tutte_polynomial(g2, x, y))))$
156
157flow_polynomial(g, x) := block(
158  [n: graph_order(g), m: graph_size(g)],
159  (-1)^(m-n+1)*ratexpand(psubst(['y=0, 'x=x], tutte_polynomial(g,'y,1-'x))))$
160
161rank_polynomial(g, x, y) := block(
162  [tp: tutte_polynomial(g, 'x, 'y), n:graph_order(g)],
163  ratexpand(x^(n-1)*psubst(['x=1+1/x, 'y=1+y], tp)))$