/*
This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as
published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version, see
<http://www.gnu.org/licenses/>.

This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty
of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
*/

/****************************************************************************************************************************

Very quick guide for Linux users
(For more details including compilation on Windows see below.)

1) Install csdp or sdpa. For example on Ubuntu you can install csdp by typing: 

 % sudo apt-get install coinor-csdp

2) Compile this file:  

 % g++ ExactDensityBounder.cpp -lm -fopenmp -O3 -o edb

3) Tell the program which SDP solver you are using: e.g. if you are using csdp then:

 % ./edb -s csdp

4) Examples of usage:
 
 % ./edb 5 {123 124 345}

 % ./edb 6 {123 124 345 156} 

 % ./edb 6 {123 124 134 234} {123 124 125 345 346} {123 124 345 156 256} {123 124 125 346 356 456} -o F_3

 % ./edb 6 {123 124 134 234} 4:{123} -i -o K4G1

 % ./edb 6 {123 124 134 234} 4:{123} -i -t + 1.1 4.2 4.3 -o ShortProof

The 1st example reproves the upper bound of 2/9 of Frankl and Furedi's result (Theorem 2.2) using 5 vertex 3-graphs.  
The proof is written in "Temp.txt". (Note that supersaturation is used so the forbidden graph K_4^- = {123 124 134} 
is automatically added.) 

The 2nd example reproves the upper bound of 2/9 in Theorem 2.3 using 6 vertex 3-graphs. The proof is written in "Temp.txt".
 
The 3rd example reproves the upper bound of 8/27 for the family F_3 in Theorem 3.4 using 6 vertex 3-graphs. The proof is
written in "F_3.txt".

The 4th example reproves the bound of 5/9 for the induced {K_4^(3),G_1} problem (Theorem 4.1 due to Razborov) using 6 vertex
3-graphs. The proof is written in "K4G1.txt".

The 5th example proves the same result as the 4th but only makes use of the single type of order 1 together with two of the
three types of order 4. This gives a shorter proof than the previous example, which is written in "ShortProof.txt".


Provided that the SDP solver is successful the program will attempt to prove that (a suitable rational approximation of the
SDP solver's output) is exact.
 
In many cases it will fail (often because the result is not sharp) however in these cases the absolute value of the SDP
output still provides a numerical upper bound for the Turan density that may be of interest.


See the comments below for more detailed usage instructions.

****************************************************************************************************************************/


//Name             : ExactDensityBounder.
//
//Author           : Rahil Baber. 
//
//Date             : 6th July 2011.
//
//Description      : An easier to use version of DensityBounder which can make approximate results exact. The program
//                   determines an upper bound for Turan densities of 3-uniform hypergraphs. It creates an SDP problem to
//                   determine the bound. The SDP problem is then solved by third party SDP solvers (such as csdp or sdpa).
//                   The program then removes rounding errors and outputs the details into a ".txt" file that can be verified
//                   by the program DensityChecker. This source code has not been written so that it is easy to understand
//                   (instead see the source code for DensityBounder which is a bit more readable).
//
//Acknowledgements : This program is largely based on the method described by Alexander Razborov in his paper
//                   "On 3-Hypergraphs with Forbidden 4-Vertex Configurations", SIAM Journal on Discrete Mathematics, 24 (3)
//                   946-963 (2010). Thanks to Emil Vaughan for pointing out that formulating the problem in the primal
//                   setting is more efficient for smaller problems.
//
//Compilation instructions : This program has been successfully compiled under Visual C++ 2010, and gcc 4.4.1. The code makes
//                   use of OpenMP to speed up the calculations when multiple processors are present, though the program will
//                   still function correctly if OpenMP is not enabled/supported by the compiler. 
//
//                   Visual C++ 2010 : The program should be compiled as a console application. It is recommended that the
//                   "Active Solution Configuration" be set to "Release" so that the compiler can fully optimize the code.
//                   OpenMP can be enabled via "Project>Properties>Configuration Properties>C/C++>Language>Open MP Support".
//
//                   gcc 4.4.1 : Simply navigate to the directory in which the source code resides and type at the terminal
//                   "g++ ExactDensityBounder.cpp -lm -O3 -o edb"
//                   this should create an optimized executable called "./edb". To enable OpenMP support instead type
//                   "g++ ExactDensityBounder.cpp -lm -fopenmp -O3 -o edb".
//
//Basic usage instructions : The program formulates the problem of finding an upper bound of the Turan density of a family
//                   of 3-graphs as an SDP (semidefinite program). The program cannot solve SDP problems itself but instead
//                   calls an external third party SDP solver such as CSDP or SDPA. As such a stand-alone binary executable
//                   of an SDP solver should be available somewhere on your system. Currently this program only produces the
//                   SDP problem in the ".dat-s" format and can only handle solutions written in a format which is consistent
//                   with that of CSDP's output or one of the SDPA family of solvers.
//
//                   The program should be run from a command-line interface as it expects arguments to be passed to
//                   it from the command-line. The very first time the program is run it should be told what the default
//                   command is for solving an SDP problem. This can be set using the "-s" switch. For example in windows we
//                   might type "ExactDensityBounder.exe -s C:\sdpa.exe" at the command-line to indicate the executable
//                   "sdpa.exe" located on the C drive should be run when the program wants to solve an SDP problem. In linux
//                   we might type "./edb -s ./csdp" to indicate that the executable "csdp" located in the current working
//                   directory should be run as the default solver.
//
//                   Once the default SDP command has been set, an upper bound for the turan density can be calculated by
//                   typing first the order of the subgraphs H which we will be searching over, followed by the edge sets of
//                   the graphs that are being forbidden. E.g. "./edb 6 {123 124 134 234}". Optionally a series of switches
//                   can be passed after the edge lists to indicate a variety of extra data. E.g. the "-o" switch is used to
//                   specify an output filename, the "-i" switch indicates that the graphs are forbidden induced subgraphs,
//                   and the "-c" switch overrides the default SDP command. See the sections below for more detailed
//                   information.
//
//                   The program analyses the solution to the SDP problem given by the SDP solvers and tries to guess what
//                   the bound should be if there were no rounding errors. It then tries to rigorously prove this bound. If
//                   it fails it should return the best upper bound that it could prove, however, this functionality has not
//                   yet been fully implemented (see the "To do" list below).
//
//Examples of basic usage :
//    ExactDensityBounder.exe -s csdp.exe
//    ./edb 6 {123 124 134 234} -o K4
//    ExactDensityBounder 5 {234} 4 : {123 124 134 234} -i -o "Induced Prob"
//    ./edb 5 {123 124 134} -c ./sdpa
//
//Usage :
//<Executable>
//    Displays (basic) usage.
//
//<Executable> Test <Filename>
//    Wipes <Filename> and writes "<Modify>" into it. Used to internally test the system() command. This is not intended for
//    users but included for completeness.
//
//<Executable> -s <SDP Command>
//    Use to set the default SDP command for solving SDP problems. See the "Arguments" section below.
//
//<Executable> <Order> <Graphs> [-i] [-c <SDP Command>] [-o <Filename>] [-t <+|-> [<type order>.<type index> ...]] [-<p|d>]
//[-r] [-a <number of decimal places>] [-b]
//    Use to calculate an upper bound for the Turan density of <Graphs>. See the "Arguments" section below. The switches can
//    be specified in any order after <Graphs>.
//
//Arguments :
//    <> : A mandatory field.
//
//    [] : An optional field.
//
//    <Executable> : The name of the compiled program. For example "ExactDensityBounder.exe" on windows, "./edb" or "./a.out"
//    in linux.
//
//    <Order> : The order of the subgraphs H. This should be an integer between 3 and 9. Increasing the order will result in
//    a better bound but the computation will be more intensive. Currently orders higher than 9 are not supported by this
//    program.
//
//    <Graphs> : A list of 3-uniform hypergraphs (separated by spaces). Graphs with orders higher than 9 are not supported.
//    Each graph should either be of the form "<Order of graph> : <Edge set>" or simply "<Edge set>". "<Order of graph>"
//    should be a single digit and if it is omitted then it is taken to be the largest label of a vertex in the edge set.
//    The vertices are assumed to be labelled by integers starting from 1. The edge set should be enclosed in curly braces
//    and contain a list of 3 digit numbers (separated by spaces) representing the edges. See the "Examples" section above.
//
//    [-i] : Indicates that we wish to consider only induced subgraphs. If "-i" is used then the program only forbids
//    <Graphs> from appearing as *induced* subgraphs. If "-i" is omitted then <Graphs> are forbidden from appearing as
//    subgraphs (as in the normal definition of Turan density).
//
//    [-c] : Used to override the default command (that runs the SDP solver) with <SDP command>. See also the
//    "<SDP command>" and "-s" sections below.
//
//    <SDP command> : The command used to run an SDP solver on the underlying semidefinite programming problem. For example
//    "./csdp" or "sdpa.exe". The program runs the SDP solver by running the command
//    "<SDP command> <Filename>.dat-s <Filename>.out". It always assumes the SDP solver completed successfully and that the
//    solution is stored in "<Filename>.out". It can currently only read the ".out" file if it adheres to the same format as
//    csdp or sdpa's output files. See also the sections on "[-c]", "-s", and "<Filename>".
//
//    [-o] : Used to indicate <Filename> as the name to be used for the files the program creates. If <Filename> is not set
//    using "-o" then the default name "Temp" is used. See also the section on "<Filename>".
//
//    <Filename> : A filename (without extension). It is used to name the temporary ".dat-s" input file to the SDP solver
//    and the ".out" output file. It is also used to name the ".txt" file which will hold the data which forms the proof.
//    The program DensityChecker can be used to verify the calculation from the ".txt" file.
//
//    -s : Used to save <SDP command> as the default command that runs the SDP solver. The command is saved in the file
//    "ConfigEDB.cfg". See the "[-c]" and "<SDP command>" sections.
//
//    [-t] : Sets which types to include or exclude. Run with "-t -" to see the "<type order>.<type index>" of each of the
//    types. "-t -" is used to exclude specific types, "-t +" is used to include specific types. For example "-t + 4.1 0.1"
//    means include only the types "4.1" and "0.1". Note that <type index> starts from 1 not 0. If this switch is not used
//    the default behaviour is to include all types (which is equivalent to the "-t -" command). 
//
//    [-p] : Forces the program to work in the primal SDP setting (maximizing Tr(CX) s.t. X>=0 and Tr(A_i X)=a_i). If both
//    -p and -d are omitted then the application will choose the form which produces the smallest problem.
//
//    [-d] : Forces the program to work in the dual SDP setting (minimizing sum(a_i y_i) s.t. sum(A_i y_i)-C>=0). If both -p
//    and -d are omitted then the application will choose the form which produces the smallest problem.
//
//    [-r] : Causes the default rounding scheme to be ignored, and instead asks the user to choose the target density, the
//    sharp inequalities, and the 0-eigenvectors.
//
//    [-a] : Overrides the accuracyProjectSolution variable (which describes how many decimal places the SDP solution should
//    be rounded too) with <number of decimal places>. <number of decimal places> should be a positive integer less than or
//    equal to 1000 (values less than 10 are recommended).
//
//    [-b] : If this switch is omitted (and "-i" is omitted) then the program will automatically add graphs (which have
//    blow-ups that are supergraphs of a graph in <Graphs>) to the forbidden graph list. This does not affect the Turan
//    density by a supersaturation argument, but can improve the bound the program is able to find. If "-b" (or "-i") is
//    specified then no additional graphs are added to the forbidden graph list.

//*** Unfinished. ***
//To do:
//- A parameter file to allow the user to tweak all the tolerances without recompilation.
//- Need to add the functionality that if making the Turan density exact fails it should output its best bound.
//  The work around of using "-r" with target density "1/1" should sometimes work.
//- Currently does not take the extremal constructions into account (which may help in making a result exact).
//- Need to neaten the code and add more comments.
//It is unlikely to get finished.



#ifdef _OPENMP
#include <omp.h>
#endif

#include <math.h>
#include <vector>
#include <list>
#include <limits>
#include <iomanip>
#include <iostream>
#include <fstream>
#include <string>
#include <sstream>
#include <cstring>
#include <stdlib.h>
using namespace std;

/***Start of the list of tolerance and threshold values which you may want to tweak.***/

//Used to allocate memory for CLargeInt objects. If an overflow error occurs it most likely
//means that LI_MAX was set too low. The larger LI_MAX is, the slower the program will run
//and the more memory it will require.
//*** You may need to change this value.***
#define LI_MAX 1000

//Used by ProjectSolution().
//The higher the value the closer the solution will be to the SDP's output (but the longer the program will take to run).
//Roughly speaking it will round the output to "accuracyProjectSolution" number of decimal places.
//Should be set to a positive integer. Can be changed via the "-a" command-line switch.
long accuracyProjectSolution = 6;

//Used by CalculateSharpIneq() to decide which inequalities are sufficiently close to the target density that they
//should be made sharp.
//Should be set to a small non-negative value.
double toleranceSharpIneq = 1e-6;

//Used in MakeMatricesPositiveDefinite() to decide how many 0-eigenvectors there are in each matrix.
//If the "Threshold for zero" value in the "_constraints.txt" file is above thresholdMaxConstraint then it is excluded as
//a choice for the number of 0-eigenvectors.
//Should be set to a small non-negative value.
double thresholdMaxConstraint = 1e-5;

//Used in MakeMatricesPositiveDefinite() to decide how many 0-eigenvectors there are in each matrix.
//If the "max error in rationalized values" value in the "_constraints.txt" file is above thresholdMaxError then it is
//excluded as a choice for the number of 0-eigenvectors.
//Should be set to a small non-negative value.
double thresholdMaxError = 1e-5;

//Used in MakeMatricesPositiveDefinite() to decide how many 0-eigenvectors there are in each matrix.
//If the "largest denominator in program's guess" value in the "_constraints.txt" file is above thresholdSimpleRational
//then it is excluded as a choice for the number of 0-eigenvectors.
//Should be set to a small positive value.
long thresholdSimpleRational = 20;

//Used in MakeMatricesPositiveDefinite() to rationalize the 0-eigenvectors.
//The denominator of the rationalized eigenvector entry is limited to "accuracySimpleRational".
//Should be set to a positive value.
long accuracySimpleRational = 100;

//Used by ProjectSolution() to identify zero entries when solving a linear system of constraints.
//Should be a small non-negative value.
double toleranceProjectSolution = 1e-6;

//Used by FindZeroEigenData() to identify entries which should be considered zero.
//Should be a small non-negative value.
double toleranceEigenConstraint = 1e-10;

//Used by CalculateSharpIneq() to decide what the target density we should be aiming for should be.
//The denominator is limited to be less than or equal to TargetDensityAccuracy.
//Should be set to a positive integer. It is unlikely that you will need to change this value.
long accuracyTargetDensity = 100;

/***End of the list of tolerances and thresholds.***/

//Set to false to manually choose the target density, sharp inequalities, and the 0-eigenvectors.
//Can be set to false via the "-r" command-line switch.
bool bDefaultRoundingScheme = true;

const char* configName      = "ConfigEDB.cfg"; //Holds the default SDP command.
const char* defaultFilename = "Temp"; 

string FilenameSDP;
string FilenameOut;
string FilenameTxt;
string SDPCommand;
string RunCommand;

//Factorial[n] = n!.
long Factorial[10] = {1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880};

//Iso[][] (together with Factorial[]) will be useful in generating isomorphic copies of
//graphs. Iso[][] is an array (orderOfH! by orderOfH) of all permutations of
//{0,1,...,orderOfH-1}. Iso[p][n] is the image of n under permutation p. Iso[][] is ordered
//in such a way that {Iso[p][] : p<=r!} is the set of permutations in which only {0,1,..,r-1}
//are permuted and {r,...,orderOfH-1} are left in place. InitIso() allocates memory for
//Iso[][] and initializes it. Its memory is deallocated by Cleanup().
vector<vector<long> > Iso;

//Exception classes which are thrown when overflow or division by zero occurs.
class ErrorOverflow{};
class ErrorDivisionByZero{};

//Used to store and manipulate very large integers. Can store any integer of magnitude less
//than 2^(8*LI_MAX). The integer's magnitude is stored in a static array to avoid slow
//dynamic memory allocation. The amount of the array being used is also stored, to increase
//the speed of operations on small integers. Functions which use this class may throw an
//ErrorOverflow or ErrorDivisionByZero object.
class CLargeInt
{
public:
    //The default constructor, sets the large integer to 0.
    CLargeInt();

    //The following constructors allow the fundamental integral data types to be
    //automatically cast into CLargeInt objects. These constructors call Assign() and so may
    //throw an ErrorOverflow object. No object will be thrown if LI_MAX >=
    //sizeof(unsigned long), (sizeof(unsigned long) = 4 usually).
    CLargeInt(unsigned int  input);
    CLargeInt(int           input);
    CLargeInt(unsigned long input);
    CLargeInt(long          input);

    //The copy constructor.
    CLargeInt(const CLargeInt& input);
    
    //The assignment operator.
    CLargeInt& operator=(const CLargeInt& input);

    //The unary operators. Used to evaluate +A or -A, when A is a CLargeInt object.
    const CLargeInt operator+() const;
    const CLargeInt operator-() const;

    //The basic arithmetic operators. The +,-,* operators may throw an ErrorOverflow object.
    //The integer division operators / and % are wrappers for the Divide() function and so
    //may throw an ErrorDivisionByZero object.
    friend const CLargeInt operator+(const CLargeInt& A, const CLargeInt& B);
    friend const CLargeInt operator-(const CLargeInt& A, const CLargeInt& B);
    friend const CLargeInt operator*(const CLargeInt& A, const CLargeInt& B);
    friend const CLargeInt operator/(const CLargeInt& A, const CLargeInt& B);
    friend const CLargeInt operator%(const CLargeInt& A, const CLargeInt& B);

    //The comparison operators. They are wrappers for the Compare() function.
    friend const bool operator==(const CLargeInt& A, const CLargeInt& B);
    friend const bool operator!=(const CLargeInt& A, const CLargeInt& B);
    friend const bool operator< (const CLargeInt& A, const CLargeInt& B);
    friend const bool operator<=(const CLargeInt& A, const CLargeInt& B);
    friend const bool operator> (const CLargeInt& A, const CLargeInt& B);
    friend const bool operator>=(const CLargeInt& A, const CLargeInt& B);

    //Used to convert the CLargeInt object into a long. If the value stored in the object is
    //larger than what can be stored in a long variable it throws an ErrorOverflow object.
    long ConvertToLong() const;

private:
    //Used by constructors to convert the fundamental integral data types to CLargeInt
    //objects. It may throw an ErrorOverflow object but only if LI_MAX<sizeof(unsigned long),
    //(sizeof(unsigned long) = 4 usually).
    void Assign(unsigned long input);

    //Used by the comparison operators to compare two large integers A, B.
    //Returns -1 if A < B,
    //         0 if A = B,
    //         1 if A > B.
    static long Compare(const CLargeInt& A, const CLargeInt& B);

    //Used by the integer division operators / and %. Define |input| to be the magnitude of
    //input, and sign(input) to be +1 if input>=0 and -1 if input<0. The function calculates
    //outputQ, and outputR as follows:
    //outputQ = sign(inputN)*sign(inputD)*(|inputN|/|inputD|), where '/' is integer division,
    //outputR = sign(inputN)*(|inputN|%|inputD|).
    //The choice of signs for outputQ and outputR were chosen so that inputN = inputD *
    //outputQ + outputR. If inputD = 0 an ErrorDivisionByZero class is thrown. The function
    //assumes outputQ resides in a different address to inputN, inputD, and outputR. It also
    //assumes outputR resides in a different address to inputN, inputD, and outputQ.
    static void Divide(const CLargeInt& inputN, const CLargeInt& inputD,
                       CLargeInt& outputQ, CLargeInt& outputR);

    //Data[], Size, and bPositive hold the data that represents the large integer. Data[] and
    //Size hold the magnitude of the large integer, bPositive holds the sign. bPositive =
    //true if the integer is greater than or equal to 0, and false otherwise. The magnitude
    //of the integer is given by:
    //Data[Size-1]*256^(Size-1)+Data[Size-2]*256^(Size-2)+...+Data[1]*256+Data[0].
    //Size holds the number of bytes the integer takes up in Data[]. To avoid leading zeros,
    //which can artificially make a small integer use up a lot of the Data[] array, we add
    //the following restriction. If the integer is 0 then Size=0, otherwise Data[Size-1] must
    //be non-zero.
    unsigned char Data[LI_MAX];
    long          Size; //Number of bytes. Size=0 <=> integer=0. Size!=0 => Data[Size-1]!=0.
    bool          bPositive; //bPositive=true if integer>=0, false otherwise.
};

//Allows us to use << operator to output the CLargeInt object L to the screen or a file.
ostream& operator<<(ostream& out, const CLargeInt& L);

//Used to store and manipulate fractions with very large numerators and denominators. It uses
//two CLargeInt member variables, named N and D, which store the numerator and denominator of
//the fraction respectively. Functions which use this class may throw an ErrorOverflow or
//ErrorDivisionByZero object.
class CLargeFraction
{
public:
    //The default constructor, sets the fraction to 0.
    CLargeFraction();
        
    //The following constructors allow the fundamental integral data types to be
    //automatically cast into CLargeFraction objects. They may throw an ErrorOverflow object.
    //No object will be thrown if LI_MAX >= sizeof(unsigned long), (sizeof(unsigned long) = 4
    //usually).
    CLargeFraction(unsigned int  input);
    CLargeFraction(int           input);
    CLargeFraction(unsigned long input);
    CLargeFraction(long          input);

    //A constructor that allows CLargeInt objects to be automatically cast into
    //CLargeFraction objects.
    CLargeFraction(const CLargeInt& input);
    
    //The unary operators. Used to evaluate +A or -A, when A is a CLargeFraction object.
    //These functions may throw an ErrorDivisionByZero object.
    const CLargeFraction operator+() const;
    const CLargeFraction operator-() const;

    //The basic arithmetic operators. The +,-,*,/ operators may throw an ErrorOverflow or
    //ErrorDivisionByZero object.
    friend const CLargeFraction operator+(const CLargeFraction& A, const CLargeFraction& B);
    friend const CLargeFraction operator-(const CLargeFraction& A, const CLargeFraction& B);
    friend const CLargeFraction operator*(const CLargeFraction& A, const CLargeFraction& B);
    friend const CLargeFraction operator/(const CLargeFraction& A, const CLargeFraction& B);

    //The comparison operators. They are wrappers for the Compare() function. They may throw
    //an ErrorOverflow or ErrorDivisionByZero object.
    friend const bool operator==(const CLargeFraction& A, const CLargeFraction& B);
    friend const bool operator!=(const CLargeFraction& A, const CLargeFraction& B);
    friend const bool operator< (const CLargeFraction& A, const CLargeFraction& B);
    friend const bool operator<=(const CLargeFraction& A, const CLargeFraction& B);
    friend const bool operator> (const CLargeFraction& A, const CLargeFraction& B);
    friend const bool operator>=(const CLargeFraction& A, const CLargeFraction& B);

    //Removes any common factors between the numerator and the denominator and makes the
    //denominator positive. It may throw an ErrorDivisionByZero object.
    void Simplify();

    CLargeInt N; //Numerator.
    CLargeInt D; //Denominator.

private:
    //Used by the comparison operators to compare two fractions A, B.
    //Returns -1 if A < B,
    //         0 if A = B,
    //         1 if A > B.
    //It may throw an ErrorOverflow or ErrorDivisionByZero object.
    static long Compare(const CLargeFraction& A, const CLargeFraction& B);
};

//Allows us to use << operator to output the CLargeFraction object F to the screen or a file.
ostream& operator<<(ostream& out, const CLargeFraction& F);

//The default constructor, sets the large integer to 0.
CLargeInt::CLargeInt()
{
    Size      = 0;
    bPositive = true;

    return;
}

//Allows an unsigned int to be automatically cast into a CLargeInt object. It may throw an
//ErrorOverflow object.
CLargeInt::CLargeInt(unsigned int input)
{
    Assign(input);
    return;
}

//Allows an int to be automatically cast into a CLargeInt object. It may throw an
//ErrorOverflow object.
CLargeInt::CLargeInt(int input)
{
    unsigned long temp;

    if(input>=0)
        Assign(input);
    else
    {
        //Call Assign() on -input.

        //-input may be 2147483648 which could cause an overflow. To avoid this situation we
        //cast -(input+1) to an unsigned long then add one.
        input++;
        temp = -input;
        temp++;

        Assign(temp);
        bPositive = false;
    }

    return;
}

//Allows an unsigned long to be automatically cast into a CLargeInt object. It may throw an
//ErrorOverflow object.
CLargeInt::CLargeInt(unsigned long input)
{
    Assign(input);
    return;
}

//Allows a long to be automatically cast into a CLargeInt object. It may throw an
//ErrorOverflow object.
CLargeInt::CLargeInt(long input)
{
    unsigned long temp;

    if(input>=0)
        Assign(input);
    else
    {
        //Call Assign() on -input.

        //-input may be 2147483648 which could cause an overflow. To avoid this situation we
        //cast -(input+1) to an unsigned long then add one.
        input++;
        temp = -input;
        temp++;

        Assign(temp);
        bPositive = false;
    }

    return;
}

//The copy constructor.
CLargeInt::CLargeInt(const CLargeInt& input)
{
    long i;
    
    //Check if we are trying to copy the object to itself.
    if(this==&input)
        return;

    Size      = input.Size;
    bPositive = input.bPositive;
    
    //We do not need to copy all of Data[].
    for(i=0;i<Size;i++)
        Data[i] = input.Data[i];

    return;
}

//The assignment operator.
CLargeInt& CLargeInt::operator=(const CLargeInt& input)
{
    long i;

    //Check if we are trying to assign the object to itself.
    if(this==&input)
        return *this;

    Size      = input.Size;
    bPositive = input.bPositive;
    
    //We do not need to copy all of Data[].
    for(i=0;i<Size;i++)
        Data[i] = input.Data[i];

    return *this;
}

//The + unary operator. Used to evaluate +A, when A is a CLargeInt object.
const CLargeInt CLargeInt::operator+() const
{
    return *this;
}

//The - unary operator. Used to evaluate -A, when A is a CLargeInt object.
const CLargeInt CLargeInt::operator-() const
{
    CLargeInt output = *this;
    
    if(output.Size!=0)
        output.bPositive = !output.bPositive;
    
    return output;
}

//The + arithmetic operator. Returns A+B. It may throw an ErrorOverflow object.
const CLargeInt operator+(const CLargeInt& A, const CLargeInt& B)
{
    long      i;
    long      sum, carryOver;
    CLargeInt output;

    //Deal with A=0 or B=0 as a special case.
    if(A.Size==0)
        return B;

    if(B.Size==0)
        return A;

    //A!=0 and B!=0.
    if(A.bPositive==B.bPositive)
    {
        //Add the magnitudes of A and B.
        output.bPositive = A.bPositive;

        //output will take up max{A.Size,B.Size}+{1 or 0} bytes. Provisionally set the size.
        if(A.Size>=B.Size)
            output.Size = A.Size;
        else
            output.Size = B.Size;

        //Do the addition byte by byte.
        carryOver = 0;
        for(i=0;i<output.Size;i++)
        {
            sum = carryOver;

            if(i<A.Size)
                sum += A.Data[i];

            if(i<B.Size)
                sum += B.Data[i];
            
            output.Data[i] = sum%256;
            carryOver      = sum/256;
        }

        if(carryOver!=0)
        {
            if(output.Size>=LI_MAX)
                throw ErrorOverflow(); //Overflow occured.

            output.Data[output.Size] = carryOver;
            output.Size++;
        }

        return output;
    }
    else
    {
        //Find the difference between the magnitudes of A and B.

        //If A.Size==B.Size assume |A|>=|B|.
        if(A.Size>=B.Size)
        {
            output.bPositive = A.bPositive;
            output.Size      = A.Size;

            //Do the subtraction of B from A, byte by byte.
            carryOver = 0;
            for(i=0;i<output.Size;i++)
            {
                sum = carryOver;

                if(i<A.Size)
                    sum += A.Data[i];

                if(i<B.Size)
                    sum -= B.Data[i];
                            
                if(sum>=0)
                {
                    output.Data[i] = sum;
                    carryOver      = 0;
                }
                else
                {
                    output.Data[i] = 256+sum;
                    carryOver      = -1;
                }
            }

            if(carryOver==0)
            {
                //Remove leading zeros.
                while(true)
                {
                    if(output.Size==0)
                    {
                        output.bPositive = true;
                        return output;
                    }

                    if(output.Data[output.Size-1]!=0)
                        return output;

                    output.Size--;
                }
            }

            //A.Size==B.Size but |A|<|B|. 
        }

        //A.Size<B.Size or (A.Size==B.Size and carryOver!=0) => |A|<|B|
        output.bPositive = B.bPositive;
        output.Size      = B.Size;

        //Do the subtraction of A from B, byte by byte.
        carryOver = 0;
        for(i=0;i<output.Size;i++)
        {
            sum = carryOver;

            if(i<A.Size)
                sum -= A.Data[i];

            if(i<B.Size)
                sum += B.Data[i];
                        
            if(sum>=0)
            {
                output.Data[i] = sum;
                carryOver      = 0;
            }
            else
            {
                output.Data[i] = 256+sum;
                carryOver      = -1;
            }
        }

        //Remove leading zeros.
        while(true)
        {
            if(output.Size==0)
            {
                output.bPositive = true;
                return output;
            }

            if(output.Data[output.Size-1]!=0)
                return output;

            output.Size--;
        }
    }
}

//The - arithmetic operator. Returns A-B. It may throw an ErrorOverflow object.
const CLargeInt operator-(const CLargeInt& A, const CLargeInt& B)
{
    return A+(-B);
}

//The * arithmetic operator. Returns A*B. It may throw an ErrorOverflow object.
const CLargeInt operator*(const CLargeInt& A, const CLargeInt& B)
{
    long      i, j, k;
    long      sum, carryOver;
    CLargeInt output;

    //Deal with A=0 or B=0 as a special case.
    if(A.Size==0)
        return A;
    
    if(B.Size==0)
        return B;
    
    //A!=0 and B!=0.

    //Calculate the sign of A*B.
    if(A.bPositive==B.bPositive)
        output.bPositive = true;
    else
        output.bPositive = false;

    //output will take up A.Size+B.Size-{1 or 0} bytes. Provisionally set the size.
    if(A.Size+B.Size-1>LI_MAX)
        throw ErrorOverflow(); //Overflow will occur.

    output.Size = A.Size+B.Size-1;

    //Initialize output to 0.
    for(i=0;i<output.Size;i++)
        output.Data[i] = 0;
        
    //Multiply byte by byte.
    for(i=0;i<A.Size;i++)
    for(j=0;j<B.Size;j++)
    {
        carryOver = A.Data[i]*B.Data[j];
        
        //Add carryOver*(256^(i+j)) to output.

        if(carryOver==0)
            continue;

        //Do the addition byte by byte.
        for(k=i+j;k<output.Size;k++)
        {
            sum = output.Data[k]+carryOver;
            
            output.Data[k] = sum%256;
            carryOver      = sum/256;

            if(carryOver==0)
                break;
        }

        if(carryOver!=0)
        {
            if(output.Size>=LI_MAX)
                throw ErrorOverflow();

            output.Data[output.Size] = carryOver;
            output.Size++;
        }
    }

    return output;
}

//The / arithmetic operator. Returns A/B where / is integer division (see Divide() for a more
//precise definition of /). If B=0 an ErrorDivisionByZero object is thrown.
const CLargeInt operator/(const CLargeInt& A, const CLargeInt& B)
{
    CLargeInt output, Temp;

    CLargeInt::Divide(A,B,output,Temp);
    
    return output;
}

//The % arithmetic operator. Returns A%B (see Divide() for a more precise definition of %).
//If B=0 an ErrorDivisionByZero object is thrown.
const CLargeInt operator%(const CLargeInt& A, const CLargeInt& B)
{
    CLargeInt output, Temp;

    CLargeInt::Divide(A,B,Temp,output);
    
    return output;
}

//Returns (A==B).
const bool operator==(const CLargeInt& A, const CLargeInt& B)
{
    return (CLargeInt::Compare(A,B)==0);
}

//Returns (A!=B).
const bool operator!=(const CLargeInt& A, const CLargeInt& B)
{
    return (CLargeInt::Compare(A,B)!=0);
}

//Returns (A<B).
const bool operator<(const CLargeInt& A, const CLargeInt& B)
{
    return (CLargeInt::Compare(A,B)==-1);
}

//Returns (A<=B).
const bool operator<=(const CLargeInt& A, const CLargeInt& B)
{
    return (CLargeInt::Compare(A,B)!=1);
}

//Returns (A>B).
const bool operator>(const CLargeInt& A, const CLargeInt& B)
{
    return (CLargeInt::Compare(A,B)==1);
}

//Returns (A>=B).
const bool operator>=(const CLargeInt& A, const CLargeInt& B)
{
    return (CLargeInt::Compare(A,B)!=-1);
}

//Used to convert the CLargeInt object into a long. If the value stored in the object is
//larger than what can be stored in a long variable it throws an ErrorOverflow object.
long CLargeInt::ConvertToLong() const
{
    long i;
    long max, min;
    long result;
    
    //Check if overflow occurs.
    if(Size>=sizeof(long))
    {
        //max =  2^{8*sizeof(long)-1} - 1,
        //min = -2^{8*sizeof(long)-1}.
        //To avoid overflow, we define them as follows:
        max  = (1L<<(8*sizeof(long)-2))-1;
        max += 1L<<(8*sizeof(long)-2);
        min  = -max;
        min -= 1;

        if(*this<min || max<*this)
            throw ErrorOverflow();
    }

    //Convert the object into a long.
    result = 0;
    for(i=0;i<sizeof(long) && i<Size;i++)
    {
        if(bPositive==true)
            result += Data[i]<<(8*i);
        else
            result -= Data[i]<<(8*i);
    }
    
    return result;
}

//Used by constructors to convert the fundamental integral data types to CLargeInt objects.
//It may throw an ErrorOverflow object but only if LI_MAX < sizeof(unsigned long),
//(sizeof(unsigned long) = 4 usually).
void CLargeInt::Assign(unsigned long input)
{
    long i;

    Size      = 0;
    bPositive = true;

    if(input==0)
        return;

    //Initialize Data[] so it equals input.
    for(i=0;i<sizeof(unsigned long);i++)
    {
        //input!=0, fill in Data[i]. 
        
        if(i>=LI_MAX)
            throw ErrorOverflow();
        
        //Remove the least significant byte of input and store it in Data[i].
        Data[i] = input%256;
        input   = input/256;

        if(input==0)
        {
            Size = i+1;
            return;
        }
    }

    return;
}

//Used by the comparison operators to compare two large integers A, B.
//Returns -1 if A < B,
//         0 if A = B,
//         1 if A > B.
long CLargeInt::Compare(const CLargeInt& A, const CLargeInt& B)
{
    long i;
    
    if(A.bPositive!=B.bPositive)
    {
        if(A.bPositive==true)
            return 1;  //A >= 0 > B.
        else
            return -1; //A < 0 <= B.
    }

    //A and B have the same sign.

    if(A.Size!=B.Size)
    {
        if((A.bPositive==true  && A.Size>B.Size)
        || (A.bPositive==false && A.Size<B.Size))
            return 1;  //A > B.
        else
            return -1; //A < B.
    }

    //A.Size = B.Size and A.bPositive = B.bPositive.
    
    //Check the magnitude of the integers, beginning with the most significant byte. 
    for(i=A.Size-1;i>=0;i--)
    {
        if(A.Data[i]>B.Data[i])
        {
            if(A.bPositive==true)
                return 1;  //A > B.
            else
                return -1; //A < B.
        }

        if(A.Data[i]<B.Data[i])
        {
            if(A.bPositive==true)
                return -1; //A < B.
            else
                return 1;  //A > B.
        }

        //A.Data[i] = B.Data[i]. Check the next most significant byte.
    }

    return 0; //A = B.
}

//Used by the integer division operators / and %. Define |input| to be the magnitude of
//input, and sign(input) to be +1 if input>=0 and -1 if input<0. The function calculates
//outputQ, and outputR as follows:
//outputQ = sign(inputN)*sign(inputD)*(|inputN|/|inputD|), where '/' is integer division,
//outputR = sign(inputN)*(|inputN|%|inputD|).
//The choice of signs for outputQ and outputR were chosen so that inputN = inputD * 
//outputQ + outputR. If inputD = 0 an ErrorDivisionByZero class is thrown. The function
//assumes outputQ resides in a different address to inputN, inputD, and outputR. It also
//assumes outputR resides in a different address to inputN, inputD, and outputQ.
void CLargeInt::Divide(const CLargeInt& inputN, const CLargeInt& inputD,
                       CLargeInt& outputQ, CLargeInt& outputR)
{
    long i, j;
    long sum, carryOver;
    bool bSubtractD;
    long bitD, bitN;
    long nBit, nByte;
    long rBit, rByte;
    long overflowByte;

    //Deal with inputN=0 and inputD=0 as a special case.
    if(inputD.Size==0)
        throw ErrorDivisionByZero();
    
    if(inputN.Size==0)
    {
        outputQ = 0;
        outputR = 0;
        return;
    }

    //We will do the division via the long division method in binary.

    //Find the leading bit of inputD.
    j = inputD.Data[inputD.Size-1];
    for(i=7;i>=0;i--)
    {
        if((j>>i)%2==1)
            break;
    }

    bitD = i;
        
    //Find the leading bit of inputN.
    j = inputN.Data[inputN.Size-1];
    for(i=7;i>=0;i--)
    {
        if((j>>i)%2==1)
            break;
    }

    bitN = i;

    //Check if |inputD|>|inputN| obviously holds.
    if(inputD.Size>inputN.Size || (inputD.Size==inputN.Size && bitD>bitN))
    {
        outputQ = 0;
        outputR = inputN;
        return;
    }
    
    //|inputN| is stored in at least as many bits as |inputD|.

    //|inputD| is stored in "(inputD.Size-1)*8+bitD+1" bits. Fill outputR with the
    //"(inputD.Size-1)*8+bitD+1" most significant bits of inputN.
    outputR.Size = inputD.Size;

    for(i=0;i<outputR.Size;i++)
        outputR.Data[i] = 0;

    nBit  = bitN;
    nByte = inputN.Size-1;
    rBit  = bitD;
    rByte = inputD.Size-1;
    while(true)
    {
        //Copy bit.
        outputR.Data[rByte] |= ((inputN.Data[nByte]>>nBit)&1)<<rBit;

        //Move to next bit.
        rBit--;

        if(rBit==-1)
        {
            rBit = 7;
            rByte--;
        }

        if(rByte==-1)
            break;

        nBit--;

        if(nBit==-1)
        {
            nBit = 7;
            nByte--;
        }
    }

    //outputR has been initialized.

    //nBit and nByte hold the position of the last bit in inputN to be added to outputR. This
    //is also the most significant bit possible in outputQ.

    //Initialize outputQ.
    outputQ.Size = nByte+1;

    for(i=0;i<outputQ.Size;i++)
        outputQ.Data[i] = 0;
    
    //Calculate the value of the bits of outputQ.
    overflowByte = 0; //Used when outputR overflows.
    while(true)
    {
        //Calculate the value of the "nByte*8+nBit" bit in outputQ.

        //Check if |outputR| >= |inputD|. If it is we set the bit in outputQ to 1, and
        //subtract |inputD| from |outputR|.
        if(overflowByte!=0)
            bSubtractD = true; //|outputR| takes up more bytes than |inputD|.
        else
        {
            for(i=outputR.Size-1;i>=0;i--)
            {
                if(outputR.Data[i]>inputD.Data[i])
                {
                    bSubtractD = true;
                    break;
                }

                if(outputR.Data[i]<inputD.Data[i])
                {
                    bSubtractD = false;
                    break;
                }
            }
            
            if(i==-1)
                bSubtractD = true; //|outputR| = |inputD|.
        }
        
        if(bSubtractD==true)
        {
            //Subtract inputD from outputR.
            carryOver = 0;
            for(i=0;i<outputR.Size;i++)
            {
                sum = outputR.Data[i]-inputD.Data[i]+carryOver;
                
                if(sum>=0)
                {
                    outputR.Data[i] = sum;
                    carryOver       = 0;
                }
                else
                {
                    outputR.Data[i] = 256+sum;
                    carryOver       = -1;
                }
            }

            //Set the bit of outputQ to 1.
            outputQ.Data[nByte] |= 1<<nBit; 
        }

        if(nByte==0 && nBit==0)
            break; //Found all the bits of outputQ.

        //Move to the next bit.
        nBit--;

        if(nBit==-1)
        {
            nBit = 7;
            nByte--;
        }

        //Shift the bits of outputR to the left (by multiplying by 2) and append the
        //"nByte*8+nbit" bit from inputN to the end.
        carryOver = (inputN.Data[nByte]>>nBit)&1;
        for(i=0;i<outputR.Size;i++)
        {
            sum             = 2*outputR.Data[i]+carryOver;
            outputR.Data[i] = sum%256;
            carryOver       = sum/256;
        }
        
        overflowByte = carryOver;

        //Calculate the next bit of outputQ.
    }

    //outputR holds the magnitude of the remainder.

    //Calculate the signs of outputQ and outputR.
    if(inputN.bPositive==inputD.bPositive)
        outputQ.bPositive = true;
    else
        outputQ.bPositive = false;

    outputR.bPositive = inputN.bPositive;
    
    //Remove leading zeros.
    while(true)
    {
        if(outputQ.Size==0)
        {
            outputQ.bPositive = true;
            break;
        }

        if(outputQ.Data[outputQ.Size-1]!=0)
            break;

        outputQ.Size--;
    }

    while(true)
    {
        if(outputR.Size==0)
        {
            outputR.bPositive = true;
            break;
        }

        if(outputR.Data[outputR.Size-1]!=0)
            break;

        outputR.Size--;
    }

    return;
}

//Allows us to use << operator to output the CLargeInt object L to the screen or a file.
ostream& operator<<(ostream& out, const CLargeInt& L)
{
    long          i, n;
    unsigned char digit[3*LI_MAX];
    CLargeInt     Q, R;
    stringstream  sout;

    if(L==0)
    {
        sout << "0";
        out << sout.str();
        return out;
    }

    //Set Q to be the magnitude of L.
    if(L<0)
    {
        sout << "-"; //Display the sign of the integer.
        Q = -L;
    }
    else
        Q = L;

    //Remove the least significant decimal digit of Q and store it in digit[].
    n = 0;
    while(Q!=0)
    {
        R = Q%10;
        Q = Q/10;
        
        digit[n] = R.ConvertToLong();
        n++;
    }
    
    //Display the digits of the integer.
    for(i=n-1;i>=0;i--)
        sout << char('0'+digit[i]);

    out << sout.str();

    return out;
}

//The default constructor, sets the fraction to 0.
CLargeFraction::CLargeFraction()
{
    N = 0;
    D = 1;
    return;
}

//Allows an unsigned int to be automatically cast into a CLargeFraction object. It may throw
//an ErrorOverflow object.
CLargeFraction::CLargeFraction(unsigned int input)
{
    N = input;
    D = 1;
    return;
}

//Allows an int to be automatically cast into a CLargeFraction object. It may throw an
//ErrorOverflow object.
CLargeFraction::CLargeFraction(int input)
{
    N = input;
    D = 1;
    return;
}

//Allows an unsigned long to be automatically cast into a CLargeFraction object. It may throw
//an ErrorOverflow object.
CLargeFraction::CLargeFraction(unsigned long input)
{
    N = input;
    D = 1;
    return;
}

//Allows a long to be automatically cast into a CLargeFraction object. It may throw an
//ErrorOverflow object.
CLargeFraction::CLargeFraction(long input)
{
    N = input;
    D = 1;
    return;
}

//Allows a CLargeInt object to be automatically cast into a CLargeFraction object.
CLargeFraction::CLargeFraction(const CLargeInt& input)
{
    N = input;
    D = 1;
    return;
}

//The + unary operator. Used to evaluate +A, when A is a CLargeFraction object. It may throw
//an ErrorDivisionByZero object.
const CLargeFraction CLargeFraction::operator+() const
{
    CLargeFraction output = *this;
    output.Simplify();
    return output;
}

//The - unary operator. Used to evaluate -A, when A is a CLargeFraction object. It may throw
//an ErrorDivisionByZero object.
const CLargeFraction CLargeFraction::operator-() const
{
    CLargeFraction output;

    output.N = -N;
    output.D = D;

    output.Simplify();
    return output;
}

//The + arithmetic operator. Returns A+B. It may throw an ErrorOverflow or
//ErrorDivisionByZero object.
const CLargeFraction operator+(const CLargeFraction& A, const CLargeFraction& B)
{
    CLargeFraction output;

    output.N = A.N*B.D+B.N*A.D;
    output.D = A.D*B.D;

    output.Simplify();
    return output;
}

//The - arithmetic operator. Returns A-B. It may throw an ErrorOverflow or
//ErrorDivisionByZero object.
const CLargeFraction operator-(const CLargeFraction& A, const CLargeFraction& B)
{
    return A+(-B);
}

//The * arithmetic operator. Returns A*B. It may throw an ErrorOverflow or
//ErrorDivisionByZero object.
const CLargeFraction operator*(const CLargeFraction& A, const CLargeFraction& B)
{
    CLargeFraction output;

    output.N = A.N*B.N;
    output.D = A.D*B.D;
    
    output.Simplify();
    return output;
}

//The / arithmetic operator. Returns A/B. It may throw an ErrorOverflow or
//ErrorDivisionByZero object.
const CLargeFraction operator/(const CLargeFraction& A, const CLargeFraction& B)
{
    CLargeFraction output;

    if(B.D==0)
        throw ErrorDivisionByZero();

    output.N = A.N*B.D;
    output.D = B.N*A.D;
    
    output.Simplify();
    return output;
}

//Returns (A==B). It may throw an ErrorOverflow or ErrorDivisionByZero object.
const bool operator==(const CLargeFraction& A, const CLargeFraction& B)
{
    return (CLargeFraction::Compare(A,B)==0);
}

//Returns (A!=B). It may throw an ErrorOverflow or ErrorDivisionByZero object.
const bool operator!=(const CLargeFraction& A, const CLargeFraction& B)
{
    return (CLargeFraction::Compare(A,B)!=0);
}

//Returns (A<B). It may throw an ErrorOverflow or ErrorDivisionByZero object.
const bool operator<(const CLargeFraction& A, const CLargeFraction& B)
{
    return (CLargeFraction::Compare(A,B)==-1);
}

//Returns (A<=B). It may throw an ErrorOverflow or ErrorDivisionByZero object.
const bool operator<=(const CLargeFraction& A, const CLargeFraction& B)
{
    return (CLargeFraction::Compare(A,B)!=1);
}

//Returns (A>B). It may throw an ErrorOverflow or ErrorDivisionByZero object.
const bool operator>(const CLargeFraction& A, const CLargeFraction& B)
{
    return (CLargeFraction::Compare(A,B)==1);
}

//Returns (A>=B). It may throw an ErrorOverflow or ErrorDivisionByZero object.
const bool operator>=(const CLargeFraction& A, const CLargeFraction& B)
{
    return (CLargeFraction::Compare(A,B)!=-1);
}

//Removes any common factors between the numerator and the denominator and makes the
//denominator positive. It may throw an ErrorDivisionByZero object.
void CLargeFraction::Simplify()
{
    CLargeInt R1, R2, R3;

    if(D==0)
        throw ErrorDivisionByZero();

    if(N==0)
    {
        D = 1;
        return;
    }
    
    //Apply the Euclidean algorithm to determine the highest common factor.
    R1 = N;
    R2 = D;
    while(true)
    {
        R3 = R1%R2;

        if(R3==0)
            break; //R2 is the highest common factor.

        R1 = R2;
        R2 = R3;
    }
    
    //Divide by the highest common factor.
    N = N/R2;
    D = D/R2;

    //Make the denominator positive.
    if(D<0)
    {
        N = -N;
        D = -D;
    }

    return;
}

//Used by the comparison operators to compare two fractions A, B.
//Returns -1 if A < B,
//         0 if A = B,
//         1 if A > B.
//It may throw an ErrorOverflow or ErrorDivisionByZero object.
long CLargeFraction::Compare(const CLargeFraction& A, const CLargeFraction& B)
{
    long      result;
    CLargeInt L;

    if(A.D==0 || B.D==0)
        throw ErrorDivisionByZero();

    L = A.N*B.D-A.D*B.N;

    if(L==0)
        return 0; //A = B.

    if(L>0)
        result = 1; //A > B, if both A.D>0 and B.D>0.
    else
        result = -1;//A < B, if both A.D>0 and B.D>0.
    
    //Take the signs of A.D and B.D into account.
    if(A.D<0)
        result = -result;

    if(B.D<0)
        result = -result;

    return result;
}

//Allows us to use << operator to output the CLargeFraction object F to the screen or a file.
ostream& operator<<(ostream& out, const CLargeFraction& F)
{
    stringstream sout;
    sout << F.N << "/" << F.D;
    out << sout.str();
    return out;
}

//A simple class to store 3-graphs that have at most nine vertices. We take our vertex set to be
//{0,1,...,nV-1}, where nV is the number of vertices in our graph. We can represent a flag as
//(g, n) a CGraph g and an integer n, by taking the labelled vertices to be 0,1,...,n-1 and to
//avoid confusion we label them 0,1,...,n-1 respectively (instead of 1,2,...,n).
class CGraph
{
public:
    //Number of vertices.
    long nV;

    //The adjacency matrix.
    //M[i][j][k] = true  if edge ijk is in the 3-graph,
    //             false if edge ijk is not in the 3-graph.
    bool M[9][9][9];
};

//A structure that stores information about flags generated from a type.
class CFlagInfo
{
public:
    //The graph of the type that the labelled vertices of the flags induce. Formally the type
    //is (Type, Type.nV).
    CGraph Type;

    //FlagList defines a list of flags (FlagList[], Type.nV) for which the labelled
    //vertices induce the graph Type, i.e. FlagList[n].M[i][j][k] should be identical to
    //Type.M[i][j][k] for all i,j,k < Type.nV. 
    vector<CGraph> FlagList;

    //An element of the basis is a linear combination of flags (with coefficients 1 or -1).
    //How many flags the element is composed of is stored in Basis[r][i].size() (where 
    //r=0 if the element is invariant, r=1 if it isn't, and 0 <= i < Basis[r].size() indexes
    //elements in the basis). The indices of where the flags, from the element, appear in
    //FlagList is stored in Basis[r][i][]. The coefficients of the flags can
    //be determined based on whether we are looking at an invariant element or not. If the
    //element is invariant we take all the coefficients to be 1. If it is not invariant we
    //take the coefficient of the first flag to be 1 and the coefficient of the second to be
    //-1.
    vector<vector<long> > Basis[2];
};

//CHashNode and CHashTable are used to define a very basic hash table.
//CHashTable::HashFunction() needs to be defined for each <Key,Value> type.
//The Key class should provide an operator== function.

//Used to create a linked list of key, value pairs in CHashTable.
template<typename Key, typename Val>
class CHashNode
{
public:
    CHashNode<Key,Val>* pNext;
    long hashedKey;
    Key  key;
    Val  val;
};

//A hash table where collisions are resolved via a linked list for each "bucket".
template<typename Key, typename Val>
class CHashTable
{
public:
    CHashTable():noOfEntries(0),sizeOfTable(0),Table(NULL),mask(0){}
    ~CHashTable(){Free();}

    //Free the memory allocated to Table[] and set sizeOfTable to 0. 
    void Free();

    //Adds an element to the table.
    void Add(const Key& key, const Val& val);

    //Gets a pointer to the value associated with the first matching key it can find.
    //Returns NULL if no matching key can found.
    Val* Get(const Key& key) const;

    //Used to map keys into buckets.
    static long HashFunction(const Key& key);

    //Displays the lengths of the lists. Useful for debugging purposes.
    //Long list lengths indicate a poor hash function or a small table size.
    void CollisionStats();

    long noOfEntries;
    long sizeOfTable;
    CHashNode<Key,Val>** Table; //An array of linked lists.

private:
    //Allocates the size and memory for the table.
    void Allocate(long sizeOfTable);

    //Used to resize Table if it is getting too full.
    void DoubleSize();

    //Used to do fast modding.
    long mask;
};

//Free the memory allocated to Table[] and set sizeOfTable to 0. 
template<typename Key, typename Val>
void CHashTable<Key,Val>::Free()
{
    long i;
    CHashNode<Key,Val>* pNode;

    if(Table!=NULL)
    {
        for(i=0;i<sizeOfTable;i++)
        {
            //Free the linked list.
            while(Table[i]!=NULL)
            {
                pNode = Table[i]->pNext;
                delete Table[i];
                Table[i] = pNode;
            }
        }

        delete[] Table;
        Table = NULL;
    }

    sizeOfTable = 0;
    noOfEntries = 0;
    mask        = 0;

    return;
}

//Allocates the size and memory for the table.
template<typename Key, typename Val>
void CHashTable<Key,Val>::Allocate(long sizeOfTable)
{
    long i;

    //Clear the table.
    Free();

    if(sizeOfTable<=0)
        return;
    
    //Allocate the table.
    Table             = new CHashNode<Key,Val>*[sizeOfTable];
    this->sizeOfTable = sizeOfTable;
    
    //Initialize the table.
    for(i=0;i<sizeOfTable;i++)
        Table[i] = NULL;
    
    return;
}

template<typename Key, typename Val>
void CHashTable<Key,Val>::DoubleSize()
{
    long i;
    long bucketIndex;
    CHashTable<Key,Val> ht;
    CHashNode<Key,Val>* pNode;
    long  newSize;

    if(sizeOfTable==0)
        newSize = 2;
    else
        newSize = sizeOfTable<<1;

    ht.Allocate(newSize);
    ht.mask = (mask<<1) | 1;

    //Move the entries accross.
    for(i=0;i<sizeOfTable;i++)
        while(Table[i]!=NULL)
        {
            pNode    = Table[i];
            Table[i] = Table[i]->pNext;
                
            bucketIndex           = (pNode->hashedKey)&ht.mask;
            pNode->pNext          = ht.Table[bucketIndex];
            ht.Table[bucketIndex] = pNode;
            ht.noOfEntries++;
        }

    //Clear the table.
    Free();

    //Swap the tables.
    Table          = ht.Table;
    sizeOfTable    = ht.sizeOfTable;
    noOfEntries    = ht.noOfEntries;
    mask           = ht.mask;
    ht.Table       = NULL;
    ht.sizeOfTable = 0;
    ht.noOfEntries = 0;
    ht.mask        = 0;
        
    return;
}

//Adds an element to the table.
template<typename Key, typename Val>
void CHashTable<Key,Val>::Add(const Key& key, const Val& val)
{
    long hashedKey;
    long bucketIndex;
    CHashNode<Key,Val>* pNew;

    //Check if the table is 3/4 full and if so create a larger Table.
    if(noOfEntries/3+1>sizeOfTable/4)
        DoubleSize();

    hashedKey   = HashFunction(key);
    bucketIndex = hashedKey & mask;
    
    //Create a new node.
    pNew = new CHashNode<Key,Val>;

    //Set the new node as the head of the list.
    pNew->pNext        = Table[bucketIndex];
    Table[bucketIndex] = pNew;
    noOfEntries++;

    //Fill in the rest of the node.
    pNew->hashedKey = hashedKey;
    pNew->key       = key;
    pNew->val       = val;

    return;
}

//Gets a pointer to the value associated with the first matching key it can find.
//Returns NULL if no matching key can be found.
template<typename Key, typename Val>
Val* CHashTable<Key,Val>::Get(const Key& key) const
{
    long hashedKey;
    CHashNode<Key,Val>* pNode;

    if(sizeOfTable==0)
        return NULL;

    hashedKey   = HashFunction(key);
    
    pNode = Table[hashedKey & mask];
    while(pNode!=NULL)
    {
        //Check if pNode hashedKey matches.
        if(pNode->hashedKey==hashedKey)
        {
            //Check the keys match.
            if(pNode->key==key)
                return &(pNode->val);
        }

        pNode = pNode->pNext;
    }

    return NULL;
}

//Displays the lengths of the lists. Useful for debugging purposes.
//Long list lengths indicate a poor hash function or a small table size.
template<typename Key, typename Val>
void CHashTable<Key,Val>::CollisionStats()
{
    long i, j;
    long lengths[6];
    CHashNode<Key,Val>* pNode; 

    cout << "Stats : " << endl;
    cout << noOfEntries << " / " << sizeOfTable << endl;

    for(i=0;i<6;i++)
        lengths[i] = 0;

    for(i=0;i<sizeOfTable;i++)
    {
        j = 0; //The length of the list.
        
        pNode = Table[i];
        while(pNode!=NULL)
        {
            j++;
            pNode = pNode->pNext;
        }

        if(j<=4)
            lengths[j]++;
        else
            lengths[5]++;
    }

    for(i=0;i<=4;i++)
        cout << i << "  : " << lengths[i] << endl;
    cout << "5+ : " << lengths[5] << endl;

    return;
}

//sw, and sz should not be zero.
//Returns an integer in [0,2^31-1]
long GetRandomInt(bool bSeed=false, unsigned long sw=123, unsigned long sz=345678)
{
    static unsigned long w = 123;
    static unsigned long z = 3456789;

    if(bSeed==true)
    {
        w = sw;
        z = sz;
    }

    z = 36969*(z&65535)+(z>>16);
    w = 18000*(w&65535)+(w>>16);
    
    return ((z<<16)+w)&0x7FFFFFFF;
}

class KeyFlag
{
public:
    long   typeOrder;
    CGraph flag;

    //Used when searching through the hash table.
    const bool operator==(const KeyFlag& in) const;

    //Used to hash the key.
    static long hashOrder[9];
    static long hashM[9][9][9];
};

long KeyFlag::hashOrder[9];
long KeyFlag::hashM[9][9][9];

bool IsCopy(const CGraph& g1, const CGraph& g2);

const bool KeyFlag::operator==(const KeyFlag& in) const
{
    if(typeOrder!=in.typeOrder)
        return false;

    if(IsCopy(flag,in.flag)==false)
        return false;

    return true;
}

class ValFlag
{
public:
    //Indices.
    long iType; //Index of the type in FInfo[typeOrder][].
    long iFlag; //Index of the flag in FlagList.
};

//Used to map keys into buckets.
template<>
long CHashTable<KeyFlag,ValFlag>::HashFunction(const KeyFlag& key)
{
    long i, j, k;
    long hashedKey;

    hashedKey = 0;
    
    hashedKey ^= key.hashOrder[key.typeOrder];
    hashedKey ^= key.hashOrder[key.flag.nV];

    for(i=0;i<key.flag.nV;i++)
    for(j=0;j<i;j++)
    for(k=0;k<j;k++)
        if(key.flag.M[i][j][k]==true)
            hashedKey ^= key.hashM[i][j][k];

    return hashedKey;
}

class KeyProd
{
public:
    unsigned char nType;
    long          iType;
    long         iFlag1;
    long         iFlag2;
    
    //Used when searching through the hash table.
    const bool operator==(const KeyProd& in) const;
};

//Used when searching through the hash table.
const bool KeyProd::operator==(const KeyProd& in) const
{
    if(in.nType==nType && in.iType==iType
    && in.iFlag1==iFlag1 && in.iFlag2==iFlag2)
        return true;

    return false;
}

//Used to map keys into buckets.
template<>
long CHashTable<KeyProd,long>::HashFunction(const KeyProd& key)
{
    long hashedKey = 0;
    
    hashedKey ^= 196723186*(key.nType);
    hashedKey ^= 26764*(key.iType&0xFFFF);
    hashedKey ^= 16223*(key.iType>>16);
    hashedKey ^=  9218*(key.iFlag1&0xFFFF);
    hashedKey ^=  8137*(key.iFlag1>>16);
    hashedKey ^= 24855*(key.iFlag2&0xFFFF);
    hashedKey ^= 29529*(key.iFlag2>>16);

    return hashedKey;
}

long              orderOfH;
vector<CGraph>    H;
vector<CGraph>    ForbiddenGraphs;
bool              bForbidInduced;
vector<CFlagInfo> FInfo[8]; //[Max order of a type+1]

bool         bDisplayTypeNumber;
bool         bIncludeTypesByDefault;
vector<long> modifyTypeList[8]; //[Max order of a type+1]

bool bForbidBlowups      = false;
long ForbiddenNonBlowups = 0;

//Indicates which setting the SDP should take.
// 1 = Primal.
// 0 = Undecided.
//-1 = Dual.
long ForceForm;

//A hash table to quickly find the position of flags in FInfo[][].FlagList[].
CHashTable<KeyFlag,ValFlag> htFlags;

//A hash table for each member of H to quickly find the product of flags.
vector<CHashTable<KeyProd,long> > htProd;

//Outputs the edge set of a CGraph object to the screen or a file.
void DisplayEdges(ostream& out, const CGraph& g)
{
    long i, j, k;
    bool bFirst;
    stringstream sout;

    sout << "{";

    //Used to help us put the commas in the right place.
    bFirst = true;

    for(i=0;i<g.nV;i++)
    for(j=0;j<i;j++)
    for(k=0;k<j;k++)
        if(g.M[k][j][i]==true)
        {
            //Display the commas that separate the edges.
            if(bFirst==false)
                sout << ", ";
            else
                bFirst = false;
                
            //Display the edge ijk.
            sout << k+1 << j+1 << i+1;
        }

    sout << "}";

    out << sout.str();

    return;
}

//Allows us to use << operator to output CGraph objects to the screen or a file.
ostream& operator<<(ostream& out, const CGraph& g)
{
    stringstream sout;

    //Display the number of vertices.
    sout << "Order " << g.nV << " : ";

    //Display the edges.
    sout << "Edges = ";
    DisplayEdges(sout,g);

    out << sout.str();

    return out;
}

//Allocates memory for Iso[][] and initializes it. The function assumes orderOfH has been set
//to a value between 3 and 9. Iso[][] is initialized to an orderOfH! by orderOfH array
//representing all permutations of {0,1,...,orderOfH-1}. Iso[p][n] will be the image of n
//under permutation p. Iso[][] will be ordered in such a way that {Iso[p][] : p<=r!} is the
//set of permutations in which only {0,1,..,r-1} are permuted and {r,...,orderOfH-1} are left
//in place.
void InitIso()
{
    long i, j, k, p, n;

    //Destroy Iso[][].
    Iso.clear();
    vector<vector<long> >().swap(Iso);

    //Allocate memory for Iso[][].
    Iso.resize(Factorial[orderOfH]);

    for(i=0;i<Factorial[orderOfH];i++)
        Iso[i].resize(orderOfH);

    //Initialize Iso[][].

    //Set Iso[0][] to be the identity map.
    for(n=0;n<orderOfH;n++)
        Iso[0][n] = n;

    //p = number of permutations filled in.
    p = 1; 

    //Inductively build up Iso[][]:
    //Assume {Iso[k][] : k<=i!} are all permutations which permute {0,1,...,i-1} but leave
    //{i,i+1,...,orderOfH-1}. Generate new permutations by swapping i and j in Iso[k][], for
    //all j<i and k<=i!.
    for(i=1;i<orderOfH;i++)
    for(j=0;j<i;j++)
    for(k=0;k<Factorial[i];k++)
    {
        for(n=0;n<orderOfH;n++)
        {
            if(Iso[k][n]==i)
                Iso[p][n] = j;
            else
            if(Iso[k][n]==j)
                Iso[p][n] = i;
            else
                Iso[p][n] = Iso[k][n];
        }

        p++;
    }
    
    return;
}

//Returns the number of edges in the graph g.
long NumberOfEdges(const CGraph& g)
{
    long i, j, k;
    long edges;

    edges = 0;

    for(i=0;i<g.nV;i++)
    for(j=0;j<i;j++)
    for(k=0;k<j;k++)
        if(g.M[i][j][k]==true)
            edges++;

    return edges;
}

//Number of digits needed to display n.
long NumberOfDigits(unsigned long n)
{
    long result = 0;

    if(n==0)
        return 1;

    while(n!=0)
    {
        result++;
        n /= 10;
    }

    return result;
}

//Returns true if g1 and g2 are of the same order and have the same edge set.
bool IsCopy(const CGraph& g1, const CGraph& g2)
{
    long i, j, k;

    if(g1.nV!=g2.nV)
        return false;

    for(i=0;i<g1.nV;i++)
    for(j=0;j<i;j++)
    for(k=0;k<j;k++)
        if(g1.M[i][j][k]!=g2.M[i][j][k])
            return false;

    return true;
}

//If bInduced = false, IsSubgraph() returns true if gSub is a subgraph of gSuper.
//If bInduced = true,  IsSubgraph() returns true if gSub is an induced subgraph of gSuper.
//Assumes Iso[][] has been initialized.
bool IsSubgraph(const CGraph& gSub, const CGraph& gSuper, bool bInduced)
{
    long i, j, k, p;
    bool bIsSub;

    if(gSuper.nV<gSub.nV)
        return false;

    if(bInduced==false)
    {
        //Permute the vertices of gSuper, to create all its isomorphic graphs, and check its edge
        //set with that of gSub. 
        for(p=0;p<Factorial[gSuper.nV];p++)
        {
            //Due to the way we constructed Iso[][],
            //{Iso[p][0], Iso[p][1], ..., Iso[p][gSuper.nV-1]} is a permutation of 
            //{0, 1, ..., gSuper.nV-1}.

            //Assume gSub is a subgraph. If at some point we've shown it isn't, then
            //try a different permutation.
            bIsSub = true;

            //Check the triple ijk (where i>j>k).
            for(i=0;i<gSub.nV;i++)
            {
                for(j=0;j<i;j++)
                {
                    for(k=0;k<j;k++)
                    {
                        //Check if the edge set of gSub is a subset.
                        if(gSub.M[i][j][k]==true
                        && gSuper.M[Iso[p][i]][Iso[p][j]][Iso[p][k]]==false)
                        {
                            //Try a different permutation p.
                            bIsSub = false;
                            break;
                        }
                    }

                    if(bIsSub==false)
                        break;
                }

                if(bIsSub==false)
                    break;
            }

            if(bIsSub==true)
                return true;
        }

        //Tested all permutations with no success.
        return false;
    }
    else
    {
        //The procedure for induced subgraphs is almost identical.
        for(p=0;p<Factorial[gSuper.nV];p++)
        {
            bIsSub = true;

            for(i=0;i<gSub.nV;i++)
            {
                for(j=0;j<i;j++)
                {
                    for(k=0;k<j;k++)
                    {
                        //Check if the edge sets are identical.
                        if(gSub.M[i][j][k]!=gSuper.M[Iso[p][i]][Iso[p][j]][Iso[p][k]])
                        {
                            bIsSub = false;
                            break;
                        }
                    }

                    if(bIsSub==false)
                        break;
                }

                if(bIsSub==false)
                    break;
            }

            if(bIsSub==true)
                return true;
        }

        return false;
    }

    return false;
}

//Each edge set of a graph can be assigned a unique integer. FindMinIso() creates gOut a
//graph with the property that the flag (gOut, orderOfType) is isomorphic to the flag 
//(gIn, orderOfType) and has the smallest integer associated with its edge set out of all
//such graphs. This function is used by CalculateFlagProducts() and GenerateBasis() to speed
//up searching through lists of flags. In order for these functions to work correctly the way
//FindMinIso() constructs the integer representing the edge set is important. The function
//assumes Iso[][] has been initialized.
void FindMinIso(CGraph gIn, CGraph& gOut, long orderOfType)
{
    long i, j, k, i2, j2, k2, p;
    bool bFoundSmaller, bBreak;

    //gOut holds the graph with the smallest integer associated with its edge set so far.
    gOut = gIn;

    //We wish to test all isomorphic flags of (gIn, orderOfType). We can construct them by
    //permutating vertices {orderOfType, orderOfType+1, ..., gIn.nV-1} of gIn, whilst keeping
    //vertices {0,1,...,orderOfType-1} fixed. Due to the way we constructed Iso[][], 
    //{gIn.nV-1-Iso[p][gIn.nV-1], gIn.nV-1-Iso[p][gIn.nV-2], ..., gIn.nV-1-Iso[p][0]}
    //describes just such a permutation (for 0 <= p < Factorial[gIn.nV-orderOfType]).
    for(p=0;p<Factorial[gIn.nV-orderOfType];p++)
    {
        bFoundSmaller = false;
        bBreak        = false;

        //We construct the unique integer associated with the edge set of a graph g by
        //considering an ordering of all triples ijk (with g.nV > i > j > k >=0).
        //Specifically ijk > rst if
        //   i>r,
        //or i==r and j>s,
        //or i==r and j==s and k>t.
        //Our integer is a "gIn.nV choose 3" bit number where each bit is associated with a
        //triple in our ordering. The largest triple is associated with the least significant
        //bit, and the second largest triple with the second least significant bit, etc. The
        //value of the bits in the integer are given by the associated triples, with a value
        //of 1 if ijk is an edge in g and 0 if it isn't.

        //We will not explicitly construct such integers, but by examining the edge set of
        //gOut and gIn, under the permutation p, we will make a comparison of which
        //associated integer would be smaller.

        //Cycle through the triples from the smallest to the largest. Note that we can take
        //i >= orderOfType because of the way we are permuting gIn.
        for(i=orderOfType;i<gIn.nV;i++)
        {
            i2 = gIn.nV-1-Iso[p][gIn.nV-1-i];

            for(j=0;j<i;j++)
            {
                j2 = gIn.nV-1-Iso[p][gIn.nV-1-j];
                
                for(k=0;k<j;k++)
                {
                    k2 = gIn.nV-1-Iso[p][gIn.nV-1-k];

                    //Compare bits

                    if(gIn.M[i2][j2][k2]==false && gOut.M[i][j][k]==true)
                    {
                        //gIn under permutation p is smaller than gOut. Store gIn under
                        //permutation p in gOut. 
                        bFoundSmaller = true;
                        bBreak        = true;
                        break;
                    }

                    if(gIn.M[i2][j2][k2]==true && gOut.M[i][j][k]==false)
                    {
                        //gIn under permutation p is larger than gOut. Try a different value
                        //of p.
                        bFoundSmaller = false;
                        bBreak        = true;
                        break;
                    }

                    //The integers are equal so far. Try the next triple ijk.
                }

                if(bBreak==true)
                    break;
            }

            if(bBreak==true)
                break;
        }

        if(bFoundSmaller==true)
        {
            //Copy the transformed gIn into gOut
            for(i=0;i<gIn.nV;i++)
            {
                i2 = gIn.nV-1-Iso[p][gIn.nV-1-i];

                for(j=0;j<gIn.nV;j++)
                {
                    j2 = gIn.nV-1-Iso[p][gIn.nV-1-j];
                    
                    for(k=0;k<gIn.nV;k++)
                    {
                        k2 = gIn.nV-1-Iso[p][gIn.nV-1-k];

                        gOut.M[i][j][k] = gIn.M[i2][j2][k2];
                    }
                }
            }
        }

        //Check if the next permutation produces a smaller integer.
    }

    return;
}

//This function returns true if the output of FindMinIso() would be the same as the input.
//This function is considerably faster than calling FindMinIso().
//This function is used by Generate().
//The function assumes Iso[][] has been initialized.
bool IsMinIso(const CGraph& gIn, long orderOfType)
{
    long i, j, k, i2, j2, k2, p;
    bool bBreak;
    
    //This function is almost identical to FindMinIso().
    //The difference is that this function returns as soon as anything smaller is found.
    for(p=0;p<Factorial[gIn.nV-orderOfType];p++)
    {
        bBreak = false;

        for(i=orderOfType;i<gIn.nV;i++)
        {
            i2 = gIn.nV-1-Iso[p][gIn.nV-1-i];

            for(j=0;j<i;j++)
            {
                j2 = gIn.nV-1-Iso[p][gIn.nV-1-j];
                
                for(k=0;k<j;k++)
                {
                    k2 = gIn.nV-1-Iso[p][gIn.nV-1-k];

                    //Compare bits

                    if(gIn.M[i2][j2][k2]==false && gIn.M[i][j][k]==true)
                        return false; //Found something smaller.

                    if(gIn.M[i2][j2][k2]==true && gIn.M[i][j][k]==false)
                    {
                        //gIn under permutation p is larger than gOut. Try a different value
                        //of p.
                        bBreak = true;
                        break;
                    }

                    //The integers are equal so far. Try the next triple ijk.
                }

                if(bBreak==true)
                    break;
            }

            if(bBreak==true)
                break;
        }
        
        //Check if the next permutation produces a smaller integer.
    }

    return true;
}

//Used in GenerateH() and GenerateFlags() to create a list of admissible graphs/flags of a
//given order not containing isomorphic copies. Takes as input the flag
//(inputFlag, orderOfType), which it uses to create admissible flags of order inputFlag.nV+1
//that have inputFlag as an induced subgraph. It stores the new admissible graphs it finds in
//List.
void Generate(CGraph inputFlag, long orderOfType, list<CGraph>& List)
{
    long    i, j, k;
    long    b, n;
    long    noOfEdgesToSet;
    bool    e;
    CGraph& g = inputFlag;

    //Create a new graph which is inputFlag with an extra vertex v.
    g.nV++;

    //We need to determine which edges containing vertex v our new graph will contain.
    noOfEdgesToSet = (g.nV-1)*(g.nV-2)/2;

    //There are 2^noOfEdgesToSet possible choices for the set of edges containing v. We let b
    //encode which of the edges we include in g via the value of its bits.
    //Bits    0,   1,   2,   3,   4,   5,   6,   7,   8, ... of b tells us whether edges 
    //edges 01v, 02v, 12v, 03v, 13v, 23v, 04v, 14v, 24v, ... respectively are in graph g.
   
    //Test all sets of edges containing v.
    for(b=0;b<(1<<noOfEdgesToSet);b++)
    {
        n = 0;      //n = position of bit in b. 
        i = g.nV-1; //i = vertex v.
        for(j=0;j<i;j++)
        for(k=0;k<j;k++)
        {
            //Check if ikj is in our graph, by inspecting bit n of b.
            if((b>>n)%2==1)
                e = true;
            else
                e = false;

            //Update the adjacency matrix.
            g.M[i][j][k] = e;
            g.M[i][k][j] = e;
            g.M[j][i][k] = e;
            g.M[j][k][i] = e;
            g.M[k][i][j] = e;
            g.M[k][j][i] = e;
            
            n++;
        }
                
        //Check if the graph is admissible.
        for(i=0;i<ForbiddenGraphs.size();i++)
            if(IsSubgraph(ForbiddenGraphs[i],g,bForbidInduced)==true)
                break;

        if(i!=ForbiddenGraphs.size())
            continue; //A forbidden graph was found.

        //Due to the way FindMinIso() chooses a minimal representive of the isomorphic copies
        //of a flag and the way we are inductively building up our list of flags from smaller
        //order flags which are minimally represented, we are guaranteed to generate all
        //flags which are minimally represented. Hence if the flag is not already the
        //minimal representation of all its isomorphic copies we can safely ignore it.
        if(IsMinIso(g,orderOfType)==false)
            continue;

        //Add to List.
        List.push_back(g);
    }
    
    return;
}

//Generates 'H' a linked list representing $\mathcal{H}$ the family of admissible graphs of
//order orderOfH. The function assumes orderOfH, ForbiddenGraphs, and bForbidInduced, have
//been initialized.
bool GenerateH()
{
    long        i, j, k, n;
    list<CGraph> OldList;
    list<CGraph> NewList;
    list<CGraph>::const_iterator it;
    CGraph      g;
    long        graphIndex, total, foundSoFar;
#ifdef _OPENMP
    long globalError;
    long loopError;
    vector<list<CGraph> > listArray;
    vector<list<CGraph>::const_iterator> listIndex;
#endif
    
    cout << "Generating H..." << endl;
    
    vector<CGraph>().swap(H);

    //We take 'OldList', a list of graphs of order n-1, and add a vertex to each graph in the
    //list to generate 'NewList' a list of graphs of order n. We then make 'NewList' the new
    //'OldList' and repeat the procedure until we get the desired order of graphs.

    cout << "\rChecking for order 0 forbidden graphs." << flush;

    //Check if order 0 graphs are forbidden.
    for(i=0;i<ForbiddenGraphs.size();i++)
        if(ForbiddenGraphs[i].nV==0)
        {
            cout << "\rFound 0 graphs.                       " << endl;
            cout << endl;
            cout << "ERROR : There are no graphs in H." << endl << endl;
            return false;
        }

    cout << "\rGenerating order 0 graphs :           " << flush;

    //Fill in the admissible graph of order 0.
    g.nV = 0;
    for(i=0;i<9;i++)
    for(j=0;j<9;j++)
    for(k=0;k<9;k++)
        g.M[i][j][k] = false;
    
    OldList.clear();
    OldList.push_back(g);
    
    NewList.clear();

    for(n=1;n<=orderOfH;n++)
    {
        //Generate NewList, a list of graphs of order n, from OldList a list of graphs of
        //order n-1.

        //graphIndex and total are used to keep track of the progress.
        graphIndex = 0;
        total      = OldList.size();
        foundSoFar = 0;

#ifdef _OPENMP
        //Create the listIndex of OldList so we can have random/direct access.
        listIndex.resize(total);
        for(i=0,it=OldList.begin();it!=OldList.end();i++,it++)
            listIndex[i] = it;

        //Create an empty list for each graph in OldList.
        listArray.resize(total);
        for(i=0;i<total;i++)
            listArray[i].clear();
                    
        cout << "\rGenerating order " << n << " graphs : ";
        cout << graphIndex << " of " << total << " graphs extended. ";
        cout << "Found " << foundSoFar << "." << flush;

        globalError = 0; //No error.
        #pragma omp parallel for schedule(dynamic) default(shared) private(i,loopError)
        for(i=0;i<total;i++)
        {
            try
            {
                #pragma omp critical (GenH_Progress)
                {
                    loopError = globalError; //No error.
                }

                if(loopError==0)
                    Generate(*(listIndex[i]),0,listArray[i]);
            }
            catch(bad_alloc& e)
            {
                loopError = 2; //bad_alloc error.
            }
            catch(...)
            {
                loopError = 1; //Unknown error
            }
        
            #pragma omp critical (GenH_Progress)
            {
                if(loopError>globalError)
                    globalError = loopError;

                if(globalError==0)
                {
                    foundSoFar += listArray[i].size();
                    graphIndex++;
                    cout << "\rGenerating order " << n << " graphs : ";
                    cout << graphIndex << " of " << total << " graphs extended. ";
                    cout << "Found " << foundSoFar << "." << flush;
                }
                //Move to the next graph in OldList.
            }
        }

        //Display any errors.
        if(globalError!=0)
        {
            cout << endl << endl;

            if(globalError==1)
                cout << "ERROR : An unknown error occurred.";

            if(globalError==2)
                cout << "ERROR : Allocating memory failed.";
            
            cout << endl << endl;
            return false;
        }

        //Append ListArray into NewList.
        for(i=0;i<total;i++)
            NewList.splice(NewList.begin(),listArray[i]);
#else
        for(it=OldList.begin();it!=OldList.end();it++)
        {
            foundSoFar = NewList.size();

            graphIndex++;
            cout << "\rGenerating order " << n << " graphs : ";
            cout << "Extending graph " << graphIndex << " of " << total << ". ";
            cout << "Found " << foundSoFar << "." << flush;


            Generate(*it,0,NewList);
            
            //Move to the next graph in OldList.
        }
#endif

        //Make NewList the new OldList.
        OldList.clear();
        NewList.swap(OldList);
        
        cout << "\rGenerating order " << n << " graphs :                       ";
        cout << "       ";
        j = 2*NumberOfDigits(total)+NumberOfDigits(foundSoFar);
        for(i=0;i<j;i++)
            cout << " ";
    }

    //Move OldList to H.
    H.clear();
    while(OldList.empty()==false)
    {
        H.push_back(OldList.front());
        OldList.pop_front();
    }
    
    //Display how large H is.
    n = H.size();

    cout << "\r                           ";
    if(n==1)
        cout << "\rFound 1 graph." << endl;
    else
        cout << "\rFound " << n << " graphs."<< endl;

    cout << endl;
    
    if(H.size()==0)
    {
        cout << "ERROR : There are no graphs in H." << endl << endl;
        return false;
    }

    return true;
}

//Generates admissible Flags of order FlagOrder. The type is given by FI.Type. 
void GenerateFlags(CFlagInfo& FI, long FlagOrder)
{
    long         i, j, k;
    list<CGraph> OldList;
    list<CGraph> NewList;
    list<CGraph>::const_iterator it;
    
    //Set the undefined parts of M[][][] to false.
    for(i=0;i<9;i++)
    for(j=0;j<9;j++)
    for(k=0;k<9;k++)
    {
        if(i==j || i==k || j==k    || i>=FI.Type.nV || j>=FI.Type.nV || k>=FI.Type.nV)
            FI.Type.M[i][j][k] = false;
    }

    FI.FlagList.clear();
    vector<CGraph>().swap(FI.FlagList);

    //Deal with the nonsensical case of creating a flag with FlagOrder < FI.Type.nV.
    if(FlagOrder<FI.Type.nV)
        return;

    //We take 'OldList', a list of graphs of order n-1, and add a vertex to each graph in the
    //list to generate 'NewList' a list of graphs of order n. We then make 'NewList' the new
    //'OldList' and repeat the procedure until we get the desired order of flags.
    
    //Initialize OldList with the only flag of order FI.Type.nV.
    OldList.clear();
    OldList.push_back(FI.Type);
    
    NewList.clear();

    for(i=FI.Type.nV+1;i<=FlagOrder;i++)
    {
        //For each flag in OldList generate flags for NewList by adding an extra vertex. 
        for(it=OldList.begin();it!=OldList.end();it++)
            Generate(*it,FI.Type.nV,NewList);

        //Swap the lists.
        OldList.clear();
        NewList.swap(OldList);
    }

    for(it=OldList.begin();it!=OldList.end();it++)
        FI.FlagList.push_back(*it);

    return;
}

//Generates a new basis of flags from the information in FI. It stores the necessary data in
//FI.Basis[].
void GenerateBasis(CFlagInfo& FI)
{
    long         i, j, k, m, n, p;
    vector<long> minFlag;
    CGraph       g;

    //Generate the orbit of each flag under the automorphism group of FI.Type, and store the
    //flag with the smallest index in minFlag[].
    minFlag.resize(FI.FlagList.size());

    //To start with we assume the unmodified flag has the smallest index out of the flags
    //that form its orbit.
    for(i=0;i<FI.FlagList.size();i++)
        minFlag[i] = i;

    //Find an element of the automorphism group, i.e. a permutation p that does not change
    //FI.Type.
    for(p=0;p<Factorial[FI.Type.nV];p++)
    {
        //Check p leaves FI.Type unchanged.
        g.nV = FI.Type.nV;
        
        for(i=0;i<FI.Type.nV;i++)
        for(j=0;j<FI.Type.nV;j++)
        for(k=0;k<FI.Type.nV;k++)
            g.M[i][j][k] = FI.Type.M[Iso[p][i]][Iso[p][j]][Iso[p][k]];

        if(IsCopy(g,FI.Type)==false)
            continue; //FI.Type under p is not FI.Type.

        //Apply p to all flags to create a member of their orbit.
        for(n=0;n<FI.FlagList.size();n++)
        {
            //Apply p to the flag.
            g.nV = FI.FlagList[n].nV;
                    
            for(i=0;i<g.nV;i++)
            for(j=0;j<g.nV;j++)
            for(k=0;k<g.nV;k++)
                g.M[i][j][k] = FI.FlagList[n].M[Iso[p][i]][Iso[p][j]][Iso[p][k]];

            //Find g in FI.FlagList[].
            FindMinIso(g,g,FI.Type.nV);

            for(m=0;m<FI.FlagList.size();m++)
                if(IsCopy(g,FI.FlagList[m])==true)
                    break;

            //Check if this new flag has a smaller index.
            if(m<minFlag[n])
                minFlag[n] = m;
        }
    }

    //Fill in the basis data.
    FI.Basis[0].clear();
    FI.Basis[1].clear();

    for(i=0;i<FI.FlagList.size();i++)
    {
        //Find the orbit in which the flag with the smallest index is i (if such an orbit
        //exists).
        
        if(minFlag[i]!=i)
            continue; //No such orbit.

        //Found a new invariant basis element i.e. the sum of the flags with minFlag[] = i.
        //The first flag in this element is FI.FlagList[i].
        FI.Basis[0].push_back(vector<long>(1,i));
        
        //Find the other flags in the orbit.
        for(j=i+1;j<FI.FlagList.size();j++)
        {
            if(minFlag[j]!=i)
                continue; //Not part of the orbit.

            //Add the flag FI.FlagList[j] to the invariant element.
            FI.Basis[0].back().push_back(j);
            
            //Add a new element which is not invariant.
            FI.Basis[1].push_back(vector<long>(2,0));
            FI.Basis[1].back()[0] = j;
            FI.Basis[1].back()[1] = i;
        }
    }

    return;
}

//Generate all the types of order less than or equal to orderOfH-2.
void FillFlagInfo()
{
    long      i, j, k, p, n;
    bool      bKeep, bFirst;
    long      m[4];
    CFlagInfo FI;

    //Set the type.
    FI.Type.nV = 0;

    for(i=0;i<9;i++)
    for(j=0;j<9;j++)
    for(k=0;k<9;k++)
        FI.Type.M[i][j][k] = false;

    cout << "\rFinding the types..." << flush;

    for(n=0;n<=orderOfH-2;n++)
    {
        GenerateFlags(FI,n);
        
        FInfo[n].resize(FI.FlagList.size());
        
        for(i=0;i<FI.FlagList.size();i++)
            FInfo[n][i].Type = FI.FlagList[i];
    }

    cout << "\rFinding the flags..." << flush;
    
    for(n=0;n<=orderOfH-2;n++)
    for(i=0;i<FInfo[n].size();i++)
        GenerateFlags(FInfo[n][i],(FInfo[n][i].Type.nV+orderOfH)/2);
    
    cout << "\rGenerating the basis elements..." << flush;

    for(n=0;n<=orderOfH-2;n++)
    for(i=0;i<FInfo[n].size();i++)
        GenerateBasis(FInfo[n][i]);

    //For display purposes work out the maximum number of digits in each field.
    {
        for(i=0;i<4;i++)
            m[i] = 0;

        for(n=0;n<=orderOfH-2;n++)
        {
            if(m[0]<FInfo[n].size())
                m[0] = FInfo[n].size();

            for(i=0;i<FInfo[n].size();i++)
            {
                if(m[1]<FInfo[n][i].FlagList.size())
                    m[1] = FInfo[n][i].FlagList.size();

                if(m[2]<FInfo[n][i].Basis[0].size())
                    m[2] = FInfo[n][i].Basis[0].size();

                if(m[3]<FInfo[n][i].Basis[1].size())
                    m[3] = FInfo[n][i].Basis[1].size();
            }
        }

        for(i=0;i<4;i++)
            m[i] = NumberOfDigits(m[i]);
    }

    cout << "\rNumber of flags (invariant, anti-invariant basis size) and Type :" << endl;
    for(n=0;n<=orderOfH-2;n++)
    for(i=0;i<FInfo[n].size();i++)
    {
        //Check if the type has been removed by the "-t" switch.
        if(bIncludeTypesByDefault==true)
            bKeep = true;
        else
            bKeep = false;

        for(j=0;j<modifyTypeList[n].size();j++)
            if(modifyTypeList[n][j]==i)
            {
                bKeep = !bKeep;
                break;
            }
        
        if(bDisplayTypeNumber==true)
        {
            if(bKeep==true)
                cout << "O ";
            else
                cout << "X ";

            cout << n << "." << i+1 << " ";

            for(j=0;j<m[0]-NumberOfDigits(i+1);j++)
                cout << " ";
        }

        if(bDisplayTypeNumber==true || bKeep==true)
        {
            for(j=0;j<m[1]-NumberOfDigits(FInfo[n][i].FlagList.size());j++)
                cout << " ";
            
            cout << FInfo[n][i].FlagList.size() << " (";
            
            for(j=0;j<m[2]-NumberOfDigits(FInfo[n][i].Basis[0].size());j++)
                cout << " ";
            
            cout << FInfo[n][i].Basis[0].size() << ", ";
            
            for(j=0;j<m[3]-NumberOfDigits(FInfo[n][i].Basis[1].size());j++)
                cout << " ";

            cout << FInfo[n][i].Basis[1].size() << ") ";
            cout << FInfo[n][i].Type.nV << ":{";

            bFirst = true;
            for(p=0;p<FInfo[n][i].Type.nV;p++)
            for(j=0;j<p;j++)
            for(k=0;k<j;k++)
                if(FInfo[n][i].Type.M[k][j][p]==true)
                {
                    if(bFirst==false)
                        cout << ", ";
                    else
                        bFirst = false;
                
                    cout << k+1 << j+1 << p+1;
                }

            cout << "}" << endl;
        }

        //Delete types.
        if(bKeep==false)
        {
            FInfo[n][i].Type.nV = 0;
            FInfo[n][i].FlagList.clear();
            FInfo[n][i].Basis[0].clear();
            FInfo[n][i].Basis[1].clear();
        }
    }

    //Shrink FInfo[n].
    for(n=0;n<=orderOfH-2;n++)
    {
        vector<CFlagInfo> Temp;

        //Count how many types there are.
        k = 0;
        for(i=0;i<FInfo[n].size();i++)
            if(FInfo[n][i].FlagList.size()>0)
                k++;

        Temp.resize(k);

        k = 0;
        for(i=0;i<FInfo[n].size();i++)
        {
            if(FInfo[n][i].FlagList.size()<=0)
                continue;

            //Store FInfo[n][i]
            Temp[k].Type = FInfo[n][i].Type;
            Temp[k].FlagList.swap(FInfo[n][i].FlagList);
            Temp[k].Basis[0].swap(FInfo[n][i].Basis[0]);
            Temp[k].Basis[1].swap(FInfo[n][i].Basis[1]);
            k++;
        }

        FInfo[n].swap(Temp);
    }
        
    cout << endl;
    return;
}

//Calculates the products of flags.
//Uses FI.FlagList[] and H[] to fill in pFProd->Data[][][].
bool CalculateFlagProducts()
{
    long    i, j, k, n, p;
    long    h, hTotal;
    long    orderOfType, orderOfFlag;
    KeyFlag key; 
    ValFlag val, *pValT, *pValF1, *pValF2;
    KeyProd keyProd;
    long*   pValProd;
    long    mapping1[9], mapping2[9];
    CGraph  f1, f2;
    long    loopError, globalError;
    long    oldOrder, oldP; ///Used to speed up displaying of the progress.

    cout << "Calculating the flag products..." << endl;

    hTotal = H.size();
    
    //Create the hash tables.
    {
        cout << "\rCreating hash tables..." << flush;

        //Set the random values used by the hash functions.
        for(i=0;i<9;i++)
            KeyFlag::hashOrder[i] = GetRandomInt();

        for(i=0;i<9;i++)
        for(j=0;j<9;j++)
        for(k=0;k<9;k++)
            KeyFlag::hashM[i][j][k] = GetRandomInt();
                
        //Add the flags to the htFlags hash table.
        for(n=0;n<=orderOfH-2;n++)
        for(i=0;i<FInfo[n].size();i++)
        {
            //Add the type.
            {
                key.typeOrder = n;
                key.flag      = FInfo[n][i].Type;
                val.iType     = i;
                val.iFlag     = -1;

                htFlags.Add(key,val);
            }

            //Add the flags.
            for(j=0;j<FInfo[n][i].FlagList.size();j++)
            {
                key.typeOrder = n;
                key.flag      = FInfo[n][i].FlagList[j];
                val.iType     = i;
                val.iFlag     = j;

                htFlags.Add(key,val);
            }
        }

        //Create hash tables to store the products of flags.
        htProd.resize(hTotal);
    }

    //Used to speed up displaying of the progress.
    oldOrder = -1;
    oldP     = -1;

    for(orderOfType=0;orderOfType<=orderOfH-2;orderOfType++)
    {
        cout << "\rChecking order " << orderOfType << " of " << orderOfH-2 << " : " << flush;

        if(FInfo[orderOfType].size()<=0)
            continue;

        orderOfFlag = (orderOfH+orderOfType)/2;

        //Use Iso[p][] to choose vertices which will be part of our flags.
        for(p=0;p<Factorial[orderOfH];p++)
        {
            //Create a mapping from Iso[p][]. Vertex v in the flag f1 is the vertex
            //"mapping1[v]" in H[].

            //The induced type of the flag are given by vertices 0, 1, ..., orderOfType-1, which
            //will be vertices Iso[p][0], Iso[p][1], ..., Iso[p][orderOfType-1], in H[].
            for(i=0;i<orderOfType;i++)
            {
                mapping1[i] = Iso[p][i];
                mapping2[i] = Iso[p][i];
            }

            //Map the remaining vertices of the flags.
            for(i=orderOfType;i<orderOfFlag;i++)
            {
                //The vertices orderOfType, orderOfType+1, ..., orderOfFlag-1 of f1 are the
                //vertices Iso[p][orderOfType], Iso[p][orderOfType+1], ..., Iso[p][orderOfFlag-1]
                //in H[].
                mapping1[i] = Iso[p][i];

                //The vertices orderOfType, orderOfType+1, ..., orderOfFlag-1 of f2 are the
                //vertices Iso[p][orderOfFlag], Iso[p][orderOfFlag+1], ...,
                //Iso[p][2*orderOfFlag-orderOfType-1] in H[].
                mapping2[i] = Iso[p][i+orderOfFlag-orderOfType];
            }

            //To avoid repetition of work we will only consider mappings where the following hold:
            //mapping1[orderOfType] < mapping2[orderOfType],
            //mapping1[orderOfType+i] < mapping1[orderOfType+j], for all i<j,
            //mapping2[orderOfType+i] < mapping2[orderOfType+j], for all i<j.
            if(mapping1[orderOfType]>mapping2[orderOfType])
                continue;

            for(i=orderOfType;i<orderOfFlag-1;i++)
            {
                if(mapping1[i]>mapping1[i+1])
                    break;

                if(mapping2[i]>mapping2[i+1])
                    break;
            }

            //Check if we broke out the for loop early.
            if(i!=orderOfFlag-1)
                continue;

            if(orderOfType>oldOrder || p>oldP+10)
            {
                oldOrder = orderOfType;
                oldP     = p;

                cout << "\rChecking order " << orderOfType << " of " << orderOfH-2 << " : ";
                cout << "Applying mapping " << p+1 << " of " << Factorial[orderOfH] << ".     " << flush;
            }

            //Apply the mappings to each H to get two flags f1, f2.
            globalError = 0; //No error
            #pragma omp parallel for schedule(dynamic) default(shared) private(h,i,j,k,f1,f2,key,pValT,pValF1,pValF2,keyProd,pValProd,loopError)
            for(h=0;h<hTotal;h++)
            {
                try
                {
                    #pragma omp critical (CalcProd_ErrorHandling)
                    {
                        loopError = globalError;
                    }

                    if(loopError==0)
                    {
                        //Create the two flags.
                        f1.nV = orderOfFlag;
                        f2.nV = orderOfFlag;

                        for(i=0;i<orderOfFlag;i++)
                        for(j=0;j<orderOfFlag;j++)
                        for(k=0;k<orderOfFlag;k++)
                        {
                            f1.M[i][j][k] = H[h].M[mapping1[i]][mapping1[j]][mapping1[k]];
                            f2.M[i][j][k] = H[h].M[mapping2[i]][mapping2[j]][mapping2[k]];
                        }
                        
                        //Check if the type matches something in FInfo[orderOfType][].
                        key.flag      = f1;
                        key.flag.nV   = orderOfType;
                        key.typeOrder = orderOfType;
                
                        pValT = htFlags.Get(key);

                        if(pValT!=NULL)
                        {
                            //Find the isomorphic copy of each flag in FI.FlagList[]
                            FindMinIso(f1,f1,orderOfType);
                            FindMinIso(f2,f2,orderOfType);

                            key.flag = f1;
                            pValF1   = htFlags.Get(key);

                            key.flag = f2;
                            pValF2   = htFlags.Get(key);

                            if(pValF1!=NULL && pValF2!=NULL)
                            {
                                keyProd.nType  = orderOfType;
                                keyProd.iType  = pValT->iType;

                                if(pValF1->iFlag<=pValF2->iFlag)
                                {
                                    keyProd.iFlag1 = pValF1->iFlag;
                                    keyProd.iFlag2 = pValF2->iFlag;
                                }
                                else
                                {
                                    keyProd.iFlag1 = pValF2->iFlag;
                                    keyProd.iFlag2 = pValF1->iFlag;
                                }

                                //Store the amount to increment by in i.
                                i = Factorial[orderOfFlag-orderOfType]*Factorial[orderOfFlag-orderOfType];
                                if(pValF1->iFlag==pValF2->iFlag)
                                    i *= 2;
                        
                                pValProd = htProd[h].Get(keyProd);

                                if(pValProd==NULL)
                                    htProd[h].Add(keyProd,i); //Add a new entry.
                                else
                                    *pValProd += i;
                            }
                        }
                    }
                }
                catch(bad_alloc&)
                {
                    loopError = 2; //bad_alloc error.
                }
                catch(...)
                {
                    loopError = 1; //Unknown error.
                }
        
                #pragma omp critical (CalcProd_ErrorHandling)
                {
                    if(loopError>globalError)
                        globalError = loopError;
                }
            }

            //Display any errors.
            if(globalError!=0)
            {
                cout << endl << endl;

                if(globalError==1)
                    cout << "ERROR : An unknown error occurred.";

                if(globalError==2)
                    cout << "ERROR : Allocating memory failed.";
            
                cout << endl << endl;
                return false;
            }
        }
    }
    
    cout << "\rDone.                                                     " << flush;
    cout << endl;
    cout << endl;

    return true;
}

//Creates a ".dat-s" file that can be used as input to a SDP solver.
//Returns false if some sort of error was encountered, and true if successful.
//This function formulates the optimization as a primal sdp problem.
bool CreatePrimalSDPFile()
{
    long     i, j, m, n, t;
    long     b, v;
    long     r, c;
    long     nVar, nBlocks, noOfH;
    long     coeff;
    long     product;
    ofstream FileOut;
    long*    pVal;
    KeyProd  key;
    
    //Used to buffer the output.
    long         nBuf;
    const long   maxBuf = 1000000;
    vector<char> memBuffer;
    char*        buffer;

    //Create the ".dat-s" file for an SDP solver.

    //Create a character array to store the output to the file.
    {
        //maxBuf    = The amount the buffer is allowed to be filled before a write operation to the file.
        //memBuffer = A vector which allocates (and automatically deallocates) the memory for the buffer.
        //buffer    = A pointer that allows us to use the memory as if it was a C-style array.
        //nBuf      = The amount of the buffer that is filled.

        //Add some extra space to account for characters which may be
        //written during a write operation which exceeds the buffer size.
        memBuffer.resize(maxBuf+1000);

        buffer = &memBuffer[0];
        nBuf   = 0;
    }
    
    //Calculate the number of blocks and the number of variables we are optimizing over.

    //We have one constraint matrix for each H.
    noOfH = H.size();
    nVar  = noOfH; 

    nBlocks = 0;
    for(n=0;n<=orderOfH-2;n++)
    for(t=0;t<FInfo[n].size();t++)
    for(j=0;j<2;j++)
        if(FInfo[n][t].Basis[j].size()>0)
            nBlocks++;

    //As well as the matrices which we are constraining to be positive semidefinite, we have
    //a diagonal matrix of slack variables and $\rho$ ($\rho$ represents the density). $\rho$
    //will be positioned in the last entry of the diagonal matrix. 
    nBlocks++;

    cout << "Opening the \".dat-s\" file..." << flush;
    
    FileOut.open(FilenameSDP.c_str(),ios::trunc);

    if(FileOut.good()==false)
    {
        cout << endl << endl;
        cout << "ERROR: Cannot create the file \"" << FilenameSDP << "\"." << endl;
        cout << endl;
        return false;
    }

    ForceForm = 1;
    cout << "\rCreating the (primal) SDP problem...   " << endl;

    //Output the number of blocks and variables we will be optimizing over.
    nBuf += sprintf(buffer+nBuf,"%li\n%li\n",nVar,nBlocks);

    //Output the dimensions of the blocks.
    for(n=0;n<=orderOfH-2;n++)
    for(t=0;t<FInfo[n].size();t++)
    for(j=0;j<2;j++)
        if(FInfo[n][t].Basis[j].size()>0)
        {
            nBuf += sprintf(buffer+nBuf,"%li ",(long)FInfo[n][t].Basis[j].size());
            if(nBuf>=maxBuf)
            {
                FileOut.write(buffer,nBuf);
                nBuf = 0;
            }
        }

    //Output the dimension of the block that contains the slack variables and $\rho$.
    //The block is diagonal and we encode this information by prefixing a minus sign.
    nBuf += sprintf(buffer+nBuf,"-%li\n",noOfH+1);

    //Set the constant variables "a"
    for(i=0;i<noOfH;i++)
    {
        //Set the value of the entry, which is the "-density of H[i]" =
        //NumberOfEdges(H[i])/(Factorial[orderOfH] choose 3) multiplied by
        //Factorial[orderOfH] (in order to keep all entries integers).
        nBuf += sprintf(buffer+nBuf,"%li",-NumberOfEdges(H[i])*Factorial[orderOfH-3]*6);

        if(i!=noOfH-1)
            nBuf += sprintf(buffer+nBuf," ");
        else
            nBuf += sprintf(buffer+nBuf,"\n");

        if(nBuf>=maxBuf)
        {
            FileOut.write(buffer,nBuf);
            nBuf = 0;
        }
    }

    //Set the matrix entries. Entries take the form:
    //"Matrix number" "Block number" row col value
    //Indexes start at 1 not 0.

    //Set the constant matrix C, so that the objective function we are trying to maximize is $-\rho$.
    nBuf += sprintf(buffer+nBuf,"0 %li %li %li -1\n",nBlocks,noOfH+1,noOfH+1);
    if(nBuf>=maxBuf)
    {
        FileOut.write(buffer,nBuf);
        nBuf = 0;
    }

    //Set the other constraint matrices.
    for(v=0;v<noOfH;v++)
    {
        cout << "\rCalculating constraint matrix " << v+1 << " of " << noOfH << "." << flush; 

        b = 1; //Current block number.
        for(n=0;n<=orderOfH-2;n++)
        for(t=0;t<FInfo[n].size();t++)
        {
            //Partially set the key (for looking up the product of flags). 
            key.nType  = n;
            key.iType  = t;

            for(m=0;m<2;m++)
            {
                if(FInfo[n][t].Basis[m].size()<=0)
                    continue; //Ignore blocks of dimension 0.

                //r = row index,
                //c = column index.
                for(r=0;r<FInfo[n][t].Basis[m].size();r++)
                for(c=r;c<FInfo[n][t].Basis[m].size();c++)
                { 
                    //Multiply basis element r with element c and store the result in product.

                    //Clear product.
                    product = 0;

                    //Multiply flag i of element r with flag j of element c. If m = 1 then we
                    //also must take into account the coefficients of the flags which are 1 when
                    //i,j = 0, and -1 when i,j = 1. 
                    for(i=0;i<FInfo[n][t].Basis[m][r].size();i++)
                    for(j=0;j<FInfo[n][t].Basis[m][c].size();j++)
                    {
                        if(m==0 || (m==1 && i==j))
                            coeff = 1;
                        else 
                            coeff = -1;

                        //Look up the product of the flags.

                        //Set the key.
                        if(FInfo[n][t].Basis[m][r][i]<=FInfo[n][t].Basis[m][c][j])
                        {
                            key.iFlag1 = FInfo[n][t].Basis[m][r][i];
                            key.iFlag2 = FInfo[n][t].Basis[m][c][j];
                        }
                        else
                        {
                            key.iFlag1 = FInfo[n][t].Basis[m][c][j];
                            key.iFlag2 = FInfo[n][t].Basis[m][r][i];
                        }
                                        
                        pVal = htProd[v].Get(key);
                                        
                        if(pVal!=NULL)
                            product += coeff*(*pVal);
                    }

                    //Set the entry.
                    if(product==0)
                        continue; //We only need to set the non-zero entries.

                    nBuf += sprintf(buffer+nBuf,"%li %li %li %li %li\n",v+1,b,r+1,c+1,product);
                    if(nBuf>=maxBuf)
                    {
                        FileOut.write(buffer,nBuf);
                        nBuf = 0;
                    }
                }

                b++;
            }
        }

        //Set the slack variable.
        nBuf += sprintf(buffer+nBuf,"%li %li %li %li %li\n",v+1,nBlocks,v+1,v+1,Factorial[orderOfH]);

        //Set the coefficient of $\rho$.
        nBuf += sprintf(buffer+nBuf,"%li %li %li %li %li\n",v+1,nBlocks,noOfH+1,noOfH+1,-Factorial[orderOfH]);

        if(nBuf>=maxBuf)
        {
            FileOut.write(buffer,nBuf);
            nBuf = 0;
        }
    }

    FileOut.write(buffer,nBuf);
    nBuf = 0;
    
    FileOut.close();

    cout << "\r                                   ";
    for(i=0;i<NumberOfDigits(noOfH);i++)
        cout << "  ";

    cout << "\rDone." << endl;
    cout << endl;
    cout << endl;

    return true;
}

//Creates a ".dat-s" file that can be used as input to a SDP solver.
//Returns false if some sort of error was encountered, and true if successful.
//This function formulates the optimization as a dual sdp problem.
bool CreateDualSDPFile()
{
    long     i, j, k, m, n, t;
    long     b, v;
    long     r, c;
    long     nVar, nBlocks, noOfH;
    long     coeff;
    long     product;
    ofstream FileOut;
    long*    pVal;
    KeyProd  key;
    
    //Used to buffer the output.
    long         nBuf;
    const long   maxBuf = 1000000;
    vector<char> memBuffer;
    char*        buffer;

    //Create the ".dat-s" file for an SDP solver.

    //Create a character array to store the output to the file.
    {
        //maxBuf    = The amount the buffer is allowed to be filled before a write operation to the file.
        //memBuffer = A vector which allocates (and automatically deallocates) the memory for the buffer.
        //buffer    = A pointer that allows us to use the memory as if it was a C-style array.
        //nBuf      = The amount of the buffer that is filled.

        //Add some extra space to account for characters which may be
        //written during a write operation which exceeds the buffer size.
        memBuffer.resize(maxBuf+1000);

        buffer = &memBuffer[0];
        nBuf   = 0;
    }
    
    //Calculate the number of blocks and the number of variables we are optimizing over.
    noOfH = H.size();

    nVar    = 0; 
    nBlocks = 0;
    for(n=0;n<=orderOfH-2;n++)
    for(t=0;t<FInfo[n].size();t++)
    for(j=0;j<2;j++)
        if(FInfo[n][t].Basis[j].size()>0)
        {
            nBlocks++;
            
            //We only need to determine the upper right half of the positive semidefinite
            //matrices as they are symmetric.
            nVar += FInfo[n][t].Basis[j].size()*(FInfo[n][t].Basis[j].size()+1)/2;
        }

    //We have an extra variable $\rho$, that represents the upper bound of the Turan density,
    //which we are trying to minimize.
    nVar++;

    //As well as the matrices which we are constraining to be positive semidefinite, we have
    //a series of lower bounds on $\rho$ that must also be satisfied. We encode this in an
    //extra block.
    nBlocks++;

    cout << "Opening the \".dat-s\" file..." << flush;
    
    FileOut.open(FilenameSDP.c_str(),ios::trunc);

    if(FileOut.good()==false)
    {
        cout << endl << endl;
        cout << "ERROR: Cannot create the file \"" << FilenameSDP << "\"." << endl;
        cout << endl;
        return false;
    }

    ForceForm = -1;
    cout << "\rCreating the (dual) SDP problem...   " << endl;

    //Output the number of blocks and variables we will be optimizing over.
    nBuf += sprintf(buffer+nBuf,"%li\n%li\n",nVar,nBlocks);

    //Output the dimensions of the blocks.
    for(n=0;n<=orderOfH-2;n++)
    for(t=0;t<FInfo[n].size();t++)
    for(j=0;j<2;j++)
        if(FInfo[n][t].Basis[j].size()>0)
        {
            nBuf += sprintf(buffer+nBuf,"%li ",(long)FInfo[n][t].Basis[j].size());
            if(nBuf>=maxBuf)
            {
                FileOut.write(buffer,nBuf);
                nBuf = 0;
            }
        }

    //Output the dimension of the block that represents the inequalities.
    //The block is diagonal and we encode this information by prefixing a minus sign.
    nBuf += sprintf(buffer+nBuf,"-%li\n",noOfH);

    //Set the objective function to $\rho$.
    
    //We take the first variable to be $\rho$.
    nBuf += sprintf(buffer+nBuf,"1");

    //Set the coefficients of all other variables to be 0.
    for(i=0;i<nVar-1;i++)
    {
        nBuf += sprintf(buffer+nBuf," 0");
        if(nBuf>=maxBuf)
        {
            FileOut.write(buffer,nBuf);
            nBuf = 0;
        }
    }
    
    nBuf += sprintf(buffer+nBuf,"\n");

    //Set the matrix entries. Entries take the form:
    //"Matrix number" "Block number" row col value
    //Indexes start at 1 not 0.

    //For each variable we have a matrix that consists of blocks down its diagonal. We start
    //by defining such a matrix associated with 0 th variable which is defined as the
    //constant 1. The only constant terms are $d(H_i)$ which appear in the constraints of
    //$\rho$. We multiply the constraints by Factorial[orderOfH] to keep all the coefficients
    //integers. 
    for(i=0;i<noOfH;i++)
        if(NumberOfEdges(H[i])!=0)
        {
            //We only need to set the non-zero entries of the matrix.
            nBuf += sprintf(buffer+nBuf,"0 %li %li %li %li\n",nBlocks,i+1,i+1,NumberOfEdges(H[i])*Factorial[orderOfH-3]*6);
            if(nBuf>=maxBuf)
            {
                FileOut.write(buffer,nBuf);
                nBuf = 0;
            }
        }

    //Set the matrix corresponding to variable 1, which we've defined to be $\rho$.
    for(i=0;i<noOfH;i++)
    {
        nBuf += sprintf(buffer+nBuf,"1 %li %li %li %li\n",nBlocks,i+1,i+1,Factorial[orderOfH]);
        if(nBuf>=maxBuf)
        {
            FileOut.write(buffer,nBuf);
            nBuf = 0;
        }
    }

    //Set the matrices for the other variables.

    b = 1; //Current block number.
    v = 2; //Current variable number.
    for(n=0;n<=orderOfH-2;n++)
    for(t=0;t<FInfo[n].size();t++)
    {
        //Partially set the key (for looking up the product of flags). 
        key.nType  = n;
        key.iType  = t;

        for(m=0;m<2;m++)
        {
            if(FInfo[n][t].Basis[m].size()<=0)
                continue; //Ignore blocks of dimension 0.

            //r = row index,
            //c = column index.
            for(r=0;r<FInfo[n][t].Basis[m].size();r++)
            for(c=r;c<FInfo[n][t].Basis[m].size();c++)
            {
                cout << "\rCalculating matrix " << v << " of " << nVar << "." << flush; 

                //Set the position of the variable v, in block b, row r, column c.
                nBuf += sprintf(buffer+nBuf,"%li %li %li %li 1\n",v,b,r+1,c+1);
                if(nBuf>=maxBuf)
                {
                    FileOut.write(buffer,nBuf);
                    nBuf = 0;
                }

                //Multiply basis element r with element c and store the result in product.
                for(k=0;k<noOfH;k++)
                {
                    //Clear product.
                    product = 0;

                    //Multiply flag i of element r with flag j of element c. If m = 1 then we
                    //also must take into account the coefficients of the flags which are 1 when
                    //i,j = 0, and -1 when i,j = 1. 
                    for(i=0;i<FInfo[n][t].Basis[m][r].size();i++)
                    for(j=0;j<FInfo[n][t].Basis[m][c].size();j++)
                    {
                        if(m==0 || (m==1 && i==j))
                            coeff = 1;
                        else 
                            coeff = -1;

                        //Look up the product of the flags.

                        //Set the key.
                        if(FInfo[n][t].Basis[m][r][i]<=FInfo[n][t].Basis[m][c][j])
                        {
                            key.iFlag1 = FInfo[n][t].Basis[m][r][i];
                            key.iFlag2 = FInfo[n][t].Basis[m][c][j];
                        }
                        else
                        {
                            key.iFlag1 = FInfo[n][t].Basis[m][c][j];
                            key.iFlag2 = FInfo[n][t].Basis[m][r][i];
                        }
                                        
                        pVal = htProd[k].Get(key);
                                        
                        if(pVal!=NULL)
                            product += coeff*(*pVal);
                    }

                    //Variables not on the diagonal make twice the contribution because the
                    //matrices are symmetric.
                    if(r!=c)
                        product *= 2;

                    //Set the entry.
                    if(product==0)
                        continue; //We only need to set the non-zero entries.

                    nBuf += sprintf(buffer+nBuf,"%li %li %li %li %li\n",v,nBlocks,k+1,k+1,-product);
                    if(nBuf>=maxBuf)
                    {
                        FileOut.write(buffer,nBuf);
                        nBuf = 0;
                    }
                }

                v++;
            }

            b++;
        }
    }

    FileOut.write(buffer,nBuf);
    nBuf = 0;
    
    FileOut.close();

    cout << "\r                        ";
    for(i=0;i<NumberOfDigits(nVar);i++)
        cout << "  ";

    cout << "\rDone." << endl;
    cout << endl;
    cout << endl;

    return true;
}

//A short description of the usage of the program.
void ShortUsage()
{
    cout << endl;
    cout << "Correct Usage :" << endl;
    cout << "    <Executable> <Order> <Graphs> [-i] [-c <SDP command>] [-o <Filename>]" << endl;
    cout << "    <Executable> -s <SDP Command>" << endl;
    cout << endl;
    cout << "Examples :" << endl;
    cout << "    ExactDensityBounder.exe -s csdp.exe" << endl;
    cout << "    ./edb 6 {123 124 134, 234} -o K4" << endl;
    cout << "    ExactDensityBounder 5 {234} 4 : {123 124 134 234} -i -o \"Induced Prob\"" << endl;
    cout << "    ./edb 5 {123 124 134} -c ./sdpa" << endl;
    cout << endl;
    cout << "Run the program with no arguments for more details." << endl;
    cout << endl;

    return;
}

//A description of the usage of the program.
void Usage()
{
    cout << "Name :" << endl;
    cout << "    ExactDensityBounder" << endl;
    cout << endl;
    cout << "Description :" << endl;
    cout << "    Finds an upper bound of the Turan density for a family of 3-graphs." << endl;
    cout << endl;
    cout << "Usage :" << endl;
    cout << "    <Executable> <Order> <Graphs> [-i] [-c <SDP command>] [-o <Filename>]" << endl;
    cout << "    <Executable> -s <SDP Command>" << endl;
    cout << endl;
    cout << "Examples :" << endl;
    cout << "    ExactDensityBounder.exe -s csdp.exe" << endl;
    cout << "    ./edb 6 {123 124 134 234} -o K4" << endl;
    cout << "    ExactDensityBounder 5 {234} 4 : {123 124 134 234} -i -o \"Induced Prob\"" << endl;
    cout << "    ./edb 5 {123 124 134} -c ./sdpa" << endl;
    cout << endl;
    cout << "Arguments :" << endl;
    cout << "    <> : A mandatory field." << endl;
    cout << endl;
    cout << "    [] : An optional field." << endl;
    cout << endl;
    cout << "    <Executable> : The name of the compiled program. For example" << endl;
    cout << "    \"ExactDensityBounder.exe\" on windows, \"./edb\" or \"./a.out\" in linux." << endl;
    cout << endl;
    cout << "    <Order> : The order of the subgraphs H. This should be an integer between" << endl;
    cout << "    3 and 9. Increasing the order will result in a better bound but the" << endl;
    cout << "    computation will be more intensive. Currently orders higher than 9 are not" << endl;
    cout << "    supported by this program." << endl;
    cout << endl;
    cout << "    <Graphs> : A list of 3-uniform hypergraphs (separated by spaces). Graphs" << endl;
    cout << "    with orders higher than 9 are not supported. Each graph should either be" << endl;
    cout << "    of the form \"<Order of graph> : <Edge set>\" or simply \"<Edge set>\"." << endl;
    cout << "    \"<Order of graph>\" should be a single digit and if it is omitted then it is" << endl;
    cout << "    taken to be the largest label of a vertex in the edge set. The vertices are" << endl;
    cout << "    assumed to be labelled by integers starting from 1. The edge set should be" << endl;
    cout << "    enclosed in curly braces and contain a list of 3 digit numbers (separated" << endl;
    cout << "    by spaces) representing the edges. See the \"Examples\" section above." << endl;
    cout << endl;
    cout << "    [-i] : Indicates that we wish to consider only induced subgraphs. If \"-i\"" << endl;
    cout << "    is used then the program only forbids <Graphs> from appearing as *induced*" << endl;
    cout << "    subgraphs. If \"-i\" is omitted then <Graphs> are forbidden from appearing as" << endl;
    cout << "    subgraphs (as in the normal definition of Turan density)." << endl;
    cout << endl;
    cout << "    [-c] : Used to override the default command (that runs the SDP solver) with" << endl;
    cout << "    <SDP command>. See also the \"<SDP command>\" and \"-s\" sections below." << endl;
    cout << endl;
    cout << "    <SDP command> : The command used to run an SDP solver on the underlying" << endl;
    cout << "    semidefinite programming problem. For example \"./csdp\" or \"sdpa.exe\". The" << endl;
    cout << "    program runs the SDP solver by running the command" << endl;
    cout << "    \"<SDP command> <Filename>.dat-s <Filename>.out\". It always assumes the SDP" << endl;
    cout << "    solver completed successfully and that the solution is stored in" << endl;
    cout << "    \"<Filename>.out\". It can currently only read the \".out\" file if it adheres" << endl;
    cout << "    to the same format as csdp or sdpa's output files. See also the sections on" << endl;
    cout << "    \"[-c]\", \"-s\", and \"<Filename>\"." << endl;
    cout << endl;
    cout << "    [-o] : Used to indicate <Filename> as the name to be used for the files the" << endl;
    cout << "    program creates. If <Filename> is not set using \"-o\" then the default name" << endl;
    cout << "    \"Temp\" is used. See also the section on \"<Filename>\"." << endl;
    cout << endl;
    cout << "    <Filename> : A filename (without extension). It is used to name the" << endl;
    cout << "    temporary \".dat-s\" input file to the SDP solver and the \".out\" output file." << endl;
    cout << "    It is also used to name the \".txt\" file which will hold the data which" << endl;
    cout << "    forms the proof. The program DensityChecker can be used to verify the" << endl;
    cout << "    calculation from the \".txt\" file." << endl;
    cout << endl;
    cout << "    -s : Used to save <SDP command> as the default command that runs the SDP" << endl;
    cout << "    solver. The command is saved in the file \"ConfigEDB.cfg\". See the \"[-c]\"" << endl;
    cout << "    and \"<SDP command>\" sections." << endl;
    cout << endl;

    return;
}

//Outputs the full filename of the config file based on the path used to execute the program.
string GetConfigName(const char* argv0)
{
    long   i;
    string Filename;

    Filename = argv0;

    for(i=Filename.size()-1;i>=0;i--)
        if(Filename[i]=='/' || Filename[i]=='\\')
            break;

    Filename.resize(i+1);
    Filename += configName;

    return Filename;
}

//Saves command[] in the config file.
//Returns true if successful.
bool SaveSDPCommand(const char* argv0, const char* command)
{
    ofstream FileConfig;
    string   Filename;

    Filename = GetConfigName(argv0);

    cout << endl;
    cout << "Opening \"" << configName << "\"."<< endl;
    FileConfig.open(Filename.c_str(),ios::trunc);

    if(FileConfig.good()==false)
    {
        cout << endl;
        cout << "ERROR : Could not open \"" << Filename << "\" to write the command." << endl << endl;
        return false;
    }

    cout << "Writing the command..." << endl;
    FileConfig << command << endl;
        
    if(FileConfig.good()==false)
    {
        cout << endl;
        cout << "ERROR : Writing to the file failed." << endl << endl;
        
        FileConfig.close();
        return false;
    }

    cout << "Closing \"" << configName << "\"."<< endl << endl;
    FileConfig.close();

    return true;
}

//Reads 'command' from the config file.
//Returns true if successful.
bool ReadSDPCommand(const char* argv0, string& command)
{
    long     i;
    ifstream FileConfig;
    string   Filename;
    string   blank;

    Filename = GetConfigName(argv0);
    
    //Create a blank string the same size as configName.
    blank.resize(strlen(configName));
    for(i=0;i<blank.size();i++)
        blank[i] = ' ';

    cout << "Reading the default command   : Opening \"" << configName << "\"." << flush;
    FileConfig.open(Filename.c_str());

    if(FileConfig.good()==false)
    {
        cout << endl << endl;
        cout << "ERROR : Could not open \"" << Filename << "\"." << endl;
        cout << "Make sure the default command to invoke the SDP solver has been set by running" << endl;
        cout << "the program with the \"-s\" option." << endl;
        ShortUsage();
        return false;
    }

    cout << "\rReading the default command   :          " << blank << "  ";
    cout << "\rReading the default command   : Reading the command..." << flush;
    getline(FileConfig,command);
        
    cout << "\rReading the default command   :                       ";
    cout << "\rReading the default command   : Closing \"" << configName << "\"." << flush;
    FileConfig.close();

    //Clear the line.
    cout << "\r                                         " << blank << "  \r" << flush;
    
    return true;
}

//Saves str[] in the file filename[].
//Used to test the command interpreter's argument format.
//Returns true if successful.
bool SaveFormatTest(const char* filename, const char* str)
{
    ofstream File;

    File.open(filename,ios::trunc);

    if(File.good()==false)
    {
        cout << endl;
        cout << "ERROR : Could not open \"" << filename << "\" for writing." << endl;
        return false;
    }

    File << str << endl;
        
    if(File.good()==false)
    {
        cout << endl;
        cout << "ERROR : Writing to the file \"" << filename << "\" failed." << endl;
        File.close();
        return false;
    }

    File.close();
    return true;
}

//Reads str in the file filename[].
//Used to test the command interpreter's argument format.
//Returns true if successful.
bool ReadFormatTest(const char* filename, string& str)
{
    ifstream File;

    File.open(filename);

    if(File.good()==false)
    {
        cout << endl;
        cout << "ERROR : Could not open \"" << filename << "\" for reading." << endl << endl;
        return false;
    }

    getline(File,str);

    File.close();
    return true;
}

//Parses the command line arguments passed to the program via main().
//Returns true if the program should continue and construct the ".dat-s" file.
//Returns false if the program should terminate.
bool ProcessArguments(long  argc, const char* const argv[])
{
    long   i, j, k, n;
    long   aN, cN;
    bool   bMidway, bReadEdgeSet, bReadColon;
    bool   bError;
    CGraph g;
    string Filename;

    //If the second argument is "Test" then the program has been called to test the format of the system() call.
    if(argc==3 && string(argv[1])=="Test")
    {
        SaveFormatTest(argv[2],"<Modify>");
        return false;
    }

    cout << endl;

    if(argc<=1)
    {
        Usage();
        return false;
    }

    //argv[0] is the application name.
    
    //argv[1] should be a single digit integer or -s.
    bError = false;

    if(strlen(argv[1])==1)
    {
        if('3'<=argv[1][0] && argv[1][0]<='9')
            orderOfH = argv[1][0]-'0';
        else
            bError = true;
    }
    else
    if(strlen(argv[1])==2)
    {
        if(argv[1][0]=='-' && (argv[1][1]=='s' || argv[1][1]=='S'))
        {
            cout << "Saving the default SDP command : ";

            if(argc<=2)
            {
                cout << endl;
                cout << endl;
                cout << "ERROR : Expected a second argument indicating the SDP command to save." << endl;
                ShortUsage();
                return false;
            }

            cout << "\"" << argv[2] << "\"" << endl;

            if(argc>3)
            {
                cout << endl;
                cout << "ERROR : Unexpected extra arguments." << endl;
                ShortUsage();
                return false;
            }

            //Save the command.
            if(SaveSDPCommand(argv[0],argv[2])==false)
                return false;

            //Even though we succeeded, return false to cause the program to terminate.
            return false;
        }
        else
            bError = true;
    }
    else
        bError = true;

    cout << "OpenMP (multiprocessing)      : ";
#ifdef _OPENMP
    cout << "Enabled (" << omp_get_num_procs() << " processors)." << endl;
#else
    cout << "Disabled" << endl;
#endif

    if(bError==true)
    {
        cout << endl;
        cout << "ERROR : Expected the first argument to be \"-s\" or an integer between 3 and 9." << endl;
        ShortUsage();
        return false;
    }
    
    cout << "Order of the subgraphs H      : " << orderOfH << endl; 
    
    //Read in the graphs.
    {
        ForbiddenGraphs.clear();

        cout << "Forbidden graphs              : " << endl;

        aN      = 1;
        cN      = 0;
        bMidway = false;
        while(true)
        {
            //Move to the next character.
            cN++;
            while(argv[aN][cN]=='\0')
            {
                aN++;
                cN = 0;

                if(aN>=argc)
                    break;
            }

            //Check if the arguments have ended.
            if(aN>=argc)
            {
                //Reached the end of the arguments.
                if(bMidway==true)
                {
                    cout << endl;
                    cout << "ERROR : Input unexpectedly ended midway through a definition of a graph." << endl;
                    ShortUsage();
                    return false;
                }

                break;
            }

            //Check if we've reached the end of the graph list.
            if(cN==0 && bMidway==false && argv[aN][cN]=='-')
                break;

            if(bMidway==false)
            {
                //Expect to read a comma, space, integer, ':' or '{'.

                if(argv[aN][cN]==',' || argv[aN][cN]==' ')
                    continue;
                else
                if('0'<=argv[aN][cN] && argv[aN][cN]<='9')
                {
                    g.nV = argv[aN][cN]-'0';

                    for(i=0;i<9;i++)
                    for(j=0;j<9;j++)
                    for(k=0;k<9;k++)
                        g.M[i][j][k] = false;

                    bMidway      = true;
                    bReadColon   = true;
                    bReadEdgeSet = false;

                    continue;
                }
                else
                if(argv[aN][cN]==':')
                {
                    g.nV = -1;

                    for(i=0;i<9;i++)
                    for(j=0;j<9;j++)
                    for(k=0;k<9;k++)
                        g.M[i][j][k] = false;
                    
                    bMidway      = true;
                    bReadColon   = false;
                    bReadEdgeSet = false;

                    continue;
                }
                else
                if(argv[aN][cN]=='{')
                {
                    g.nV = -1;

                    for(i=0;i<9;i++)
                    for(j=0;j<9;j++)
                    for(k=0;k<9;k++)
                        g.M[i][j][k] = false;
                    
                    bMidway      = true;
                    bReadEdgeSet = true;

                    continue;
                }
                else
                {
                    bError = true;
                    break;
                }
            }

            if(bMidway==true)
            {
                if(bReadEdgeSet==false)
                {
                    //Read the separating colon or the starting edge set brace '{'.

                    if(bReadColon==true)
                    {
                        //Read the separating colon.

                        //Ignore spaces
                        if(argv[aN][cN]==' ')
                            continue;
                        else
                        if(argv[aN][cN]==':')
                        {
                            bReadColon = false;
                            continue;
                        }
                        else
                        {
                            bError = true;
                            break;
                        }
                    }
                    else
                    {
                        //Read the starting brace.
                                                
                        //Ignore spaces
                        if(argv[aN][cN]==' ')
                            continue;
                        else
                        if(argv[aN][cN]=='{')
                        {
                            bReadEdgeSet = true;
                            continue;
                        }
                        else
                        {
                            bError = true;
                            break;
                        }
                    }
                }

                if(bReadEdgeSet==true)
                {
                    //Read the edges.

                    //Expects to read an integer, space, comma, or a closing brace.

                    if(argv[aN][cN]=='}')
                    {
                        bMidway = false;

                        if(g.nV==-1)
                        {
                            //Find the largest index
                            for(i=0;i<9;i++)
                            for(j=i+1;j<9;j++)
                            for(k=j+1;k<9;k++)
                            {
                                if(g.M[i][j][k]==false)
                                    continue;

                                if(k>g.nV)
                                    g.nV = k;
                            }

                            g.nV++;
                        }

                        cout << g << endl;

                        //Save the graph
                        ForbiddenGraphs.push_back(g);

                        continue;
                    }
                    else
                    if('1'<=argv[aN][cN] && argv[aN][cN]<='9')
                    {
                        //Read the edge "ijk".
                        i = argv[aN][cN]-'1';
                        
                        cN++;
                        if(argv[aN][cN]<'1' || '9'<argv[aN][cN])
                        {
                            bError = true;
                            break;
                        }

                        j = argv[aN][cN] -'1';

                        cN++;
                        if(argv[aN][cN]<'1' || '9'<argv[aN][cN])
                        {
                            bError = true;
                            break;
                        }

                        k = argv[aN][cN] -'1';

                        //The next character shouldn't be a digit.
                        if('1'<=argv[aN][cN+1] && argv[aN][cN+1]<='9')
                        {
                            bError = true;
                            break;
                        }

                        //ijk should be distinct.
                        if(i==j || i==k || j==k)
                        {
                            cout << endl;
                            cout << "ERROR : The vertices in an edge are not all distinct." << endl;
                            ShortUsage();
                            return false;
                        }

                        //Check the vertices aren't too large.
                        if(g.nV!=-1 && (g.nV<=i || g.nV<=j || g.nV<=k))
                        {
                            cout << endl;
                            cout << "ERROR : The edge set and the order of the graph are inconsistent." << endl;
                            ShortUsage();
                            return false;
                        }

                        //Store the edge.
                        g.M[i][j][k] = true;
                        g.M[i][k][j] = true;
                        g.M[j][i][k] = true;
                        g.M[j][k][i] = true;
                        g.M[k][i][j] = true;
                        g.M[k][j][i] = true;

                        continue;
                    }
                    else
                    if(argv[aN][cN]==' ' || argv[aN][cN]==',')
                        continue;
                    else
                    {
                        bError = true;
                        break;
                    }
                }
            }
        }

        if(bError==true)
        {
            cout << endl;
            cout << "ERROR : Unexpected character in the graph list." << endl;
            ShortUsage();
            return false;
        }
    }

    //Read in the switches.
    bool bI  = false;
    bool bC  = false;
    bool bO  = false;
    bool bT  = false;
    bool bPD = false;
    bool bR  = false;
    bool bA  = false;
    bool bB  = false;

    try
    {
        while(true)
        {
            if(aN>=argc)
                break;

            if(strlen(argv[aN])!=2 || argv[aN][0]!='-'
            ||(argv[aN][1]!='i' && argv[aN][1]!='I'
            && argv[aN][1]!='c' && argv[aN][1]!='C'
            && argv[aN][1]!='o' && argv[aN][1]!='O'
            && argv[aN][1]!='t' && argv[aN][1]!='T'
            && argv[aN][1]!='p' && argv[aN][1]!='P'
            && argv[aN][1]!='d' && argv[aN][1]!='D'
            && argv[aN][1]!='r' && argv[aN][1]!='R'
            && argv[aN][1]!='a' && argv[aN][1]!='A'
            && argv[aN][1]!='b' && argv[aN][1]!='B'))
                throw 0;

            if(argv[aN][1]=='i' || argv[aN][1]=='I')
            {
                if(bI==true)
                    throw 0;
                else
                    bI = true;

                bForbidInduced = true;
                aN++;

                continue;
            }

            if(argv[aN][1]=='c' || argv[aN][1]=='C')
            {
                if(bC==true)
                    throw 0;
                else
                    bC = true;

                aN++;
                if(aN>=argc)
                {
                    cout << endl;
                    cout << endl;
                    cout << "ERROR : Expected another argument indicating the command to run the SDP solver." << endl;
                    ShortUsage();
                    return false;
                }

                SDPCommand = argv[aN];
                aN++;

                continue;
            }
        
            if(argv[aN][1]=='o' || argv[aN][1]=='O')
            {
                if(bO==true)
                    throw 0;
                else
                    bO = true;

                aN++;
                if(aN>=argc)
                {
                    cout << endl;
                    cout << endl;
                    cout << "ERROR : Expected another argument indicating the filename." << endl;
                    ShortUsage();
                    return false;
                }

                Filename = argv[aN];
                aN++;

                continue;
            }
        
            //Set the type list.
            //expects arguments of the form:
            //+ 4.5 7.13
            //- 0.1
            //+   = Remove all types except the ones specified.
            //-   = Include all types except the ones specified.
            //n.x = n is the order of the type (a single digit), x is the index (starting from 1).
            if(argv[aN][1]=='t' || argv[aN][1]=='T')
            {
                if(bT==true)
                    throw 0;
                else
                    bT = true;

                //Read the types to add/remove.
                bool bReadFirstArg, bReadSecondArg;

                bDisplayTypeNumber = true;

                //Clear the modifiedTypeList array.
                for(i=0;i<=9-2;i++)
                    modifyTypeList[i].clear();

                bReadFirstArg = false;
                while(true)
                {
                    aN++;

                    if(aN>=argc)
                    {
                        if(bReadFirstArg==false)
                        {
                            cout << endl << endl;
                            cout << "ERROR : Expected another argument indicating the default treatment of types." << endl;
                            ShortUsage();
                            return false;
                        }
                        else
                            break; //Successfully read all the arguments.
                    }

                    if(bReadFirstArg==false)
                    {
                        if(strlen(argv[aN])!=1 || (argv[aN][0]!='+' && argv[aN][0]!='-'))
                        {
                            cout << endl << endl;
                            cout << "ERROR : Expected to read \"+\" or \"-\"." << endl;
                            ShortUsage();
                            return false;
                        }

                        if(argv[aN][0]=='+')
                            bIncludeTypesByDefault = false;
                        else
                            bIncludeTypesByDefault = true;

                        bReadFirstArg  = true;
                        bReadSecondArg = false;
                        continue;
                    }

                    //Read in a type.
                    //"order.index"
                    {
                        //Check if the next switch has been encountered.
                        if(strlen(argv[aN])==2 && argv[aN][0]=='-')
                            break; //Finished reading the types.

                        if(strlen(argv[aN])<3)
                            throw 0;

                        //Read the order.
                        n = argv[aN][0]-'0';
                        if(n<0 || orderOfH-2<n)
                            throw 0;
                
                        //Read '.'
                        if(argv[aN][1]!='.')
                            throw 0;

                        //Read the index.

                        cN = 2;
                        j  = 0;
                        while(argv[aN][cN]!='\0')
                        {
                            i = argv[aN][cN]-'0';
                            if(i<0 || 9<i)
                                throw 0;

                            j = j*10+i;

                            cN++;
                        }

                        if(j<=0)
                            throw 0;

                        modifyTypeList[n].push_back(j-1);

                        if(bReadSecondArg==false)
                            bReadSecondArg = true;
                    }
                }

                continue;
            }
        
            //Set the form of the SDP problem.
            if(argv[aN][1]=='p' || argv[aN][1]=='P' || argv[aN][1]=='d' || argv[aN][1]=='D')
            {
                if(bPD==true)
                    throw 0;
                else
                    bPD = true;

                if(argv[aN][1]=='p' || argv[aN][1]=='P')
                    ForceForm = 1;
                else
                    ForceForm = -1;

                aN++;

                continue;
            }

            //Set the rounding scheme.
            if(argv[aN][1]=='r' || argv[aN][1]=='R')
            {
                if(bR==true)
                    throw 0;
                else
                    bR = true;

                bDefaultRoundingScheme = false;

                aN++;

                continue;
            }

            //Set the accuracy for ProjectSolution().
            if(argv[aN][1]=='a' || argv[aN][1]=='A')
            {
                if(bA==true)
                    throw 0;
                else
                    bA = true;

                aN++;

                if(aN>=argc)
                {
                    cout << endl;
                    cout << endl;
                    cout << "ERROR : Expected another argument specifying the number of decimal places." << endl;
                    ShortUsage();
                    return false;
                }

                //Read the number of decimal places which should be an integer between 1 and 1000 inclusive.
                {
                    //Remove white spaces.
                    cN = 0;
                    while(true)
                    {
                        if(argv[aN][cN]==' ' || argv[aN][cN]=='\t' || argv[aN][cN]=='\r' || argv[aN][cN]=='\n')
                        {
                            cN++;
                            continue;
                        }

                        break;
                    }

                    //Read the sign.
                    if(argv[aN][cN]=='-')
                        throw 0;

                    if(argv[aN][cN]=='+')
                        cN++;

                    //Check that the first character is a digit.
                    if(argv[aN][cN]-'0'<0 || 9<argv[aN][cN]-'0')
                        throw 0;

                    //Read the digits into n, until a non-digit is encountered (or n>1000).
                    n = 0;
                    while(true)
                    {
                        if(0<=argv[aN][cN]-'0' && argv[aN][cN]-'0'<=9)
                        {
                            n = 10*n+argv[aN][cN]-'0';

                            if(n>1000)
                                throw 0;

                            cN++;
                            continue;
                        }

                        break;
                    }

                    //Check n is not 0.
                    if(n<=0)
                        throw 0;

                    //Check that the rest of the characters are white spaces.
                    while(true)
                    {
                        if(argv[aN][cN]=='\0')
                            break;

                        if(argv[aN][cN]==' ' || argv[aN][cN]=='\t' || argv[aN][cN]=='\r' || argv[aN][cN]=='\n')
                        {
                            cN++;
                            continue;
                        }

                        throw 0;
                    }
                }

                accuracyProjectSolution = n;
                aN++;

                continue;
            }

            if(argv[aN][1]=='b' || argv[aN][1]=='B')
            {
                if(bB==true)
                    throw 0;
                else
                    bB = true;

                aN++;

                continue;
            }
        }

        //Set the defaults.
        {
            if(bI==false)
                bForbidInduced = false;

            if(bC==false)
            {
                if(ReadSDPCommand(argv[0],SDPCommand)==false)
                    return false;
            }

            if(bO==false)
                Filename = defaultFilename;

            FilenameSDP = Filename+".dat-s";
            FilenameOut = Filename+".out";
            FilenameTxt = Filename+".txt";

            if(bT==false)
            {
                //Do not remove any types.
                bDisplayTypeNumber     = false;
                bIncludeTypesByDefault = true;

                //Clear the modifiedTypeList array.
                for(i=0;i<=9-2;i++)
                    modifyTypeList[i].clear();
            }

            if(bPD==false)
                ForceForm = 0;

            if(bB==false && bForbidInduced==false)
                bForbidBlowups = true;
            else
                bForbidBlowups = false;
        }

        //Display the Results.
        {
            cout << "Forbid only induced subgraphs : " << (bForbidInduced?"True":"False") << endl;
            cout << "SDP command                   : \"" << SDPCommand << "\"" << endl;
            cout << "Filename                      : \"" << Filename << "\"" << endl;

            if(bT==true)
            {
                cout << "Types specified               : " << (bIncludeTypesByDefault?"All":"None");

                for(n=0;n<=9-2;n++)
                {
                    if(modifyTypeList[n].size()<=1)
                        continue;

                    //Order modifyTypeList[n].
                    for(i=0;i<modifyTypeList[n].size();i++)
                    for(j=i+1;j<modifyTypeList[n].size();j++)
                        if(modifyTypeList[n][i]>modifyTypeList[n][j])
                        {
                            //Swap.
                            k                    = modifyTypeList[n][i];
                            modifyTypeList[n][i] = modifyTypeList[n][j];
                            modifyTypeList[n][j] = k;
                        }
                }
                                
                bool bFirstType = false;
                for(n=0;n<=9-2;n++)
                {
                    if(bFirstType==false && modifyTypeList[n].size()>0)
                    {
                        cout << " except";
                        bFirstType = true;
                    }

                    for(j=0;j<modifyTypeList[n].size();j++)
                        cout << " " << n << "." << modifyTypeList[n][j]+1; 
                }

                cout << endl;
            }

            if(bPD==true)
                cout << "Problem formulation           : " << ((ForceForm==1)?"Primal":"Dual") << endl;

            if(bR==true)
                cout << "Rounding scheme               : " << ((bDefaultRoundingScheme==true)?"Default":"Manual") << endl;

            if(bA==true)
                cout << "Decimal place accuracy        : " << accuracyProjectSolution << endl;

            cout << "Forbid blowups                : " << (bForbidBlowups?"True":"False") << endl;
        }
    }
    catch(int)
    {
        cout << endl;
        cout << "ERROR : Unexpected argument." << endl;
        ShortUsage();
        return false;
    }
    
    return true;
}

//Tries to determine how arguments should be passed to the command interpreter via system().
//'cmd.exe' on windows expects the arguments to be enclosed by an extra pair of quotation marks,
//this is incompatible with 'sh' on linux.
//To determine which format works the program tries to call itself through system() with the argument "Test".
//If the call is successful the file FilenameOut will be modified.
//Returns true if the program should continue and construct the ".dat-s" file.
//Returns false if the program should terminate.
bool FindSystemFormat(const char* argv0)
{
    long   ret;
    string data;
    string toTest;
    
    cout << endl << endl;
    cout << "Determine the command interperter's argument format..." << endl;

    if(system(NULL)==0)
    {
        RunCommand = "";
        cout << "No command interpreter is available." << endl;
        cout << "The SDP solver will not be run." << endl;
        cout << endl << endl;
        throw 0;
        return true;
    }

    //Create an "<Empty>" .dat-s file.
    if(SaveFormatTest(FilenameSDP.c_str(),"<Empty>")==false)
    {
        ShortUsage();
        return false;
    }

    //Check it was created successfully. 
    if(ReadFormatTest(FilenameSDP.c_str(),data)==false)
        return false;

    if(data!="<Empty>")
    {
        cout << endl;
        cout << "ERROR : Retrieving data from the file \"" << FilenameSDP.c_str() << "\" failed." << endl << endl;
        return false;
    }

    //Try the 'sh' style format argument format.
    {
        //"<executable>" "Test" "<FilenameSDP>"
        toTest = string("\"")+argv0+"\" \"Test\" \""+FilenameSDP+"\"";

        //"<SDPCommand>" "<FilenameSDP>" "<FilenameOut>"
        RunCommand = "\""+SDPCommand+"\" \""+FilenameSDP+"\" \""+FilenameOut+"\"";

        //Test the command interperter.
        ret = system(toTest.c_str());

        //Check if the file was modified. 
        if(ReadFormatTest(FilenameSDP.c_str(),data)==false)
            return false;

        if(data!="<Modify>" && data!="<Empty>")
        {
            cout << endl;
            cout << "ERROR : Retrieving data from the file \"" << FilenameSDP.c_str() << "\" failed." << endl << endl;
            return false;
        }

        if(data=="<Modify>")
        {
            cout << "Format 1 succeeded." << endl;
            cout << endl << endl;
            return true;
        }
    }
    
    //'sh' style format failed.
    //Try the 'cmd.exe' style argument format.
    {
        cout << "Format 1 failed. Trying format 2..." << endl;

        //""<executable>" "Test" "<FilenameSDP>""
        toTest = "\""+toTest+"\"";

        //""<SDPCommand>" "<FilenameSDP>" "<FilenameOut>""
        RunCommand = "\""+RunCommand+"\""; //cmd.exe style format.
        
        //Test the command interperter.
        ret = system(toTest.c_str());

        //Check if the file was modified. 
        if(ReadFormatTest(FilenameSDP.c_str(),data)==false)
            return false;

        if(data!="<Modify>" && data!="<Empty>")
        {
            cout << endl;
            cout << "ERROR : Retrieving data from the file \"" << FilenameSDP.c_str() << "\" failed." << endl << endl;
            return false;
        }

        if(data=="<Modify>")
        {
            cout << "Format 2 succeeded." << endl;
            cout << endl << endl;
            return true;
        }
    }

    RunCommand = "";
    cout << "Format 2 failed." << endl;
    cout << endl;
    cout << "Could not get the command interpreter to run successfully." << endl;
    cout << "The SDP solver will not be run." << endl;
    cout << endl << endl;

    return true;
}

bool ChooseSDPFormat()
{
    long j, n, t;
    long nVar;

    if(ForceForm==0)
    {
        //Determine which formulation is more efficient.

        //Calculate the number of constraint matrices in the dual formulation.
        nVar = 1;
        for(n=0;n<=orderOfH-2;n++)
        for(t=0;t<FInfo[n].size();t++)
        for(j=0;j<2;j++)
            if(FInfo[n][t].Basis[j].size()>0)
                nVar += FInfo[n][t].Basis[j].size()*(FInfo[n][t].Basis[j].size()+1)/2;
        
        if(H.size()<nVar)
            ForceForm = 1;  //Primal.
        else
            ForceForm = -1; //Dual.
    }

    if(ForceForm==-1)
        return CreateDualSDPFile();

    return CreatePrimalSDPFile();
}

//Run the SDPCommand on FilenameSDP to create FilenameOut.
void RunSDPSolver()
{
    long ret;

    cout << "Running the SDP solver command : ";
    
    if(RunCommand=="")
        cout << "<None>" << endl << endl << endl;
    else
    {
        cout << RunCommand << endl << endl << endl;
        ret = system(RunCommand.c_str());
        cout << endl;
    }

    cout << endl;

    cout << "Control returned to ExactDensityBounder." << endl;
    cout << endl;

    return;
}

template<typename T>
vector<vector<T> > Inverse(vector<vector<T> > Matrix)
{
    long i, j, n;
    long dim;
    long maxIndex;
    T    maxSoFar;
    T    temp;
    vector<vector<T> > Inv; //Will hold the inverse of Matrix.

    dim = Matrix.size();

    //Allocate Inv matrix.
    Inv.resize(dim);
    for(i=0;i<dim;i++)
        Inv[i].resize(dim);

    //Set Inv to the identity.
    for(i=0;i<dim;i++)
    for(j=0;j<dim;j++)
    {
        if(i==j)
            Inv[i][j] = 1;
        else
            Inv[i][j] = 0;
    }

    //Do row operations to Matrix and Inv.
    for(n=0;n<dim;n++)
    {
        //Deal with row n.

        //Find the largest magnitude element in column n.
        {
            //Initially set the max element to occur in row n.
            maxIndex = n;

            if(Matrix[maxIndex][n]>=0)
                maxSoFar = Matrix[maxIndex][n];
            else
                maxSoFar = -Matrix[maxIndex][n];

            for(i=n+1;i<dim;i++)
            {
                if(maxSoFar<Matrix[i][n])
                {
                    maxSoFar = Matrix[i][n];
                    maxIndex = i;
                    continue;
                }

                if(maxSoFar<-Matrix[i][n])
                {
                    maxSoFar = -Matrix[i][n];
                    maxIndex = i;
                    continue;
                }
            }
        }

        //swap row n and row maxIndex.
        if(maxIndex!=n)
        {
            for(i=0;i<dim;i++)
            {
                temp                = Matrix[n][i];
                Matrix[n][i]        = Matrix[maxIndex][i];
                Matrix[maxIndex][i] = temp;
            }

            //Repeat row operations for Inv.
            for(i=0;i<dim;i++)
            {
                temp             = Inv[n][i];
                Inv[n][i]        = Inv[maxIndex][i];
                Inv[maxIndex][i] = temp;
            }
        }

        //Divide row n by Matrix[n][n].
        {
            temp = Matrix[n][n];

            for(i=0;i<dim;i++)
                Matrix[n][i] = Matrix[n][i]/temp;

            //Repeat row operations for Inv.
            for(i=0;i<dim;i++)
                Inv[n][i] = Inv[n][i]/temp;
        }

        //Subtract row n from other rows.
        for(i=0;i<dim;i++)
        {
            if(i==n)
                continue;

            temp = Matrix[i][n];

            //Subtract temp * row n from row i
            for(j=0;j<dim;j++)
                Matrix[i][j] = Matrix[i][j]-temp*Matrix[n][j];

            //Repeat row operations for Inv.
            for(j=0;j<dim;j++)
                Inv[i][j] = Inv[i][j]-temp*Inv[n][j];
        }
    }

    //Matrix is now the identity, and Inv holds the inverse of the original Matrix.
    return Inv;
}

//Prob represents the linear constraint problem (Prob[variable number][constraint number]), i.e.
//Prob[0][0]*x[0]+Prob[1][0]*x[1]+Prob[2][0]*x[2]+... = Prob[Prob.size()-1][0]
//Prob[0][1]*x[0]+Prob[1][1]*x[1]+Prob[2][1]*x[2]+... = Prob[Prob.size()-1][1]
//....
//The solutions are stored in Soln[variable number][kernel element], i.e.
//x[0] = Soln[0][Soln[0].size()-1]+c_0*Sol[0][0]+c_1*Sol[0][1]+c_2*Sol[0][2]+...
//x[1] = Soln[1][Soln[1].size()-1]+c_0*Sol[1][0]+c_1*Sol[1][1]+c_2*Sol[1][2]+...
//...
//where c_0, c_1, c_2, ... are real numbers.
//If no solutions exist Soln[0].size() will be 0.
//tolerance is used to decide if entries are zero (specifically they are treated as zero
//if their magnitude is <= tolerance).
template<typename T>
void LinearConstraintSolver(vector<vector<T> > Prob, vector<vector<T> >& Soln, T tolerance)
{
    long i, j, k;
    long nVar, nConstraints;
    long maxV, maxC;
    T    maxSoFar;
    T    temp;
    vector<long> leadingVariable, leadingConstraint;

    vector<vector<T> >().swap(Soln);

    //The number of variables we are trying to find.
    nVar = Prob.size()-1;

    if(nVar<=0)
    {
        //No variables to find.
        vector<vector<T> >().swap(Soln);
        return;
    }

    Soln.resize(nVar);

    //The number of constraints.
    nConstraints = Prob[0].size();

    if(nConstraints<=0)
    {
        //There are no constraints.
        for(i=0;i<nVar;i++)
        {
            Soln[i].resize(nVar+1,0);
            Soln[i][i] = 1;
        }

        return;
    }

    //For each constraint we need to associate a unique "leading" variable.
    
    //To begin with we will set it to -1 to represent no variable.
    leadingVariable.resize(nConstraints,-1);

    //Each variable is currently the leading variable for none of the constraints.
    leadingConstraint.resize(nVar,-1);

    for(k=0;k<nConstraints;k++)
    {
        //Find the largest element that is not part of a constraint which has a leading variable.
        //Make this element a new leading variable.
        maxSoFar =  0;
        maxV     = -1;
        maxC     = -1;
        for(i=0;i<nVar;i++)
        {
            //Check if the variable is already a leading variable for a constraint.
            if(leadingConstraint[i]!=-1)
                continue;
            
            if(maxSoFar<Prob[i][k])
            {
                maxSoFar = Prob[i][k];
                maxV     = i;
                maxC     = k;
                continue;
            }
            
            if(maxSoFar<-Prob[i][k])
            {
                maxSoFar = -Prob[i][k];
                maxV     = i;
                maxC     = k;
                continue;
            }
        }

        if(maxV==-1 || maxC==-1 || maxSoFar<=tolerance)
        {
            //For the system of constraints to be consistent, the constant must be zero.
            if(Prob[nVar][k]<-tolerance || tolerance<Prob[nVar][k])
                return; //Not consistent.

            continue;
        }

        //Set maxV as the leading variable for constraint maxC
        leadingVariable[maxC]   = maxV;
        leadingConstraint[maxV] = maxC;

        //Normalize the constraint.
        temp = Prob[maxV][maxC];
        for(i=0;i<=nVar;i++)
            Prob[i][maxC] = Prob[i][maxC]/temp;

        //Subtract the constraint maxC from other constraints.
        for(i=0;i<nConstraints;i++)
        {
            if(i==maxC)
                continue;

            temp = Prob[maxV][i];
            if(-tolerance<=temp && temp<=tolerance)
                continue;

            for(j=0;j<=nVar;j++)
                Prob[j][i] = Prob[j][i]-temp*Prob[j][maxC];
        }
    }
    
    //Find the elements of the kernel.
    for(i=0;i<nVar;i++)
    {
        if(leadingConstraint[i]!=-1)
            continue;

        for(j=0;j<nVar;j++)
        {
            if(j==i)
                Soln[j].push_back(-1);
            else
            if(leadingConstraint[j]==-1)
                Soln[j].push_back(0);
            else
                Soln[j].push_back(Prob[i][leadingConstraint[j]]);

        }
    }

    //Find a specific solution.
    for(i=0;i<nVar;i++)
    {
        if(leadingConstraint[i]==-1)
            Soln[i].push_back(0);
        else
            Soln[i].push_back(Prob[nVar][leadingConstraint[i]]);
    }

    return;
}

template<typename T>
class SparseComponent
{
public:
    T    coeff;
    long index;
};

//Zero entries are (mostly) omitted from Soln to save space.
//Soln[0] to Soln[Soln.size()-2] describe the kernel vectors.
//Soln[Soln.size()-1] describes a solution.
template<typename T>
void SparseLinearConstraintSolver(vector<vector<T> >& Prob, vector<vector<SparseComponent<T> > >& Soln, T tolerance)
{
    long i, j;
    long nVar, nConstraints;
    long maxV, maxC;
    T    maxSoFar;
    T    temp;
    vector<long> leadingVariable, leadingConstraint;

    vector<vector<SparseComponent<T> > >().swap(Soln);

    //The number of variables we are trying to find.
    nVar = Prob.size()-1;

    if(nVar<=0)
    {
        //No variables to find.
        vector<vector<SparseComponent<T> > >().swap(Soln);
        return;
    }
    
    //The number of constraints.
    nConstraints = Prob[0].size();

    if(nConstraints<=0)
    {
        //There are no constraints.
        Soln.resize(nVar+1);
        for(i=0;i<nVar;i++)
        {
            Soln[i].resize(1);
            Soln[i][0].index = i;
            Soln[i][0].coeff = 1;
        }

        return;
    }

    //For each constraint we need to associate a unique "leading" variable.
    
    //To begin with we will set it to -1 to represent no variable.
    leadingVariable.resize(nConstraints,-1);

    //Each variable is currently the leading variable for none of the constraints.
    leadingConstraint.resize(nVar,-1);

    while(true)
    {
        //Find the largest element that is not part of a constraint which has a leading variable.
        //Make this element a new leading variable.
        maxSoFar =  0;
        maxV     = -1;
        maxC     = -1;
        for(i=0;i<nVar;i++)
        {
            //Check if the variable is already a leading variable for a constraint.
            if(leadingConstraint[i]!=-1)
                continue;

            for(j=0;j<nConstraints;j++)
            {
                //Check if the constraint already has a leading variable.
                if(leadingVariable[j]!=-1)
                    continue;

                if(maxSoFar<Prob[i][j])
                {
                    maxSoFar = Prob[i][j];
                    maxV     = i;
                    maxC     = j;
                    continue;
                }

                if(maxSoFar<-Prob[i][j])
                {
                    maxSoFar = -Prob[i][j];
                    maxV     = i;
                    maxC     = j;
                    continue;
                }
            }
        }

        if(maxV==-1 || maxC==-1 || maxSoFar<=tolerance)
            break; //No more non-zero elements.

        //Set maxV as the leading variable for constraint maxC
        leadingVariable[maxC]   = maxV;
        leadingConstraint[maxV] = maxC;

        //Normalize the constraint.
        temp = Prob[maxV][maxC];
        for(i=0;i<=nVar;i++)
            Prob[i][maxC] = Prob[i][maxC]/temp;

        //Subtract the constraint maxC from other constraints.
        for(i=0;i<nConstraints;i++)
        {
            if(i==maxC)
                continue;

            temp = Prob[maxV][i];
            for(j=0;j<=nVar;j++)
                Prob[j][i] = Prob[j][i]-temp*Prob[j][maxC];
        }
    }

    //Check if there are any constraints which have no leading variable.
    for(i=0;i<nConstraints;i++)
    {
        if(leadingVariable[i]==-1)
        {
            //For the system of constraints to be consistent, the constant must be zero.
            if(Prob[nVar][i]<-tolerance || tolerance<Prob[nVar][i])
                return; //Not consistent.
        }
    }

    //Find the elements of the kernel.
    for(i=0;i<nVar;i++)
    {
        if(leadingConstraint[i]!=-1)
            continue;

        Soln.push_back(vector<SparseComponent<T> >());

        for(j=0;j<nVar;j++)
        {
            if(j==i)
            {
                Soln.back().push_back(SparseComponent<T>());
                Soln.back().back().index = j;
                Soln.back().back().coeff = -1;
            }
            else
            if(leadingConstraint[j]==-1)
                continue;
            else
            {
                Soln.back().push_back(SparseComponent<T>());
                Soln.back().back().index = j;
                Soln.back().back().coeff = Prob[i][leadingConstraint[j]];
            }
        }
    }

    //Find a specific solution.
    Soln.push_back(vector<SparseComponent<T> >());

    for(i=0;i<nVar;i++)
    {
        if(leadingConstraint[i]==-1)
            continue;
        else
        {
            Soln.back().push_back(SparseComponent<T>());
            Soln.back().back().index = i;
            Soln.back().back().coeff = Prob[nVar][leadingConstraint[i]];
        }
    }

    return;
}

double ConvertFracToDouble(CLargeFraction F)
{
    long      i;
    double    N, D;
    bool      bPositive;
    long      eN, eD;
    CLargeInt temp;

    if(F.N>=0)
        bPositive = true;
    else
    {
        bPositive = false;
        F.N = -F.N;
    }

    N  = 0;
    eN = 0;
    while(F.N!=0)
    {
        temp = F.N%256;
        F.N  = F.N/256;

        N += temp.ConvertToLong();
        N /= 256;
        eN++;
    }

    if(bPositive==false)
        N = -N;

    if(F.D>=0)
        bPositive = true;
    else
    {
        bPositive = false;
        F.D = -F.D;
    }

    D  = 0;
    eD = 0;
    while(F.D!=0)
    {
        temp = F.D%256;
        F.D  = F.D/256;

        D += temp.ConvertToLong();
        D /= 256;
        eD++;
    }

    if(bPositive==false)
        D = -D;

    //Floating value is N/D*256^(eN-eD).
    N  = N/D;
    eN = eN-eD;

    if(eN>0)
    {
        for(i=0;i<eN;i++)
            N *= 256;
    }
    else
    {
        eN = -eN;
        for(i=0;i<eN;i++)
            N /= 256;
    }

    return N;
}

CLargeFraction ConvertDoubleToFrac(double d, const CLargeInt& Denominator)
{
    long n;
    long temp;
    bool bPositive;
    CLargeFraction Result;
    CLargeInt      Pow256;

    d *= ConvertFracToDouble(Denominator);

    if(d<0)
    {
        d = -d;
        bPositive = false;
    }
    else
        bPositive = true;

    Result.D = Denominator;

    Result.N = 0;

    n      = 0;
    Pow256 = 1;
    while(d>=1)
    {
        d /= 256;

        n++;
        Pow256 = Pow256*256;
    }

    while(true)
    {
        temp = (int)d;
        Result.N = Result.N+temp*Pow256;

        d -= temp;
        d *= 256;
        Pow256 = Pow256/256;
        n--;

        if(n==-1)
            break;
    }

    if(bPositive==false)
        Result.N = -Result.N;

    return Result;
}

class CTerm
{
public:
    vector<vector<double> >         dMatrix;
    vector<vector<CLargeFraction> > fMatrix;

    //Flag information.
    long nType;
    long iType;

    vector<vector<long> >           BasisFlag;
    vector<vector<CLargeFraction> > BasisCoeff;
};

vector<CTerm>  Soln;
vector<long>   sharpIneq;
CLargeFraction TargetTuranDensity;

CLargeInt GreatestCommonDivisor(CLargeInt A, CLargeInt B)
{
    CLargeInt& R1 = A;
    CLargeInt& R2 = B;
    CLargeInt  R3;

    if(R1<0)
        R1 = -R1;

    if(R2<0)
        R2 = -R2;

    if(R1==0)
        return R2;

    if(R2==0)
        return R1;
  
    //Apply the Euclidean algorithm to determine the highest common factor.
    while(true)
    {
        R3 = R1%R2;

        if(R3==0)
            break; //R2 is the highest common factor.

        R1 = R2;
        R2 = R3;
    }

    return R2;
}

CLargeInt LowestCommonMultiple(CLargeInt A, CLargeInt B)
{
    CLargeInt& R1 = A;
    CLargeInt& R2 = B;
    CLargeInt  R3;
    CLargeInt  Product;

    if(R1<0)
        R1 = -R1;

    if(R2<0)
        R2 = -R2;

    if(R1==0)
        return R2;

    if(R2==0)
        return R1;
  
    return R1*R2/GreatestCommonDivisor(R1,R2);
}

CLargeFraction NearestRational(double toMakeRational, long maxDenominator)
{
    long   i;
    double diff;
    double bestDiff;
    CLargeFraction F, Best;

    Best     = 0;
    bestDiff = toMakeRational - ConvertFracToDouble(Best);
    if(bestDiff<0)
        bestDiff = -bestDiff;

    for(i=1;i<=maxDenominator;i++)
    {
        F    = ConvertDoubleToFrac(toMakeRational,i);
        diff = toMakeRational-ConvertFracToDouble(F);
        
        if(diff<0)
            diff = -diff;

        if(diff<bestDiff)
        {
            bestDiff = diff;
            Best = F;
        }

        F.N  = F.N-1;
        diff = toMakeRational-ConvertFracToDouble(F);
        
        if(diff<0)
            diff = -diff;

        if(diff<bestDiff)
        {
            bestDiff = diff;
            Best = F;
        }

        F.N  = F.N+2;
        diff = toMakeRational-ConvertFracToDouble(F);
        
        if(diff<0)
            diff = -diff;

        if(diff<bestDiff)
        {
            bestDiff = diff;
            Best = F;
        }
    }

    Best.Simplify();
    return Best;
}

bool CalculateSharpIneq()
{
    double  tolerance = toleranceSharpIneq;
    long    i, j, r, c, h, n;
    long    dim;
    double  product;
    long    maxIneq;
    long*   pVal;
    KeyProd key;
    vector<double> Ineq;
    vector<long>   sortIndex; 
    ofstream       File;
    string         filename;

    cout << "Finding the sharp inequalities..." << endl; 

    Ineq.resize(H.size());

    for(i=0;i<H.size();i++)
        Ineq[i] = double(NumberOfEdges(H[i]))*6/(orderOfH*(orderOfH-1)*(orderOfH-2));

    //Calculate each term's contribution.
    for(n=0;n<Soln.size();n++)
    {
        dim = Soln[n].dMatrix.size();

        //Calculate the term's contribution to Ineq[h].
        for(h=0;h<H.size();h++)
        {
            key.nType  = Soln[n].nType;
            key.iType  = Soln[n].iType;

            //Calculate the entry dMatrix[r][c]'s contribution.
            for(r=0;r<dim;r++)
            for(c=r;c<dim;c++)
            {
                product = 0;

                //Multiply the flags in the basis.
                for(i=0;i<Soln[n].BasisFlag[r].size();i++)
                for(j=0;j<Soln[n].BasisFlag[c].size();j++)
                {
                    //Look up the product of the flags.

                    //Set the key.
                    if(Soln[n].BasisFlag[r][i]<=Soln[n].BasisFlag[c][j])
                    {
                        key.iFlag1 = Soln[n].BasisFlag[r][i];
                        key.iFlag2 = Soln[n].BasisFlag[c][j];
                    }
                    else
                    {
                        key.iFlag1 = Soln[n].BasisFlag[c][j];
                        key.iFlag2 = Soln[n].BasisFlag[r][i];
                    }
                                        
                    pVal = htProd[h].Get(key);
                                        
                    if(pVal!=NULL)
                        product += ConvertFracToDouble(Soln[n].BasisCoeff[r][i]*Soln[n].BasisCoeff[c][j])
                                   *(*pVal)/Factorial[orderOfH];
                }

                Ineq[h] += product*Soln[n].dMatrix[r][c]*((r==c)?1:2);
            }
        }
    }

    //Sort the inequalities
    sortIndex.resize(H.size());
    for(h=0;h<H.size();h++)
        sortIndex[h] = h;

    for(h=0;h<H.size();h++)
    {
        maxIneq = h;
        for(i=h;i<H.size();i++)
            if(Ineq[sortIndex[i]]>Ineq[sortIndex[maxIneq]])
                maxIneq = i;
        
        //Swap h and maxIneq
        i                  = sortIndex[h];
        sortIndex[h]       = sortIndex[maxIneq];
        sortIndex[maxIneq] = i;
    }

    TargetTuranDensity = NearestRational(Ineq[sortIndex[0]],accuracyTargetDensity);

    if(bDefaultRoundingScheme==false)
    {
        cout << setprecision(20);
        cout << "Density     : " << Ineq[sortIndex[0]] << " ("<< TargetTuranDensity << ")" << endl;
        cout << "Numerator   : ";
        cin >> i;
        cout << "Denominator : ";
        cin >> j;

        TargetTuranDensity.N = i;
        TargetTuranDensity.D = j;
        TargetTuranDensity.Simplify();

        cout << endl;
    }
    else
    {
        cout << "Target density : " << TargetTuranDensity << endl;
    }

    double dTarget;
    dTarget = ConvertFracToDouble(TargetTuranDensity);
    for(n=0;n<H.size();n++)
    {
        if(Ineq[sortIndex[n]]<dTarget-tolerance)
            break;
    }
    
    filename = FilenameOut;
    filename.resize(filename.size()-4);
    filename += "_ineq.txt";
    File.open(filename.c_str(),ios::trunc);

    if(File.good()==false)
    {
        cout << endl;
        cout << "ERROR : Could not open \"" << filename << "\" for writing." << endl;
        return false;
    }

    File << "(" << n << ")" << endl;

    File << scientific << setprecision(15);
    for(h=0;h<H.size();h++)
    {
        File << setw(4) << h+1 << " : " << showpos << Ineq[sortIndex[h]] << noshowpos;// << endl;
        File << " : ";
        DisplayEdges(File,H[sortIndex[h]]);
        File << endl;
    }

    File.close();

    if(bDefaultRoundingScheme==false)
    {
        cout << "Number of sharp inequalities (" << n << ") : ";
        cin >> n;
    }

    vector<long>().swap(sharpIneq);
    sharpIneq.resize(n);

    for(i=0;i<n;i++)
        sharpIneq[i] = sortIndex[i];

    //for(i=0;i<sharpIneq.size();i++)
    //    cout << sharpIneq[i] << " " << Ineq[sharpIneq[i]] << endl; 

    return true;
}

//Scale each Basis element so that all coefficients are integers.
//Also ensure the first non-zero coefficient in each basis element is positive.
void ScaleBasis()
{
    long      i, j, n;
    CLargeInt LCM, GCD;
    long      firstSign;
    
    for(n=0;n<Soln.size();n++)
    {
        for(i=0;i<Soln[n].BasisCoeff.size();i++)
        {
            firstSign = 0;
            LCM       = 1;
            for(j=0;j<Soln[n].BasisCoeff[i].size();j++)
            {
                Soln[n].BasisCoeff[i][j].Simplify();
                LCM = LowestCommonMultiple(LCM,Soln[n].BasisCoeff[i][j].D);

                if(firstSign==0 && Soln[n].BasisCoeff[i][j]!=0)
                {
                    if(Soln[n].BasisCoeff[i][j]>0)
                        firstSign = 1;
                    else
                        firstSign = -1;
                }
            }

            LCM = firstSign*LCM;

            //Scale the basis.
            for(j=0;j<Soln[n].BasisCoeff[i].size();j++)
                Soln[n].BasisCoeff[i][j] = Soln[n].BasisCoeff[i][j]*LCM;

            //Scale the fMatrix.
            for(j=0;j<Soln[n].fMatrix.size();j++)
            {
                Soln[n].fMatrix[i][j] = Soln[n].fMatrix[i][j]/LCM;
                Soln[n].fMatrix[j][i] = Soln[n].fMatrix[j][i]/LCM;
            }
        }
    }

    //Repeat this time dividing by the greatest common divisor.
    for(n=0;n<Soln.size();n++)
    {
        for(i=0;i<Soln[n].BasisCoeff.size();i++)
        {
            GCD = 0;
            for(j=0;j<Soln[n].BasisCoeff[i].size();j++)
                GCD = GreatestCommonDivisor(GCD,Soln[n].BasisCoeff[i][j].N);
            
            //Scale the basis.
            for(j=0;j<Soln[n].BasisCoeff[i].size();j++)
                Soln[n].BasisCoeff[i][j].N = Soln[n].BasisCoeff[i][j].N/GCD;

            //Scale the fMatrix.
            for(j=0;j<Soln[n].fMatrix.size();j++)
            {
                Soln[n].fMatrix[i][j] = Soln[n].fMatrix[i][j]*GCD;
                Soln[n].fMatrix[j][i] = Soln[n].fMatrix[j][i]*GCD;
            }
        }
    }

    return;
}

//Used by CalculateDensityBound().
CLargeFraction CalculateCoeff(long h)
{
    long    i, j, n, r, c;
    long    dim;
    long*   pVal;
    KeyProd key;
    CLargeInt      product;
    CLargeFraction Coeff;
        
    Coeff = CLargeFraction(NumberOfEdges(H[h]))*6/(orderOfH*(orderOfH-1)*(orderOfH-2));
    
    //Calculate each term's contribution.
    for(n=0;n<Soln.size();n++)
    {
        dim = Soln[n].fMatrix.size();
        
        key.nType  = Soln[n].nType;
        key.iType  = Soln[n].iType;

        //Calculate the entry dMatrix[r][c]'s contribution.
        for(r=0;r<dim;r++)
        for(c=r;c<dim;c++)
        {
            product = 0;

            //Multiply the flags in the basis.
            for(i=0;i<Soln[n].BasisFlag[r].size();i++)
            for(j=0;j<Soln[n].BasisFlag[c].size();j++)
            {
                //Look up the product of the flags.

                //Set the key.
                if(Soln[n].BasisFlag[r][i]<=Soln[n].BasisFlag[c][j])
                {
                    key.iFlag1 = Soln[n].BasisFlag[r][i];
                    key.iFlag2 = Soln[n].BasisFlag[c][j];
                }
                else
                {
                    key.iFlag1 = Soln[n].BasisFlag[c][j];
                    key.iFlag2 = Soln[n].BasisFlag[r][i];
                }
                                        
                pVal = htProd[h].Get(key);
                                        
                if(pVal!=NULL)
                    product = product+Soln[n].BasisCoeff[r][i].N*Soln[n].BasisCoeff[c][j].N*(*pVal);
            }

            Coeff = Coeff+product*Soln[n].fMatrix[r][c]*((r==c)?1:2)/Factorial[orderOfH];
        }
    }

    return Coeff;
}

//This algorithm calls ScaleBasis(). 
bool CalculateDensityBound(CLargeFraction& Density)
{
    long i, h;
    long noOfH = H.size();
    long globalError;
    long loopError;
    long doneSoFar;
    vector<CLargeFraction> Ineq;

    cout << "Calculating the upper bound..." << endl;

    ScaleBasis();

    //Calculate the coefficient of each H.
    Ineq.resize(noOfH);

    doneSoFar   = 0;
    globalError = 0; //No error.
    #pragma omp parallel for schedule(dynamic) default(shared) private(h,loopError)
    for(h=0;h<noOfH;h++)
    {
        try
        {
            #pragma omp critical (CalculateDensityBound_Progress)
            {
                loopError = globalError;
            }

            if(loopError==0)
                Ineq[h] = CalculateCoeff(h);
        }
        catch(ErrorDivisionByZero&)
        {
            loopError = 4;
        }
        catch(ErrorOverflow&)
        {
            loopError = 3;
        }
        catch(bad_alloc&)
        {
            loopError = 2;
        }
        catch(...)
        {
            loopError = 1;
        }
        
        #pragma omp critical (CalculateDensityBound_Progress)
        {
            if(loopError>globalError)
                globalError = loopError;

            if(globalError==0)
            {
                doneSoFar++;
                cout << "\r" << doneSoFar << " of " << noOfH << " coefficients calculated." << flush; 
            }
        }
    }

    //Display any errors.
    if(globalError!=0)
    {
        cout << endl << endl;
        
        if(globalError==4)
            cout << "ERROR : Division by zero.";

        if(globalError==3)
            cout << "ERROR : Numeric overflow occurred. Try increasing LI_MAX in the source code.";

        if(globalError==2)
            cout << "ERROR : Allocating memory failed.";

        if(globalError==1)
            cout << "ERROR : An unknown error occurred.";

        cout << endl << endl;
        return false;
    }

    cout << "\r                             ";
    for(i=0;i<2*NumberOfDigits(noOfH);i++)
        cout << " ";

    cout << "\rDone."<< endl;

    //Find the largest coefficient.
    Density = 0;
    for(h=0;h<noOfH;h++)
        if(Density<Ineq[h])
            Density = Ineq[h];

    return true;
}

//Matrix[dim][dim].
//vecSystem[system number][vector no][element of vector index].
//maxVal[constraint row].
void FindZeroEigenData(vector<vector<double> > Matrix, vector<vector<vector<double> > >& vecSystem, vector<double>& maxVal)
{
    long   i, j, k;
    long   dim;
    long   maxV, maxC;
    double maxSoFar;
    double tolerance = toleranceEigenConstraint;
    double temp;
    vector<long> leadingVariable, leadingConstraint;
    
    //Works in much the same way as LinearConstraintSolver(), but the number of zero rows is varied.
    //The constants of the constraints are all zero and so have been omitted.

    vector<double>().swap(maxVal);
    vector<vector<vector<double> > >().swap(vecSystem);

    dim = Matrix.size();

    if(dim<=0)
    {
        //No vectors to find.
        vecSystem.push_back(vector<vector<double> >());

        return;
    }

    //For each constraint record its maximum value.
    maxVal.resize(dim,0);

    //Add the system when all rows are considered zero.
    vecSystem.push_back(vector<vector<double> >());

    vecSystem.back().resize(dim);
    for(i=0;i<dim;i++)
    {
        vecSystem.back()[i].resize(dim,0);
        vecSystem.back()[i][i] = 1;
    }

    //For each constraint we need to associate a unique "leading" variable.
    
    //To begin with we will set it to -1 to represent no variable.
    leadingVariable.resize(dim,-1);

    //Each variable is currently the leading variable for none of the constraints.
    leadingConstraint.resize(dim,-1);

    for(k=0;k<dim;k++)
    {
        //Find the largest element that is not part of a constraint which has a leading variable.
        //Make this element a new leading variable.
        maxSoFar =  0;
        maxV     = -1;
        maxC     = -1;
        for(i=0;i<dim;i++)
        {
            //Check if the constraint already has a leading variable.
            if(leadingVariable[i]!=-1)
                continue;

            for(j=0;j<dim;j++)
            {
                //Check if the variable is already a leading variable for a constraint.
                if(leadingConstraint[j]!=-1)
                    continue;
            
                if(maxSoFar<Matrix[i][j])
                {
                    maxSoFar = Matrix[i][j];
                    maxV     = j;
                    maxC     = i;
                    continue;
                }
            
                if(maxSoFar<-Matrix[i][j])
                {
                    maxSoFar = -Matrix[i][j];
                    maxV     = j;
                    maxC     = i;
                    continue;
                }
            }
        }

        if(maxV==-1 || maxC==-1 || maxSoFar<=tolerance)
        {
            //No more variables were found or they are all below tolerance.
            //Set the rest of maxVal to maxSoFar.

            for(i=k;i<dim;i++)
                maxVal[i] = maxSoFar;

            return;
        }

        maxVal[k] = maxSoFar;

        //Set maxV as the leading variable for constraint maxC
        leadingVariable[maxC]   = maxV;
        leadingConstraint[maxV] = maxC;

        //Normalize the constraint.
        temp = Matrix[maxC][maxV];
        for(i=0;i<dim;i++)
            Matrix[maxC][i] /= temp;

        //Subtract the constraint maxC from other constraints.
        for(i=0;i<dim;i++)
        {
            if(i==maxC)
                continue;

            temp = Matrix[i][maxV];
            for(j=0;j<dim;j++)
                Matrix[i][j] -= temp*Matrix[maxC][j];
        }

        //Find the elements of the kernel, i.e. the 0-eigenvectors.
        vecSystem.push_back(vector<vector<double> >());

        for(i=0;i<dim;i++)
        {
            if(leadingConstraint[i]!=-1)
                continue;

            //Determine the new kernel vector.
            vecSystem.back().push_back(vector<double>());

            for(j=0;j<dim;j++)
            {
                if(j==i)
                    vecSystem.back().back().push_back(-1);
                else
                if(leadingConstraint[j]==-1)
                    vecSystem.back().back().push_back(0);
                else
                    vecSystem.back().back().push_back(Matrix[leadingConstraint[j]][i]);

            }
        }
    }
    
    return;
}

//*** Not used ***
//Used to separate the 0-eigenvectors into vectors which have approximately rational components.
//The vectors should appear as columns in vectors[][].
void SimplifyVectors(vector<vector<double> >& vectors)
{
    long   i, j, n;
    long   nVec, dim;
    long   indexOfMax;
    double maxSoFar, temp;

    //Store the number of components in a vector.
    dim = vectors.size();

    //Store the number of vectors.
    nVec = vectors[0].size();

    for(n=0;n<nVec;n++)
    {
        //Find the largest element in vector n.
        indexOfMax = 0;
        maxSoFar   = vectors[indexOfMax][n];
        if(maxSoFar<0)
            maxSoFar = -maxSoFar;

        for(i=0;i<dim;i++)
        {
            if(maxSoFar<vectors[i][n])
            {
                maxSoFar   = vectors[i][n];
                indexOfMax = i;
            }

            if(maxSoFar<-vectors[i][n])
            {
                maxSoFar   = -vectors[i][n];
                indexOfMax = i;
            }
        }

        //Divide through by the max element.
        temp = vectors[indexOfMax][n];
        for(i=0;i<dim;i++)
            vectors[i][n] /= temp;

        //Subtract from other vectors.
        for(i=0;i<nVec;i++)
        {
            if(i==n)
                continue;

            temp = vectors[indexOfMax][i];
            for(j=0;j<dim;j++)
                vectors[j][i] -= temp*vectors[j][n];
        }
    }
    
    return;
}

//This function is used in parsing an input file. It reads one character from the file. If
//the end of the file is reached it returns '\0'. If the character is '#' it interprets this
//as the start of a comment. It will then skip over all characters until '\n','\r', or '\0'
//is encountered and returns that instead.
char ReadChar(istream& in)
{
    char c;

    in.get(c);

    //Check if the end of the file has been reached.
    if(in.eof()==true)
        return '\0';

    //Check if a comment has been encountered.
    if(c=='#')
    {
        //Remove the comment.
        while(true)
        {
            in.get(c);

            if(in.eof()==true)
                return '\0';

            if(c=='\n' || c=='\r' || c=='\0')
                return c;
        }
    }

    return c;
}

//This function is used in parsing an input file. It discards whitespace characters until it
//reaches a non-whitespace character. The last character that was read (but not yet
//interpreted) should be passed to the function via 'c'. When the function returns 'c' will
//hold the last character read by the function.
void RemoveWhitespace(istream& in, char& c)
{
    while(c==' ' || c=='\t' || c=='\n' || c=='\r')
        c = ReadChar(in);

    return;
}

//This function is used in parsing an input file. It reads in a CLargeInt value and returns
//it via the variable 'L'. It may throw an ErrorOverflow object if the integer is too large.
//The first character must be +,-,0-9. If the first character is + or - the next character
//must be one of 0-9 or a series of whitespaces followed by 0-9. No whitespaces are allowed
//between the numeric digits. If the function reads in an integer of this form successfully
//it returns true, otherwise it returns false. The last character that was read (but not yet
//interpreted) should be passed to the function via 'c'. When the function returns 'c' will
//hold the last character read by the function.
bool ReadLargeInteger(istream& in, CLargeInt& L, char& c)
{
    bool bPositive;
    bool bReadDigit;

    L = 0;

    //Check if the first character is +,-,0-9.
    if(c=='+')
    {
        bReadDigit = false;
        bPositive  = true;
    }
    else
    if(c=='-')
    {
        bReadDigit = false;
        bPositive  = false;
    }
    else
    if('0'<=c && c<='9')
    {
        L = c-'0';
        bReadDigit = true;
        bPositive  = true;
    }
    else
        return false;

    //If the first character was + or -, remove any whitespaces and read in the first digit.
    if(bReadDigit==false)
    {
        c = ReadChar(in);
        RemoveWhitespace(in,c);

        if('0'<=c && c<='9')
            L = c-'0';
        else
            return false;
    }
    
    //Read the remaining digits.
    while(true)
    {
        c = ReadChar(in);

        if('0'<=c && c<='9')
            L = 10*L+c-'0'; //May throw ErrorOverflow object.
        else
            break; //Reached the end of the integer.
    }

    if(bPositive==false)
        L = -L;

    return true;
}

bool MakeMatricesPositiveDefinite()
{
    //double   thresholdMaxConstraint  = 1e-5;
    //double   thresholdMaxError       = 1e-5;
    //long     thresholdSimpleRational = 20;
    //long     accuracy                = 100; //Used to rationalize the vector elements.
    long accuracy = accuracySimpleRational; 
    long i, j, k, n, p, q;
    vector<long>                              noOfZeroVectors;
    vector<vector<double> >                   zeroVectors;
    vector<vector<vector<CLargeFraction> > >  ChangeOfBasis;
    vector<vector<double> >                   eigenvalues; //Holds the max constraint element not the eigenvalues.
    vector<vector<vector<double> > >          eigenvectors;
    vector<vector<vector<vector<double> > > > eSystem;

    noOfZeroVectors.resize(Soln.size());

    cout << "Computing the eigenvector constraints..." << endl;
    {
        ofstream File;
        string   filename;

        filename = FilenameOut;
        filename.resize(filename.size()-4);
        filename += "_constraints.txt";
        File.open(filename.c_str(),ios::trunc);

        if(File.good()==false)
        {
            cout << endl;
            cout << "ERROR : Could not open \"" << filename << "\" for writing." << endl;
            return false;
        }

        File << "Matrix [matrix number] : ([program's guess at the number of 0-eigenvectors])" << endl;
        File << "[assumed number of 0-eigenvectors] :  [Threshold for zero] : ";
        File << "[largest denominator in program's guess] : [max error in rationalized values]" << endl;
        File << endl << endl;

        eigenvalues.resize(Soln.size());
        eigenvectors.resize(Soln.size());
        eSystem.resize(Soln.size());

        File << setprecision(10) << scientific;
        for(n=0;n<Soln.size();n++)
        {
            FindZeroEigenData(Soln[n].dMatrix,eSystem[n],eigenvalues[n]);

            //Acumulate data to help us decide how many 0-eigenvectors there are. 
            double            error;
            vector<double>    maxError;
            vector<CLargeInt> D;
            CLargeFraction    F;

            maxError.resize(eigenvalues[n].size());
            D.resize(eigenvalues[n].size());

            for(i=0;i<eigenvalues[n].size();i++)
            {
                //Find the largest denominator in the eigensystem and the maximum error in the rationalization.
                //(This gives us a measure of how complicated the system is.)
                D[i]        = 0;
                maxError[i] = 0;
                if(i<eSystem[n].size())
                {
                    for(j=0;j<eSystem[n][i].size();j++)
                    for(k=0;k<eSystem[n][i][j].size();k++)
                    {
                        F = NearestRational(eSystem[n][i][j][k],accuracy);

                        if(D[i]<F.D)
                            D[i] = F.D;

                        error = eSystem[n][i][j][k] - ConvertFracToDouble(F);
                        if(error<0)
                            error = -error;

                        if(error>maxError[i])
                            maxError[i] = error;
                    }
                }
            }

            //Decide how many 0-eigenvectors there are.
            noOfZeroVectors[n] = 0;
            for(i=0;i<eigenvalues[n].size();i++)
            {
                p = eigenvalues[n].size()-1-i;

                if(eigenvalues[n][p]<=thresholdMaxConstraint
                && D[p]<=thresholdSimpleRational
                && D[p]!=0
                && maxError[p]<=thresholdMaxError)
                    noOfZeroVectors[n] = i+1;
            }

            File << "Matrix " << n+1 << " : ";
            File << "(" << noOfZeroVectors[n] << ")" << endl;

            for(i=0;i<eigenvalues[n].size();i++)
            {
                p = eigenvalues[n].size()-1-i;
                
                //File << setw(NumberOfDigits(eigenvalues[n].size())) << i+1 << " : ";
                File << setw(3) << i+1 << " : ";
                File << eigenvalues[n][p] << " : ";
                File << setw(3) << D[p] << " : ";
                File << maxError[p];
                File << endl;
            }

            File << endl;
        }

        File.close();
    }
        
    if(bDefaultRoundingScheme==false)
    {
        cout << endl;
        cout << "Number of 0-eigenvectors : " << endl;

        for(n=0;n<Soln.size();n++)
        {
            cout << "Matrix " << setw(NumberOfDigits(Soln.size())) << n+1 << " (" << noOfZeroVectors[n] << ") : ";
            cin >> noOfZeroVectors[n];
        }

        cout << endl;
    }

    cout << "Computing the eigenvectors..." << endl;
    {
        ofstream File;
        string   filename;

        filename = FilenameOut;
        filename.resize(filename.size()-4);
        filename += "_eigenvectors.txt";
        File.open(filename.c_str(),ios::trunc);

        if(File.good()==false)
        {
            cout << endl;
            cout << "ERROR : Could not open \"" << filename << "\" for writing." << endl;
            return false;
        }

        File << "Matrix [matrix number] :" << endl;
        File << "[0-eigenvector number] :  [floating point entry] ";
        File << "([guessed rational] [error in guess]), ..." << endl;
        File << endl << endl;
        File << scientific;

        ChangeOfBasis.resize(Soln.size());

        for(n=0;n<Soln.size();n++)
        {
            //Index into eSystem (may need to increase the number of 0-eigenvectors).
            p = eigenvalues[n].size()-noOfZeroVectors[n];
            if(eSystem[n].size()<=p)
            {
                p                  = eSystem[n].size()-1;
                noOfZeroVectors[n] = eigenvalues[n].size()-p;
            }
            
            zeroVectors.resize(Soln[n].dMatrix.size());
            for(i=0;i<zeroVectors.size();i++)
                zeroVectors[i].resize(noOfZeroVectors[n]);

            ChangeOfBasis[n].resize(Soln[n].dMatrix.size());
            for(i=0;i<Soln[n].dMatrix.size();i++)
                ChangeOfBasis[n][i].resize(Soln[n].dMatrix.size());
            
            for(k=0;k<noOfZeroVectors[n];k++)
            for(i=0;i<zeroVectors.size();i++)
                zeroVectors[i][k] = eSystem[n][p][k][i];

            File << "Matrix " << n+1 << " : " << endl;

            for(k=0;k<zeroVectors[0].size();k++)
            {
                //File << " " << setw(NumberOfDigits(zeroVectors[0].size())) << k+1 << "  :";
                File << setw(3) << k+1 << " :";

                for(i=0;i<zeroVectors.size();i++)
                {
                    double error;

                    ChangeOfBasis[n][i][k] = NearestRational(zeroVectors[i][k],accuracy);
                    error = zeroVectors[i][k] - ConvertFracToDouble(ChangeOfBasis[n][i][k]);
                    if(error<0)
                        error = -error;

                    File << "  " << showpos << setprecision(10) << zeroVectors[i][k] << noshowpos;
                    File << " (" << setw(7) << ChangeOfBasis[n][i][k];
                    File << " " << setprecision(2) << error << ")";

                    if(i!=zeroVectors.size()-1)
                        File << ",";
                }

                File << endl;
            }

            File << endl;
        }

        //Write out template input file for user defined rounding.
        File << "\n\nTemplate input file :" << endl << endl;
        for(n=0;n<Soln.size();n++)
        {
            File << "#Matrix " << n+1 << " : " << endl;

            for(k=0;k<noOfZeroVectors[n];k++)
            {
                for(i=0;i<Soln[n].dMatrix.size();i++)
                    File << setw(7) << ChangeOfBasis[n][i][k] << " ";
                File << endl;
            }

            File << endl;
        }

        File.close();
    }

    if(bDefaultRoundingScheme==false)
    {
        char     c;
        bool     bUseFile;
        string   filename;
        ifstream FileIn;
        CLargeInt      L;
        CLargeFraction F;

        while(true)
        {
            cout << "Use program's guess for eigenvectors [y/n] : ";
            cin.clear();
            cin.ignore(numeric_limits<streamsize>::max(),'\n');
            cin >> c;
        
            if(c=='y' || c=='Y')
            {
                bUseFile = false;
                break;
            }

            if(c=='n' || c=='N')
            {
                bUseFile = true;
                break;
            }
        }

        if(bUseFile==true)
        {
            cout << "Name of vector file : " << endl;
            cin.clear();
            cin.ignore(numeric_limits<streamsize>::max(),'\n');
            getline(cin,filename);
                
            FileIn.open(filename.c_str());
            if(FileIn.good()==false)
            {
                cout << endl;
                cout << "ERROR : Could not open \"" << filename << "\" for reading." << endl;
                return false;
            }

            c = ReadChar(FileIn);
        
            for(n=0;n<Soln.size();n++)
            for(k=0;k<noOfZeroVectors[n];k++)
            for(i=0;i<Soln[n].dMatrix.size();i++)
            {
                RemoveWhitespace(FileIn,c);

                if(ReadLargeInteger(FileIn,L,c)==false)
                {
                    cout << endl;
                    cout << "ERROR : Expected to read a fraction." << endl;
                    return false;
                }

                F = L;

                RemoveWhitespace(FileIn,c);

                //Check if it is a fraction.
                if(c=='/')
                {
                    c = ReadChar(FileIn);
                    RemoveWhitespace(FileIn,c);

                    if(ReadLargeInteger(FileIn,L,c)==false)
                    {
                        cout << endl;
                        cout << "ERROR : Expected to read a fraction." << endl;
                        return false;
                    }

                    if(L==0)
                    {
                        cout << endl;
                        cout << "ERROR : Denominator is zero." << endl;
                        return false;
                    }

                    F = F/L;
                }

                //Store F
                F.Simplify();
                ChangeOfBasis[n][i][k] = F;
            }
        }
    }

    cout << "Applying a change of basis..." << endl;

    //Fill in the rest of ChangeOfBasis.
    {
        vector<vector<CLargeFraction> > LinearProb;
        vector<vector<CLargeFraction> > LinearSoln;
        vector<vector<double> >         dChangeOfBasis;
        CTerm                           newTerm;

        for(n=0;n<Soln.size();n++)
        {
            //Find all vectors which evaluate to zero when dot producted with the 0-eigenvectors.

            LinearProb.resize(Soln[n].dMatrix.size()+1);
            for(i=0;i<LinearProb.size();i++)
                LinearProb[i].resize(noOfZeroVectors[n],0);

            for(i=0;i<LinearProb.size()-1;i++)
            for(j=0;j<LinearProb[0].size();j++)
                LinearProb[i][j] = ChangeOfBasis[n][i][j];

            LinearConstraintSolver<CLargeFraction>(LinearProb,LinearSoln,0);

            //Store the kernel in ChangeOfBasis.
            for(i=0;i<LinearSoln.size();i++)
            for(j=0;j<LinearSoln[0].size()-1;j++)
                ChangeOfBasis[n][i][noOfZeroVectors[n]+j] = LinearSoln[i][j];

            //Convert the entries of ChangeOfBasis into doubles.
            dChangeOfBasis.resize(ChangeOfBasis[n].size());
            for(i=0;i<dChangeOfBasis.size();i++)
                dChangeOfBasis[i].resize(ChangeOfBasis[n].size());

            for(i=0;i<dChangeOfBasis.size();i++)
            for(j=0;j<dChangeOfBasis.size();j++)
                dChangeOfBasis[i][j] = ConvertFracToDouble(ChangeOfBasis[n][i][j]);

            //Transform Soln[n].dMatrix by dChangeOfBasis.

            //Clear newTerm.
            newTerm.nType = Soln[n].nType;
            newTerm.iType = Soln[n].iType;
            vector<vector<CLargeFraction> >().swap(newTerm.fMatrix);
            vector<vector<double> >().swap(newTerm.dMatrix);
            vector<vector<long> >().swap(newTerm.BasisFlag);
            vector<vector<CLargeFraction> >().swap(newTerm.BasisCoeff);

            newTerm.dMatrix.resize(Soln[n].dMatrix.size()-noOfZeroVectors[n]);
            for(i=0;i<newTerm.dMatrix.size();i++)
                newTerm.dMatrix[i].resize(newTerm.dMatrix.size());

            for(i=0;i<newTerm.dMatrix.size();i++)
            for(j=0;j<newTerm.dMatrix.size();j++)
            {
                newTerm.dMatrix[i][j] = 0;

                for(p=0;p<Soln[n].dMatrix.size();p++)
                for(q=0;q<Soln[n].dMatrix.size();q++)
                    newTerm.dMatrix[i][j] += 
                        dChangeOfBasis[p][i+noOfZeroVectors[n]]
                       *Soln[n].dMatrix[p][q]
                       *dChangeOfBasis[q][j+noOfZeroVectors[n]];
            }

            //Calculate the new basis.
            ChangeOfBasis[n] = Inverse(ChangeOfBasis[n]);

            newTerm.BasisFlag.resize(newTerm.dMatrix.size());
            newTerm.BasisCoeff.resize(newTerm.dMatrix.size());

            for(i=0;i<newTerm.dMatrix.size();i++)
            {
                vector<CLargeFraction> coeff;

                coeff.resize(FInfo[newTerm.nType][newTerm.iType].FlagList.size(),0);

                for(j=0;j<ChangeOfBasis[n].size();j++)
                for(k=0;k<Soln[n].BasisFlag[j].size();k++)
                    coeff[Soln[n].BasisFlag[j][k]] = coeff[Soln[n].BasisFlag[j][k]]
                                                +ChangeOfBasis[n][i+noOfZeroVectors[n]][j]*Soln[n].BasisCoeff[j][k];

                for(j=0;j<coeff.size();j++)
                    if(coeff[j]!=0)
                    {
                        newTerm.BasisFlag[i].push_back(j);
                        newTerm.BasisCoeff[i].push_back(coeff[j]);
                    }
            }

            //Replace Soln[n].
            Soln[n] = newTerm;
        }
    }

    return true;
}

bool ProjectSolution()
{
    long    i, j, k, n, r, c, h, var, cons;
    long*   pVal;
    KeyProd key;
    CLargeFraction          product;
    vector<vector<long> >   vectorToMatrixMap;
    vector<vector<long> >   unusedVarMap;
    vector<vector<SparseComponent<CLargeFraction> > > SharpSoln;
    vector<vector<double> > ProjSoln;
    CLargeInt               Accuracy;

    Accuracy = 1;
    for(i=0;i<accuracyProjectSolution;i++)
        Accuracy = Accuracy*10;
    
    //Create a mapping from matrix elements to a vector.
    vectorToMatrixMap.resize(0);
    for(n=0;n<Soln.size();n++)
    for(i=0;i<Soln[n].dMatrix.size();i++)
    for(j=i;j<Soln[n].dMatrix.size();j++)
    {
        vectorToMatrixMap.push_back(vector<long>());
        vectorToMatrixMap.back().resize(3);

        vectorToMatrixMap.back()[0] = n;
        vectorToMatrixMap.back()[1] = i;
        vectorToMatrixMap.back()[2] = j;
    }

    if(vectorToMatrixMap.size()==0)
        return true; //No matrix entries.

    cout << "Solving the sharp inequality problem..." << endl;
    //Find all solution vectors which cause the sharp inequalties to equal the TargetTuranDensity.
    {
        vector<vector<CLargeFraction> > SharpProb;

        SharpProb.resize(vectorToMatrixMap.size()+1);
        for(i=0;i<SharpProb.size();i++)
            SharpProb[i].resize(sharpIneq.size());

        //Set the constants.
        for(i=0;i<SharpProb[0].size();i++)
            SharpProb[SharpProb.size()-1][i] = TargetTuranDensity
                        -CLargeFraction(NumberOfEdges(H[sharpIneq[i]])*6)/(orderOfH*(orderOfH-1)*(orderOfH-2));

        //Set the rest of the variables.
        for(var=0;var<SharpProb.size()-1;var++)
        for(cons=0;cons<SharpProb[0].size();cons++)
        {
            h = sharpIneq[cons];
            n = vectorToMatrixMap[var][0];
            r = vectorToMatrixMap[var][1];
            c = vectorToMatrixMap[var][2];

            key.nType = Soln[n].nType;
            key.iType = Soln[n].iType;
            
            //Calculate the entry's contribution.
            SharpProb[var][cons] = 0;

            //Multiply the flags in the basis.
            for(i=0;i<Soln[n].BasisFlag[r].size();i++)
            for(j=0;j<Soln[n].BasisFlag[c].size();j++)
            {
                //Look up the product of the flags.

                //Set the key.
                if(Soln[n].BasisFlag[r][i]<=Soln[n].BasisFlag[c][j])
                {
                    key.iFlag1 = Soln[n].BasisFlag[r][i];
                    key.iFlag2 = Soln[n].BasisFlag[c][j];
                }
                else
                {
                    key.iFlag1 = Soln[n].BasisFlag[c][j];
                    key.iFlag2 = Soln[n].BasisFlag[r][i];
                }
                                        
                pVal = htProd[h].Get(key);
                                        
                if(pVal!=NULL)
                    SharpProb[var][cons] = SharpProb[var][cons]
                                          +Soln[n].BasisCoeff[r][i]*Soln[n].BasisCoeff[c][j]*(*pVal)/Factorial[orderOfH];
            }

            if(r!=c)
                SharpProb[var][cons] = SharpProb[var][cons]*2;
        }

        //Get rid of variables which do not make a contribution to the SharpProb inequalities.
        {
            vector<vector<CLargeFraction> > SharpProbShrink;
            vector<vector<long> >           usedVarMap;
            bool isUsed;

            vector<vector<long> >().swap(unusedVarMap);

            for(var=0;var<SharpProb.size()-1;var++)
            {
                isUsed = false;
                for(cons=0;cons<SharpProb[0].size();cons++)
                    if(SharpProb[var][cons]!=0)
                    {
                        isUsed = true;
                        break;
                    }

                if(isUsed==true)
                {
                    //Add to SharpProbShrink.
                    SharpProbShrink.push_back(SharpProb[var]);

                    //Add to usedVarMap.
                    usedVarMap.push_back(vectorToMatrixMap[var]);
                }
                else
                    unusedVarMap.push_back(vectorToMatrixMap[var]); //Add to unusedVarMap.
            }

            //Add constants.
            SharpProbShrink.push_back(SharpProb[SharpProb.size()-1]);

            //Swap.
            vectorToMatrixMap.swap(usedVarMap);
            SharpProb.swap(SharpProbShrink);
        }

        //Solve the linear equations.

        //Check if there are no variables to find in SharpProb.
        if(SharpProb.size()==1)
        {
            //Check all the constants are 0.
            for(i=0;i<SharpProb[0].size();i++)
                if(SharpProb[0][i]!=0)
                {
                    cout << endl;
                    cout << "ERROR : The matrices could not satisfy the sharp inequalities." << endl;
                    return false;
                }

            //Rationalize Soln and return.

            //Allocate the fraction version of dMatrix. 
            for(n=0;n<Soln.size();n++)
            {
                Soln[n].fMatrix.resize(Soln[n].dMatrix.size());
                for(i=0;i<Soln[n].fMatrix.size();i++)
                    Soln[n].fMatrix[i].resize(Soln[n].fMatrix.size());
            }

            //Rationalize dMatrix.

            //Accuracy = Accuracy*Accuracy;
            for(n=0;n<Soln.size();n++)
            {
                for(i=0;i<Soln[n].fMatrix.size();i++)
                for(j=0;j<Soln[n].fMatrix.size();j++)
                    Soln[n].fMatrix[i][j] = ConvertDoubleToFrac(Soln[n].dMatrix[i][j],Accuracy);
            }

            //Make symmetric.
            CLargeFraction F;
            for(n=0;n<Soln.size();n++)
            {
                for(i=0;i<Soln[n].fMatrix.size();i++)
                for(j=i+1;j<Soln[n].fMatrix.size();j++)
                {
                    F = (Soln[n].fMatrix[i][j]+Soln[n].fMatrix[j][i])/2;
                    Soln[n].fMatrix[i][j] = F;
                    Soln[n].fMatrix[j][i] = F;
                }
            }

            return true;
        }

        SparseLinearConstraintSolver<CLargeFraction>(SharpProb,SharpSoln,0);

        if(SharpSoln.size()<=0)
        {
            cout << endl;
            cout << "ERROR : The matrices could not satisfy the sharp inequalities." << endl;
            return false;
        }
    }
            
    cout << "Solving the projection problem..." << endl;
    //Project the dMatrix solution onto the SharpSoln.
    {
        vector<vector<double> > ProjProb;
        vector<double>          offset;
        vector<vector<long> >   commonEntry;      //[component index][list index] = kernel number.
        vector<vector<double> > commonEntryCoeff; //[component index][list index] = coefficient of kernel entry.

        offset.resize(vectorToMatrixMap.size());
        for(var=0;var<offset.size();var++)
        {
            n = vectorToMatrixMap[var][0];
            r = vectorToMatrixMap[var][1];
            c = vectorToMatrixMap[var][2];
            
            offset[var] = Soln[n].dMatrix[r][c];
        }

        for(i=0;i<SharpSoln.back().size();i++)
            offset[SharpSoln.back()[i].index] -= ConvertFracToDouble(SharpSoln.back()[i].coeff);

        //We need to find a coefficient for each member of the kernel in SharpSoln.
        ProjProb.resize(SharpSoln.size());

        //We get one constraint per kernel vector.
        for(i=0;i<ProjProb.size();i++)
            ProjProb[i].resize(SharpSoln.size()-1);
        
        //Set the constants.
        for(cons=0;cons<SharpSoln.size()-1;cons++)
        {
            //The constant is the offset dot producted with the kernel vector.
            ProjProb[ProjProb.size()-1][cons] = 0;
            
            for(i=0;i<SharpSoln[cons].size();i++)
                ProjProb[ProjProb.size()-1][cons] += offset[SharpSoln[cons][i].index]
                                                    *ConvertFracToDouble(SharpSoln[cons][i].coeff);
        }

        //Calculate the common entries of the kernel elements.
        commonEntry.resize(offset.size());
        commonEntryCoeff.resize(offset.size());
        for(i=0;i<ProjProb.size()-1;i++)
        for(j=0;j<SharpSoln[i].size();j++)
        {
            commonEntry[SharpSoln[i][j].index].push_back(i);
            commonEntryCoeff[SharpSoln[i][j].index].push_back(ConvertFracToDouble(SharpSoln[i][j].coeff));
        }

        //Set the rest of the entries.
        //ProjProb[var][cons] = dot product of SharpSoln[var] and SharpSoln[cons].
        for(var=0;var<ProjProb.size()-1;var++)
        for(cons=0;cons<ProjProb[0].size();cons++)
            ProjProb[var][cons] = 0;

        for(i=0;i<commonEntry.size();i++)
        for(j=0;j<commonEntry[i].size();j++)
        for(k=0;k<commonEntry[i].size();k++)
            ProjProb[commonEntry[i][j]][commonEntry[i][k]] += commonEntryCoeff[i][j]*commonEntryCoeff[i][k];

        //Solve the linear equations.
        LinearConstraintSolver<double>(ProjProb,ProjSoln,toleranceProjectSolution);
    }
    
    cout << "Converting the solution into matrices..." << endl;
    //Map the projected solution back into matrices.
    {
        //Convert ProjSoln into fractions.
        vector<CLargeFraction> fProjSoln;

        fProjSoln.resize(ProjSoln.size()); 
        for(i=0;i<fProjSoln.size();i++)
            fProjSoln[i] = ConvertDoubleToFrac(ProjSoln[i][ProjSoln[i].size()-1],Accuracy);

        //Allocate the fraction version of dMatrix. 
        for(n=0;n<Soln.size();n++)
        {
            Soln[n].fMatrix.resize(Soln[n].dMatrix.size());
            for(i=0;i<Soln[n].fMatrix.size();i++)
                Soln[n].fMatrix[i].resize(Soln[n].fMatrix.size());
        }

        //Zero the matrix.
        for(n=0;n<Soln.size();n++)
        {
            for(i=0;i<Soln[n].fMatrix.size();i++)
            for(j=0;j<Soln[n].fMatrix.size();j++)
                Soln[n].fMatrix[i][j] = 0;
        }

        //Apply the unused variables.

        //Accuracy = Accuracy*Accuracy;
        for(k=0;k<unusedVarMap.size();k++)
        {
            n = unusedVarMap[k][0];
            r = unusedVarMap[k][1];
            c = unusedVarMap[k][2];

            Soln[n].fMatrix[r][c] = ConvertDoubleToFrac(Soln[n].dMatrix[r][c],Accuracy);
            Soln[n].fMatrix[c][r] = ConvertDoubleToFrac(Soln[n].dMatrix[c][r],Accuracy);
        }
        
        //Apply the SharpSoln.back() solution.
        for(k=0;k<SharpSoln.back().size();k++)
        {
            n = vectorToMatrixMap[SharpSoln.back()[k].index][0];
            r = vectorToMatrixMap[SharpSoln.back()[k].index][1];
            c = vectorToMatrixMap[SharpSoln.back()[k].index][2];

            Soln[n].fMatrix[r][c] = Soln[n].fMatrix[r][c]+SharpSoln.back()[k].coeff;
        }
        
        //Apply the kernel vectors.
        for(i=0;i<SharpSoln.size()-1;i++)
        for(k=0;k<SharpSoln[i].size();k++)
        {
            n = vectorToMatrixMap[SharpSoln[i][k].index][0];
            r = vectorToMatrixMap[SharpSoln[i][k].index][1];
            c = vectorToMatrixMap[SharpSoln[i][k].index][2];

            Soln[n].fMatrix[r][c] = Soln[n].fMatrix[r][c]+SharpSoln[i][k].coeff*fProjSoln[i];
        }

        //Make symmetric.
        for(k=0;k<vectorToMatrixMap.size();k++)
        {
            n = vectorToMatrixMap[k][0];
            r = vectorToMatrixMap[k][1];
            c = vectorToMatrixMap[k][2];
            Soln[n].fMatrix[c][r] = Soln[n].fMatrix[r][c];
        }

        CLargeFraction F;
        for(n=0;n<Soln.size();n++)
        {
            for(i=0;i<Soln[n].fMatrix.size();i++)
            for(j=i+1;j<Soln[n].fMatrix.size();j++)
            {
                F = (Soln[n].fMatrix[i][j]+Soln[n].fMatrix[j][i])/2;
                Soln[n].fMatrix[i][j] = F;
                Soln[n].fMatrix[j][i] = F;
            }
        }
    }

    return true;
}

//Returns true if all matrices in Soln[] are positive semidefinite (psd). A bad_alloc
//or ErrorOverflow object may be thrown.
bool IsPSD()
{
    long i, j, k, s, n;
    long nTerm;
    long noOfTerms;
    long dim;
    long globalError;
    long loopError;
    CLargeFraction   F;
    vector<vector<CLargeFraction> > M;

    cout << endl;
    cout << "Checking the matrices are positive semidefinite..." << endl;
    
    noOfTerms = Soln.size();
    for(nTerm=1;nTerm<=noOfTerms;nTerm++)
    {
        if(Soln[nTerm-1].fMatrix.size()<=0)
            continue;

        cout << "\rChecking matrix " << nTerm << " of " << noOfTerms << " : ";
        cout << "Copying matrix.   " << flush;

        M   = Soln[nTerm-1].fMatrix;
        dim = M.size();

        //Make M[][] symmetric.
        for(i=0;i<dim;i++)
        for(j=i+1;j<dim;j++)
        {
            F       = (M[i][j]+M[j][i])/2;
            M[i][j] = F;
            M[j][i] = F;
        }

        //The matrix M is psd if and only if R^T M R is psd, where R is any invertible
        //matrix. The transformation of adding lambda times column i to column j of matrix M
        //can be represented by a matrix C such that M C gives the desired result. The matrix
        //C is invertible when i and j are distinct. Furthermore C^T M is equivalent to
        //adding lambda times row i to row j of M. Hence by doing symmetric row and column
        //operations on M we can transform it into another symmetric matrix which retains
        //whether it is psd or not.
        for(n=0;n<dim;n++)
        {
            cout << "\rChecking matrix " << nTerm << " of " << noOfTerms << " : ";
            cout << "Checking row " << n+1 << " of " << dim << "." << flush;

            if(M[n][n]<0)
            {
                cout << endl;
                cout << endl;
                cout << "*** ERROR: Matrix is not positive semidefinite. ***" << endl;
                return false;
            }

            if(M[n][n]==0)
            {
                //Check all other entries are 0.
                for(i=n+1;i<dim;i++)
                    if(M[n][i]!=0)
                    {
                        //The submatrix {M[n][n] M[n][i]}, {M[i][n] M[i][i]}, is not positive
                        //semidefinite and hence neither is M[][].
                        cout << endl;
                        cout << endl;
                        cout << "*** ERROR: Matrix is not positive semidefinite. ***" << endl;
                        return false;
                    }

                continue;
            }

            //Subtract a multiple of row n from row i, such that M[i][n] = 0 (for i in
            //[n+1,dim-1]). Then repeat the operations for columns instead of rows (to keep
            //the matrix in the form R^T M R). 

            //To speed up the calculation temporarily divide row n by M[n][n].
            for(i=n+1;i<dim;i++)
                M[n][i] = M[n][i]/M[n][n];

            //Subtract row n from other rows.
            
            //Was:
            //for(i=n+1;i<dim;i++)
            //for(j=n+1;j<dim;j++)
            //   M[i][j] = M[i][j]-M[i][n]*M[n][j];
        

            globalError = 0; //No error.
            s           = (dim-n-1)*(dim-n-1);
            #pragma omp parallel for schedule(dynamic) default(shared) private(i,j,k,loopError)
            for(k=0;k<s;k++)
            {
                i = k/(dim-n-1)+n+1;
                j = k%(dim-n-1)+n+1;
                
                try
                {
                    #pragma omp critical (IsPSD_Progress)
                    {
                        loopError = globalError;
                    }

                    if(loopError==0)
                        M[i][j] = M[i][j]-M[i][n]*M[n][j];
                }
                catch(ErrorDivisionByZero&)
                {
                    loopError = 4;
                }
                catch(ErrorOverflow&)
                {
                    loopError = 3;
                }
                catch(bad_alloc&)
                {
                    loopError = 2;
                }
                catch(...)
                {
                    loopError = 1;
                }
        
                #pragma omp critical (IsPSD_Progress)
                {
                    if(loopError>globalError)
                        globalError = loopError;
                }
            }

            //Display any errors.
            if(globalError!=0)
            {
                cout << endl << endl;
        
                if(globalError==4)
                    cout << "ERROR : Division by zero.";

                if(globalError==3)
                    cout << "ERROR : Numeric overflow occurred. Try increasing LI_MAX in the source code.";

                if(globalError==2)
                    cout << "ERROR : Allocating memory failed.";

                if(globalError==1)
                    cout << "ERROR : An unknown error occurred.";

                cout << endl << endl;
                return false;
            }

            //Now we need to restore row n, and do the column operations. It is easy to show
            //that this will not effect the entries M[i][j] where i,j>=n+1. Furthermore row n
            //and column n will end up consisting entirely of zeros, with the exception of
            //M[n][n] which will be some positive value. Hence to save time we will skip
            //these operations. It is enough to check the submatrix consisting of the entries
            //M[i][j] where i,j>=n+1. The entire matrix is psd if and only if the submatrix
            //is psd.
        }

        cout << "\rChecking matrix " << nTerm << " of " << noOfTerms << " : ";
        cout << "                  ";
        for(i=0;i<2*NumberOfDigits(dim);i++)
            cout << " ";
    }

    cout << "\r                                         ";
    for(i=0;i<2*NumberOfDigits(noOfTerms);i++)
        cout << " ";

    cout << "\rDone." << endl;
    cout << endl;

    return true;
}

void InitSoln()
{
    long i, j, n, t, m;

    vector<CTerm>().swap(Soln);
            
    //Determine the blocks we are optimizing over.
    for(n=0;n<=orderOfH-2;n++)
    for(t=0;t<FInfo[n].size();t++)
    for(m=0;m<2;m++)
        if(FInfo[n][t].Basis[m].size()>0)
        {
            Soln.push_back(CTerm());

            Soln.back().nType = n;
            Soln.back().iType = t;

            //Copy the basis.
            Soln.back().BasisFlag.resize(FInfo[n][t].Basis[m].size());
            Soln.back().BasisCoeff.resize(FInfo[n][t].Basis[m].size());

            for(i=0;i<FInfo[n][t].Basis[m].size();i++)
            {
                Soln.back().BasisFlag[i].resize(FInfo[n][t].Basis[m][i].size());
                Soln.back().BasisCoeff[i].resize(FInfo[n][t].Basis[m][i].size());

                for(j=0;j<FInfo[n][t].Basis[m][i].size();j++)
                {
                    Soln.back().BasisFlag[i][j] = FInfo[n][t].Basis[m][i][j];

                    if(m==0 || (m==1 && j==0))
                        Soln.back().BasisCoeff[i][j] = 1;
                    else
                        Soln.back().BasisCoeff[i][j] = -1;
                }
            }

            //Allocate memory for the dMatrix.
            Soln.back().dMatrix.resize(Soln.back().BasisFlag.size());
            for(i=0;i<Soln.back().dMatrix.size();i++)
                Soln.back().dMatrix[i].resize(Soln.back().dMatrix.size(),0);
        }

    return;
}

bool ReadOutFile_PrimalCSDP()
{
    long     i, j, m, n;
    double   d;
    ifstream File;
    
    InitSoln();
    
    File.open(FilenameOut.c_str());

    if(File.good()==false)
    {
        cout << endl;
        cout << "ERROR : Could not open \"" << FilenameOut << "\" for reading." << endl << endl;
        return false;
    }

    //Read in the variables y (and skip them).
    for(i=0;i<H.size();i++)
        File >> d;

    while(true)
    {
        //Read "Matrix number" "Block number" row column value.
        File >> m >> n >> i >> j >> d;

        if(File.eof()==true)
            break;

        if(File.good()==false)
        {
            cout << endl;
            cout << "ERROR : Reading \"" << FilenameOut << "\" failed." << endl << endl;
            return false;
        }

        if(m!=2 || n>Soln.size())
            continue;

        if(i>Soln[n-1].dMatrix.size() || j>Soln[n-1].dMatrix.size())
        {
            cout << endl;
            cout << "ERROR : Reading \"" << FilenameOut << "\" failed." << endl << endl;
            return false;
        }

        Soln[n-1].dMatrix[i-1][j-1] = d;
        Soln[n-1].dMatrix[j-1][i-1] = d;
    }

    File.close();
    return true;
}

bool ReadOutFile_DualCSDP()
{
    long     n, r, c;
    double   d;
    ifstream File;
    
    InitSoln();
    
    File.open(FilenameOut.c_str());

    if(File.good()==false)
    {
        cout << endl;
        cout << "ERROR : Could not open \"" << FilenameOut << "\" for reading." << endl << endl;
        return false;
    }

    //Read in the density bound $\rho$ (and ignore it).
    File >> d;

    for(n=0;n<Soln.size();n++)
    for(r=0;r<Soln[n].dMatrix.size();r++)
    for(c=r;c<Soln[n].dMatrix.size();c++)
    {
        File >> d;
                
        if(File.eof()==true || File.good()==false)
        {
            cout << endl;
            cout << "ERROR : Reading \"" << FilenameOut << "\" failed." << endl << endl;
            return false;
        }

        Soln[n].dMatrix[r][c] = d;
        Soln[n].dMatrix[c][r] = d;
    }

    File.close();
    return true;
}

//Reads through "in" until it finds "search".
//If successful it returns true.
//If "search" cannot be found it returns false.
bool SearchForHeader(istream& in, string search)
{
    long i;
    char c;
    vector<long> found;

    if(search.size()==0)
        return true;

    found.push_back(0);
    while(true)
    {
        //Check if there are any members of found which indicate they've completed.
        for(i=0;i<found.size();i++)
        {
            if(found[i]==search.size())
                return true;
        }

        c = ReadChar(in);
                
        if(c=='\0')
            return false; //Reached the end of the stream.

        //Update found[].
        for(i=0;i<found.size();i++)
        {
            //Check if entry is unused.
            if(found[i]==-1)
                continue;

            if(search[found[i]]==c)
                found[i]++; //The next character of "search" has been found.
            else
                found[i] = -1; //Restart.
        }

        //Find an unused entry to become 0.
        for(i=0;i<found.size();i++)
            if(found[i]==-1)
                break;

        if(i==found.size())
            found.push_back(0); //No unused entries, so add an extra one.
        else
            found[i] = 0;
    }

    return false;
}

//A rough method to handle doubles in SDPA out files.
//Returns false if it reaches the end of the stream before a double could be read.
bool ReadSDPADouble(istream& in, double& d)
{
    char c;
    stringstream sstr;

    d = 0;

    //Ignore "{}, \t\r\n".
    while(true)
    {
        c = ReadChar(in);

        if(c=='\0')
            return false;

        if(c=='{' || c=='}' || c==',' || c==' ' || c=='\t' || c=='\r' || c=='\n')
            continue;

        break;
    }

    //Read in the characters forming the double.
    sstr << c;
    while(true)
    {
        c = ReadChar(in);

        if(c=='\0')
            break;

        if(c=='{' || c=='}' || c==',' || c==' ' || c=='\t' || c=='\r' || c=='\n')
            break;

        sstr << c;
    }

    sstr >> d;

    if(sstr.fail()==true)
        return false;

    return true;
}

bool ReadOutFile_PrimalSDPA()
{
    long     n, r, c;
    double   d;
    ifstream File;
    
    InitSoln();
    
    File.open(FilenameOut.c_str());

    if(File.good()==false)
    {
        cout << endl;
        cout << "ERROR : Could not open \"" << FilenameOut << "\" for reading." << endl << endl;
        return false;
    }

    //Find "yMat ="
    if(SearchForHeader(File,"yMat =")==false)
    {
        cout << endl;
        cout << "ERROR : Reading \"" << FilenameOut << "\" failed." << endl << endl;
        return false;
    }

    //Read in the matrix entries.
    for(n=0;n<Soln.size();n++)
    {
        for(r=0;r<Soln[n].dMatrix.size();r++)
        for(c=0;c<Soln[n].dMatrix.size();c++)
            if(ReadSDPADouble(File,Soln[n].dMatrix[r][c])==false)
            {
                cout << endl;
                cout << "ERROR : Reading \"" << FilenameOut << "\" failed." << endl << endl;
                return false;
            }

        //Make symmetric.
        for(r=0;r<Soln[n].dMatrix.size();r++)
        for(c=r+1;c<Soln[n].dMatrix.size();c++)
        {
            d = (Soln[n].dMatrix[r][c]+Soln[n].dMatrix[c][r])/2;
            Soln[n].dMatrix[r][c] = d;
            Soln[n].dMatrix[c][r] = d;
        }
    }

    File.close();
    return true;
}

bool ReadOutFile_DualSDPA()
{
    long     n, r, c;
    double   d;
    ifstream File;
    
    InitSoln();
    
    File.open(FilenameOut.c_str());

    if(File.good()==false)
    {
        cout << endl;
        cout << "ERROR : Could not open \"" << FilenameOut << "\" for reading." << endl << endl;
        return false;
    }

    //Find "xVec ="
    if(SearchForHeader(File,"xVec =")==false)
    {
        cout << endl;
        cout << "ERROR : Reading \"" << FilenameOut << "\" failed." << endl << endl;
        return false;
    }
    
    //Read in the density bound $\rho$ (and ignore it).
    if(ReadSDPADouble(File,d)==false)
    {
        cout << endl;
        cout << "ERROR : Reading \"" << FilenameOut << "\" failed." << endl << endl;
        return false;
    }

    //Read in the matrix entries.
    for(n=0;n<Soln.size();n++)
    for(r=0;r<Soln[n].dMatrix.size();r++)
    for(c=r;c<Soln[n].dMatrix.size();c++)
    {    
        if(ReadSDPADouble(File,d)==false)
        {
            cout << endl;
            cout << "ERROR : Reading \"" << FilenameOut << "\" failed." << endl << endl;
            return false;
        }

        Soln[n].dMatrix[r][c] = d;
        Soln[n].dMatrix[c][r] = d;
    }

    File.close();
    return true;
}

bool ReadOutFile()
{
    bool     returnCode;
    bool     bIsSDPA;
    ifstream File;

    cout << "Checking .out file..." << flush;
    File.open(FilenameOut.c_str());

    if(File.good()==false)
    {
        cout << endl;
        cout << endl;
        cout << "ERROR : Could not open \"" << FilenameOut << "\" for reading." << endl << endl;
        return false;
    }

    //Find "xVec ="
    if(SearchForHeader(File,"xVec =")==false)
        bIsSDPA = false;
    else
        bIsSDPA = true;

    File.close();

    if(bIsSDPA==true)
        cout << "\rReading SDPA .out file..." << endl;
    else
        cout << "\rReading CSDP .out file..." << endl;

    if(bIsSDPA==true)
    {
        if(ForceForm==1)
            returnCode = ReadOutFile_PrimalSDPA();
        else
            returnCode = ReadOutFile_DualSDPA();
    }
    else
    {
        if(ForceForm==1)
            returnCode = ReadOutFile_PrimalCSDP();
        else
            returnCode = ReadOutFile_DualCSDP();
    }

    if(returnCode==false)
        return false;

    cout << "Done." << endl;
    cout << endl;

    return true;
}

//Displays 'Density' as both a fraction and a decimal. It may throw an
//ErrorOverflow object.
void DisplayDensity(CLargeFraction Density, long noOfDecimalPlaces)
{
    long      i;
    bool      bNegative;
    bool      bRounded;
    CLargeInt Digits, PowerOf10;
    
    Density.Simplify();

    cout << "\rUpper bound of the Turan density" << endl;
    
    if(Density.D==1)
    {
        cout << " = " << Density.N << endl;
        cout << endl;
        return;
    }

    cout << " = " << Density;
    
    //Display Density as a decimal to 'noOfDecimalPlaces' decimal places.

    //Check if Density is negative.
    if(Density<0)
    {
        Density   = -Density;
        bNegative = true;
    }
    else
        bNegative = false;

    //Calculate the digits that will make up the decimal.
    
    //Calculate 10^noOfDecimalPlaces.
    PowerOf10 = 1;
    for(i=0;i<noOfDecimalPlaces;i++)
        PowerOf10 = PowerOf10*10;

    Density = Density*PowerOf10;

    //'Digits' holds Density*10^noOfDecimalPlaces rounded down. We wish to round up.
    Digits = Density.N/Density.D;
    
    //Check if rounding up is necessary.
    if(Density.N%Density.D==0)
        bRounded = false;
    else
    {
        //Rounding is necessary, currently the magnitude has been rounded down. If the result
        //is negative then this is fine. If the result is positive then we add one so that
        //the magnitude is rounded up.
        bRounded = true;

        if(bNegative==false)
            Digits = Digits+1;
    }

    //Start displaying the digits.
    cout << endl;
    cout << " = ";
    
    if(bNegative==true)
        cout << "-";

    //Display the digits that lie to the left of the decimal point.
    cout << Digits/PowerOf10 << ".";

    //Display the digits that lie to the right of the decimal point, one at a time.
    Digits = Digits%PowerOf10;
    for(i=0;i<noOfDecimalPlaces;i++)
    {
        //Display the digit (10*Digits)/PowerOf10 rounded down.
        Digits = Digits*10;

        cout << Digits/PowerOf10;

        Digits = Digits%PowerOf10;
    }

    if(bRounded==true)
        cout << " (rounded up)";
    
    cout << endl;
    cout << endl;

    return;
}

bool WriteProofFile()
{
    long     i, j, k, m, n;
    ofstream FileOut;

    ScaleBasis();

    cout << "Writing the data into \"" << FilenameTxt << "\"." << endl;
    
    FileOut.open(FilenameTxt.c_str(),ios::trunc);

    if(FileOut.good()==false)
    {
        cout << endl << endl;
        cout << "ERROR: Cannot create the file." << endl;
        return false;
    }
    
    FileOut << "#Density : ";
    if(TargetTuranDensity.D==1)
        FileOut << TargetTuranDensity.N << endl;
    else
        FileOut << TargetTuranDensity << endl;

    FileOut << endl;
    FileOut << "Program                       : DensityChecker" <<endl;
    FileOut << endl;
    FileOut << "Order of subgraphs H          : " << orderOfH << endl;
    FileOut << "Number of forbidden graphs    : " << ForbiddenGraphs.size() << endl;
    FileOut << "Forbidden graphs              :" << endl;

    for(n=0;n<ForbiddenGraphs.size();n++)
    {
        if(bForbidBlowups==true && n==ForbiddenNonBlowups)
            FileOut << "#Forbidden graphs added by supersaturation argument :" << endl;

        FileOut << ForbiddenGraphs[n].nV << " : ";
        DisplayEdges(FileOut,ForbiddenGraphs[n]);
        FileOut << endl;
    }
    
    FileOut << "Forbid only induced subgraphs : " << (bForbidInduced==true?"yes":"no") << endl;
    FileOut << "Number of terms               : ";

    m = 0;
    for(n=0;n<Soln.size();n++)
        if(Soln[n].fMatrix.size()>0)
            m++;

    FileOut << m << endl;
    FileOut << endl;
    
    //Write out the Terms.
    m = 0;
    for(n=0;n<Soln.size();n++)
    {
        CTerm&       Term = Soln[n];
        CFlagInfo&   FI   = FInfo[Soln[n].nType][Soln[n].iType];
        vector<long> FlagMap;

        //Ignore zero dimensional matrices.
        if(Term.fMatrix.size()<=0)
            continue;

        m++;

        FileOut << endl;
        FileOut << "--- Term " << m << " ---" << endl;
        FileOut << endl;
        FileOut << "Order of type    : " << FI.Type.nV << endl;
        FileOut << "Order of flags   : " << FI.FlagList[0].nV << endl;

        //Relabel flags to eliminate flags which do not appear in the basis.
        {
            FlagMap.resize(FI.FlagList.size(),-1);

            for(i=0;i<Term.fMatrix.size();i++)
            for(j=0;j<Term.BasisFlag[i].size();j++)
                if(Term.BasisCoeff[i][j]!=0)
                    FlagMap[Term.BasisFlag[i][j]] = 0;

            //Relabel.
            k = 0;
            for(i=0;i<FlagMap.size();i++)
                if(FlagMap[i]==0)
                {
                    k++;
                    FlagMap[i] = k;
                }
        }

        FileOut << "Number of flags  : " << k << endl;
        FileOut << "Matrix dimension : " << Term.fMatrix.size() << endl;
        FileOut << endl;
        FileOut << "Flags :" << endl;

        for(i=0;i<FI.FlagList.size();i++)
        {
            if(FlagMap[i]==-1)
                continue;

            FileOut << "F" << FlagMap[i];
            
            for(j=0;j<NumberOfDigits(k)-NumberOfDigits(FlagMap[i]);j++)
                FileOut << " ";
            FileOut << " : ";
            
            DisplayEdges(FileOut,FI.FlagList[i]);
            FileOut << endl;
        }
        FileOut << endl;

        FileOut << "Basis :" << endl;
        for(i=0;i<Term.BasisFlag.size();i++)
        {
            bool      bFirst;
            CLargeInt bCoeff;

            FileOut << "B" << i+1;
    
            for(j=0;j<NumberOfDigits(Term.BasisFlag.size())-NumberOfDigits(i+1);j++)
                FileOut << " ";
            FileOut << " : ";

            bFirst = true;
            for(j=0;j<Term.BasisFlag[i].size();j++)
            {
                bCoeff = Term.BasisCoeff[i][j].N;

                if(bCoeff!=0)
                {
                    if(bCoeff>0)
                    {
                        if(bFirst==false)
                            FileOut << "+";
                    }
                    else
                    {
                        FileOut << "-";
                        bCoeff = -bCoeff;
                    }

                    if(bCoeff!=1)
                        FileOut << bCoeff;

                    FileOut << "(F";
                    FileOut << FlagMap[Term.BasisFlag[i][j]] << ")";

                    bFirst = false;
                }
            }

            FileOut << endl;
        }
        FileOut << endl;
            
        FileOut << "Matrix :" << endl;
        for(i=0;i<Term.fMatrix.size();i++)
        {
            for(j=0;j<Term.fMatrix.size();j++)
            {
                FileOut << Term.fMatrix[i][j];

                if(j!=Term.fMatrix.size()-1)
                    FileOut << " ";
            }

            FileOut << endl;
        }
        FileOut << endl;
    }

    FileOut.close();
    cout << endl;

    return true;
}

//Generate an un-blownup version of the graph gIn. Also removes isolated vertices.
void ReduceGraph(CGraph gIn, CGraph& gOut)
{
    long i, j, k, r;
    long noInMap;
    long map[9];

    for(i=0;i<9;i++)
    for(j=0;j<9;j++)
    for(k=0;k<9;k++)
        gOut.M[i][j][k] = false;

    for(i=0;i<9;i++)
    for(j=0;j<9;j++)
    for(k=0;k<9;k++)
        if(i==j || i==k || j==k)
            gIn.M[i][j][k] = false;

    noInMap = 0;
    for(r=0;r<gIn.nV;r++)
    {
        //Check if vertex r is isolated.
        for(i=0;i<gIn.nV;i++)
        {
            for(j=0;j<gIn.nV;j++)
                if(gIn.M[r][i][j]==true)
                    break;

            if(j!=gIn.nV)
                break;
        }

        if(i==gIn.nV)
            continue; //Vertex r is isolated.

        //Check if vertex r is a clone of another vertex.

        for(i=r+1;i<gIn.nV;i++)
        {
            //Check if r and i appear in an edge together.
            for(k=0;k<gIn.nV;k++)
                if(gIn.M[r][i][k]==true)
                    break;

            if(k!=gIn.nV)
                continue; //Found an edge.

            for(j=0;j<gIn.nV;j++)
            {
                for(k=0;k<gIn.nV;k++)
                    if(gIn.M[r][j][k]!=gIn.M[i][j][k])
                        break;

                if(k!=gIn.nV)
                    break; //Found a difference.
            }

            if(j==gIn.nV)
                break; //No differences found.
        }

        if(i!=gIn.nV)
            continue; //r is a clone of i.

        map[noInMap] = r;
        noInMap++;
    }

    //Create gOut the reduced graph.
    gOut.nV = noInMap;

    for(i=0;i<noInMap;i++)
    for(j=0;j<noInMap;j++)
    for(k=0;k<noInMap;k++)
        gOut.M[i][j][k] = gIn.M[map[i]][map[j]][map[k]];

    FindMinIso(gOut,gOut,0);

    return;
}

//Creates a family of graphs with |V|="order", which are blow-ups of S.
void Blowup(vector<CGraph>& Family, long order, const CGraph& S)
{
    long     i, j, k, n;
    long     c[10];
    CGraph   g, gMin;
    KeyFlag  key;
    ValFlag  val;
    CHashTable<KeyFlag,ValFlag> htList;
    vector<CGraph>              classFamily[10];
    vector<long>                index;
    vector<vector<long> >       classVertices;
        
    vector<CGraph>().swap(Family);
    
    key.typeOrder = 0;
    val.iFlag     = 0;
    val.iType     = 0;
    
    for(i=0;i<9;i++)
    for(j=0;j<9;j++)
    for(k=0;k<9;k++)
        KeyFlag::hashM[i][j][k] = GetRandomInt();

    for(i=0;i<9;i++)
        KeyFlag::hashOrder[i] = GetRandomInt();

    //Fill in the first few classFamily members.
    {
        for(i=0;i<9;i++)
        for(j=0;j<9;j++)
        for(k=0;k<9;k++)
            g.M[i][j][k] = false;

        g.nV = 0;
        classFamily[0].push_back(g);

        g.nV = 1;
        classFamily[1].push_back(g);

        g.nV = 2;
        classFamily[2].push_back(g);

        g.nV = 3;
        classFamily[3].push_back(g);

        if(NumberOfEdges(S)>0)
        {
            g.M[0][1][2] = true;
            g.M[0][2][1] = true;
            g.M[1][0][2] = true;
            g.M[1][2][0] = true;
            g.M[2][0][1] = true;
            g.M[2][1][0] = true;

            classFamily[3].push_back(g);
        }
    }

    if(order<=3)
    {
        Family = classFamily[order];
        return;
    }

    n = order;
    //for(n=4;n<=order;n++)
    {
        htList.Free();

        g.nV = n;

        for(i=0;i<g.nV+1;i++)
            c[i] = 0;

        c[0] = -1; 
        while(true)
        {
            //Increment. 
            i = 0;
            while(true)
            {
                c[i]++;

                if(c[i]==S.nV)
                {
                    c[i] = 0;
                    i++;
                }
                else
                    break;
            }

            //End condition.
            if(c[g.nV]!=0)
                break;

            //Check the labels are decreasing.
            for(i=0;i<g.nV-1;i++)
                if(c[i]<c[i+1])
                    break;

            if(i!=g.nV-1)
                continue;

            //Main body.

            for(i=0;i<9;i++)
            for(j=0;j<9;j++)
            for(k=0;k<9;k++)
            {
                g.M[i][j][k] = false;

                if(i==j || i==k || j==k || i>=g.nV || j>=g.nV || k>=g.nV)
                    continue;

                if(S.M[c[i]][c[j]][c[k]]==true)
                    g.M[i][j][k] = true;
            }

            FindMinIso(g,gMin,0);
        
            key.flag = gMin;
            if(htList.Get(key)==NULL)
            {
                htList.Add(key,val);
                Family.push_back(gMin);
            }
        }
    }
    
    return;
}

void FindForbiddenBlowups()
{
    long   i, j, k, p, q;
    CGraph g;
    list<CGraph> NewForbiddenGraphs;
    list<CGraph> TempBlowups;
    list<CGraph>::iterator it, jt, kt, endt;
    list<vector<CGraph> >::const_iterator itF;

    ForbiddenNonBlowups = ForbiddenGraphs.size();

    if(ForbiddenNonBlowups<=0)
        return;
    
    cout << "Extending the forbidden graph list..." << endl;

    //Copy ForbiddenGraphs vector into a list.
    NewForbiddenGraphs.clear();
    for(i=0;i<ForbiddenGraphs.size();i++)
        NewForbiddenGraphs.push_back(ForbiddenGraphs[i]);

    //Keep track of last element in NewForbiddenGraphs to ensure we don't delete any of the graphs
    //that were in the user defined ForbddenGraphs list.
    endt = NewForbiddenGraphs.end();
    endt--;

    //Go through each forbidden graph and find a non-covering pair of vertices to collapse into a new smaller forbidden
    //graph.
    it = NewForbiddenGraphs.begin();
    while(it!=NewForbiddenGraphs.end())
    {
        //TempBlowups will temporarily hold our new forbidden graphs before we add them to the list.
        TempBlowups.clear();

        for(i=0;i<it->nV;i++)
        for(j=i+1;j<it->nV;j++)
        {
            //Check if i,j are a non-covering pair of vertices.
            for(k=0;k<it->nV;k++)
            {
                if(k==i || k==j)
                    continue;

                if(it->M[i][j][k]==true)
                    break;
            }

            if(k==it->nV)
            {
                //Found a non-covering pair of vertices.

                //Union vertex j into i.
                g = *it;
                for(p=0;p<g.nV;p++)
                for(q=0;q<g.nV;q++)
                    if(it->M[j][p][q]==true)
                    {
                        g.M[i][p][q] = true;
                        g.M[p][i][q] = true;
                        g.M[p][q][i] = true;
                    }

                //Isolate vertex j.
                for(p=0;p<g.nV;p++)
                for(q=0;q<g.nV;q++)
                {
                    g.M[j][p][q] = false;
                    g.M[p][j][q] = false;
                    g.M[p][q][j] = false;
                }

                //ReduceGraph will remove vertex j, and reduce the graph further if it is a blowup.
                ReduceGraph(g,g);

                TempBlowups.push_back(g);
            }
        }

        //Remove subgraphs within TempBlowups.
        for(jt=TempBlowups.begin();jt!=TempBlowups.end();jt++)
        for(kt=TempBlowups.begin();kt!=TempBlowups.end();)
        {
            if(jt==kt)
            {
                kt++;
                continue;
            }

            if(IsSubgraph(*jt,*kt,false)==true)
                kt = TempBlowups.erase(kt);
            else
                kt++;
        }

        for(jt=NewForbiddenGraphs.begin();jt!=NewForbiddenGraphs.end();jt++)
        for(kt=TempBlowups.begin();kt!=TempBlowups.end();)
        {
            if(IsSubgraph(*jt,*kt,false)==true)
                kt = TempBlowups.erase(kt);
            else
                kt++;
        }

        //Remove subgraphs from NewForbiddnGraphs (making sure not to delete the original user specified graphs).
        kt = endt;
        kt++;
        while(kt!=NewForbiddenGraphs.end())
        {
            for(jt=TempBlowups.begin();jt!=TempBlowups.end();jt++)
                if(IsSubgraph(*jt,*kt,false)==true)
                    break;

            if(jt==TempBlowups.end())
            {
                kt++;
                continue;
            }

            //Erase.

            //Make sure "it" doesn't become invalidated.
            if(kt==it)
                it--; 

            kt = NewForbiddenGraphs.erase(kt);
        }

        //Add TempBlowups.
        NewForbiddenGraphs.splice(NewForbiddenGraphs.end(),TempBlowups);

        //Move to the next element.
        it++;
    }

    ForbiddenGraphs.clear();
    for(it=NewForbiddenGraphs.begin();it!=NewForbiddenGraphs.end();it++)
        ForbiddenGraphs.push_back(*it);

    cout << "Done." << endl << endl;

    if(ForbiddenNonBlowups<ForbiddenGraphs.size())
    {
        cout << "Additional forbidden graphs : " << endl;
        for(i=ForbiddenNonBlowups;i<ForbiddenGraphs.size();i++)
        {
            DisplayEdges(cout,ForbiddenGraphs[i]);
            cout << endl;
        }

        cout << endl;
    }

    return;
}

//Expicitly free the memory associated with the global variables.
//Note that the destructors of the objects will automatically free the memory anyway.
void Cleanup()
{
    long i;

    cout << "\rFreeing memory..." << flush;

    //Clear strings.
    {
        string().swap(FilenameSDP);
        string().swap(FilenameOut);
        string().swap(FilenameTxt);
        string().swap(SDPCommand);
        string().swap(RunCommand);
    }
    
    //Clear Iso.
    vector<vector<long> >().swap(Iso);

    //Clear H.
    vector<CGraph>().swap(H);

    //Clear the forbidden graph list.
    vector<CGraph>().swap(ForbiddenGraphs);

    //Clear FInfo[] and the type list.
    for(i=0;i<8;i++)
    {
        vector<CFlagInfo>().swap(FInfo[i]);
        vector<long>().swap(modifyTypeList[i]);
    }

    //Clear the flag hash table.
    htFlags.Free();

    //Clear the product hash tables.
    vector<CHashTable<KeyProd,long> >().swap(htProd);
    
    //Clear Soln and sharpIneq.
    vector<CTerm>().swap(Soln);
    vector<long>().swap(sharpIneq);

    cout << "\r                 \r" << flush;

    return;
}

int main(int argc, char* argv[])
{
    try
    {
        if(ProcessArguments(argc,argv)==false)
            return 0;
    
        if(FindSystemFormat(argv[0])==false)
            return 0;

        InitIso();

        if(bForbidBlowups==true)
            FindForbiddenBlowups();

        FillFlagInfo();
    
        if(GenerateH()==false)
            return 0;

        if(CalculateFlagProducts()==false)
            return 0;
        
        if(ChooseSDPFormat()==false)
            return 0;

        RunSDPSolver();

        if(ReadOutFile()==false)
            return 0;

        //To do: Rationalize the solution if it is not exact.

        if(CalculateSharpIneq()==false)
            return 0;

        if(MakeMatricesPositiveDefinite()==false)
            return 0;

        if(ProjectSolution()==false)
            return 0;

        //Check matrices are positive definite.
        if(IsPSD()==false)
            return 0;

        CLargeFraction Density;

        if(CalculateDensityBound(Density)==false)
            return 0;

        cout << endl;
        DisplayDensity(Density,10);

        if(TargetTuranDensity!=Density)
        {
            cout << endl << endl;
            cout << "ERROR : Failed to reach the target density." << endl << endl;
            //return 0;

            //Output best bound anyway.
            TargetTuranDensity = Density;
        }

        if(WriteProofFile()==false)
            return 0;
    }
    catch(ErrorOverflow&)
    {
        cout << endl << endl;
        cout << "ERROR : Numeric overflow occurred. Try increasing LI_MAX in the source code.";
        cout << endl << endl;
        return 0;
    }
    catch(ErrorDivisionByZero&)
    {
        cout << endl << endl;
        cout << "ERROR : Division by zero.";
        cout << endl << endl;
        return 0;
    }
    catch(bad_alloc&)
    {
        cout << endl << endl;
        cout << "ERROR : Allocating memory failed.";
        cout << endl << endl;
        return 0;
    }
    catch(exception& e)
    {
        cout << endl << endl;
        cout << "ERROR : An exception occurred." << endl;
        cout << e.what();
        cout << endl << endl;
        return 0;
    }
    catch(...)
    {
        cout << endl << endl;
        cout << "ERROR : An unknown error occurred.";
        cout << endl << endl;
        return 0;
    }

    Cleanup();

    cout << "--- Finished ---" << endl << endl;

    return 0;
}