(* :Title: Combinatorica *)

(* :Author: Steven S. Skiena *)

(* :Summary:
This package contains all the programs from the book, Implementing
Discrete Mathematics: Combinatorics and Graph Theory with Mathematica,
by Steven S. Skiena, Addison-Wesley Publishing Co., Advanced Book Program,
350 Bridge Parkway, Redwood City CA 94065.  ISBN 0-201-50943-1.
For ordering information, call 1-800-447-2226.
*)



(* :Discussion:
The programs from the book "Implementing Discrete Mathematics: Combinatorics
and Graph Theory with Mathematica" are available by anonymous ftp.cs.sunysb.edu
in the pub/Combinatorica directory.  



Any comments, bug reports, or requests to get on the Combinatorica
mailing list should be forwarded to:

	Steven Skiena
	Department of Computer Science
	State University of New York
	Stony Brook, NY 11794

	skiena@sbcs.sunysb.edu

	(516)-632-9026 / 8470
*)
(* :Context: DiscreteMath`Combinatorica`  
*)
(* :Package Version: 1.2
*)
(* :Copyright: Copyright 1990-1999 by Steven S. Skiena

This package may be copied in its entirety for nonprofit purposes only.
Sale, other than for the direct cost of the media, is prohibited.  This
copyright notice must accompany all copies.

The author, Wolfram Research, and Addison-Wesley Publishing Company,
Inc. make no representations, express or implied, with respect to this
documentation, or the software it describes and contains, including
without limitations, any implied warranties of mechantability or fitness
for a particular purpose, all of which are expressly disclaimed.  The
author, Wolfram Research, or Addison-Wesley, their licensees,
distributors and dealers shall in no event be liable for any indirect,
incidental, or consequential damages.
*)

(* :History:
	Version 1.2 by Steven S. Skiena, March 1997. 
	Version 1.1 modification by ECM, March 1996.  
		Replaced K with CompleteGraph because K is now the
			default generic name for the summation index in
			symbolic sum.
		Added CombinatorialFunctions.m and Permutations.m to
			BeginPackage, and commented out CatalanNumber,
			PermutationQ, ToCycles, FromCycles, and
			RandomPermutation so there would be no shadowing of
			symbols among the DiscreteMath packages.	
		Replaced old BinarySearch with new code by Paul Abbott
			correctly implementing binary search.
        Version 1.0 by Steven S. Skiena, April 1995.
        Version .9 by Steven S. Skiena, February 1992.
	Version .8 by Steven S. Skiena, July 1991.
	Version .7 by Steven S. Skiena, January 1991. 
	Version .6 by Steven S. Skiena, June 1990.
*)

(* :Keywords:
	adjacency, automorphism, chromatic, clique, coloring,
	combination, composition, connected components, connectivity, cycle,
	de Bruijn, degree, derangement, Dijkstra, Durfee,
	embedding, equivalence, Eulerian, Ferrers,
	geodesic, graph, Gray code, group, Hamiltonian cycle, Harary, Hasse,
	heap, hypercube, interval, inversion, involution, isomorphism,
	Josephus, network,
	partition, perfect, permutation, planar graph, Polya, pseudograph,
	self-loop, sequence, signature, simple, spanning tree,
	stable marriage, star, Stirling,
	transitive closure, traveling salesman tour, tree, Turan,
	vertex cover, wheel, Young tableau
*)
(* :Source:
	Steven Skiena: Implementing Discrete Mathematics: Combinatorics
			and Graph Theory with Mathematica, 
			Addison-Wesley Publishing Co.
*)
(* :Mathematica Version: 3.0
*)

BeginPackage["DiscreteMath`Combinatorica`", 
		"DiscreteMath`CombinatorialFunctions`",
 (* needed for CatalanNumber *)
		"DiscreteMath`Permutations`"
 (* needed for PermutationQ, ToCycles, FromCycles, and RandomPermutation *)
		]



Unprotect[
AcyclicQ, AddEdge, AddVertex, AllPairsShortestPath, ArticulationVertices,
Automorphisms, Backtrack, BiconnectedComponents, BiconnectedQ,
BinarySearch, BinarySubsets, BipartiteMatching, BipartiteQ,
BreadthFirstTraversal, Bridges, CartesianProduct, (* CatalanNumber, *)
ChangeEdges, ChangeVertices, ChromaticNumber, ChromaticPolynomial,
CirculantGraph, CircularVertices, CliqueQ, CodeToLabeledTree, Cofactor,
CompleteGraph,
CompleteQ, Compositions, ConnectedComponents, ConnectedQ, ConstructTableau,
Contract, CostOfPath, Cycle, DeBruijnSequence, DegreeSequence,
DeleteCycle, DeleteEdge, DeleteFromTableau, DeleteVertex, DepthFirstTraversal,
DerangementQ, Derangements, Diameter, Dijkstra, DilateVertices,
DistinctPermutations, Distribution, DurfeeSquare, Eccentricity,
EdgeChromaticNumber, EdgeColoring, EdgeConnectivity, Edges, Element,
EmptyGraph, EmptyQ, EncroachingListSet, EquivalenceClasses, 
EquivalenceRelationQ, Equivalences, EulerianCycle, EulerianQ, Eulerian,
ExactRandomGraph, ExpandGraph, ExtractCycles, FerrersDiagram, FindCycle,
FindSet, FirstLexicographicTableau, FromAdjacencyLists, (* FromCycles, *)
FromInversionVector, FromOrderedPairs, FromUnorderedPairs, FunctionalGraph,
Girth, GraphCenter, GraphComplement, GraphDifference, GraphIntersection,
GraphJoin, GraphPower, GraphProduct, GraphSum, GraphUnion, GraphicQ,
GrayCode, GridGraph, HamiltonianCycle, HamiltonianQ, Harary,
HasseDiagram, HeapSort, Heapify, HideCycles, Hypercube, IdenticalQ,
IncidenceMatrix, IndependentSetQ, Index, InduceSubgraph, InitializeUnionFind,
InsertIntoTableau, IntervalGraph, InversePermutation, Inversions,
InvolutionQ, IsomorphicQ, IsomorphismQ, Isomorphism, Josephus, KSubsets,
LabeledTreeToCode, LastLexicographicTableau,
LexicographicPermutations, LexicographicSubsets, LineGraph,
LongestIncreasingSubsequence, M, MakeGraph, MakeSimple, MakeUndirected,
MaximalMatching, MaximumAntichain, MaximumClique, MaximumIndependentSet,
MaximumSpanningTree, MinimumChainPartition, MinimumChangePermutations,
MinimumSpanningTree, MinimumVertexCover, MultiplicationTable,
NetworkFlowEdges, NetworkFlow, NextComposition, NextKSubset,
NextPartition, NextPermutation, NextSubset, NextTableau, NormalizeVertices,
NthPair, NthPermutation, NthSubset, NumberOfCompositions,
NumberOfDerangements, NumberOfInvolutions, NumberOfPartitions,
NumberOfPermutationsByCycles, NumberOfSpanningTrees, NumberOfTableaux,
OrientGraph, PartialOrderQ, PartitionQ, Partitions, PathConditionGraph,
Path, PerfectQ, PermutationGroupQ, (* PermutationQ, *) Permute, PlanarQ,
PointsAndLines, Polya, PseudographQ, RadialEmbedding, Radius,
RandomComposition, RandomGraph, RandomHeap, RandomKSubset,
RandomPartition, RandomPermutation1, RandomPermutation2,
(* RandomPermutation, *)
RandomSubset, RandomTableau, RandomTree, RandomVertices, RankGraph,
RankPermutation, RankSubset, RankedEmbedding, ReadGraph, 
RealizeDegreeSequence, RegularGraph, RegularQ, RemoveSelfLoops,
RevealCycles, RootedEmbedding, RotateVertices, Runs, SamenessRelation,
SelectionSort, SelfComplementaryQ, ShakeGraph, ShortestPathSpanningTree,
ShortestPath, ShowGraph, ShowLabeledGraph, SignaturePermutation,
SimpleQ, Spectrum, SpringEmbedding, StableMarriage, Star,
StirlingFirst, StirlingSecond, Strings, StronglyConnectedComponents,
Subsets, TableauClasses, TableauQ, TableauxToPermutation, Tableaux,
ToAdjacencyLists, (* ToCycles, *) ToInversionVector, ToOrderedPairs,
ToUnorderedPairs, TopologicalSort, TransitiveClosure, TransitiveQ,
TransitiveReduction, TranslateVertices, TransposePartition, TransposeTableau,
TravelingSalesmanBounds, TravelingSalesman, TreeQ, TriangleInequalityQ,
Turan, TwoColoring, UndirectedQ, UnionSet, UnweightedQ,
V, VertexColoring, VertexConnectivity, VertexCoverQ, Vertices,
WeaklyConnectedComponents, Wheel, WriteGraph, DilworthGraph
]



$NewMessage[ K, "usage"] (* reset the usage of K to the System usage *)
K::usage = StringJoin[ K::usage, "  The use of K to create a complete
graph is obsolete.  Use CompleteGraph to create a complete graph."]

$NewMessage[ All, "usage"] (* reset the usage of All to the System usage *)
All::usage = StringJoin[ All::usage, "   All is also an option to 
certain Combinatorica functions specifying that all solutions should be
returned, instead of just the first one."]

AcyclicQ::usage = "AcyclicQ[g] yields True if graph g is acyclic. AcyclicQ[g,
Directed] yields True if g is a directed acyclic graph."

AddEdge::usage = "AddEdge[g, {x,y}] gives graph g with a new undirected edge {x,y}, while AddEdge[g, {x,y}, Directed] gives graph g with a new directed edge {x,y}."

AddVertex::usage = "AddVertex[g] adds a disconnected vertex to graph g."

AllPairsShortestPath::usage = "AllPairsShortestPath[g] gives a matrix, where the (i,j)th entry is the length of the shortest path in g between vertices i and j."

ArticulationVertices::usage = "ArticulationVertices[g] gives a list of all articulation vertices in graph g, vertices whose removal will disconnect the graph."

Automorphisms::usage = "Automorphisms[g] finds the automorphism group of a graph g, the set of isomorphisms of g with itself."

Backtrack::usage = "Backtrack[s, partialQ, solutionQ] performs a backtrack search of the state space s, expanding a partial solution so long as partialQ is True and returning the first complete solution, as identified by solutionQ."

BiconnectedComponents::usage = "BiconnectedComponents[g] gives a list of the biconnected components of graph g."

BiconnectedQ::usage = "BiconnectedQ[g] yields True if graph g is biconnected."

BinarySearch::usage =
"BinarySearch[l, k] searches sorted list l for key k and gives the position
of l containing k.  BinarySearch[l, k, f] gives the position using a 
function f to extract the key from an element of l."

BinarySubsets::usage = "BinarySubsets[l] gives all subsets of l ordered according to the binary string defining each subset."

BipartiteMatching::usage = "BipartiteMatching[g] gives the list of edges associated with a maximum matching in bipartite graph g."

BipartiteQ::usage = "BipartiteQ[g] yields True if graph g is bipartite."

BreadthFirstTraversal::usage = "BreadthFirstTraversal[g,v] performs a breadth-first traversal of graph g starting from vertex v, and gives a list of vertices in the order in which they were encountered."

Bridges::usage = "Bridges[g] gives a list of the bridges of graph g, the edges whose removal disconnects the graph."

CartesianProduct::usage = "CartesianProduct[l1,l2] gives the Cartesian product of lists l1 and l2."

(*
CatalanNumber::usage = "CatalanNumber[n] gives the nth Catalan number, for a positive integer n."
*)

ChangeEdges::usage =
"ChangeEdges[g,e] gives a graph with the adjacency matrix e and the embedding of
graph g."

ChangeVertices::usage =
"ChangeVertices[g,v] gives a graph with the adjacency matrix of graph g and the
list v as its embedding."

ChromaticNumber::usage = "ChromaticNumber[g] gives the chromatic number of the graph, the fewest number of colors necessary to color the graph."

ChromaticPolynomial::usage = "ChromaticPolynomial[g,z] gives the chromatic polynomial P(z) of graph g, which counts the number of ways to color g with at most z colors."

CirculantGraph::usage =
"CirculantGraph[n,l] gives a circulant graph on n vertices, meaning the ith
vertex is adjacent to the (i+j)th and (i-j)th vertex, for each j in list l."

CircularVertices::usage =
"CircularVertices[n] gives a list of n points equally spaced on a circle."

CliqueQ::usage = "CliqueQ[g,c] yields True if the list of vertices c defines a clique in graph g."

CodeToLabeledTree::usage =
"CodeToLabeledTree[l] gives the unique labeled tree on n vertices from the
Prufer code l, which consists of a list of n-2 integers from 1 to n."

Cofactor::usage = "Cofactor[m,{i,j}] gives the (i,j)th cofactor of matrix m."

CompleteGraph::usage =
"CompleteGraph[n] creates a complete graph on n vertices. CompleteGraph[a,b,c,...,k] creates a complete k-partite graph of the prescribed shape."

CompleteQ::usage = "CompleteQ[g] yields True if graph g is complete."

Compositions::usage = "Compositions[n,k] gives a list of all compositions of integer n into k parts."

ConnectedComponents::usage = "ConnectedComponents[g] gives the vertices of graph g partitioned into connected components."

ConnectedQ::usage = "ConnectedQ[g] yields True if undirected graph g is connected. ConnectedQ[g, Directed] and ConnectedQ[g, Undirected] yields True if g is strongly or weakly connected, respectively."

ConstructTableau::usage = "ConstructTableau[p] performs the bumping algorithm repeatedly on each element of permutation p, resulting in a distinct Young tableau."

Contract::usage = "Contract[g, {x,y}] gives the graph resulting from contracting edge {x,y} of graph g."

CostOfPath::usage = "CostOfPath[g, p] sums up the weights of the edges in graph g defined by the path p."

Cycle::usage =
"Cycle[n] gives the cycle on n vertices, a 2-regular connected graph."

DeBruijnSequence::usage =
"DeBruijnSequence[a, n] gives a de Bruijn sequence on the alphabet described by
list a, the shortest sequence such that every string of length n on a occurs as
a contiguous subrange of the sequence."

DegreeSequence::usage = "DegreeSequence[g] gives the sorted degree sequence of graph g."

DeleteCycle::usage = "DeleteCycle[g, c] deletes undirected cycle c from graph g. DeleteCycle[g,c,Directed] deletes directed cycle c from graph g."

DeleteEdge::usage = "DeleteEdge[g, {x,y}] gives graph g minus undirected edge {x,y}, while DeleteEdge[g, {x,y}, Directed] gives graph g minus directed edge {x,y}."

DeleteFromTableau::usage = "DeleteFromTableau[t, r] deletes the last element of row r from Young tableaux t."

DeleteVertex::usage = "DeleteVertex[g, v] deletes vertex v from graph g."

DepthFirstTraversal::usage = "DepthFirstTraversal[g, v] performs a depth-first traversal of graph g starting from vertex v, and gives a list of vertices in the order in which they were encountered."

DerangementQ::usage = "DerangementQ[p] tests whether permutation p is a derangement, a permutation without a fixed point."

Derangements::usage =
"Derangements[p] gives all derangements of permutation p."

Diameter::usage = "Diameter[g] gives the diameter of graph g, the length of the longest shortest path between two vertices of g."

Dijkstra::usage = "Dijkstra[g, v] gives the shortest path spanning tree and associated distances from vertex v of graph g."

DilateVertices::usage = "DilateVertices[v, d] multiplies each coordinate of each vertex position in list l by d, thus dilating the embedding."

Directed::usage = "Directed is an option to inform certain functions that the graph is directed."

DistinctPermutations::usage = "DistinctPermutations[l] gives all permutations of the multiset described by list l."

Distribution::usage = "Distribution[l, set] lists the frequency of occurrence of each element of set in list l."

DurfeeSquare::usage = "DurfeeSquare[p] gives the number of rows involved in the Durfee square of partition p, the side of the largest-sized square contained within the Ferrers diagram of p."

Eccentricity::usage = "Eccentricity[g] gives the eccentricity of each vertex v of graph g, the length of the longest shortest path from v."

Edge::usage = "Edge is an option to inform certain functions to work with edges instead of vertices."

EdgeChromaticNumber::usage = "EdgeChromaticNumber[g] gives the fewest number of colors necessary to color each edge of graph g, so that no two edges incident on the same vertex have the same color."

EdgeColoring::usage = "EdgeColoring[g] uses Brelaz's heuristic to find a good, but not necessarily minimal, edge coloring of graph g."

EdgeConnectivity::usage = "EdgeConnectivity[g] gives the minimum number of edges whose deletion from graph g disconnects it."

Edges::usage = "Edges[g] gives the adjacency matrix of graph g."

Element::usage = "Element[a, l] gives the lth element of nested list a, where l is a list of indices"

EmptyGraph::usage = "EmptyGraph[n] generates an empty graph on n vertices."

EmptyQ::usage = "EmptyQ[g] yields True if graph g contains no edges."

EncroachingListSet::usage =
"EncroachingListSet[p] gives the encroaching list set associated with
permutation p."

EquivalenceClasses::usage = "EquivalenceClasses[r] identifies the equivalence classes among the elements of matrix r."

EquivalenceRelationQ::usage = "EquivalenceRelationQ[r] yields True if the matrix r defines an equivalence relation. EquivalenceRelationQ[g] tests whether the adjacency matrix of graph g defines an equivalence relation."

Equivalences::usage = "Equivalences[g, h] lists the vertex equivalence classes between graphs g and h defined by the all-pairs shortest path heuristic."

EulerianCycle::usage = "EulerianCycle[g] finds an Eulerian circuit of undirected graph g if one exists. EulerianCycle[g, Directed] finds an Eulerian circuit of directed graph g if one exists."

EulerianQ::usage = "EulerianQ[g] yields True if the graph g is Eulerian, meaning there exists a tour which includes each edge exactly once. EulerianQ[g, Directed] yields True if the directed graph g is Eulerian."

Eulerian::usage = "Eulerian[n, k] gives the number of permutations of length n with k runs."

ExactRandomGraph::usage =
"ExactRandomGraph[n, e] gives a random labeled graph of exactly e edges and n
vertices."

ExpandGraph::usage = "ExpandGraph[g, n] expands graph g to n vertices by adding disconnected vertices."

ExtractCycles::usage = "ExtractCycles[g] gives a list of edge disjoint cycles in graph g."

FerrersDiagram::usage = "FerrersDiagram[p] draws a Ferrers diagram of integer partition p."

FindCycle::usage = "FindCycle[g] finds a list of vertices which define an undirected cycle in graph g. FindCycle[g, Directed] finds a directed cycle in graph g."

FindSet::usage = "FindSet[n, s] gives the root of the set containing n in union-find data structure s."

FirstLexicographicTableau::usage =
"FirstLexicographicTableau[p] gives the first Young tableau with shape
described by partition p."

FromAdjacencyLists::usage =
"FromAdjacencyLists[l] gives an adjacency matrix representation for a graph with
adjacency lists l, using a circular embedding. FromAdjacencyLists[l, v] uses
v as the embedding for the resulting graph."

(*
FromCycles::usage = "FromCycles[c] restores a cycle structure c to the original permutation."
*)

FromInversionVector::usage =
"FromInversionVector[v] gives the unique permutation with inversion vector v."

FromOrderedPairs::usage =
"FromOrderedPairs[l] gives an adjacency matrix representation from a list of
ordered pairs l, using a circular embedding. FromOrderedPairs[l, v] uses v as
the embedding for the resulting graph."

FromUnorderedPairs::usage =
"FromUnorderedPairs[l] gives an adjacency matrix representation from a list of
unordered pairs l, using a circular embedding. FromUnorderedPairs[l, v] uses
v as the embedding for the resulting graph."

FunctionalGraph::usage =
"FunctionalGraph[f,n] gives the functional digraph on n vertices defined by
integer function f."

Girth::usage = "Girth[g] gives the length of the shortest cycle in unweighted graph g."

Graph::usage = "Graph[g,v] represents a graph object where g is an adjacency matrix and v is a list of vertices."

GraphCenter::usage = "GraphCenter[g] gives a list of the vertices of graph g with minimum eccentricity."

GraphComplement::usage = "GraphComplement[g] gives the complement of graph g."

GraphDifference::usage =
"GraphDifference[g, h] gives the graph resulting from subtracting the adjacency
 matrix of graph g from that of graph h."

GraphIntersection::usage = "GraphIntersection[g, h] gives the graph defined by
the edges which are in both graph g and graph h."

GraphJoin::usage = "GraphJoin[g, h] gives the join of graphs g and h."

GraphPower::usage = "GraphPower[g, k] gives the kth power of graph g, meaning there is an edge between any pair of vertices of g with a path between them of length at most k."

GraphProduct::usage = "GraphProduct[g, h] gives the product of graphs g and h."

GraphSum::usage = "GraphSum[g, h] gives the graph resulting from adding the
adjacency matrices of graphs g and h."

GraphUnion::usage = "GraphUnion[g, h] gives the union of graphs g and h.
GraphUnion[n, g] gives n copies of graph g, where n is an integer."

GraphicQ::usage = "GraphicQ[s] yields True if the list of integers s is graphic, and thus represents a degree sequence of some graph."

GrayCode::usage = "GrayCode[l] gives a binary reflected Gray code on set l."

GridGraph::usage = "GridGraph[n, m] gives an n*m grid graph, the product of
paths on n and m vertices."

HamiltonianCycle::usage = "HamiltonianCycle[g] finds a Hamiltonian cycle in graph g if one exists. HamiltonianCycle[g, All] gives all Hamiltonian cycles of graph g."

HamiltonianQ::usage = "HamiltonianQ[g] yields True if there exists a Hamiltonian cycle in graph g, in other words, if there exists a cycle which visits each vertex exactly once."

Harary::usage = 
"Harary[k, n] gives the minimal k-connected graph on n vertices."

HasseDiagram::usage = 
"HasseDiagram[g] gives a Hasse diagram of the relation defined by directed
acyclic graph g."

HeapSort::usage = "HeapSort[l] performs a heap sort on the items of list l."

Heapify::usage = "Heapify[p] builds a heap from permutation p."

HideCycles::usage = "HideCycles[c] canonically encodes the cycle structure c into a unique permutation."

Hypercube::usage = "Hypercube[n] gives an n-dimensional hypercube."

IdenticalQ::usage = "IdenticalQ[g, h] yields True if graphs g and h have identical adjacency matrices."

IncidenceMatrix::usage = "IncidenceMatrix[g] gives the (0,1) incidence matrix of graph g, which has a row for each vertex and column for each edge and (v,e)=1 if and only if vertex v is incident upon edge e."

IndependentSetQ::usage = "IndependentSetQ[g, i] yields True if the vertices in list i define an independent set in graph g."

Index::usage = "Index[p] gives the index of permutation p, the sum of all subscripts j such that p[j] is greater than p[j+1]."

InduceSubgraph::usage = 
"InduceSubgraph[g, s] gives the subgraph of graph g induced by the list of
vertices s."

InitializeUnionFind::usage = "InitializeUnionFind[n] initializes a union-find data structure for n elements."

InsertIntoTableau::usage = "InsertIntoTableau[e, t] inserts integer e into Young tableau t using the bumping algorithm."

IntervalGraph::usage =
"IntervalGraph[l] gives the interval graph defined by the list of intervals l."

InversePermutation::usage = "InversePermutation[p] yields the multiplicative inverse of permutation p."

Inversions::usage = "Inversions[p] counts the number of inversions in permutation p."

InvolutionQ::usage = "InvolutionQ[p] yields True if permutation p is its own inverse."

IsomorphicQ::usage = "IsomorphicQ[g, h] yields True if graphs g and h are isomorphic."

IsomorphismQ::usage = "IsomorphismQ[g, h, p] tests if permutation p defines an isomorphism between graphs g and h."

Isomorphism::usage = "Isomorphism[g, h] gives an isomorphism between graphs g and h if one exists."

Josephus::usage = "Josephus[n, m] generates the inverse of the permutation defined by executing every mth member in a circle of n members."

KSubsets::usage = "KSubsets[l, k] gives all subsets of set l containing exactly k elements, ordered lexicographically."

LabeledTreeToCode::usage = "LabeledTreeToCode[g] reduces the tree g to its Prufer code."

LastLexicographicTableau::usage =
"LastLexicographicTableau[p] gives the last Young tableau with shape described
by partition p."

LexicographicPermutations::usage = 
"LexicographicPermutations[l] gives all permutations of list l in lexicographic
order."

LexicographicSubsets::usage = "LexicographicSubsets[l] gives all subsets of set l in lexicographic order."

LineGraph::usage = "LineGraph[g] gives the line graph of graph g."

LongestIncreasingSubsequence::usage = "LongestIncreasingSubsequence[p] finds the longest increasing scattered subsequence of permutation p."

M::usage = "M[g] gives the number of edges in undirected graph g."

MakeGraph::usage = 
"MakeGraph[v, f] gives the binary relation defined by function f on all pairs of
elements of list v."

MakeSimple::usage = "MakeSimple[g] gives an undirected, unweighted graph derived from directed graph g."

MakeUndirected::usage = "MakeUndirected[g] gives a graph with an undirected edge for each directed edge of graph g."

MaximalMatching::usage = "MaximalMatching[g] gives the list of edges associated with a maximal matching of graph g."

MaximumAntichain::usage = "MaximumAntichain[g] gives a largest set of unrelated vertices in partial order g."

MaximumClique::usage = "MaximumClique[g] finds the largest clique in graph g."

MaximumIndependentSet::usage = "MaximumIndependentSet[g] finds the largest independent set of graph g."

MaximumSpanningTree::usage = "MaximumSpanningTree[g] uses Kruskal's algorithm to find a maximum spanning tree of graph g."

MinimumChainPartition::usage = "MinimumChainPartition[g] partitions partial order g into a minimum number of chains."

MinimumChangePermutations::usage = "MinimumChangePermutations[l] gives all
permutations of list l such that adjacent permutations differ by only one
transposition."

MinimumSpanningTree::usage = "MinimumSpanningTree[g] uses Kruskal's algorithm to find a minimum spanning tree of graph g."

MinimumVertexCover::usage = "MinimumVertexCover[g] finds the minimum vertex cover of graph g."

MultiplicationTable::usage = "MultiplicationTable[l, f] gives the complete
transition table defined by the binary relation function f on the elements of
list l."

NetworkFlowEdges::usage = "NetworkFlowEdges[g, source, sink] gives the adjacency matrix showing the distribution of the maximum flow from source to sink in graph g."

NetworkFlow::usage = "NetworkFlow[g, source, sink] finds the maximum flow through directed graph g from source to sink."

NextComposition::usage = "NextComposition[l] gives the integer composition
which follows l in a canonical order."

NextKSubset::usage = "NextKSubset[l, s] gives the k-subset of list l which appears after k-subsets s in lexicographic order."

NextPartition::usage = "NextPartition[p] gives the integer partition following p in reverse lexicographic order."

NextPermutation::usage = "NextPermutation[p] gives the permutation following p in lexicographic order."

NextSubset::usage = "NextSubset[l, s] gives the subset of l following subset
s in canonical order."

NextTableau::usage = "NextTableau[t] gives the tableau of shape t which follows t in lexicographic order."

NormalizeVertices::usage = "NormalizeVertices[v] gives a list of vertices with the same structure as v but with all coordinates of all points between 0 and 1."

NthPair::usage = "NthPair[n] gives the nth unordered pair of positive integers, when sequenced to minimize the size of the larger integer."

NthPermutation::usage = "NthPermutation[n, l] gives the nth lexicographic permutation of list l."

NthSubset::usage = "NthSubset[n, l] gives the nth subset of list l in canonical order."

NumberOfCompositions::usage = "NumberOfCompositions[n, k] counts the number of distinct compositions of integer n into k parts."

NumberOfDerangements::usage = "NumberOfDerangements[n] counts the derangements on n elements, the permutations without any fixed points."

NumberOfInvolutions::usage = "NumberOfInvolutions[n] counts the number of involutions on n elements."

NumberOfPartitions::usage = "NumberOfPartitions[n] counts the number of distinct integer partitions of n."

NumberOfPermutationsByCycles::usage = "NumberOfPermutationsByCycles[n, m] gives the number of permutations of length n with exactly m cycles."

NumberOfSpanningTrees::usage = "NumberOfSpanningTrees[g] gives the number of distinct labeled spanning trees of graph g."

NumberOfTableaux::usage = "NumberOfTableaux[p] uses the hook length formula to count the number of Young tableaux with shape defined by partition p."

OrientGraph::usage = "OrientGraph[g] assigns a direction to each edge of a bridgeless, undirected graph g, so that the graph is strongly connected."

PartialOrderQ::usage = "PartialOrderQ[g] yields True if the binary relation defined by the adjacency matrix of graph g is a partial order, meaning it is transitive, reflexive, and anti-symmetric."

PartitionQ::usage = "PartitionQ[p] yields True if p is an integer partition."

Partitions::usage =
"Partitions[n] gives all partitions of integer n in reverse lexicographic order."

PathConditionGraph::usage = "PathConditionGraph[g] replaces each non-edge of a graph by an infinite cost, so shortest path algorithms work correctly."

Path::usage = "Path[n] gives a tree consisting only of a path on n vertices."

PerfectQ::usage = "PerfectQ[g] yields True if g is a perfect graph, meaning that for every induced subgraph of g the size of the largest clique equals the chromatic number."

PermutationGroupQ::usage = "PermutationGroupQ[l] yields True if the list of permutations l forms a permutation group."

(*
PermutationQ::usage = "PermutationQ[p] yields True if p represents a permutation and False otherwise."
*)

Permute::usage = "Permute[l, p] permutes list l according to permutation p."

PlanarQ::usage = "PlanarQ[g] yields True if graph g is planar, meaning it can be drawn in the plane so no two edges cross."

PointsAndLines::usage = 
"PointsAndLines[g] gives a partial graphics representation of a graph g."

Polya::usage = "Polya[g, m] gives the polynomial giving the number of colorings, with m colors, of a structure defined by the permutation group g."

PseudographQ::usage = "PseudographQ[g] yields True if graph g is a pseudograph, meaning it contains self-loops."

RadialEmbedding::usage = 
"RadialEmbedding[g] gives a radial embedding of graph g, radiating from the
center of the graph."

Radius::usage = "Radius[g] gives the radius of graph g, the minimum eccentricity of any vertex of g."

RandomComposition::usage = "RandomComposition[n, k] gives a random composition
of integer n into k parts."

RandomGraph::usage =
"RandomGraph[n, p] gives a random labeled graph on n vertices with an edge
probability of p and unit edge weights.  RandomGraph[n, p, {l,h}] gives
a random labeled graph with edge weights of integers drawn uniformly at random
from the range (l,h).  RandomGraph[n, p, {l,h}, Directed] similarly gives
a random directed graph."  

RandomHeap::usage = "RandomHeap[n] gives a random heap on n elements."

RandomKSubset::usage = "RandomKSubset[l, k] gives a random subset of set l with exactly k elements."

RandomPartition::usage =
"RandomPartition[n] gives a random partition of integer n."

RandomPermutation1::usage = "RandomPermutation1[n] sorts random numbers to generate a random permutation."

RandomPermutation2::usage = "RandomPermutation2[n] uses random transpositions to generate random permutations."

(*
RandomPermutation::usage = "RandomPermutation[n] gives a random permutation of length n."
*)

RandomSubset::usage = "RandomSubset[l] creates a random subset of set l."

RandomTableau::usage = 
"RandomTableau[p] gives a random Young tableau of shape p."

RandomTree::usage = "RandomTree[n] gives a random labeled tree on n vertices."

RandomVertices::usage = "RandomVertices[g] assigns a random embedding to graph g."

RankGraph::usage = "RankGraph[g, l] partitions the vertices into classes based on the shortest geodesic distance to a member of list l."

RankPermutation::usage = "RankPermutation[p] gives the rank of permutation p in lexicographic order."

RankSubset::usage = "RankSubset[l, s] gives the rank, in canonical order, of subset s of set l."

RankedEmbedding::usage = "RankedEmbedding[g, l] performs a ranked embedding of graph g, with the vertices ranked in terms of geodesic distance from a member of list l."

ReadGraph::usage = "ReadGraph[f] reads a graph represented as edge lists from file f, and returns a graph object."

RealizeDegreeSequence::usage = 
"RealizeDegreeSequence[s] gives a semirandom graph with degree sequence s."

RegularGraph::usage = 
"RegularGraph[k, n] gives a semirandom k-regular graph on n vertices, if such
a graph exists."

RegularQ::usage = "RegularQ[g] yields True if g is a regular graph."

RemoveSelfLoops::usage = 
"RemoveSelfLoops[g] gives a graph g with the same edges except for any
self-loops."

RevealCycles::usage = "RevealCycles[p] gives the canonical hidden cycle structure of permutation p."

RootedEmbedding::usage = "RootedEmbedding[g, v] gives a rooted embedding of graph g with vertex v as the root."

RotateVertices::usage = "RotateVertices[v, theta] rotates each vertex position in list v by theta radians around the origin (0,0)."

Runs::usage = "Runs[p] partitions p into contiguous increasing subsequences."

SamenessRelation::usage = "SamenessRelation[l] gives a binary relation from a list of permutations l which is an equivalence relation if l is a permutation group."

SelectionSort::usage = "SelectionSort[l, f] sorts list l using ordering function f."

SelfComplementaryQ::usage = "SelfComplementaryQ[g] yields True if graph g is self-complementary, meaning it is isomorphic to its complement."

ShakeGraph::usage = "ShakeGraph[g, d] performs a random perturbation of the vertices of graph g, with each vertex moving at most a distance d from its original position."

ShortestPathSpanningTree::usage = "ShortestPathSpanningTree[g, v] gives the
shortest-path spanning tree originating from v, so that the shortest path in
graph g from v to any other vertex is the path in the tree."

ShortestPath::usage = "ShortestPath[g, start, end] finds the shortest path between vertices start and end in graph g."

ShowGraph::usage = "ShowGraph[g] displays graph g according to its embedding. ShowGraph[g, Directed] displays directed graph g according to its embedding, with arrows illustrating the orientation of each edge."

ShowLabeledGraph::usage = "ShowLabeledGraph[g] displays graph g according to its embedding, with each vertex labeled with its vertex number. ShowLabeledGraph[g, l] uses the ith element of list l as the label for vertex i."

SignaturePermutation::usage = "SignaturePermutation[p] gives the signature of permutation p."

SimpleQ::usage = "SimpleQ[g] yields True if g is a simple graph, meaning it is unweighted and contains no self-loops."

Spectrum::usage = "Spectrum[g] gives the eigenvalues of graph g."

SpringEmbedding::usage = "SpringEmbedding[g] beautifies the embedding of graph g by modeling the embedding as a system of springs."

StableMarriage::usage = "StableMarriage[mpref, fpref] finds the male optimal stable marriage defined by lists of permutations describing male and female preferences."

Star::usage =
"Star[n] gives a star on n vertices, which is a tree with one vertex of
degree n-1."

StirlingFirst::usage = "StirlingFirst[n, k] gives the Stirling numbers of the first kind."

StirlingSecond::usage = "StirlingSecond[n, k] gives the Stirling numbers of the second kind."

Strings::usage =
"Strings[l, n] gives all possible strings of length n from the elements of
list l."

StronglyConnectedComponents::usage = "StronglyConnectedComponents[g] gives the strongly connected components of directed graph g."

Subsets::usage = "Subsets[l] gives all subsets of set l."

TableauClasses::usage = "TableauClasses[p] partitions the elements of permutation p into classes according to their initial columns during Young tableaux construction."

TableauQ::usage = "TableauQ[t] yields True if and only if t represents a Young tableau."

TableauxToPermutation::usage = "TableauxToPermutation[t1, t2] gives the unique
permutation associated with Young tableaux t1 and t2, where both tableaux have
the same shape. "

Tableaux::usage = "Tableaux[p] gives all tableaux whose shape is given by integer partition p."

ToAdjacencyLists::usage = "ToAdjacencyLists[g] gives an adjacency list representation for graph g."

(*
ToCycles::usage = "ToCycles[p] gives the cycle structure of permutation p."
*)

ToInversionVector::usage = "ToInversionVector[p] gives the inversion vector associated with permutation p."

ToOrderedPairs::usage = "ToOrderedPairs[g] gives a list of ordered pairs representing the edges of undirected graph g."

ToUnorderedPairs::usage = "ToUnorderedPairs[g] gives a list of vertex pairs representing graph g, with one pair per undirected edge."

TopologicalSort::usage = "TopologicalSort[g] gives a permutation of the vertices of directed acyclic graph g such that an edge (i,j) implies vertex i appears before vertex j."

TransitiveClosure::usage = "TransitiveClosure[g] finds the transitive closure of graph g, the superset of g which contains edge {x,y} iff there is a path from x to y."

TransitiveQ::usage = "TransitiveQ[g] yields True if graph g defines a transitive relation."

TransitiveReduction::usage = "TransitiveReduction[g] finds the smallest graph which has the same transitive closure as g."

TranslateVertices::usage = "TranslateVertices[v, {x,y}] adds the vector {x,y} to each vertex in list v."

TransposePartition::usage = "TransposePartition[p] reflects a partition p of k parts along the main diagonal, creating a partition with maximum part k."

TransposeTableau::usage = "TransposeTableau[t] reflects a Young tableau t along the main diagonal, creating a different tableau."

TravelingSalesmanBounds::usage = "TravelingSalesmanBounds[g] gives upper and lower bounds on the minimum cost traveling salesman tour of graph g."

TravelingSalesman::usage = "TravelingSalesman[g] finds the optimal traveling salesman tour in graph g."

TreeQ::usage = "TreeQ[g] yields True if graph g is a tree."

TriangleInequalityQ::usage = "TriangleInequalityQ[g] yields True if the weight function defined by the adjacency matrix of graph g satisfies the triangle inequality."

Turan::usage = "Turan[n,p] gives the Turan graph, the extremal graph on n vertices which does not contain CompleteGraph[p]."

TwoColoring::usage = "TwoColoring[g] finds a two-coloring of graph g if g is bipartite."

Undirected::usage = "Undirected is an option to inform certain functions that the graph is undirected."

UndirectedQ::usage = "UndirectedQ[g] yields True if graph g is undirected."

UnionSet::usage = "UnionSet[a, b, s] merges the sets containing a and b in union-find data structure s."

UnweightedQ::usage = "UnweightedQ[g] yields True if all entries in the adjacency matrix of graph g are zero or one."

V::usage = "V[g] gives the order or number of vertices of graph g."

VertexColoring::usage = "VertexColoring[g] uses Brelaz's heuristic to find a good, but not necessarily minimal, vertex coloring of graph g."

VertexConnectivity::usage = "VertexConnectivity[g] gives the minimum number of vertices whose deletion from graph g disconnects it."

VertexCoverQ::usage = "VertexCoverQ[g, c] yields True if the vertices in list c define a vertex cover of graph g."

Vertices::usage = "Vertices[g] gives the embedding of graph g."

WeaklyConnectedComponents::usage = "WeaklyConnectedComponents[g] gives the weakly connected components of directed graph g."

Wheel::usage = "Wheel[n] gives a wheel on n vertices, which is the join of CompleteGraph[1] and Cycle[n-1]."

WriteGraph::usage = "WriteGraph[g, f] writes graph g to file f using an edge list representation."

Begin["`Private`"]

(*
PermutationQ[p_List] := (Sort[p] == Range[Length[p]])
*)

Permute[l_List,p_?PermutationQ] := l [[ p ]]

LexicographicPermutations[{}] := {{}}

LexicographicPermutations[{l_}] := {{l}}

LexicographicPermutations[{a_,b_}] := {{a,b},{b,a}}

LexicographicPermutations[l_List] :=
	Module[{i,n=Length[l]},
		Apply[
			Join,
			Table[ 
				Map[
					(Prepend[#,l[[i]]])&,
					LexicographicPermutations[
						Complement[l,{l[[i]]}]
					]
				],
				{i,n}
			]
		]
	]

RankPermutation[{1}] = 0

RankPermutation[p_?PermutationQ] := (p[[1]]-1) (Length[Rest[p]]!) +
	RankPermutation[ Map[(If[#>p[[1]], #-1, #])&, Rest[p]] ]

NthPermutation[n1_Integer,l_List] :=
	Module[{k, n=n1, s=l, i},
		Table[
			n = Mod[n,(i+1)!];
			k = s [[Quotient[n,i!]+1]];
			s = Complement[s,{k}];
			k,
			{i,Length[l]-1,0,-1}
		]
	]

NextPermutation[p_?PermutationQ] :=
	NthPermutation[ RankPermutation[p]+1, Sort[p] ]

RandomPermutation1[n_Integer?Positive] :=
	Map[ Last, Sort[ Map[({Random[],#})&,Range[n]] ] ]

RandomPermutation2[n_Integer?Positive] :=
	Module[{p = Range[n],i,x},
		Do [
			x = Random[Integer,{1,i}];
			{p[[i]],p[[x]]} = {p[[x]],p[[i]]},
			{i,n,2,-1}
		];
		p
	]

(*
RandomPermutation[n_Integer?Positive] := RandomPermutation1[n]
*)

MinimumChangePermutations[l_List] :=
	Module[{i=1,c,p=l,n=Length[l],k},
		c = Table[1,{n}];
		Join[
			{l},
			Table[
				While [ c[[i]] >= i, c[[i]] = 1; i++];
				If[OddQ[i], k=1, k=c[[i]] ];
				{p[[i]],p[[k]]} = {p[[k]],p[[i]]};
				c[[i]]++;
				i = 2;
				p,
				{n!-1}
			]
		]
	]

Backtrack[space_List,partialQ_,solutionQ_,flag_:One] :=
	Module[{n=Length[space],all={},done,index,v=2,solution},
		index=Prepend[ Table[0,{n-1}],1];
		While[v > 0,
			done = False;
			While[!done && (index[[v]] < Length[space[[v]]]),
				index[[v]]++;
				done = Apply[partialQ,{Solution[space,index,v]}];
			];
			If [done, v++, index[[v--]]=0 ];
			If [v > n,
				solution = Solution[space,index,n];
				If [Apply[solutionQ,{solution}],
					If [SameQ[flag,All],
						AppendTo[all,solution],
						all = solution; v=0
					]
				];
				v--
			]
		];
		all
	]

Solution[space_List,index_List,count_Integer] :=
	Module[{i}, Table[space[[ i,index[[i]] ]], {i,count}] ]

DistinctPermutations[s_List] := Permutations[s] /; (Length[s] == 1)

DistinctPermutations[s_List] :=
	Module[{freq,alph=Union[s],n=Length[s]},
		freq = Map[ (Count[s,#])&, alph];
		Map[
			(alph[[#]])&,
			Backtrack[
				Table[Range[Length[alph]],{n}],
				(Count[#,Last[#]] <= freq[[Last[#]]])&,
				(Count[#,Last[#]] <= freq[[Last[#]]])&,
				All
			]
		]
	]

MinOp[l_List,f_] :=
	Module[{min=First[l]},
		Scan[ (If[ Apply[f,{#,min}], min = #])&, l];
		Return[min];
	]

SelectionSort[l_List,f_] :=
	Module[{where,item,unsorted=l},
		Table[
			item = MinOp[unsorted, f];
			{where} = First[ Position[unsorted,item] ];
			unsorted = Drop[unsorted,{where,where}];
			item,
			{Length[l]}
		]
	]

(* BinarySearch code contributed by Paul Abbott, August 1995 *)

SetAttributes[BinarySearch,HoldFirst];

BinarySearch[l_, k_Integer, f_:Identity]:=
	Module[{lo = 1, mid, hi, el, ls = l},
		(
		    While[lo <= hi,
			If[(el=f[ls[[mid =
				    Floor[ (lo + hi)/2 ]]]])==k,
			   Return[mid]
			];
			If[el > k, hi = mid-1, lo = mid+1]
		    ];
		    Return[lo-1/2]
		) /; (hi = Length[l]; hi > 0)
 	];
 

MultiplicationTable[elems_List,op_] :=
	Module[{i,j,n=Length[elems],p},
		Table[
			p = Position[elems, Apply[op,{elems[[i]],elems[[j]]}]];
			If [p === {}, 0, p[[1,1]]],
			{i,n},{j,n}
		]
	]

InversePermutation[p_?PermutationQ] :=
	Module[{inverse=p},
		inverse[[p]] = Range[Length[p]];
		inverse
	]

EquivalenceRelationQ[r_?SquareMatrixQ] :=
	ReflexiveQ[r] && SymmetricQ[r] && TransitiveQ[r]
EquivalenceRelationQ[g_Graph] := EquivalenceRelationQ[ Edges[g] ]

SquareMatrixQ[{}] = True
SquareMatrixQ[r_] := MatrixQ[r] && (Length[r] == Length[r[[1]]])

ReflexiveQ[r_?SquareMatrixQ] := 
	Module[{i}, Apply[And, Table[(r[[i,i]]!=0),{i,Length[r]}] ] ]

TransitiveQ[r_?SquareMatrixQ] := TransitiveQ[ Graph[r,RandomVertices[Length[r]]] ]
TransitiveQ[r_Graph] := IdenticalQ[r,TransitiveClosure[r]]

SymmetricQ[r_?SquareMatrixQ] := (r === Transpose[r])

EquivalenceClasses[r_List?EquivalenceRelationQ] :=
	ConnectedComponents[ Graph[r,RandomVertices[Length[r]]] ]
EquivalenceClasses[g_Graph?EquivalenceRelationQ] := ConnectedComponents[g]

PermutationGroupQ[perms_List] :=
	FreeQ[ MultiplicationTable[perms,Permute], 0] &&
		EquivalenceRelationQ[SamenessRelation[perms]]

SamenessRelation[perms_List] :=
	Module[{positions = Transpose[perms], i, j, n=Length[First[perms]]},
		Table[
			If[ MemberQ[positions[[i]],j], 1, 0],
			{i,n}, {j,n}
		]
	] /; perms != {}

ToCycles[p1_?PermutationQ] :=
	Module[{p=p1,m,n,cycle,i},
		Select[
			Table[
				m = n = p[[i]];
				cycle = {};
				While[ p[[n]] != 0,
					AppendTo[cycle,m=n];
					n = p[[n]];
					p[[m]] = 0 
				];
				cycle,
				{i,Length[p]}
			],
			(# =!= {})&
		]
	]

(*
FromCycles[cyc_List] := 
	Module[{p=Table[0,{Length[Flatten[cyc]]}], pos},
		Scan[
			(pos = Last[#];
			 Scan[ Function[c, pos = p[[pos]] = c], #])&,
			cyc
		];
		p
	]
*)

HideCycles[c_List] := 
	Flatten[
		Sort[
			Map[(RotateLeft[#,Position[#,Min[#]] [[1,1]] - 1])&, c],
			(#1[[1]] > #2[[1]])&
		]
	]

RevealCycles[p_?PermutationQ] :=
	Module[{start=end=1, cycles={}},
		While [end <= Length[p],
			If [p[[start]] > p[[end]],
				AppendTo[ cycles, Take[p,{start,end-1}] ];
				start = end,
				end++
			]
		];
		Append[cycles,Take[p,{start,end-1}]]
	]

NumberOfPermutationsByCycles[n_Integer,m_Integer] := (-1)^(n-m) StirlingS1[n,m]

StirlingFirst[n_Integer,m_Integer] := StirlingFirst1[n,m]

StirlingFirst1[n_Integer,0] := If [n == 0, 1, 0] 
StirlingFirst1[0,m_Integer] := If [m == 0, 1, 0]

StirlingFirst1[n_Integer,m_Integer] := StirlingFirst1[n,m] =
	(n-1) StirlingFirst1[n-1,m] + StirlingFirst1[n-1, m-1]

StirlingSecond[n_Integer,m_Integer] := StirlingSecond1[n,m]

StirlingSecond1[n_Integer,0] := If [n == 0, 1, 0]
StirlingSecond1[0,m_Integer] := If [m == 0, 1, 0]

StirlingSecond1[n_Integer,m_Integer] := StirlingSecond1[n,m] =
	m StirlingSecond1[n-1,m] + StirlingSecond1[n-1,m-1]

SignaturePermutation[p_?PermutationQ] := (-1) ^ (Length[p]-Length[ToCycles[p]])

Polya[g_List,m_] := Apply[ Plus, Map[(m^Length[ToCycles[#]])&,g] ] / Length[g]

ToInversionVector[p_?PermutationQ] :=
	Module[{i,inverse=InversePermutation[p]},
		Table[
			Length[ Select[Take[p,inverse[[i]]], (# > i)&] ],
			{i,Length[p]-1}
		]
	]

FromInversionVector[vec_List] :=
	Module[{n=Length[vec]+1,i,p},
		p={n};
		Do [
			p = Insert[p, i, vec[[i]]+1],
			{i,n-1,1,-1}
		];
		p
	]

Inversions[p_?PermutationQ] := Apply[Plus,ToInversionVector[p]]

Index[p_?PermutationQ]:=
	Module[{i},
		Sum[ If [p[[i]] > p[[i+1]], i, 0], {i,Length[p]-1} ]
	]

Runs[p_?PermutationQ] :=
	Map[
		(Apply[Take,{p,{#[[1]]+1,#[[2]]}}])&,
		Partition[
			Join[
				{0},
				Select[Range[Length[p]-1], (p[[#]]>p[[#+1]])&],
				{Length[p]}
			],
			2,
			1
		]
	]

Eulerian[n_Integer,k_Integer] := Eulerian1[n,k]

Eulerian1[0,k_Integer] := If [k==1, 1, 0]
Eulerian1[n_Integer,k_Integer] := Eulerian1[n,k] =
	k Eulerian1[n-1,k] + (n-k+1) Eulerian1[n-1,k-1]

InvolutionQ[p_?PermutationQ] := p[[p]] == Range[Length[p]]

NumberOfInvolutions[n_Integer] :=
	Module[{k},
		n! Sum[1/((n - 2k)! 2^k k!), {k, 0, Quotient[n, 2]}]
	]

DerangementQ[p_?PermutationQ] :=
	!(Apply[ Or, Map[( # == p[[#]] )&, Range[Length[p]]] ])

NumberOfDerangements[0] = 1;
NumberOfDerangements[n_] := n * NumberOfDerangements[n-1] + (-1)^n

Derangements[n_Integer] := Derangements[Range[n]]
Derangements[p_?PermutationQ] := Select[ Permutations[p], DerangementQ ]

Josephus[n_Integer,m_Integer] :=
	Module[{live=Range[n],next},
		InversePermutation[
			Table[
				next = RotateLeft[live,m-1];
				live = Rest[next];
				First[next],
				{n}
			]
		]
	]

Heapify[p_List] :=
	Module[{j,heap=p},
		Do [
			heap = Heapify[heap,j],
			{j,Quotient[Length[p],2],1,-1}
		];
		heap
	]

Heapify[p_List, k_Integer] :=
	Module[{hp=p, i=k, l, n=Length[p]},
		While[ (l = 2 i) <= n,
			If[ (l < n) && (hp[[l]] > hp[[l+1]]), l++ ];
			If[ hp[[i]] > hp[[l]],
				{hp[[i]],hp[[l]]}={hp[[l]],hp[[i]]};
				i = l,
				i = n+1
			];
		];
		hp
	]

RandomHeap[n_Integer] := Heapify[RandomPermutation[n]]

HeapSort[p_List] :=
	Module[{heap=Heapify[p],min},
		Append[
			Table[
				min = First[heap];
				heap[[1]] = heap[[n]];
				heap = Heapify[Drop[heap,-1],1];
				min,
				{n,Length[p],2,-1}
			],
			Max[heap]
		]
	]

Strings[l_List,0] := { {} }

Strings[l_List,k_Integer?Positive] :=
	Module[{oneless = Strings[l,k-1],i,n=Length[l]},
		Apply[Join, Table[ Map[(Prepend[#,l[[i]]])&, oneless], {i,n}] ]
	]

NthSubset[n_Integer,m_Integer] := NthSubset[n,Range[m]]
NthSubset[n_Integer,l_List] :=
	l[[ Flatten[ Position[Reverse[IntegerDigits[ Mod[n,2^Length[l]],2]],1] ] ]]

BinarySubsets[l_List] :=
	Module[{pos=Reverse[Range[Length[l]]], n=Length[l]},
		Map[(l[[ Reverse[Select[pos*#, Positive]] ]])&, Strings[{0,1},n] ]
	]

NextSubset[set_List,subset_List] := NthSubset[ RankSubset[set,subset], set  ]

RankSubset[set_List,subset_List] :=
	Module[{i,n=Length[set]},
		Sum[ 2^(i-1) * If[ MemberQ[subset,set[[i]]], 1, 0], {i,n}]
	]

RandomSubset[set_List] := NthSubset[Random[Integer,2^(Length[set])-1],set]

GrayCode[l_List] := GrayCode[l,{{}}]

GrayCode[{},prev_List] := prev

GrayCode[l_List,prev_List] :=
	GrayCode[
		Rest[l],
		Join[ prev, Map[(Append[#,First[l]])&,Reverse[prev]] ]
	]

Subsets[l_List] := GrayCode[l]
Subsets[n_Integer] := GrayCode[Range[n]]

LexicographicSubsets[l_List] := LexicographicSubsets[l,{{}}]

LexicographicSubsets[{},s_List] := s

LexicographicSubsets[l_List,subsets_List] :=
	LexicographicSubsets[
		Rest[l],
		Join[
			subsets,
			Map[(Prepend[#,First[l]])&,LexicographicSubsets[Rest[l],{{}}] ]
		]
	]

KSubsets[l_List,0] := { {} }
KSubsets[l_List,1] := Partition[l,1]
KSubsets[l_List,k_Integer?Positive] := {l} /; (k == Length[l])
KSubsets[l_List,k_Integer?Positive] := {}  /; (k > Length[l])

KSubsets[l_List,k_Integer?Positive] :=
	Join[
		Map[(Prepend[#,First[l]])&, KSubsets[Rest[l],k-1]],
		KSubsets[Rest[l],k]
	]

NextKSubset[set_List,subset_List] :=
	Take[set,Length[subset]] /; (Take[set,-Length[subset]] === subset)

NextKSubset[set_List,subset_List] :=
	Module[{h=1, x=1},
		While [set[[-h]] == subset[[-h]], h++];
		While [set[[x]] =!= subset[[-h]], x++];
		Join[ Drop[subset,-h], Take[set, {x+1,x+h}] ]
	]

RandomKSubset[n_Integer,k_Integer] := RandomKSubset[Range[n],k]

RandomKSubset[set_List,k_Integer] := 
	Module[{s=Range[Length[set]],i,n=Length[set],x},
		set [[
			Sort[
				Table[
					x=Random[Integer,{1,i}];
					{s[[i]],s[[x]]} = {s[[x]],s[[i]]};
					s[[i]],
					{i,n,n-k+1,-1}
				]
			]
		]]
	] /; (k <= Length[set])

PartitionQ[p_List] := (Min[p]>0) && Apply[And, Map[IntegerQ,p]]

Partitions[n_Integer] := Partitions[n,n]

Partitions[n_Integer,_] := {} /; (n<0)
Partitions[0,_] := { {} }
Partitions[n_Integer,1] := { Table[1,{n}] }
Partitions[_,0] := {}

Partitions[n_Integer,maxpart_Integer] :=
	Join[
		Map[(Prepend[#,maxpart])&, Partitions[n-maxpart,maxpart]],
		Partitions[n,maxpart-1]
	]

NextPartition[p_List] := Join[Drop[p,-1],{Last[p]-1,1}]  /; (Last[p] > 1)

NextPartition[p_List] := {Apply[Plus,p]}  /; (Max[p] == 1)

NextPartition[p_List] :=
	Module[{index,k,m},
		{index} = First[ Position[p,1] ];
		k = p[[index-1]] - 1;
		m = Apply[Plus,Drop[p,index-1]] + k + 1;
		Join[
			Take[p,index-2],
			Table[k,{Quotient[m,k]}],
			If [Mod[m,k] == 0, {}, {Mod[m,k]}]
		]
	]

FerrersDiagram[p1_List] :=
	Module[{i,j,n=Length[p1],p=Sort[p1]},
		Show[
			Graphics[
				Join[
					{PointSize[ Min[0.05,1/(2 Max[p])] ]},
					Table[Point[{i,j}], {j,n}, {i,p[[j]]}]
				],
				{AspectRatio -> 1, PlotRange -> All}
			]
		]
	]

TransposePartition[{}] := {}

TransposePartition[p_List] :=
	Module[{s=Select[p,(#>0)&], i, row, r},
		row = Length[s];
		Table [
			r = row;
			While [s[[row]]<=i, row--];
			r,
			{i,First[s]}
		]
	]

DurfeeSquare[s_List] :=
	Module[{i,max=1},
		Do [
			If [s[[i]] >= i, max=i],
			{i,2,Min[Length[s],First[s]]}
		];
		max
	]

DurfeeSquare[{}] := 0

NumberOfPartitions[n_Integer] := NumberOfPartitions1[n]

NumberOfPartitions1[n_Integer] := 0  /; (n < 0)
NumberOfPartitions1[n_Integer] := 1  /; (n == 0)

NumberOfPartitions1[n_Integer] := NumberOfPartitions1[n] =
	Module[{m},
		Sum[ (-1)^(m+1) NumberOfPartitions1[n - m (3m-1)/2] +
			(-1)^(m+1) NumberOfPartitions1[n - m (3m+1)/2],
			{m, Ceiling[ (1+Sqrt[1.0 + 24n])/6 ], 1, -1}
		]
	]

(* from Paul Chase *)

RandomPartition[n_Integer?Positive] :=
  Module[{mult = Table[0, {n}], j, d, r=n, z},
    While[ (r > 0),
      d = 1;  j = 0;
      z = Random[] r PartitionsP[r];
      While [z >= 0,
         j++;
         If [r-j*d < 0, {j=1; d++;}];
         z -= j*PartitionsP[r-j*d];
      ];
      r -= j d;
      mult[[j]] += d;
    ];
    Reverse[Flatten[Table[Table[j, {mult[[j]]}], {j, Length[mult]}]]]
  ]

NumberOfCompositions[n_,k_] := Binomial[ n+k-1, n ]

RandomComposition[n_Integer,k_Integer] :=
	Map[
		(#[[2]] - #[[1]] - 1)&,
		Partition[Join[{0},RandomKSubset[Range[n+k-1],k-1],{n+k}], 2, 1]
	]

Compositions[n_Integer,k_Integer] :=
	Map[
		(Map[(#[[2]]-#[[1]]-1)&, Partition[Join[{0},#,{n+k}],2,1] ])&,
		KSubsets[Range[n+k-1],k-1]
	]

NextComposition[l_List] := 
	Module[{c=l, h=1, t},
		While[c[[h]] == 0, h++];
		{t,c[[h]]} = {c[[h]],0};
		c[[1]] = t - 1;
		c[[h+1]]++;
		c
	]

NextComposition[l_List] :=
	Join[{Apply[Plus,l]},Table[0,{Length[l]-1}]] /; Last[l]==Apply[Plus,l]

TableauQ[{}] = True
TableauQ[t_List] :=
	And [
		Apply[ And, Map[(Apply[LessEqual,#])&,t] ],
		Apply[ And, Map[(Apply[LessEqual,#])&,TransposeTableau[t]] ],
		Apply[ GreaterEqual, Map[Length,t] ],
		Apply[ GreaterEqual, Map[Length,TransposeTableau[t]] ]
	]

TransposeTableau[tb_List] :=
	Module[{t=Select[tb,(Length[#]>=1)&],row},
		Table[
			row = Map[First,t];
			t = Map[ Rest, Select[t,(Length[#]>1)&] ];
			row,
			{Length[First[tb]]}
		]
	]

ShapeOfTableau[t_List] := Map[Length,t]

InsertIntoTableau[e_Integer,{}] := { {e} }

InsertIntoTableau[e_Integer, t1_?TableauQ] :=
	Module[{item=e,row=0,col,t=t1},
		While [row < Length[t],
			row++;
			If [Last[t[[row]]] <= item,
				AppendTo[t[[row]],item];
				Return[t]
			];
			col = Ceiling[ BinarySearch[t[[row]],item] ];
			{item, t[[row,col]]} = {t[[row,col]], item};
		];
		Append[t, {item}]
	]

ConstructTableau[p_List] := ConstructTableau[p,{}]

ConstructTableau[{},t_List] := t

ConstructTableau[p_List,t_List] :=
	ConstructTableau[Rest[p], InsertIntoTableau[First[p],t]]

DeleteFromTableau[t1_?TableauQ,r_Integer]:=
	Module [{t=t1, col, row, item=Last[t1[[r]]]},
		col = Length[t[[r]]];
		If[col == 1, t = Drop[t,-1], t[[r]] = Drop[t[[r]],-1]];
		Do [
			While [t[[row,col]]<=item && Length[t[[row]]]>col, col++];
			If [item < t[[row,col]], col--];
			{item,t[[row,col]]} = {t[[row,col]],item},
			{row,r-1,1,-1}
		];
		t
	]

TableauxToPermutation[p1_?TableauQ,q1_?TableauQ] :=
	Module[{p=p1, q=q1, row, firstrow},
		Reverse[
			Table[
				firstrow = First[p];
				row = Position[q, Max[q]] [[1,1]];
				p = DeleteFromTableau[p,row];
				q[[row]] = Drop[ q[[row]], -1];
				If[ p == {},
					First[firstrow],
					First[Complement[firstrow,First[p]]]
				],
				{Apply[Plus,ShapeOfTableau[p1]]}
			]
		]
	] /; ShapeOfTableau[p1] === ShapeOfTableau[q1]

LastLexicographicTableau[s_List] :=
	Module[{c=0},
		Map[(c+=#; Range[c-#+1,c])&, s]
	]

FirstLexicographicTableau[s_List] :=
	TransposeTableau[ LastLexicographicTableau[ TransposePartition[s] ] ]

NextTableau[t_?TableauQ] :=
	Module[{s,y,row,j,count=0,tj,i,n=Max[t]},
		y = TableauToYVector[t];
		For [j=2, (j<n)  && (y[[j]]>=y[[j-1]]), j++, Null];
		If [y[[j]] >= y[[j-1]],
			Return[ FirstLexicographicTableau[ ShapeOfTableau[t] ] ]
		];
		s = ShapeOfTableau[ Table[Select[t[[i]],(#<=j)&], {i,Length[t]}] ];
		{row} = Last[ Position[ s, s[[ Position[t,j] [[1,1]] + 1 ]] ] ];
		s[[row]] --;
		tj = FirstLexicographicTableau[s];
		If[ Length[tj] < row,
			tj = Append[tj,{j}],
			tj[[row]] = Append[tj[[row]],j]
		];
		Join[
			Table[
				Join[tj[[i]],Select[t[[i]],(#>j)&]],
				{i,Length[tj]}
			],
			Table[t[[i]],{i,Length[tj]+1,Length[t]}]
		]
	]

Tableaux[s_List] :=
	Module[{t = LastLexicographicTableau[s]},
		Table[ t = NextTableau[t], {NumberOfTableaux[s]} ]
	]

Tableaux[n_Integer?Positive] := Apply[ Join, Map[ Tableaux, Partitions[n] ] ]

YVectorToTableau[y_List] :=
	Module[{k},
		Table[ Flatten[Position[y,k]], {k,Length[Union[y]]}]
	]

TableauToYVector[t_?TableauQ] :=
	Module[{i,y=Table[1,{Length[Flatten[t]]}]},
		Do [ Scan[ (y[[#]]=i)&, t[[i]] ], {i,2,Length[t]} ];
		y
	]

NumberOfTableaux[{}] := 1
NumberOfTableaux[s_List] := 
	Module[{row,col,transpose=TransposePartition[s]},
		(Apply[Plus,s])! /
		Product [
			(transpose[[col]]-row+s[[row]]-col+1),
			{row,Length[s]}, {col,s[[row]]}
		]
	]

NumberOfTableaux[n_Integer] := Apply[Plus, Map[NumberOfTableaux, Partitions[n]]]

(*
CatalanNumber[n_Integer] := Binomial[2n,n]/(n+1)	/; (n>=0)
*)

RandomTableau[shape_List] :=
	Module[{i=j=n=Apply[Plus,shape],done,l,m,h=1,k,y,p=shape},
		y= Join[TransposePartition[shape],Table[0,{n - Max[shape]}]];
		Do[
			{i,j} = RandomSquare[y,p]; done = False;
			While [!done,
				h = y[[j]] + p[[i]] - i - j;
				If[ h != 0,
					If[ Random[] < 0.5,
						j = Random[Integer,{j,p[[i]]}],
						i = Random[Integer,{i,y[[j]]}]
					],
					done = True
				];
			];
			p[[i]]--; y[[j]]--;
			y[[m]] = i,
			{m,n,1,-1}
		];
		YVectorToTableau[y]
	]

RandomSquare[y_List,p_List] :=
	Module[{i=Random[Integer,{1,First[y]}], j=Random[Integer,{1,First[p]}]},
		While[(i > y[[j]]) || (j > p[[i]]), 
			i = Random[Integer,{1,First[y]}];
			j = Random[Integer,{1,First[p]}]
		];
		{i,j}
	]

TableauClasses[p_?PermutationQ] :=
	Module[{classes=Table[{},{Length[p]}],t={}},
		Scan [
			(t = InsertIntoTableau[#,t];
			 PrependTo[classes[[Position[First[t],#] [[1,1]] ]], #])&,
			p
		];
		Select[classes, (# != {})&]
	]

LongestIncreasingSubsequence[p_?PermutationQ] :=
	Module[{c,x,xlast},
		c = TableauClasses[p];
		xlast = x = First[ Last[c] ];
		Append[
			Reverse[
				Map[
					(x = First[ Intersection[#,
					       Take[p, Position[p,x][[1,1]] ] ] ])&,
					Reverse[ Drop[c,-1] ]
				]
			],
			xlast
		]
	]

LongestIncreasingSubsequence[{}] := {}

AddToEncroachingLists[k_Integer,{}] := {{k}}

AddToEncroachingLists[k_Integer,l_List] :=
	Append[l,{k}]  /; (k > First[Last[l]]) && (k < Last[Last[l]])

AddToEncroachingLists[k_Integer,l1_List] :=
	Module[{i,l=l1},
		If [k <= First[Last[l]],
			i = Ceiling[ BinarySearch[l,k,First] ];
			PrependTo[l[[i]],k],
			i = Ceiling[ BinarySearch[l,-k,(-Last[#])&] ];
			AppendTo[l[[i]],k]
		];
		l
	]

EncroachingListSet[l_List] := EncroachingListSet[l,{}]
EncroachingListSet[{},e_List] := e

EncroachingListSet[l_List,e_List] :=
	EncroachingListSet[Rest[l], AddToEncroachingLists[First[l],e] ]

Edges[Graph[e_,_]] := e

Vertices[Graph[_,v_]] := v

V[Graph[e_,_]] := Length[e]

M[Graph[g_,_],___] := Apply[Plus, Map[(Apply[Plus,#])&,g] ] / 2
M[Graph[g_,_],Directed] := Apply[Plus, Map[(Apply[Plus,#])&,g] ]

ChangeVertices[g_Graph,v_List] := Graph[ Edges[g], v ]

ChangeEdges[g_Graph,e_List] := Graph[ e, Vertices[g] ]

AddEdge[Graph[g_,v_],{x_,y_},Directed] :=
	Module[ {gnew=g},
		gnew[[x,y]] ++;
		Graph[gnew,v]
	]

AddEdge[g_Graph,{x_,y_},flag_:Undirected] :=
	AddEdge[ AddEdge[g, {x,y}, Directed], {y,x}, Directed]

DeleteEdge[Graph[g_,v_],{x_,y_},Directed] :=
	Module[ {gnew=g},
		If [ g[[x,y]] > 1, gnew[[x,y]]--, gnew[[x,y]] = 0];
		Graph[gnew,v]
	]

DeleteEdge[g_Graph,{x_,y_},flag_:Undirected] :=
	DeleteEdge[ DeleteEdge[g, {x,y}, Directed], {y,x}, Directed]

AddVertex[g_Graph] := GraphUnion[g, CompleteGraph[1]]

DeleteVertex[g_Graph,v_Integer] := InduceSubgraph[g,Complement[Range[V[g]],{v}]]

Spectrum[Graph[g_,_]] := Eigenvalues[g]

ToAdjacencyLists[Graph[g_,_]] :=
	Map[ (Flatten[ Position[ #, _?(Function[n, n!=0])] ])&, g ]

FromAdjacencyLists[e_List] :=
	Module[{blanks = Table[0,{Length[e]}] },
		Graph[
			Map [ (MapAt[ 1&,blanks,Partition[#,1]])&, e ],
			CircularVertices[Length[e]]
		]
	]

FromAdjacencyLists[e_List,v_List] := ChangeVertices[FromAdjacencyLists[e], v]

ToOrderedPairs[g_Graph] := Position[ Edges[g], _?(Function[n,n != 0]) ]

ToUnorderedPairs[g_Graph] := Select[ ToOrderedPairs[g], (#[[1]] < #[[2]])& ]

FromOrderedPairs[l_List] := 
	Module[{n=Max[l]},
		Graph[
			MapAt[1&, Table[0,{n},{n}],l],
			CircularVertices[n]
		]
	]

FromOrderedPairs[{}] := Graph[{},{}]
FromOrderedPairs[l_List,v_List] := 
	Graph[ MapAt[1&, Table[0,{Length[v]},{Length[v]}], l], v]

FromUnorderedPairs[l_List] := MakeUndirected[ FromOrderedPairs[l] ]
FromUnorderedPairs[l_List,v_List] := MakeUndirected[ FromOrderedPairs[l,v] ]

PseudographQ[Graph[g_,_]] :=
	Module[{i},
		Apply[Or, Table[ g[[i,i]]!=0, {i,Length[g]} ]]
	]

UnweightedQ[Graph[g_,_]] := Apply[ And, Map[(#==0 || #==1)&, Flatten[g] ] ]

SimpleQ[g_Graph] := (!PseudographQ[g]) && (UnweightedQ[g])

RemoveSelfLoops[g_Graph] :=
	Module[{i,e=Edges[g]},
		Do [ e[[i,i]]=0, {i,V[g]} ];
		Graph[e, Vertices[g]]
	]	

EmptyQ[g_Graph] := Edges[g] == Table[0, {V[g]}, {V[g]}]

CompleteQ[g_Graph] := 
	(
	Edges[RemoveSelfLoops[g]] == Edges[ CompleteGraph[V[g]] ]
	)

InduceSubgraph[g_Graph,{}] := Graph[{},{}]

InduceSubgraph[Graph[g_,v_],s_List] :=
	Graph[Transpose[Transpose[g[[s]]] [[s]] ],v[[s]]] /; (Length[s]<=Length[g])

Contract[g_Graph,{u_Integer,v_Integer}] :=
	Module[{o,e,i,n=V[g],newg,range=Complement[Range[V[g]],{u,v}]},
		newg = InduceSubgraph[g,range];
		e = Edges[newg]; o = Edges[g];
		Graph[
			Append[
				Table[
					Append[e[[i]],
						If[o[[range[[i]],u]]>0 ||
							o[[range[[i]],v]]>0,1,0] ],
					{i,n-2}
				],
				Append[
					Map[(If[o[[u,#]]>0||o[[v,#]]>0,1,0])&,range],
					0
				]
			],
			Join[Vertices[newg], {(Vertices[g][[u]]+Vertices[g][[v]])/2}]
		]
	] /; V[g] > 2

Contract[g_Graph,_] := CompleteGraph[1]	/; V[g] == 2

GraphComplement[Graph[g_,v_]] :=
	(
	RemoveSelfLoops[ Graph[ Map[ (Map[ (If [#==0,1,0])&, #])&, g], v ] ]
	)

MakeUndirected[Graph[g_,v_]] :=
	Module[{i,j,n=Length[g]},
		Graph[ Table[If [g[[i,j]]!=0 || g[[j,i]]!=0,1,0],{i,n},{j,n}], v ]
	]

UndirectedQ[Graph[g_,_]] := (Apply[Plus,Apply[Plus,Abs[g-Transpose[g]]]] == 0)

MakeSimple[g_Graph] := 
	(
	MakeUndirected[RemoveSelfLoops[g]]
	)

BFS[g_Graph,start_Integer] :=
	Module[{e,bfi=Table[0,{V[g]}],cnt=1,edges={},queue={start}},
		e = ToAdjacencyLists[g];
		bfi[[start]] = cnt++;
		While[ queue != {},
			{v,queue} = {First[queue],Rest[queue]};
			Scan[
				(If[ bfi[[#]] == 0,
					bfi[[#]] = cnt++;
					AppendTo[edges,{v,#}];
					AppendTo[queue,#]
				])&,
				e[[v]]
			];
		];
		{edges,bfi}
	]

BreadthFirstTraversal[g_Graph,s_Integer,Edge] := First[BFS[g,s]]

BreadthFirstTraversal[g_Graph,s_Integer,___] := InversePermutation[Last[BFS[g,s]]]

DFS[v_Integer] :=
	(
	  dfi[[v]] = cnt++;
	  AppendTo[visit,v];
	  Scan[ (If[dfi[[#]]==0,AppendTo[edges,{v,#}];DFS[#] ])&, e[[v]] ]
	)

DepthFirstTraversal[g_Graph,start_Integer,flag_:Vertex] :=
   Block[{visit={},e=ToAdjacencyLists[g],edges={},dfi=Table[0,{V[g]}],cnt=1},
		DFS[start];
		If[ flag===Edge,
		     edges,
		     visit]
   ]

ShowGraph[g1_Graph,type_:Undirected] :=
	Module[{g=NormalizeVertices[g1]},
		Show[
			Graphics[
				Join[
					PointsAndLines[g],
					If[SameQ[type,Directed],Arrows[g],{}]
				]
			], 
			{AspectRatio->1, PlotRange->FindPlotRange[Vertices[g]]}
		]
	]

MinimumEdgeLength[v_List,pairs_List] :=
	Max[ Select[
		Chop[ Map[(Sqrt[ N[(v[[#[[1]]]]-v[[#[[2]]]]) . 
			(v[[#[[1]]]]-v[[#[[2]]]])] ])&,pairs] ],
		(# > 0)&
	], 0.001 ]

FindPlotRange[v_List] :=
	Module[{xmin=Min[Map[First,v]], xmax=Max[Map[First,v]],
			ymin=Min[Map[Last,v]], ymax=Max[Map[Last,v]]},
		{ {xmin - 0.05 Max[1,xmax-xmin], xmax + 0.05 Max[1,xmax-xmin]},
		  {ymin - 0.05 Max[1,ymax-ymin], ymax + 0.05 Max[1,ymax-ymin]} }
	]

PointsAndLines[Graph[e_List,v_List]] :=
	Module[{pairs=ToOrderedPairs[Graph[e,v]]},
		Join[
			{PointSize[ 0.025 ]},
			Map[Point,Chop[v]],
			Map[(Line[Chop[ v[[#]] ]])&,pairs]
		]
	]

Arrows[Graph[e_,v_]] :=
	Module[{pairs=ToOrderedPairs[Graph[e,v]], size, triangle},
		size = Min[0.05, MinimumEdgeLength[v,pairs]/3];
		triangle={ {0,0}, {-size,size/2}, {-size,-size/2} };
		Map[
			(Polygon[
				TranslateVertices[
					RotateVertices[
						triangle,
						Arctan[Apply[Subtract,v[[#]]]]+Pi
					],
					v[[ #[[2]] ]]
				]
			])&,
			pairs
		]
	]

ShowLabeledGraph[g_Graph] := ShowLabeledGraph[g,Range[V[g]]]
ShowLabeledGraph[g1_Graph,labels_List] :=
	Module[{pairs=ToOrderedPairs[g1], g=NormalizeVertices[g1], v},
		v = Vertices[g];
		Show[
			Graphics[
				Join[
					PointsAndLines[g],
					Map[(Line[Chop[ v[[#]] ]])&, pairs],
					GraphLabels[v,labels]
				]

			],
			{AspectRatio->1, PlotRange->FindPlotRange[v]} 
		]
	]

GraphLabels[v_List,l_List] :=
	Module[{i},
		Table[ Text[ l[[i]],v[[i]]-{0.03,0.03},{0,1} ],{i,Length[v]}]
	]

CircularVertices[0] := {}

CircularVertices[n_Integer] :=
	Module[{i,x = N[2 Pi / n]},
		Chop[ Table[ N[{ (Cos[x i]), (Sin[x i]) }], {i,n} ] ]
	]

CircularVertices[Graph[g_,_]] := Graph[ g, CircularVertices[ Length[g] ] ]

RankGraph[g_Graph, start_List] :=
	Module[ {rank = Table[0,{V[g]}],edges = ToAdjacencyLists[g],v,queue,new},
		Scan[ (rank[[#]] = 1)&, start];
		queue = start;
		While [queue != {},
			v = First[queue];
			new = Select[ edges[[v]], (rank[[#]] == 0)&];
			Scan[ (rank[[#]] = rank[[v]]+1)&, new];
			queue = Join[ Rest[queue], new];
		];
		rank
	]

RankedEmbedding[g_Graph,start_List] := Graph[ Edges[g],RankedVertices[g,start] ]

RankedVertices[g_Graph,start_List] :=
	Module[{i,m,stages,rank,freq = Table[0,{V[g]}]},
		rank = RankGraph[g,start];
		rank = ReplacePart[rank,Max[{rank}]+1,Position[rank,0]];
		stages = Distribution[ rank ];
		Table[
			m = ++ freq[[ rank[[i]] ]];
			{rank[[i]], (m-1) + (1 - stages[[ rank[[i]] ]])/2 },
			{i,V[g]}
		]
	]

Distribution[l_List] := Distribution[l, Union[l]]
Distribution[l_List, set_List] := Map[(Count[l,#])&, set]

Eccentricity[g_Graph] := Map[ Max, AllPairsShortestPath[g] ]
Eccentricity[g_Graph,start_Integer] := Map[ Max, Last[Dijkstra[g,start]] ]

Diameter[g_Graph] := Max[ Eccentricity[g] ]

Radius[g_Graph] := Min[ Eccentricity[g] ]

GraphCenter[g_Graph] := 
	Module[{eccentricity = Eccentricity[g]},
		Flatten[ Position[eccentricity, Min[eccentricity]] ]
	]

RadialEmbedding[g_Graph,ct_Integer] :=
	Module[{center=ct,ang,i,da,theta,n,v,positioned,done,next,e=ToAdjacencyLists[g]},
		ang = Table[{0,2 Pi},{n=V[g]}];
		v = Table[{0,0},{n}];
		positioned = next = done = {center};
		While [next != {},
			center = First[next];
			new = Complement[e[[center]], positioned];
			Do [
				da = (ang[[center,2]]-ang[[center,1]])/Length[new];
				ang[[ new[[i]] ]] = {ang[[center,1]] + (i-1)*da,
					ang[[center,1]] + i*da};
				theta = Apply[Plus,ang[[ new[[i]] ]] ]/2;
				v[[ new[[i]] ]] = v[[center]] +
					N[{Cos[theta],Sin[theta]}],
				{i,Length[new]}
			];
			next = Join[Rest[next],new];
			positioned = Union[positioned,new];
			AppendTo[done,center]
		];
		Graph[Edges[g],v]
	]

RadialEmbedding[g_Graph] := RadialEmbedding[g,First[GraphCenter[g]]];

RootedEmbedding[g_Graph,rt_Integer] :=
	Module[{root=rt,pos,i,x,dx,new,n=V[g],v,done,next,e=ToAdjacencyLists[g]},
		pos = Table[{-Ceiling[Sqrt[n]],Ceiling[Sqrt[n]]},{n}];
		v = Table[{0,0},{n}];
		next = done = {root};
		While [next != {},
			root = First[next];
			new = Complement[e[[root]], done];
			Do [
				dx = (pos[[root,2]]-pos[[root,1]])/Length[new];
				pos[[ new[[i]] ]] = {pos[[root,1]] + (i-1)*dx,
					pos[[root,1]] + i*dx};
				x = Apply[Plus,pos[[ new[[i]] ]] ]/2;
				v[[ new[[i]] ]] = {x,v[[root,2]]-1},
				{i,Length[new]}
			];
			next = Join[Rest[next],new];
			done = Join[done,new]
		];
		Graph[Edges[g],N[v]]
	]

TranslateVertices[v_List,{x_,y_}] := Map[ (# + {x,y})&, v ]
TranslateVertices[Graph[g_,v_],{x_,y_}] := Graph[g, TranslateVertices[v,{x,y}] ]

DilateVertices[v_List,d_] := (d * v)
DilateVertices[Graph[e_,v_],d_] := Graph[e, DilateVertices[v,d]]

RotateVertices[v_List,t_] := 
	Module[{d,theta},
		Map[
			(If[# == {0,0}, {0,0},
				d=Sqrt[#[[1]]^2 + #[[2]]^2];
			 	theta = t + Arctan[#];
				If[!FreeQ[theta, Complex],
				(* hack to get around new ArcTan behavior *)
					theta = Chop[theta]];
			 	N[{d Cos[theta], d Sin[theta]}]
			])&,
			v
		]
	]
RotateVertices[Graph[g_,v_],t_] := Graph[g, RotateVertices[v,t]]

Arctan[{x_,y_}] := Arctan1[Chop[{x,y}]]
Arctan1[{0,0}] := 0
Arctan1[{x_,y_}] := ArcTan[x,y]

NormalizeVertices[v_List] := 
	Module[{v1},
		v1 = TranslateVertices[v,{-Min[v],-Min[v]}];
		DilateVertices[v1, 1/Max[v1,0.01]]
	]

NormalizeVertices[Graph[g_,v_]] := Graph[g, NormalizeVertices[v]]

ShakeGraph[Graph[e_List,v_List], fract_:0.1] :=
	Module[{i,d,a},
		Graph[
			e,
			Table[ 
				d = Random[Real,{0,fract}];
				a = Random[Real,{0, 2 N[Pi]}];
				{N[v[[i,1]] + d Cos[a]], N[v[[i,2]] + d Sin[a]]},
				{i,Length[e]}
			]
		]
	]

CalculateForce[u_Integer,g_Graph,em_List] :=
	Module[{n=V[g],stc=0.25,gr=10.0,e=Edges[g],f={0.0,0.0},spl=1.0,v,dsquared},
		Do [
			dsquared = Max[0.001, Apply[Plus,(em[[u]]-em[[v]])^2] ];
			f += (1-e[[u,v]]) (gr/dsquared) (em[[u]]-em[[v]])
				- e[[u,v]] stc Log[dsquared/spl] (em[[u]]-em[[v]]),
			{v,n}
		];
		f
	]

SpringEmbedding[g_Graph,step_:10,inc_:0.15] :=
	Module[{new=old=Vertices[g],n=V[g],i,u,g1=MakeUndirected[g]},
		Do [
			Do [
				new[[u]] = old[[u]]+inc*CalculateForce[u,g1,old],
				{u,n}
			];
			old = new,
			{i,step}
		];
		Graph[Edges[g],new]
	]

(*	Rewritten for Version 2.0	*)

ReadGraph[file_] :=
	Module[{edgelist={}, v={},x},
		OpenRead[file];
		While[!SameQ[(x = Read[file,Number]), EndOfFile],
			AppendTo[v,Read[file,{Number,Number}]];
			AppendTo[edgelist,
				Convert[Characters[Read[file,String]]]
			];
		];
		Close[file];
		FromAdjacencyLists[edgelist,v]
	]

Toascii[s_String] := First[ ToCharacterCode[s] ]

Convert[l_List] := 
	Module[{ch,num,edge={},i=1},
		While[i <= Length[l],
			If[ DigitQ[ l[[i]] ], 
				num = 0;
				While[ ((i <= Length[l]) && (DigitQ[l[[i]]])),
					num = 10 num + Toascii[l[[i++]]] - Toascii["0"]
				];
				AppendTo[edge,num],
				i++
			];
		];
		edge
	]

WriteGraph[g_Graph,file_] := 
	Module[{edges=ToAdjacencyLists[g],v=N[NormalizeVertices[Vertices[g]]],i,x,y},
		OpenWrite[file];
		Do[
			WriteString[file,"	",ToString[i]];
			{x,y} = Chop[ v [[i]] ];
			WriteString[file,"	",ToString[x],"	",ToString[y]];
			Scan[
				(WriteString[file,"	",ToString[ # ]])&,
				edges[[i]]
			];
			Write[file],
			{i,V[g]}
		];
		Close[file];
	]

GraphUnion[g_Graph,h_Graph] :=
	Module[{maxg=Max[ Map[First,Vertices[g]] ], minh=Min[ Map[First,Vertices[h]] ]},
		FromOrderedPairs[
			Join[ ToOrderedPairs[g], (ToOrderedPairs[h] + V[g])],
			Join[ Vertices[g], Map[({maxg-minh+1,0}+#)&, Vertices[h] ] ]
		]
	]

GraphUnion[1,g_Graph] := g
GraphUnion[0,g_Graph] := EmptyGraph[0];
GraphUnion[k_Integer,g_Graph] := GraphUnion[ GraphUnion[k-1,g], g]

ExpandGraph[g_Graph,n_] := GraphUnion[ g, EmptyGraph[n - V[g]] ] /; V[g] <= n

GraphIntersection[g_Graph,h_Graph] :=
	FromOrderedPairs[
		Intersection[ToOrderedPairs[g],ToOrderedPairs[h]],
		Vertices[g]
	] /; (V[g] == V[h])

GraphDifference[g1_Graph,g2_Graph] :=
	Graph[Edges[g1] - Edges[g2], Vertices[g1]] /; V[g1]==V[g2]

GraphSum[g1_Graph,g2_Graph] :=
	Graph[Edges[g1] + Edges[g2], Vertices[g1]] /; V[g1]==V[g2]

GraphJoin[g_Graph,h_Graph] :=
	Module[{maxg=Max[ Abs[ Map[First,Vertices[g]] ] ]},
		FromUnorderedPairs[
			Join[
				ToUnorderedPairs[g],
				ToUnorderedPairs[h] + V[g],
				CartesianProduct[Range[V[g]],Range[V[h]]+V[g]]
			],
			Join[ Vertices[g], Map[({maxg+1,0}+#)&, Vertices[h]]]
		]
	]

CartesianProduct[a_List,b_List] :=
	Module[{i,j},
		Flatten[ Table[{a[[i]],b[[j]]},{i,Length[a]},{j,Length[b]}], 1]
	]

GraphProduct[g_Graph,h_Graph] :=
	Module[{k,eg=ToOrderedPairs[g],eh=ToOrderedPairs[h],leng=V[g],lenh=V[h]},
		FromOrderedPairs[
			Flatten[
				Join[
					Table[eg+(i-1)*leng, {i,lenh}],
					Map[ (Table[
						{leng*(#[[1]]-1)+k, leng*(#[[2]]-1)+k},
						{k,1,leng}
					      ])&,
					      eh
					]
				],
				1
			],
			ProductVertices[Vertices[g],Vertices[h]]
		]
	]

ProductVertices[vg_,vh_] :=
	Flatten[
		Map[
			(TranslateVertices[
				DilateVertices[vg, 1/(Max[Length[vg],Length[vh]])],
			#])&,
			 RotateVertices[vh,Pi/2]
		],
		1
	]

IncidenceMatrix[g_Graph] :=
	Map[
		( Join[
			Table[0,{First[#]-1}], {1},
			Table[0,{Last[#]-First[#]-1}], {1},
			Table[0,{V[g]-Last[#]}]
		] )&,
		ToUnorderedPairs[g]
	]

LineGraph[g_Graph] :=
	Module[{b=IncidenceMatrix[g], edges=ToUnorderedPairs[g], v=Vertices[g]},
		Graph[
			b . Transpose[b] - 2 IdentityMatrix[Length[edges]],
			Map[ ( (v[[ #[[1]] ]] + v[[ #[[2]] ]]) / 2 )&, edges]
		]
	]

CompleteGraph[0] := Graph[{},{}]
CompleteGraph[1] := Graph[{{0}},{{0,0}}]

CompleteGraph[n_Integer?Positive] := CirculantGraph[n,Range[1,Floor[(n+1)/2]]]

CirculantGraph[n_Integer?Positive,l_List] :=
	Module[{i,r},
		r = Prepend[MapAt[1&,Table[0,{n-1}], Map[List,Join[l,n-l]]], 0];
		Graph[ Table[RotateRight[r,i], {i,0,n-1}], CircularVertices[n] ]
	]

EmptyGraph[n_Integer?Positive] :=
	Module[{i},
		Graph[ Table[0,{n},{n}], Table[{0,i},{i,(1-n)/2,(n-1)/2}] ]
	]

CompleteGraph[l__] :=
	Module[{ll=List[l],t,i,x,row,stages=Length[List[l]]},
		t = FoldList[Plus,0,ll];
		Graph[
			Apply[
				Join,
				Table [
					row = Join[
						Table[1, {t[[i-1]]}],
						Table[0, {t[[i]]-t[[i-1]]}],
						Table[1, {t[[stages+1]]-t[[i]]}]
					];
					Table[row, {ll[[i-1]]}],
					{i,2,stages+1}
				]
			
			],
			Apply [
				Join,
				Table[
					Table[{x,i-1+(1-ll[[x]])/2},{i,ll[[x]]}],
					{x,stages}
				]
			]
		]
	] /; TrueQ[Apply[And, Map[Positive,List[l]]]] && (Length[List[l]]>1)

Turan[n_Integer,p_Integer] :=
	Module[{k = Floor[ n / (p-1) ], r},
		r = n - k (p-1);
		Apply[CompleteGraph, Join[ Table[k,{p-1-r}], Table[k+1,{r}] ] ]
	] /; (n > 0 && p > 1)

Cycle[n_Integer] := CirculantGraph[n,{1}]  /; n>=3

Star[n_Integer?Positive] :=
	Module[{g},
		g = Append [ Table[0,{n-1},{n}], Append[ Table[1,{n-1}], 0] ];
		Graph[
			g + Transpose[g],
			Append[ CircularVertices[n-1], {0,0}]
		]
	]

Wheel[n_Integer] :=
	Module[{i,row = Join[{0,1}, Table[0,{n-4}], {1}]},
		Graph[
			Append[
				Table[ Append[RotateRight[row,i-1],1], {i,n-1}],
				Append[ Table[1,{n-1}], 0]
			],
			Append[ CircularVertices[n-1], {0,0} ]
		]
	] /; If[n >= 4, 
		True,
		Message[Wheel::nvert];
		False]

Wheel::nvert =
"The number of vertices in a wheel must be greater than or equal to four."


Path[1] := CompleteGraph[1]
Path[n_Integer?Positive] :=
	FromUnorderedPairs[ Partition[Range[n],2,1], Map[({#,0})&,Range[n]] ]

GridGraph[n_Integer?Positive,m_Integer?Positive] :=
	GraphProduct[
		ChangeVertices[Path[n], Map[({Max[n,m]*#,0})&,Range[n]]],
		Path[m]
	]

Hypercube[n_Integer] := Hypercube1[n]

Hypercube1[0] := CompleteGraph[1]
Hypercube1[1] := Path[2]
Hypercube1[2] := Cycle[4]

Hypercube1[n_Integer] := Hypercube1[n] =
	GraphProduct[
		RotateVertices[ Hypercube1[Floor[n/2]], 2Pi/5],
		Hypercube1[Ceiling[n/2]]
	]

LabeledTreeToCode[g_Graph] :=
	Module[{e=ToAdjacencyLists[g],i,code},
		Table [
			{i} = First[ Position[ Map[Length,e], 1 ] ];
			code = e[[i,1]];
			e[[code]] = Complement[ e[[code]], {i} ];
			e[[i]] = {};
			code,
			{V[g]-2}
		]
	]

CodeToLabeledTree[l_List] :=
	Module[{m=Range[Length[l]+2],x,i},
		FromUnorderedPairs[
			Append[
				Table[
					x = Min[Complement[m,Drop[l,i-1]]];
					m = Complement[m,{x}];
					{x,l[[i]]},
					{i,Length[l]}
				],
				m
			]
		]
	]

RandomTree[n_Integer?Positive] :=
	RadialEmbedding[CodeToLabeledTree[ Table[Random[Integer,{1,n}],{n-2}] ], 1]

RandomGraph[n_Integer,p_] := RandomGraph[n,p,{1,1}]

RandomGraph[n_Integer,p_,range_List] :=
	Module[{i,g},
		g = Table[ 
			Join[
				Table[0,{i}],
				Table[ 
					If[Random[Real]<p, Random[Integer,range], 0],
					{n-i}
				]
			],
			{i,n}
		];
		Graph[ g + Transpose[g], CircularVertices[n] ]
	]

ExactRandomGraph[n_Integer,e_Integer] :=
	FromUnorderedPairs[
		Map[ NthPair, Take[ RandomPermutation[n(n-1)/2], e] ],
		CircularVertices[n]
	]

NthPair[0] := {}
NthPair[n_Integer] :=
	Module[{i=2},
		While[ Binomial[i,2] < n, i++];
		{n - Binomial[i-1,2], i}
	]

RandomVertices[n_Integer] := Table[{Random[], Random[]}, {n}]
RandomVertices[g_Graph] := Graph[ Edges[g], RandomVertices[V[g]] ]

RandomGraph[n_Integer,p_,range_List,Directed] :=
	RemoveSelfLoops[
		Graph[
			Table[If[Random[Real]<p,Random[Integer,range],0],{n},{n}],
			CircularVertices[n]
		]
	]

RandomGraph[n_Integer,p_,Directed] := RandomGraph[n,p,{1,1},Directed]

DegreeSequence[g_Graph] := Reverse[ Sort[ Degrees[g] ] ]

Degrees[Graph[g_,_]] := Map[(Apply[Plus,#])&, g]

GraphicQ[s_List] := False /; (Min[s] < 0) || (Max[s] >= Length[s])
GraphicQ[s_List] := (First[s] == 0) /; (Length[s] == 1)
GraphicQ[s_List] :=
	Module[{m,sorted = Reverse[Sort[s]]},
		m = First[sorted];
		GraphicQ[ Join[ Take[sorted,{2,m+1}]-1, Drop[sorted,m+1] ] ]
	]

RealizeDegreeSequence[d_List] :=
	Module[{i,j,v,set,seq,n=Length[d],e},
		seq = Reverse[ Sort[ Table[{d[[i]],i},{i,n}]] ];
		FromUnorderedPairs[
			Flatten[ Table[
				{{k,v},seq} = {First[seq],Rest[seq]};
				While[ !GraphicQ[
					MapAt[
						(# - 1)&,
						Map[First,seq],
						set = RandomKSubset[Table[{i},{i,n-j}],k] 
					] ],
				        Null
				];
				e = Map[(Prepend[seq[[#,2]],v])&,set];
				seq = Reverse[ Sort[
					MapAt[({#[[1]]-1,#[[2]]})&,seq,set]
				] ];
				e,
				{j,Length[d]-1}
			], 1],
			CircularVertices[n]
		]
	] /; GraphicQ[d]

RealizeDegreeSequence[d_List,seed_Integer] :=
	(SeedRandom[seed]; RealizeDegreeSequence[d])

RegularQ[Graph[g_,_]] := Apply[ Equal, Map[(Apply[Plus,#])& , g] ]

RegularGraph[k_Integer,n_Integer] := RealizeDegreeSequence[Table[k,{n}]]

MakeGraph[v_List,f_] :=
	Module[{n=Length[v],i,j},
		Graph [
			Table[If [Apply[f,{v[[i]],v[[j]]}], 1, 0],{i,n},{j,n}],
			CircularVertices[n]
		]
	]

IntervalGraph[l_List] :=
	MakeGraph[
		l,
		( ((First[#1] <= First[#2]) && (Last[#1] >= First[#2])) ||
		  ((First[#2] <= First[#1]) && (Last[#2] >= First[#1])) )&
	]

FunctionalGraph[f_,n_] :=
	Module[{i,x},
		FromOrderedPairs[
			Table[{i, x=Mod[Apply[f,{i}],n]; If[x!=0,x,n]}, {i,n} ],
			CircularVertices[n]
		]
	]

ConnectedComponents[g_Graph] :=
	Module[{untraversed=Range[V[g]],traversed,comps={}},
		While[untraversed != {},
			traversed = DepthFirstTraversal[g,First[untraversed]];
			AppendTo[comps,traversed];
			untraversed = Complement[untraversed,traversed]
		];
		comps
	]

ConnectedQ[g_Graph] := 
	(
	Length[ DepthFirstTraversal[g,1] ] == V[g]
	)

WeaklyConnectedComponents[g_Graph] := ConnectedComponents[ MakeUndirected[g] ]

ConnectedQ[g_Graph,Undirected] := Length[ WeaklyConnectedComponents[g] ] == 1

StronglyConnectedComponents[g_Graph] :=
	Block[{e=ToAdjacencyLists[g],s,c=1,i,cur={},low=dfs=Table[0,{V[g]}],scc={}},
		While[(s=Select[Range[V[g]],(dfs[[#]]==0)&]) != {},
			SearchStrongComp[First[s]];
		];
		scc
	]

SearchStrongComp[v_Integer] :=
	Block[{r},
		low[[v]]=dfs[[v]]=c++;
		PrependTo[cur,v];
		Scan[
			(If[dfs[[#]] == 0,
				SearchStrongComp[#];
				low[[v]]=Min[low[[v]],low[[#]]],
				If[(dfs[[#]] < dfs[[v]]) && MemberQ[cur,#],
					low[[v]]=Min[low[[v]],dfs[[#]] ]
				];
			])&,
			e[[v]]
		];
		If[low[[v]] == dfs[[v]],
			{r} = Flatten[Position[cur,v]];
			AppendTo[scc,Take[cur,r]];
			cur = Drop[cur,r];
		];
	]

ConnectedQ[g_Graph,Directed] := Length[ StronglyConnectedComponents[g] ] == 1

OrientGraph[g_Graph] :=
	Module[{pairs,newg,rest,cc,c,i,e},
		pairs = Flatten[Map[(Partition[#,2,1])&,ExtractCycles[g]],1];
		newg = FromUnorderedPairs[pairs,Vertices[g]];
		rest = ToOrderedPairs[ GraphDifference[ g, newg ] ];
		cc = Sort[ConnectedComponents[newg], (Length[#1]>=Length[#2])&];
		c = First[cc];
		Do[
			e = Select[rest,(MemberQ[c,#[[1]]] &&
					 MemberQ[cc[[i]],#[[2]]])&];
			rest = Complement[rest,e,Map[Reverse,e]];
			c = Union[c,cc[[i]]];
			pairs = Join[pairs, Prepend[ Rest[e],Reverse[e[[1]]] ] ],
			{i,2,Length[cc]}
		];
		FromOrderedPairs[
			Join[pairs, Select[rest,(#[[1]] > #[[2]])&] ],
			Vertices[g]
		]
	] /; SameQ[Bridges[g],{}]

FindBiconnectedComponents[g_Graph] :=
	Block[{e=ToAdjacencyLists[g],n=V[g],par,c=0,act={},back,dfs,ap=bcc={}},
		back=dfs=Table[0,{n}];
		par = Table[n+1,{n}]; 
		Map[(SearchBiConComp[First[#]])&, ConnectedComponents[g]];
		{bcc,Drop[ap, -1]}
	]

SearchBiConComp[v_Integer] :=
	Block[{r},
		back[[v]]=dfs[[v]]=++c;
		Scan[
			(If[ dfs[[#]] == 0, 
				If[!MemberQ[act,{v,#}], PrependTo[act,{v,#}]];
				par[[#]] = v;
				SearchBiConComp[#];
				If[ back[[#]] >= dfs[[v]],
					{r} = Flatten[Position[act,{v,#}]];
					AppendTo[bcc,Union[Flatten[Take[act,r]]]];
					AppendTo[ap,v];
					act = Drop[act,r]
				];
				back[[v]] = Min[ back[[v]],back[[#]] ],
				If[# != par[[v]],back[[v]]=Min[dfs[[#]],back[[v]]]]
			])&,
			e[[v]]
		];
	]

ArticulationVertices[g_Graph]  := Union[Last[FindBiconnectedComponents[g]]];

Bridges[g_Graph] := Select[BiconnectedComponents[g],(Length[#] == 2)&]

BiconnectedComponents[g_Graph] := First[FindBiconnectedComponents[g]];

BiconnectedQ[g_Graph] := Length[ BiconnectedComponents[g] ] == 1

EdgeConnectivity[g_Graph] :=
	Module[{i},
		Apply[Min, Table[NetworkFlow[g,1,i], {i,2,V[g]}]]
	]

VertexConnectivityGraph[g_Graph] :=
	Module[{n=V[g],e},
		e=Table[0,{2 n},{2 n}];
		Scan[ (e[[#-1,#]] = 1)&, 2 Range[n] ];
		Scan[
			(e[[#[[1]], #[[2]]-1]] = e[[#[[2]],#[[1]]-1]] = Infinity)&,
			2 ToUnorderedPairs[g]
		];
		Graph[e,Apply[Join,Map[({#,#})&,Vertices[g]]]]
	]

VertexConnectivity[g_Graph] :=
	Module[{p=VertexConnectivityGraph[g],k=V[g],i=0,notedges},
		notedges = ToUnorderedPairs[ GraphComplement[g] ];
		While[ i++ <= k,
			k = Min[
				Map[
					(NetworkFlow[p,2 #[[1]],2 #[[2]]-1])&,
					Select[notedges,(First[#]==i)&]
				],
				k
			]
		];
		k
	]

Harary[k_?EvenQ, n_Integer] := CirculantGraph[n,Range[k/2]]

Harary[k_?OddQ, n_?EvenQ] := CirculantGraph[n,Append[Range[k/2],n/2]]

Harary[k_?OddQ, n_?OddQ] :=
	Module[{g=Harary[k-1,n],i},
		FromUnorderedPairs[
			Join[
				ToUnorderedPairs[g],
				{ {1,(n+1)/2}, {1,(n+3)/2} },
				Table [ {i,i+(n+1)/2}, {i,2,(n-1)/2} ]
			],
			Vertices[g]
		]
	]

IdenticalQ[g_Graph,h_Graph] := Edges[g] === Edges[h]

IsomorphismQ[g_Graph,h_Graph,p_List] := False	/;
		(V[g]!=V[h]) || !PermutationQ[p] || (Length[p] != V[g])

IsomorphismQ[g_Graph,h_Graph,p_List] := IdenticalQ[g, InduceSubgraph[h,p] ]

Isomorphism[g_Graph,h_Graph,flag_:One] := {}	/; (V[g] != V[h]) 

Isomorphism[g_Graph,h_Graph,flag_:One] :=
	Module[{eg=Edges[g],eh=Edges[h],equiv=Equivalences[g,h]},
		If [!MemberQ[equiv,{}],
			Backtrack[
				equiv,
				(IdenticalQ[InduceSubgraph[g,Range[Length[#]]],
				    	InduceSubgraph[h,#] ] &&
			 	!MemberQ[Drop[#,-1],Last[#]])&,
				(IsomorphismQ[g,h,#])&,
				flag
			],
			{}
		]
	]

IsomorphicQ[g_Graph,h_Graph] := True /; IdenticalQ[g,h]
IsomorphicQ[g_Graph,h_Graph] := ! SameQ[ Isomorphism[g,h], {}]

Equivalences[g_Graph,h_Graph] :=
	Equivalences[ AllPairsShortestPath[g], AllPairsShortestPath[h]]

Equivalences[g_List,h_List] :=
	Module[{dg=Map[Sort,g],dh=Map[Sort,h],s,i},
		Table[
			Flatten[Position[dh,_?(Function[s,SameQ[s,dg[[i]] ]])]],
			{i,Length[dg]}
		]
	] /; Length[g] == Length[h]

Automorphisms[g_Graph,flag_:All] :=
	Module[{s=AllPairsShortestPath[g]},
		Backtrack[
			Equivalences[s,s],
			(IdenticalQ[InduceSubgraph[g,Range[Length[#]]],
				    InduceSubgraph[g,#] ] &&
			 !MemberQ[Drop[#,-1],Last[#]])&,
			(IsomorphismQ[g,g,#])&,
			flag
		]
	]

SelfComplementaryQ[g_Graph] := IsomorphicQ[g, GraphComplement[g]]

FindCycle[g_Graph,flag_:Undirected] :=
     Module[{edge,n=V[g],x,queue,v,seen,parent},
       edge=ToAdjacencyLists[g];
       For[ v = 1, v <= n, v++,
           parent=Table[n+1,{n}]; parent[[v]] = 0;
           seen = {}; queue = {v};
           While[ queue != {},
               {x,queue} = {First[queue], Rest[queue]};
               AppendTo[seen,x];
               If[ SameQ[ flag, Undirected],
                   Scan[ (If[ parent[[x]] != #, parent[[#]]=x])&, edge[[x]] ],
                   Scan[ (parent[[#]]=x)&, edge[[x]]]
               ];
               If[ SameQ[flag,Undirected],
                   If[ MemberQ[ edge[[x]],v ] && parent[[x]] != v,
                       Return[ FromParent[parent,x] ]
                   ],
                   If[ MemberQ[ edge[[x]],v ],
                       Return[ FromParent[parent,x] ]
                   ]
               ];
               queue = Join[ Complement[ edge[[x]], seen], queue]
           ]
       ];
     {}
     ]

FromParent[parent_List,s_Integer] :=
	Module[{i=s,lst={s}},
		While[!MemberQ[lst,(i=parent[[i]])], PrependTo[lst,i] ];
		PrependTo[lst,i];
		Take[lst, Flatten[Position[lst,i]]]
	]

AcyclicQ[g_Graph,flag_:Undirected] := SameQ[FindCycle[g,flag],{}]

TreeQ[g_Graph] := ConnectedQ[g] && (M[g] == V[g]-1)

ExtractCycles[gi_Graph,flag_:Undirected] := 
	Module[{g=gi,cycles={},c},
		While[!SameQ[{}, c=FindCycle[g,flag]],
			PrependTo[cycles,c];
			g = DeleteCycle[g,c,flag];
		];
		cycles
	]

DeleteCycle[g_Graph,cycle_List,flag_:Undirected] :=
	Module[{newg=g},
		Scan[(newg=DeleteEdge[newg,#,flag])&, Partition[cycle,2,1] ];
		newg
	]

Girth[g_Graph] := 
	Module[{v,dist,queue,n=V[g],girth=Infinity,parent,e=ToAdjacencyLists[g],x},
		Do [
			dist = parent = Table[Infinity, {n}];
			dist[[v]] = parent[[v]] = 0;
			queue = {v};
			While [queue != {},
				{x,queue} = {First[queue],Rest[queue]};
				Scan[
					(If [ (dist[[#]]+dist[[x]]<girth) &&
				     	      (parent[[x]] != #),
						girth=dist[[#]]+dist[[x]] + 1,
				 	 If [dist[[#]]==Infinity,
						dist[[#]] = dist[[x]] + 1;
						parent[[#]] = x;
						If [2 dist[[#]] < girth-1,
							AppendTo[queue,#] ]
					]])&,
					e[[ x ]]
				];
			],
			{v,n}
		];
		girth
	] /; SimpleQ[g]

EulerianQ[g_Graph,Directed] :=
	ConnectedQ[g,Undirected] && (InDegree[g] === OutDegree[g])

EulerianQ[g_Graph,flag_:Undirected] := ConnectedQ[g,Undirected] && 
	UndirectedQ[g] && Apply[And,Map[EvenQ,DegreeSequence[g]]]

OutDegree[Graph[e_List,_],n_Integer] := Length[ Select[ e[[n]], (# != 0)& ] ]
OutDegree[g_Graph] := Map[ (OutDegree[g,#])&, Range[V[g]] ]

InDegree[g_Graph,n_Integer] := OutDegree[ TransposeGraph[g], n ];
InDegree[g_Graph] := Map[ (InDegree[g,#])&, Range[V[g]] ]

TransposeGraph[Graph[g_List,v_List]] := Graph[ Transpose[g], v ]

EulerianCycle[g_Graph,flag_:Undirected] :=
	Module[{euler,c,cycles,v},
		cycles = Map[(Drop[#,-1])&, ExtractCycles[g,flag]];
		{euler, cycles} = {First[cycles], Rest[cycles]};
		Do [
			c = First[ Select[cycles, (Intersection[euler,#]=!={})&] ];
			v = First[Intersection[euler,c]];
			euler = Join[
				RotateLeft[c, Position[c,v] [[1,1]] ],
				RotateLeft[euler, Position[euler,v] [[1,1]] ]
			];
			cycles = Complement[cycles,{c}],
			{Length[cycles]}
		];
		Append[euler, First[euler]]
	] /; EulerianQ[g,flag]

DeBruijnSequence[alph_List,n_Integer] :=
        Module[{states = Strings[alph,n-1]},
                Rest[ Map[
                        (First[ states[[#]] ])&,
                        EulerianCycle[
                                MakeGraph[
                                        states,
                                        (Module[{i},
                                         MemberQ[
                                                Table[
                                                        Append[Rest[#1],alph[[i]]],
                                                        {i,Length[alph]}
                                                ],
                                                #2
                                         ]
                                        ])&
                                ],
                                Directed
                        ]
                ] ]
        ] /; n>=2

DeBruijnSequence[alph_List,n_Integer] := alph /; n==1

HamiltonianQ[g_Graph] := False /; !BiconnectedQ[g]
HamiltonianQ[g_Graph] := HamiltonianCycle[g] != {}

HamiltonianCycle[g_Graph,flag_:One] :=
	Module[{s={1},all={},done,adj=Edges[g],e=ToAdjacencyLists[g],x,v,ind,n=V[g]},
		ind=Table[1,{n}];
		While[ Length[s] > 0,
			v = Last[s];
			done = False;
			While[ ind[[v]] <= Length[e[[v]]] && !done,
				If[!MemberQ[s,(x = e[[v,ind[[v]]++]])], done=True]
			];
			If[ done, AppendTo[s,x], s=Drop[s,-1]; ind[[v]] = 1];
			If[(Length[s] == n),
				If [(adj[[x,1]]>0),
					AppendTo[all,Append[s,First[s]]];
					If [SameQ[flag,All],
						s=Drop[s,-1],
						all = Flatten[all]; s={}
					],
					s = Drop[s,-1]
				]
			]
		];
		all
	]

TravelingSalesman[g_Graph] :=
	Module[{v,s={1},sol={},done,cost,g1,e=ToAdjacencyLists[g],x,ind,best,n=V[g]},
		ind=Table[1,{n}];
		g1 = PathConditionGraph[g];
		best = Infinity;
		While[ Length[s] > 0,
			v = Last[s];
			done = False;
			While[ ind[[v]] <= Length[e[[v]]] && !done,
				x = e[[v,ind[[v]]++]];
				done = (best > CostOfPath[g1,Append[s,x]]) &&
					!MemberQ[s,x]
			];
			If[done, AppendTo[s,x], s=Drop[s,-1]; ind[[v]] = 1];
			If[(Length[s] == n),
				cost = CostOfPath[g1, Append[s,First[s]]];
				If [(cost < best), sol = s; best = cost ];
				s = Drop[s,-1]
			]
		];
		Append[sol,First[sol]]
	]

CostOfPath[Graph[g_,_],p_List] := Apply[Plus, Map[(Element[g,#])&,Partition[p,2,1]] ]

Element[a_List,{index___}] := a[[ index ]]

TriangleInequalityQ[e_?SquareMatrixQ] :=
	Module[{i,j,k,n=Length[e],flag=True},
		Do [

			If[(e[[i,k]]!=0) && (e[[k,j]]!=0) && (e[[i,j]]!=0),
				If[e[[i,k]]+e[[k,j]]<e[[i,j]],
					flag = False;
				]
			],
			{i,n},{j,n},{k,n}
		];
		flag
	]

TriangleInequalityQ[g_Graph] := TriangleInequalityQ[Edges[g]]

TravelingSalesmanBounds[g_Graph] := {LowerBoundTSP[g], UpperBoundTSP[g]}

UpperBoundTSP[g_Graph] :=
	CostOfPath[g, Append[DepthFirstTraversal[MinimumSpanningTree[g],1],1]]

LowerBoundTSP[g_Graph] := Apply[Plus, Map[Min,ReplaceAll[Edges[g],0->Infinity]]]

PartialOrderQ[g_Graph] := ReflexiveQ[g] && AntiSymmetricQ[g] && TransitiveQ[g]

TransitiveQ[g_Graph] := IdenticalQ[g,TransitiveClosure[g]]

ReflexiveQ[Graph[g_List,_]] := 
	Module[{i},
		Apply[And, Table[(g[[i,i]]!=0),{i,Length[g]}] ]
	]

AntiSymmetricQ[g_Graph] := 
	Module[{e = Edges[g], g1 = RemoveSelfLoops[g]},
		Apply[And, Map[(Element[e,Reverse[#]]==0)&,ToOrderedPairs[g1]] ]
	]

TransitiveClosure[g_Graph] :=
	Module[{i,j,k,e=Edges[g],n=V[g]},
		Do [
			If[ e[[j,i]] != 0,
				Do [
					If[ e[[i,k]] != 0, e[[j,k]]=1],
					{k,n}
				]
			],
			{i,n},{j,n}
		];
		Graph[e,Vertices[g]]
	]

TransitiveReduction[g_Graph] :=
	Module[{closure=reduction=Edges[g],i,j,k,n=V[g]},
		Do[
			If[ closure[[i,j]]!=0 && closure[[j,k]]!=0 &&
				 reduction[[i,k]]!=0 && (i!=j) && (j!=k) && (i!=k),
					reduction[[i,k]] = 0
			],
			{i,n},{j,n},{k,n}
		];
		Graph[reduction,Vertices[g]]
	] /; AcyclicQ[RemoveSelfLoops[g],Directed] 

TransitiveReduction[g_Graph] :=
	Module[{reduction=Edges[g],i,j,k,n=V[g]},
		Do[
			If[ reduction[[i,j]]!=0 && reduction[[j,k]]!=0 &&
				 reduction[[i,k]]!=0 && (i!=j) && (j!=k) && (i!=k),
					reduction[[i,k]] = 0
			],
			{i,n},{j,n},{k,n}
		];
		Graph[reduction,Vertices[g]]
	] 

(*thanks Christoph Strnadl*)

HasseDiagram[g_,fak_:1] :=
	Module[{r, rank, m, stages, freq=Table[0,{V[g]}], 
             adjm, first},
		r = TransitiveReduction[ RemoveSelfLoops[g]  ];
            adjm = ToAdjacencyLists[r];
            rank = Table[ 0,{ V[g]} ];                              
            first = Select[ Range[ V[g]], InDegree[r,#]==0& ];
            rank = MakeLevel[ first, 1, adjm, rank];          
            first = Max[rank];
            stages = Distribution[ rank ];                       
            Graph[
		Edges[r],
		Table[
			m = ++ freq[[ rank[[i]] ]];
			{ ((m-1) + (1-stages[[rank[[i]] ]])/2) fak^(first-rank[[i]]),
                          rank[[i]] },
				{i, V[g]}
			]
		]
	] /; AcyclicQ[ RemoveSelfLoops[g],Directed ]

(*
 * SetLevel[{p1,p2,...},lvl,rank] sets the positions p1, p2,.. of
 * list rank to the level lvl, if the old entry at that position
 * is less than level.
 *)
SetLevel[l_List,lvl_,rank_List] :=  
    Module[ {r=rank},
            If[ r[[#]] < lvl, r[[#]] = lvl ] & /@ l;
            r
    ]

(*
 * MakeLevel[l,level,adjm,rank] constructs recursively the ranks of
 * each vertex according to the adjacency matrix adjm of the graph.
 * rank is the current ranking, level the new level to assign and 
 * l = {v1,v2,..} the list of vertices to be set to level.
 *)
MakeLevel[{},_,_,rank_] := rank

MakeLevel[l_List,lvl_,adjm_List,r_List] :=
  Module[ {rank=r, v, lst=l }, 
          rank = SetLevel[lst,lvl,rank];  (* make this level ready *)
          While[ lst != {},
                 v = First[lst];
                 rank = MakeLevel[adjm[[v]], lvl+1,adjm,rank]; 
                 lst = Rest[lst];
          ];                            
          rank
  ]

(*
 * HasseDiagram[g] renders a graph corresponding to the HasseDiagram of
 * the partial order induced by the directed graph g. 
 * HasseDiagram[g,fac] renders the HasseDiagram in which each vertex'
 * position is stretched by factor fac. In each stage that factor
 * is taken to the power of the distance to the 1 element.
 *)                                 
 

TopologicalSort[g_Graph] :=
	Module[{g1 = RemoveSelfLoops[g],e,indeg,zeros,v},
		e=ToAdjacencyLists[g1];
		indeg=InDegree[g1];
		zeros = Flatten[ Position[indeg, 0] ];
		Table [
			{v,zeros}={First[zeros],Rest[zeros]};
			Scan[
				( indeg[[#]]--;
				  If[indeg[[#]]==0, AppendTo[zeros,#]] )&,
				e[[ v ]]
			];
			v,
			{V[g]}
		]
	] /; AcyclicQ[RemoveSelfLoops[g],Directed]

ChromaticPolynomial[g_Graph,z_]:=
	Module[{simpleGraph},
		simpleGraph = MakeSimple[g];
		SimpleChromaticPolynomial[simpleGraph,z]
	]
 
SimpleChromaticPolynomial[g_Graph,z_] := 
	(
	0
	) /; Identical[g,CompleteGraph[0]]

SimpleChromaticPolynomial[g_Graph,z_] :=
	Module[{i},
	  Product[z-i, {i,0,V[g]-1}]
	] /; CompleteQ[g]

SimpleChromaticPolynomial[g_Graph,z_] :=
	(
	 z ( z - 1 ) ^ (V[g]-1)
	) /; TreeQ[g]

SimpleChromaticPolynomial[g_Graph,z_] :=
	(
	If [M[g]>Binomial[V[g],2]/2, ChromaticDense[g,z], ChromaticSparse[g,z]]
	)

ChromaticSparse[g_Graph,z_] := 
	(
	z^V[g] 
	) /; EmptyQ[g]

ChromaticSparse[g_Graph,z_] :=
	Module[{i=1, v, e=Edges[g], none=Table[0,{V[g]}]},
        	While[e[[i]] === none, i++];
        	v = Position[e[[i]],1] [[1,1]];
		ChromaticSparse[ DeleteEdge[g,{i,v}], z ] -
			ChromaticSparse[ Contract[g,{i,v}], z ]
	]

ChromaticDense[g_Graph,z_] := 
	(
	ChromaticPolynomial[g,z]
	) /; CompleteQ[g]

ChromaticDense[g_Graph,z_] :=
	Module[{i=1, v, e=Edges[g], all=Join[Table[1,{V[g]-1}],{0}] },
		While[e[[i]] === RotateRight[all,i], i++];
		v = Last[ Position[e[[i]],0] ] [[1]];
		ChromaticDense[ AddEdge[g,{i,v}], z ] +
			ChromaticDense[ Contract[g,{i,v}], z ]
	]

ChromaticNumber[g_Graph] :=
	Block[{ways, z},
		ways[z_] = ChromaticPolynomial[g,z];
		For [z=0, z<=V[g], z++,
			If [ways[z] > 0, Return[z]]
		]
	]

TwoColoring[g_Graph] := 
	Module[{queue,elem,edges,col,flag=True,colored=Table[0,{V[g]}]},
		edges = ToAdjacencyLists[g];
		While[ MemberQ[colored,0],
			queue = First[ Position[colored,0] ];
			colored[[ First[queue] ]] = 1;
			While[ queue != {},
				elem = First[queue];
				col = colored[[elem]];
				Scan[
					(Switch[colored[[ # ]],
						col, flag = False,
						0, AppendTo[queue, # ];
						   colored[[#]] = Mod[col,2]+1
					])&,
					edges[[elem]]
				];
				queue = Rest[queue];
			]
		];
		If [!flag, colored[[1]] = 0];
		colored
	]

BipartiteQ[g_Graph] := ! MemberQ[ TwoColoring[g], 0 ]

VertexColoring[g_Graph] :=
	Module[{v,l,n=V[g],e=ToAdjacencyLists[g],x,color=Table[0,{V[g]}]},
		v = Map[(Apply[Plus,#])&, Edges[g]];
		Do[
			l = MaximumColorDegreeVertices[e,color];
			x = First[l];
			Scan[(If[ v[[#]] > v[[x]], x = #])&, l];
			color[[x]] = Min[
				Complement[ Range[n], color[[ e[[x]] ]] ]
			],
			{V[g]}
		];
		color
	]

MaximumColorDegreeVertices[e_List,color_List] :=
	Module[{n=Length[color],l,i,x},
		l = Table[ Count[e[[i]], _?(Function[x,color[[x]]!=0])], {i,n}];
		Do [ 
			If [color[[i]]!=0, l[[i]] = -1],
			{i,n}
		];
		Flatten[ Position[ l, Max[l] ] ]
	]

EdgeColoring[g_Graph] := VertexColoring[ LineGraph[g] ]

EdgeChromaticNumber[g_Graph] := ChromaticNumber[ LineGraph[g] ]

CliqueQ[g_Graph,clique_List] :=
	IdenticalQ[ CompleteGraph[Length[clique]], InduceSubgraph[g,clique] ] /; SimpleQ[g]

MaximumClique[g_Graph] := {} /; g === CompleteGraph[0]

MaximumClique[g_Graph] :=
	Module[{d = Degrees[g],i,clique=Null,k},
		i = Max[d];
		While[(SameQ[clique,Null]),
			k = CompleteGraph[i+1];
			clique = FirstExample[
				KSubsets[Flatten[Position[d,_?((#>=i)&)]], i+1],
				(IdenticalQ[k,InduceSubgraph[g,#]])&
			];
			i--;
		];
		clique
	]

FirstExample[list_List, predicate_] := Scan[(If [predicate[#],Return[#]])&,list]

VertexCoverQ[g_Graph,vc_List] :=
	CliqueQ[ GraphComplement[g], Complement[Range[V[g]], vc] ]

MinimumVertexCover[g_Graph] :=
	Complement[ Range[V[g]], MaximumClique[ GraphComplement[g] ] ]

IndependentSetQ[g_Graph,indep_List] :=
	VertexCoverQ[ g, Complement[ Range[V[g]], indep] ]

MaximumIndependentSet[g_Graph] := Complement[Range[V[g]], MinimumVertexCover[g]]

PerfectQ[g_Graph] :=
	Apply[
		And,
		Map[(ChromaticNumber[#] == Length[MaximumClique[#]])&,
			Map[(InduceSubgraph[g,#])&, Subsets[Range[V[g]]] ] ]
	]

Dijkstra[g_Graph,start_Integer] := First[ Dijkstra[g,{start}] ]

Dijkstra[g_Graph, l_List] :=
	Module[{x,start,e=ToAdjacencyLists[g],i,p,parent,untraversed},
		p=Edges[PathConditionGraph[g]];
		Table[
			start = l[[i]];
			parent=untraversed=Range[V[g]];
			dist = p[[start]]; dist[[start]] = 0;
			Scan[ (parent[[#]] = start)&, e[[start]] ];
			While[ untraversed != {} ,
				x = First[untraversed];
				Scan[(If [dist[[#]]<dist[[x]],x=#])&, untraversed];
				untraversed = Complement[untraversed,{x}];
				Scan[
					(If[dist[[#]] > dist[[x]]+p[[x,#]],
						dist[[#]] = dist[[x]]+p[[x,#]];
						parent[[#]] = x ])&,
					e[[x]]
				];
			];
			{parent, dist},
			{i,Length[l]}
		]
	]

ShortestPath[g_Graph,s_Integer,e_Integer] := 
	Module[{parent=First[Dijkstra[g,s]],i=e,lst={e}},
		While[ (i != s) && (i != parent[[i]]),
			PrependTo[lst,parent[[i]]];
			i = parent[[i]]
		];
		If[ i == s, lst, {}]
	]

ShortestPathSpanningTree[g_Graph,s_Integer] :=
	Module[{parent=First[Dijkstra[g,s]],i},
		FromUnorderedPairs[
			Map[({#,parent[[#]]})&, Complement[Range[V[g]],{s}]],
			Vertices[g]
		]
	]

AllPairsShortestPath[g_Graph] :=
	Module[{p=Edges[ PathConditionGraph[g] ],i,j,k,n=V[g]},
		Do [
			p = Table[Min[p[[i,k]]+p[[k,j]],p[[i,j]]],{i,n},{j,n}],
			{k,n}
		];
		p
	] /; Min[Edges[g]] < 0

AllPairsShortestPath[g_Graph] := Map[ Last, Dijkstra[g, Range[V[g]]]]

PathConditionGraph[Graph[e_,v_]] := RemoveSelfLoops[Graph[ReplaceAll[e,0->Infinity],v]]

GraphPower[g_Graph,1] := g

GraphPower[g_Graph,n_Integer] :=
	Module[{prod=power=p=Edges[g]},
		Do [
			prod = prod . p;
			power = prod + power,
			{n-1}
		];
		Graph[power, Vertices[g]]
	]

InitializeUnionFind[n_Integer] := Module[{i}, Table[{i,1},{i,n}] ]

FindSet[n_Integer,s_List] := If [n == s[[n,1]], n, FindSet[s[[n,1]],s] ]

UnionSet[a_Integer,b_Integer,s_List] :=
	Module[{sa=FindSet[a,s], sb=FindSet[b,s], set=s},
		If[ set[[sa,2]] < set[[sb,2]], {sa,sb} = {sb,sa} ];
		set[[sa]] = {sa, Max[ set[[sa,2]], set[[sb,2]]+1 ]};
		set[[sb]] = {sa, set[[sb,2]]};
		set
	]

MinimumSpanningTree[g_Graph] :=
	Module[{edges=Edges[g],set=InitializeUnionFind[V[g]]},
		FromUnorderedPairs[
			Select [
				Sort[
					ToUnorderedPairs[g],
					(Element[edges,#1]<=Element[edges,#2])&
				],
				(If [FindSet[#[[1]],set] != FindSet[#[[2]],set],
					set=UnionSet[#[[1]],#[[2]],set]; True,
					False
				])&
			],
			Vertices[g]
		]
	] /; UndirectedQ[g]

MaximumSpanningTree[g_Graph] := MinimumSpanningTree[Graph[-Edges[g],Vertices[g]]]

Cofactor[m_List,{i_Integer,j_Integer}] :=
	(-1)^(i+j) * Det[ Drop[ Transpose[ Drop[Transpose[m],{j,j}] ], {i,i}] ]

NumberOfSpanningTrees[Graph[g_List,_]] :=
	Cofactor[ DiagonalMatrix[Map[(Apply[Plus,#])&,g]] - g, {1,1}]

NetworkFlow[g_Graph,source_Integer,sink_Integer] :=
	Block[{flow=NetworkFlowEdges[g,source,sink], i},
		Sum[flow[[i,sink]], {i,V[g]}]
	]

NetworkFlowEdges[g_Graph,source_Integer,sink_Integer] :=
	Block[{e=Edges[g], x, y, flow=Table[0,{V[g]},{V[g]}], p, m},
		While[ !SameQ[p=AugmentingPath[g,source,sink], {}],
			m = Min[Map[({x,y}=#[[1]]; 
				 If[SameQ[#[[2]],f],e[[x,y]]-flow[[x,y]],
					flow[[x,y]]])&,p]];
			Scan[	
				({x,y}=#[[1]];
				 If[ SameQ[#[[2]],f],
					flow[[x,y]]+=m,flow[[x,y]]-=m])&,
				 p
			]
		];
		flow
	]

AugmentingPath[g_Graph,src_Integer,sink_Integer] :=
	Block[{l={src},lab=Table[0,{V[g]}],v,c=Edges[g],e=ToAdjacencyLists[g]},
		lab[[src]] = start;
		While[l != {} && (lab[[sink]]==0),
			{v,l} = {First[l],Rest[l]};
			Scan[ (If[ c[[v,#]] - flow[[v,#]] > 0 && lab[[#]] == 0,
				lab[[#]] = {v,f}; AppendTo[l,#]])&,
				e[[v]]
			];
			Scan[ (If[ flow[[#,v]] > 0 && lab[[#]] == 0,
				lab[[#]] = {v,b}; AppendTo[l,#]] )&,
				Select[Range[V[g]],(c[[#,v]] > 0)&]
			];
		];
		FindPath[lab,src,sink]
	]

FindPath[l_List,v1_Integer,v2_Integer] :=
	Block[{x=l[[v2]],y,z=v2,lst={}},
		If[SameQ[x,0], Return[{}]];
		While[!SameQ[x, start],
			If[ SameQ[x[[2]],f],
				PrependTo[lst,{{ x[[1]], z }, f}],
				PrependTo[lst,{{ z, x[[1]] }, b}]
			];
			z = x[[1]]; x = l[[z]];
		];
		lst
	]

BipartiteMatching[g_Graph] :=
	Module[{p,v1,v2,coloring=TwoColoring[g],n=V[g]},
		v1 = Flatten[Position[coloring,1]];
		v2 = Flatten[Position[coloring,2]];
		p = BipartiteMatchingFlowGraph[g,v1,v2];
		flow = NetworkFlowEdges[p,V[g]+1,V[g]+2];
		Select[ToOrderedPairs[Graph[flow,Vertices[p]]], (Max[#]<=n)&]
	] /; BipartiteQ[g]

BipartiteMatchingFlowGraph[g_Graph,v1_List,v2_List] :=
	Module[{edges = Table[0,{V[g]+2},{V[g]+2}],i,e=ToAdjacencyLists[g]},
		Do[ 
	    		Scan[ (edges[[v1[[i]],#]] = 1)&, e[[ v1[[i]] ]] ],
			{i,Length[v1]}
		];
		Scan[(edges[[V[g] + 1, #]] = 1)&, v1];
		Scan[(edges[[#, V[g] + 2]] = 1)&, v2];
		Graph[edges,RandomVertices[V[g] + 2] ]
	]

MinimumChainPartition[g_Graph] :=
	ConnectedComponents[
		FromUnorderedPairs[
			Map[(#-{0,V[g]})&, BipartiteMatching[DilworthGraph[g]]],
			Vertices[g]
		]
	]

MaximumAntichain[g_Graph] := MaximumIndependentSet[TransitiveClosure[g]]

DilworthGraph[g_Graph] :=
	FromUnorderedPairs[
		Map[
			(#+{0,V[g]})&,
			ToOrderedPairs[RemoveSelfLoops[TransitiveReduction[g]]]
		]
	]

MaximalMatching[g_Graph] :=
	Module[{match={}},
		Scan[
			(If [Intersection[#,match]=={}, match=Join[match,#]])&,
			ToUnorderedPairs[g]
		];
		Partition[match,2]
	]

StableMarriage[mpref_List,fpref_List] :=
	Module[{n=Length[mpref],freemen,cur,i,w,husband},
		freemen = Range[n];
		cur = Table[1,{n}];
		husband = Table[n+1,{n}];
		While[ freemen != {},
			{i,freemen}={First[freemen],Rest[freemen]};
			w = mpref[[ i,cur[[i]] ]];
			If[BeforeQ[ fpref[[w]], i, husband[[w]] ], 
				If[husband[[w]] != n+1,
					AppendTo[freemen,husband[[w]] ]
				];
				husband[[w]] = i,
				cur[[i]]++;
				AppendTo[freemen,i]
			];
		];
		InversePermutation[ husband ]
	] /; Length[mpref] == Length[fpref]

BeforeQ[l_List,a_,b_] :=
	If [First[l]==a, True, If [First[l]==b, False, BeforeQ[Rest[l],a,b] ] ]
 
PlanarQ[g_Graph] :=
	Module [{simpleGraph = MakeSimple[g]},
	  SimplePlanarQ[simpleGraph]
	]

SimplePlanarQ[g_Graph] :=
  Module[{components},
	components = ConnectedComponents[g];
	Apply[  And,
		Map[(
		   PlanarQ[InduceSubgraph[g,#]])&, components]
	]
  ] /; !ConnectedQ[g]

SimplePlanarQ[g_Graph] := False /;  (M[g] > 3 V[g]-6) && (V[g] > 2)
SimplePlanarQ[g_Graph] := True /;   (M[g] < V[g] + 3)
SimplePlanarQ[g_Graph] := (
			   PlanarGivenCycleQ[ g, Rest[FindCycle[g]] ]
			  )

PlanarGivenCycleQ[g_Graph, cycle_List] :=
	Module[{b, j, i},
	   {b, j} = FindBridge[g, cycle];
	   If[ InterlockQ[j, cycle],
		False,
		Apply[And, Table[SingleBridgeQ[b[[i]],j[[i]]], {i,Length[b]}]]
	   ]
	]

SingleBridgeQ[b_Graph, {_}] := PlanarQ[b]

SingleBridgeQ[b_Graph, j_List] :=
	PlanarGivenCycleQ[ JoinCycle[b,j],
		Join[ ShortestPath[b,j[[1]],j[[2]]], Drop[j,2]] ]

JoinCycle[g1_Graph, cycle_List] :=
	Module[{g=g1},
		Scan[(g = AddEdge[g,#])&, Partition[cycle,2,1] ];
		AddEdge[g,{First[cycle],Last[cycle]}]
	]

FindBridge[g_Graph, cycle_List] :=
    Module[{rg = RemoveCycleEdges[g, cycle], b, bridge, j},
	b = Map[
		(IsolateSubgraph[rg,g,cycle,#])&,
		Select[ConnectedComponents[rg], (Intersection[#,cycle]=={})&]
	];
	b = Select[b, (!EmptyQ[#])&];
	j = Join[
		Map[Function[bridge,Select[cycle, MemberQ[Edges[bridge][[#]],1]&] ], b],
		Complement[
			Select[ToOrderedPairs[g],
				(Length[Intersection[#,cycle]] == 2)&],
			Partition[Append[cycle,First[cycle]],2,1]
		]
	];
	{b, j}
    ]

RemoveCycleEdges[g_Graph, c_List] :=
	FromOrderedPairs[
		Select[ ToOrderedPairs[g], (Intersection[c,#] === {})&],
		Vertices[g]
	]

IsolateSubgraph[g_Graph,orig_Graph,cycle_List,cc_List] :=
	Module[{eg=ToOrderedPairs[g], og=ToOrderedPairs[orig]},
		FromOrderedPairs[
			Join[
				Select[eg, (Length[Intersection[cc,#]] == 2)&],
				Select[og, (Intersection[#,cycle]!={} &&
					Intersection[#,cc]!={})&]
			],
			Vertices[g]
		]
	]

InterlockQ[ bl_List, c_List ] :=
	Module[{in = out = {}, code, jp, bridgelist = bl },
		While [ bridgelist != {},
			{jp, bridgelist} = {First[bridgelist],Rest[bridgelist]};
			code = Sort[ Map[(Position[c, #][[1,1]])&, jp] ];
			If[ Apply[ Or, Map[(LockQ[#,code])&, in] ],
				If [ Apply[Or, Map[(LockQ[#,code])&, out] ],
					Return[True],
					AppendTo[out,code]
				],
				AppendTo[in,code]
			]
		];
		False
	]

LockQ[a_List,b_List] := Lock1Q[a,b] || Lock1Q[b,a]

Lock1Q[a_List,b_List] :=
	Module[{bk, aj},
		bk = Min[ Select[Drop[b,-1], (#>First[a])&] ];
		aj = Min[ Select[a, (# > bk)&] ];
		(aj < Max[b])
	]

End[]



Protect[
AcyclicQ,
AddEdge,
AddVertex,
AllPairsShortestPath,
ArticulationVertices,
Automorphisms,
Backtrack,
BiconnectedComponents,
BiconnectedQ,
BinarySearch,
BinarySubsets,
BipartiteMatching,
BipartiteQ,
BreadthFirstTraversal,
Bridges,
CartesianProduct,
(* CatalanNumber, *)
ChangeEdges,
ChangeVertices,
ChromaticNumber,
ChromaticPolynomial,
CirculantGraph,
CircularVertices,
CliqueQ,
CodeToLabeledTree,
Cofactor,
CompleteGraph,
CompleteQ,
Compositions,
ConnectedComponents,
ConnectedQ,
ConstructTableau,
Contract,
CostOfPath,
Cycle,
DeBruijnSequence,
DegreeSequence,
DeleteCycle,
DeleteEdge,
DeleteFromTableau,
DeleteVertex,
DepthFirstTraversal,
DerangementQ,
Derangements,
Diameter,
Dijkstra,
DilateVertices,
DistinctPermutations,
Distribution,
DurfeeSquare,
Eccentricity,
EdgeChromaticNumber,
EdgeColoring,
EdgeConnectivity,
Edges,
Element,
EmptyGraph,
EmptyQ,
EncroachingListSet,
EquivalenceClasses,
EquivalenceRelationQ,
Equivalences,
EulerianCycle,
EulerianQ,
Eulerian,
ExactRandomGraph,
ExpandGraph,
ExtractCycles,
FerrersDiagram,
FindCycle,
FindSet,
FirstLexicographicTableau,
FromAdjacencyLists,
(* FromCycles, *)
FromInversionVector,
FromOrderedPairs,
FromUnorderedPairs,
FunctionalGraph,
Girth,
GraphCenter,
GraphComplement,
GraphDifference,
GraphIntersection,
GraphJoin,
GraphPower,
GraphProduct,
GraphSum,
GraphUnion,
GraphicQ,
GrayCode,
GridGraph,
HamiltonianCycle,
HamiltonianQ,
Harary,
HasseDiagram,
HeapSort,
Heapify,
HideCycles,
Hypercube,
IdenticalQ,
IncidenceMatrix,
IndependentSetQ,
Index,
InduceSubgraph,
InitializeUnionFind,
InsertIntoTableau,
IntervalGraph,
InversePermutation,
Inversions,
InvolutionQ,
IsomorphicQ,
IsomorphismQ,
Isomorphism,
Josephus,
KSubsets,
LabeledTreeToCode,
LastLexicographicTableau,
LexicographicPermutations,
LexicographicSubsets,
LineGraph,
LongestIncreasingSubsequence,
M,
MakeGraph,
MakeSimple,
MakeUndirected,
MaximalMatching,
MaximumAntichain,
MaximumClique,
MaximumIndependentSet,
MaximumSpanningTree,
MinimumChainPartition,
MinimumChangePermutations,
MinimumSpanningTree,
MinimumVertexCover,
MultiplicationTable,
NetworkFlowEdges,
NetworkFlow,
NextComposition,
NextKSubset,
NextPartition,
NextPermutation,
NextSubset,
NextTableau,
NormalizeVertices,
NthPair,
NthPermutation,
NthSubset,
NumberOfCompositions,
NumberOfDerangements,
NumberOfInvolutions,
NumberOfPartitions,
NumberOfPermutationsByCycles,
NumberOfSpanningTrees,
NumberOfTableaux,
OrientGraph,
PartialOrderQ,
PartitionQ,
Partitions,
PathConditionGraph,
Path,
PerfectQ,
PermutationGroupQ,
(* PermutationQ, *)
Permute,
PlanarQ,
PointsAndLines,
Polya,
PseudographQ,
RadialEmbedding,
Radius,
RandomComposition,
RandomGraph,
RandomHeap,
RandomKSubset,
RandomPartition,
RandomPermutation1,
RandomPermutation2,
(* RandomPermutation, *)
RandomSubset,
RandomTableau,
RandomTree,
RandomVertices,
RankGraph,
RankPermutation,
RankSubset,
RankedEmbedding,
ReadGraph,
RealizeDegreeSequence,
RegularGraph,
RegularQ,
RemoveSelfLoops,
RevealCycles,
RootedEmbedding,
RotateVertices,
Runs,
SamenessRelation,
SelectionSort,
SelfComplementaryQ,
ShakeGraph,
ShortestPathSpanningTree,
ShortestPath,
ShowGraph,
ShowLabeledGraph,
SignaturePermutation,
SimpleQ,
Spectrum,
SpringEmbedding,
StableMarriage,
Star,
StirlingFirst,
StirlingSecond,
Strings,
StronglyConnectedComponents,
Subsets,
TableauClasses,
TableauQ,
TableauxToPermutation,
Tableaux,
ToAdjacencyLists,
(* ToCycles, *)
ToInversionVector,
ToOrderedPairs,
ToUnorderedPairs,
TopologicalSort,
TransitiveClosure,
TransitiveQ,
TransitiveReduction,
TranslateVertices,
TransposePartition,
TransposeTableau,
TravelingSalesmanBounds,
TravelingSalesman,
TreeQ,
TriangleInequalityQ,
Turan,
TwoColoring,
UndirectedQ,
UnionSet,
UnweightedQ,
V,
VertexColoring,
VertexConnectivity,
VertexCoverQ,
Vertices,
WeaklyConnectedComponents,
Wheel,
WriteGraph,
DilworthGraph ]

EndPackage[ ]




