# This is the calculation of the global rigidity of a graph G for a
# generic configuration in d-dimensional space.  p is a prime number
# used to replace the generic calculations.
> with(networks):
> G:=complete(5,5):
> delete(edges({5,6},G),G):
> addvertex({11},G):
> addedge({{11,5},{11,6},{11,7},{11,4}},G):
> d:=3:
> p:=1009:
# This is the graph that will be used for the calculation.  Various
# graph operations are available to define the graph desired.  The
# vertices of G are labeled 1, .. , nops(vertices(G)), and the edges
# are labled e1, .., nops(edges(G)).  
> with(linalg):
# Warning, the assigned names charpoly and rank now have a global
# binding
# Warning, the protected names norm and trace have been redefined and
# unprotected
# The following procedure 'Stress' forms an n-by-n matrix with 4
# entries of 1 corresponding to the columns i and j.
> Stress := proc( i,j,n )
>         description "Forms an n-by-n stress matrix with a stress of 1
> between node i and node j";
>         A:=Matrix(n) : A[i,i]:=1 : A[i,j]:=-1 : A[j,i]:=-1 :
> A[j,j]:=1 : A
> end proc:
# Warning, `A` is implicitly declared local to procedure `Stress`
# The following defines a list W of variables that will correspond to
# the stresses on the edges of the graph G.  The default is a generic
# sequence of variables, but this is changed later with the subs
# command. 
> W:=[seq(w[i],i=1..nops(edges(G)))]:
# The following 'edgs' computes the vertex pairs that comprise the
# edges of the graph G.
> edgs := proc( G )
>         description"Computes the node pairs of edges of the graph G";
>         seq(ends(G)[i],i=1..nops(ends(G)));
> end proc:
# The following 'StrsMat' is the stress matrix corresponding to the
# stresses given by W. (nops(W) must be nops(edges(G)).) 
> StrsMat:= proc( G,w )
>         description"Computes the general sress matrix of the graph G
> with stresses w[1]..w[nops(edges(G))";
>        
> simplify(add(w[i]*Stress(edgs(G)[i][1],edgs(G)[i][2],nops(vertices(G))
> ),i=1..nops(edges(G)))) ;
> end proc: 
# The following is a configuration matrix P consisting of vertices(G)
# column vectors in d-dimensional space. The default is that P is
# random with variables x[i,j] mod p.  
> with(LinearAlgebra:-Modular): P :=
> Random(p,d,nops(vertices(G)),integer[]):
# The following matrix X is the product of the configuration matrix P
# and the stress matrix.
> with(LinearAlgebra): X:=simplify(Multiply(P,StrsMat(G,W))):
# Warning, these names have been rebound: Adjoint, BackwardSubstitute,
# Basis, CharacteristicPolynomial, Copy, Determinant,
# ForwardSubstitute, GramSchmidt, LUDecomposition, Multiply, Transpose
# The following is a list of equations corresponding to the conditions
# that each entry of  X be 0.  These are the equilibrium equations.
> eq:=[seq(seq(X[i,j], i=1..d), j=1..nops(vertices(G)))]:
# The following is the set of equations that are the consequence of the
# equilibrium equations on the stress.  This solution is then
# substituted back into the stress vector W, with a choice of one (or
# more) stress being 1 to eliminate the generic variables.  You the
# user must see what variable(s) Maple uses as the independent
# variable(s) and then substitute a convenient (random?) value for it.
> S:=solve({op(eq)},{op(W)}) mod p;
S := {w[6] = w[6], w[1] = 1004 w[6], w[2] = 192 w[6], w[3] = 177 w[6], 

  w[15] = 779 w[6], w[14] = 68 w[6], w[27] = 347 w[6], w[13] = 792 w[6], 

  w[20] = 296 w[6], w[23] = 281 w[6], w[19] = 160 w[6], w[18] = 146 w[6], 

  w[11] = 274 w[6], w[9] = 746 w[6], w[12] = 954 w[6], w[28] = 130 w[6], 

  w[25] = 931 w[6], w[26] = 465 w[6], w[24] = 265 w[6], w[4] = 126 w[6], 

  w[7] = 197 w[6], w[5] = 441 w[6], w[16] = 947 w[6], w[17] = 626 w[6], 

  w[22] = 632 w[6], w[10] = 995 w[6], w[8] = 750 w[6], w[21] = 179 w[6]}
> W:=subs(S,W);
W := [1004 w[6], 192 w[6], 177 w[6], 126 w[6], 441 w[6], w[6], 197 w[6], 

  750 w[6], 746 w[6], 995 w[6], 274 w[6], 954 w[6], 792 w[6], 68 w[6], 

  779 w[6], 947 w[6], 626 w[6], 146 w[6], 160 w[6], 296 w[6], 179 w[6], 

  632 w[6], 281 w[6], 265 w[6], 931 w[6], 465 w[6], 347 w[6], 130 w[6]]
# Here is where the actual substitution, which you must choose, is made.
> W:=subs(w[6]=1,W);
W := [1004, 192, 177, 126, 441, 1, 197, 750, 746, 995, 274, 954, 792, 68, 779, 

  947, 626, 146, 160, 296, 179, 632, 281, 265, 931, 465, 347, 130]
# This new value for W is put back into the stress matrix to get the
# stress matrix SM mod p.  
> SM:=StrsMat(G,W);
                       [ 11 x 11 Matrix       ]
                       [ Data Type: anything  ]
                 SM := [ Storage: rectangular ]
                       [ Order: Fortran_order ]
# This converts the stress matrix mod p to all integer form which is
# friendly for the mod p rank calculation.
> B := Mod(p,SM, integer[]);

                      [ 11 x 11 Matrix        ]
                      [ Data Type: integer[4] ]
                 B := [ Storage: rectangular  ]
                      [ Order: C_order        ]
# These next commands put the machine in a mode to do the mod p rank
# calculation.   Then RowReduce does the rank calculation.
> with(LinearAlgebra:-Modular):
# Warning, these names have been rebound: Adjoint, BackwardSubstitute,
# Basis, CharacteristicPolynomial, Copy, Determinant,
# ForwardSubstitute, LUDecomposition, Multiply, Transpose
> unprotect(det):
> unprotect(linalg:-rank)# ;
# Warning, inserted missing semicolon at end of statement,
# ...linalg:-rank);
> C := Copy(p,B):
> RowReduce(p,C,nops(vertices(G)),nops(vertices(G)),nops(vertices(G)),'d
> et',0,'rank',0,0,true):
> rank;
                               7
> restart;