share/graphs/graph_polynomials.mac

   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., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA
  19
  20 */
  21
  22
  23 /*********************
  24  *
  25  * Chromatic polynomial
  26  *
  27  *********************/
  28
  29 chromatic_polynomial(gr, x) :=
  30 if graph_size(gr)=0 then x^graph_order(gr)
  31 else  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
  53 c_poly_cycle[1](x) := x$
  54 c_poly_cycle[3](x) := x*(x-1)*(x-2)$
  55 c_poly_cycle[n](x) := x*(x-1)^(n-1)-c_poly_cycle[n-1](x)$
  56 c_poly_complete(n,x) := apply("*", makelist(x-i, i, 0, n-1))$
  57
  58 /*******************
  59  *
  60  *  Matching polynomial
  61  *
  62  *******************/
  63
  64 matching_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
  78 matching_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
  88 cycle_poly(n, x) := path_poly[n](x) - path_poly[n-2](x)$
  89 path_poly[1](x) := x$
  90 path_poly[2](x) := x^2-1$
  91 path_poly[n](x) := x*path_poly[n-1](x) - path_poly[n-2](x)$
  92
  93 /*******************
  94  *
  95  *  Tutte polynomial
  96  *
  97  *******************/
  98
  99 tutte_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
 157 flow_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
 161 rank_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)))$