/*
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.
*/

//Name        : VertexInducedHypercubeDensityChecker.
//
//Author      : Rahil Baber.
//
//Date        : 10th June 2013.
//
//Usage       : <executable> <input file>
//              E.g. VertexInducedHypercubeDensityChecker.exe 3-1.txt 
//
//Description : Using Razborov's method (extended to hypercubes) it calculates a bound on the
//              maximum density of induced hypercubes. The program expects the relevant data
//              (flags, basis, matrices, etc.) to be passed to it through the input file
//              (specified through the command line). This program is designed only to verify
//              that the input file gives the claimed bound, it does not create such input
//              files.
//
//Notes       : Currently the program can only handle hypercubes of dimension at most 5.
//              Users may wish/need to change the hard coded values 'LI_MAX' and
//              'noOfDecimalPlaces' specified below.

#include <iostream>
#include <fstream>
#include <stdlib.h>
using namespace std;

//*** Hard coded values that a user may wish to change. ***

//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.
const long LI_MAX = 1000; //You may need to change this value.

//Used by DisplayDensity() to calculate the upper bound of the Turan density to
//'noOfDecimalPlaces' decimal places.
const long noOfDecimalPlaces = 6;

//*** Class definitions. ***

//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);

//A simple class to store hypercubes that have at most 32 vertices. We take our vertex set to
//be {0,1,...,2^Dim-1}, where Dim is the dimension of the hypercube. The edge set is
//implicitly defined as follows: ij is an edge if and only if the binary representations of i
//and j differ by precisely one digit. The colours of the vetices are stored in the array
//C[]. We can represent a flag as (g, n) a CGraph g and an integer n, by
//taking the labelled vertices to be 0,1,...,2^n-1.
class CGraph
{
public:
    //Dimension of the hypercube.
    long Dim;

    //The colours of the vertices.
    //C[i] = true  if vertex i is blue,
    //       false if vertex i is red.
    bool C[32];
};

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

//Used by CLinkedList to create a dynamically allocated list of graphs.
class CListNode
{
public:
    //The graph which is part of the list.
    CGraph g;

    //CLinkedList is primarily used to store the family $\mathcal{H}$ (the admissible
    //hypercubes of dimension "dimOfH"). We will use coeff to keep track of the coefficient
    //of $p(H;G)$.
    CLargeFraction coeff;

    //Pointer to the next CListNode in the list.
    //If there is no next CListNode, pNext = NULL.
    CListNode* pNext;
};

//A very basic linked list of graphs. CLinkedList is primarily used to store the family
//$\mathcal{H}$ (the admissible hypercubes of dimension "dimOfH"). It may throw a bad_alloc
//object.
class CLinkedList
{
public:
    //The constructor. Sets pHead=NULL to indicate the list is empty.
    CLinkedList();

    //The destructor. Calls Free().
    ~CLinkedList();
    
    //Used to destroy the list, and free any memory it has allocated.
    void Free();

    //Used to add a new graph to the start of the list. It will throw a bad_alloc object if
    //there is not enough memory to add the new graph.
    void Add(CGraph& g);
    
    //Returns the number of nodes in the list.
    long Size();

    //Pointer to the first element in the list. If the list is empty pHead = NULL.
    CListNode* pHead;
};

//A very basic class that holds data used to calculate the upper bound of the Turan density.
//The calculation involves summing positive semidefinite matrices multiplied by flags. Each
//CTermInfo object holds one such matrix, the basis it is written in, the flags which are
//used to describe the basis, and the common induced type of the flags. Memory is dynamically
//allocated for the matrix, basis, and flags by ReadDataFile(). Any dynamically allocated
//memory is automatically freed by the CTermInfo's destructor. The class uses CLargeFraction
//objects and so ErrorOverflow and ErrorDivisionByZero objects may be thrown.
class CTermInfo
{
public:
    //The default constructor, initializes pointers to NULL.
    CTermInfo();

    //The destructor. A wrapper for Free().
    ~CTermInfo();

    //Deallocates memory.
    void Free();

    //The member variables are all initialized (and where appropriate allocated) by the
    //function ReadDataFile().

    //The common induced type of the flags. 
    CGraph Type;

    //FlagList[] holds a list of flags, and is of size noOfFlags. 
    long    noOfFlags;
    CGraph* FlagList;

    //Matrix[][] holds the positive semidefinite matrix and is of size MatrixDim by
    //MatrixDim.
    long             MatrixDim;
    CLargeFraction** Matrix;

    //Basis[][] represents a list of linear combinations of flags by storing the coeffcients
    //of the flags. The n-th linear combination is:
    //Basis[n-1][0](FlagList[0])+Basis[n-1][1](FlagList[1])+...
    //There are MatrixDim linear combinations (hence the first index varies from 0 to
    //MatrixDim-1). The second index varies from 0 to noOfFlags-1.
    CLargeFraction** Basis;
};

//*** Global variables. ***

//The dimension of the hypercubes $H_i$. Due to the way the function GenerateH() has been
//implemented dimOfH can be at most 5. In order for the density calculation to be valid
//dimOfH must be at least that of the induced hypercubes. It is initialized by ReadDataFile().
long dimOfH;

//H is a linked list representing $\mathcal{H}$ the family of admissible hypercubes of
//dimension dimOfH. H is initialized by a call to GenerateH().
CLinkedList H;

//InducedGraph[] holds the family $\mathcal{F}$. Its size is held in
//noOfInducedGraphs. noOfInducedGraphs is initialized in ReadDataFile(). ReadDataFile()
//also allocates the memory for InducedGraph[] then initializes it. The memory is
//deallocated in Cleanup().
long    noOfInducedGraphs;
CGraph* InducedGraph = NULL;

//Term[] is an array of size noOfTerms that holds data used to calculate the upper bound of
//the Turan density. noOfTerms is initialized in ReadDataFile(). ReadDataFile() also
//allocates the memory for Term[] then initializes it. The memory is deallocated in
//Cleanup().
long       noOfTerms;
CTermInfo* Term = NULL;

//Factorial[n] = n!.
long Factorial[6] = {1, 1, 2, 6, 24, 120};

//Iso[][] (together with Factorial[]) will be useful in generating isomorphic copies of
//graphs. Iso[][] is an array (dimOfH! by dimOfH) of all permutations of {0,1,...,dimOfH-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,...,dimOfH-1} are left in place. InitIso() allocates memory for Iso[][] and initializes
//it. Its memory is deallocated by Cleanup().
long** Iso = NULL;

//*** Function declarations. ***

//Deallocates the memory allocated to H, InducedGraph[], Term[], and Iso[][]. To deallocate
//the memory associated with Iso[][] it requires dimOfH to be set.
void Cleanup();

//Allocates memory for Iso[][] and initializes it. The function assumes dimOfH has been set
//to a value between 0 and 5. Iso[][] is initialized to a dimOfH! by dimOfH array
//representing all permutations of {0,1,...,dimOfH-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,...,dimOfH-1} are left in
//place. A bad_alloc object is thrown if memory allocation fails.
void InitIso();

//Calculates the number of characters needed to display 'n'. The function is used by
//ReadDataFile().
long NumberOfDigits(unsigned long n);

//This function is used in parsing the 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);

//This function is used in parsing the input file. It checks whether the character is a
//whitespace character (i.e. an end of line character, tab, or space).
bool IsWhitespace(char c);

//This function is used in parsing the 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);

//This function is used in parsing the input file. It continually removes the character
//'repeated' until a different character is encountered. 'repeated' must appear at least once
//otherwise the function 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 ReadRepeatingChar(istream& in, char repeated, char& c);

//This function is used in parsing the input file. Checks whether 'str[]' matches the
//characters in the istream. Ignores capitalization. Allows whitespaces in 'str[]' to be
//represented by one or more whitespaces in the istream. 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 CheckString(istream& in, const char* str, char& c);

//This function is used in parsing the input file. It removes any whitespaces then reads in a
//boolean value and returns it via the variable 'b'. The function accepts "yes", "no", "y",
//"n", "true", or "false" as the boolean value (capitalization is ignored). If one of the
//listed boolean values is not encountered the function returns false, otherwise it returns
//true. 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 ReadBool(istream& in, bool& b, char& c);

//This function is used in parsing the 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);

//This function is used in parsing the input file. It removes any whitespaces then tries to
//read in "<str[]>:<integer>", e.g. "Dimension:3". The string 'str[]' is read in using
//CheckString(). The integer is read using ReadLargeInteger() and returned via 'output'. It
//allows optional whitespaces either side of the colon. It may throw an ErrorOverflow object
//if the integer is too large.
//The function returns 0 on success,
//                     1 if the input is unexpected,
//                     2 if the integer is too large to be stored in a long.
//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.
long ReadLong(istream& in, const char* str, long& output, char& c);

//This function is used in parsing the input file. It removes any whitespaces then tries to
//read in the blue coloured vertices of a hypercube. The function assumes g.Dim has been
//set and uses g.Dim to verify the vertices are valid. The vertex set should take the form
//"{<vertex>,<vertex>,...,<vertex>}", where <vertex> is an integer representing the vertex.
//For example "{0, 2, 7}" or the empty set "{}". The integer representing the vertex should
//be between 0 and 2^g.Dim-1 (inclusive). Whitespaces either side of <vertex> are ignored. 
//The function returns true if it succeeds, and false otherwise. 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 ReadHypercube(istream& in, CGraph& g, char& c);

//This function is used in parsing the input file. It removes any whitespaces then tries to
//read in a linear combination of flags and store it in Term.Basis[basisIndex][]. The
//function assumes memory has already been allocated for Term.Basis[][]. It also assumes
//Term.noOfFlags has been set and uses this to verify the flags read in are valid. The linear
//combination of flags should take the form
//"<Fraction>(<Flag>)<+ or -><Fraction>(<Flag>)<+ or ->...".
//For example "1/2(F1)-4/3(F2)+2/3(F5)+1/7(F6)". An ErrorOverflow object may be thrown if an
//integer that is read is too large to be stored.
//The function returns 0 if it succeeds,
//                     1 if the input is incorrectly formatted,
//                     2 if a flag index does not lie in [1,Term.noOfFlags],
//                     3 if a coefficient causes a division by zero.
//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.
long ReadBasis(istream& in, CTermInfo& Term, long basisIndex, char& c);

//This function parses the input file given by filename. It initializes, and where
//appropriate allocates memory, for the global variables dimOfH, noOfInducedGraphs,
//InducedGraph[], noOfTerms, and Term[]. It returns true if it succeeds and false if an error
//occurred. It makes use of the functions NumberOfDigits(), ReadChar(), IsWhitespace(),
//RemoveWhitespace(), ReadRepeatingChar(), CheckString(), ReadBool(), ReadLargeInteger(),
//ReadLong(), ReadHypercube(), and ReadBasis().
bool ReadDataFile(const char* filename);

//Returns true if all matrices in Term[] are positive semidefinite (psd). The function
//assumes Term[] and noOfTerms have been initialized. A bad_alloc or ErrorOverflow object may
//be thrown.
bool IsPSD();

//Returns the number of blue vertices in the hypercube g.
long NumberOfVertices(CGraph& g);

//Returns true if g1 and g2 are of the same dimension and have identically coloured vertices.
bool IsCopy(CGraph& g1, CGraph& g2);

//If bInduced = false, IsSubcube() returns true if gSub is a subcube of gSuper.
//If bInduced = true,  IsSubcube() returns true if gSub is an induced subcube of gSuper.
//Assumes Iso[][] has been initialized.
bool IsSubcube(CGraph& gSub, CGraph& gSuper, bool bInduced);

//Each blue coloured vertex set can be assigned a unique integer. FindMinIso() creates gOut a
//graph with the property that the flag (gOut, dimOfType) is isomorphic to the flag 
//(gIn, dimOfType) and has the smallest integer associated with its blue coloured vertex set
//out of all such graphs. This function is used by GenerateH() to avoid storing isomorphic
//copies of graphs, and by CalculateDensity() to speed up searching through lists of flags.
//In order for GenerateH() to function correctly the way FindMinIso() constructs the integer
//representing the blue vertex set is important. The function assumes Iso[][] has been
//initialized.
void FindMinIso(CGraph gIn, CGraph& gOut, long dimOfType);

//Generates 'H' a linked list representing $\mathcal{H}$ the family of admissible hypercubes
//of dimension dimOfH. The function assumes dimOfH has been initialized. A bad_alloc object
//may be thrown.
void GenerateH();

//Displays 'Density' as both a fraction and a decimal. This function is used by
//CalculateDensity() to display an upper bound of the induced density. It may throw an
//ErrorOverflow object.
void DisplayDensity(CLargeFraction Density);

//Calculates then displays the upper bound of the induced density. The function assumes
//Iso[][], Term[], noOfTerms, noOfInducedGraphs, InducedGraph[], H, and dimOfH have been
//initialized. To save memory the function modifies Term[]. A bad_alloc or ErrorOverflow
//object may be thrown.
void CalculateDensity();

//The main function. Expects argc to be 2, argv[0] to be the name of the executable and
//argv[1] to be the name of the input file.
int main(int argc, char* argv[]);

//*** Function definitions. ***

//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  = (1<<(8*sizeof(long)-2))-1;
        max += 1<<(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;

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

    //Set Q to be the magnitude of L.
    if(L<0)
    {
        out << "-"; //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--)
        out << char('0'+digit[i]);

    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)
{
    out << F.N << "/" << F.D;
    return out;
}

//Allows us to use << operator to output the CGraph object g to the screen or a file.
ostream& operator<<(ostream& out, const CGraph& g)
{
    long i;
    bool bFirst;

    //Display the number of vertices.
    out << "Dimension " << g.Dim << " : ";

    //Display the blue vertices.
    out << "Blue vertices = {";

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

    for(i=0;i<(1<<g.Dim);i++)
        if(g.C[i]==true)
        {
            //Display the commas that separate the vertices.
            if(bFirst==false)
                out << ", ";
            else
                bFirst = false;
                
            //Display the vertex i.
            out << i;
        }

    out << "}";

    return out;
}

//The constructor. Sets pHead=NULL to indicate the list is empty.
CLinkedList::CLinkedList()
{
    pHead = NULL;
    return;
}

//The destructor. Calls Free().
CLinkedList::~CLinkedList()
{
    Free();
    return;
}

//Used to destroy the list, and free any memory it has allocated.
void CLinkedList::Free()
{
    CListNode* pTemp;

    while(pHead!=NULL)
    {
        //Delete the first element in the list and set the second element as the new first
        //element.
        pTemp = pHead->pNext;
        delete pHead;
        pHead = pTemp;
    }

    return;
}

//Used to add a new graph to the start of the list. It will throw a bad_alloc object if there
//is not enough memory to add the new graph.
void CLinkedList::Add(CGraph& g)
{
    CListNode* pNew;

    pNew = new CListNode;

    pNew->g     = g;
    pNew->coeff = 0;
    pNew->pNext = pHead;
    pHead       = pNew; 

    return;
}

//Returns the number of nodes in the list.
long CLinkedList::Size()
{
    long       size;
    CListNode* pNode;

    size  = 0;
    pNode = pHead;
    while(pNode!=NULL)
    {
        size++;
        pNode = pNode->pNext;
    }

    return size;
}

//The default constructor, initializes pointers to NULL.
CTermInfo::CTermInfo()
{
    FlagList = NULL;
    Matrix   = NULL;
    Basis    = NULL;

    return;
}

//The destructor. A wrapper for Free().
CTermInfo::~CTermInfo()
{
    Free();
    return;
}

//Deallocates memory.
void CTermInfo::Free()
{
    long i;

    delete[] FlagList;
    FlagList = NULL;

    if(Matrix!=NULL)
    {
        for(i=0;i<MatrixDim;i++)
            delete[] Matrix[i];

        delete[] Matrix;
        Matrix = NULL;
    }

    if(Basis!=NULL)
    {
        for(i=0;i<MatrixDim;i++)
            delete[] Basis[i];

        delete[] Basis;
        Basis = NULL;
    }

    return;
}

//Deallocates the memory allocated to H, InducedGraph[], Term[], and Iso[][]. To deallocate
//the memory associated with Iso[][] it requires dimOfH to be set.
void Cleanup()
{
    long i;

    H.Free();

    delete[] InducedGraph;
    InducedGraph = NULL;

    delete[] Term;
    Term = NULL;

    if(Iso!=NULL)
    {
        for(i=0;i<Factorial[dimOfH];i++)
            delete[] Iso[i];

        delete[] Iso;
        Iso = NULL;
    }
    
    return;
}

//Allocates memory for Iso[][] and initializes it. The function assumes dimOfH has been set
//to a value between 0 and 5. Iso[][] is initialized to a dimOfH! by dimOfH array
//representing all permutations of {0,1,...,dimOfH-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,...,dimOfH-1} are left in
//place. A bad_alloc object is thrown if memory allocation fails.
void InitIso()
{
    long i, j, k, p, n;

    //Allocate memory for Iso[][].
    Iso = new long*[Factorial[dimOfH]];

    for(i=0;i<Factorial[dimOfH];i++)
        Iso[i] = NULL;

    if(dimOfH==0)
        return;

    for(i=0;i<Factorial[dimOfH];i++)
        Iso[i] = new long[dimOfH];

    //Initialize Iso[][].

    //Set Iso[0][] to be the identity map.
    for(n=0;n<dimOfH;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,...,dimOfH-1}. Generate new permutations by swapping i and j in Iso[k][], for
    //all j<i and k<=i!.
    for(i=1;i<dimOfH;i++)
    for(j=0;j<i;j++)
    for(k=0;k<Factorial[i];k++)
    {
        for(n=0;n<dimOfH;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;
}

//Calculates the number of characters needed to display 'n'. The function is used by
//ReadDataFile().
long NumberOfDigits(unsigned long n)
{
    long result = 0;

    if(n==0)
        return 1;

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

    return result;
}

//This function is used in parsing the 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 the input file. It checks whether the character is a
//whitespace character (i.e. an end of line character, tab, or space).
bool IsWhitespace(char c)
{
    if(c==' ' || c=='\t' || c=='\n' || c=='\r')
        return true;

    return false;
}

//This function is used in parsing the 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(IsWhitespace(c)==true)
        c = ReadChar(in);

    return;
}

//This function is used in parsing the input file. It continually removes the character
//'repeated' until a different character is encountered. 'repeated' must appear at least once
//otherwise the function 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 ReadRepeatingChar(istream& in, char repeated, char& c)
{
    if(c!=repeated)
        return false;

    while(c==repeated)
        c = ReadChar(in);

    return true;
}

//This function is used in parsing the input file. Checks whether 'str[]' matches the
//characters in the istream. Ignores capitalization. Allows whitespaces in 'str[]' to be
//represented by one or more whitespaces in the istream. 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 CheckString(istream& in, const char* str, char& c)
{
    long i;
    char s1, s2;

    i = 0;
    while(str[i]!='\0')
    {
        if(IsWhitespace(str[i])==true)
        {
            if(IsWhitespace(c)==false)
                return false;

            //Move to a non-whitespace character in the istream and str[].
            RemoveWhitespace(in,c);

            while(IsWhitespace(str[i])==true)
                i++;
        }
        else
        {
            //Ignore capitilizations.
            s1 = str[i];
            s2 = c;

            if('A'<=s1 && s1<='Z')
                s1 += 'a'-'A';

            if('A'<=s2 && s2<='Z')
                s2 += 'a'-'A';

            //Check if str[] matches the istream.
            if(s1!=s2)
                return false;

            c = ReadChar(in);
            i++;
        }
    }

    return true;
}

//This function is used in parsing the input file. It removes any whitespaces then reads in a
//boolean value and returns it via the variable 'b'. The function accepts "yes", "no", "y",
//"n", "true", or "false" as the boolean value (capitalization is ignored). If one of the
//listed boolean values is not encountered the function returns false, otherwise it returns
//true. 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 ReadBool(istream& in, bool& b, char& c)
{
    b = false;

    RemoveWhitespace(in,c);

    //Check if input is "true".
    if(c=='t' || c=='T')
    {
        b = true;
        return CheckString(in,"true",c);
    }

    //Check if input is "false".
    if(c=='f' || c=='F')
    {
        b = false;
        return CheckString(in,"false",c);
    }

    //Check if input is "y" or "yes".
    if(c=='y' || c=='Y')
    {
        b = true;

        //Check if it is "yes".
        c = ReadChar(in);

        if(c=='e' || c=='E')
            return CheckString(in,"es",c);

        return true;
    }

    //Check if input is "n" or "no".
    if(c=='n' || c=='N')
    {
        b = false;

        //Check if it is "no".
        c = ReadChar(in);

        if(c=='o' || c=='O')
            c = ReadChar(in);

        return true;
    }

    //An unexpected character was encountered.
    return false;
}

//This function is used in parsing the 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;
}

//This function is used in parsing the input file. It removes any whitespaces then tries to
//read in "<str[]>:<integer>", e.g. "Dimension:3". The string 'str[]' is read in using
//CheckString(). The integer is read using ReadLargeInteger() and returned via 'output'. It
//allows optional whitespaces either side of the colon. It may throw an ErrorOverflow object
//if the integer is too large.
//The function returns 0 on success,
//                     1 if the input is unexpected,
//                     2 if the integer is too large to be stored in a long.
//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.
long ReadLong(istream& in, const char* str, long& output, char& c)
{
    CLargeInt L;

    output = 0;
    
    RemoveWhitespace(in,c);

    if(CheckString(in,str,c)==false)
        return 1; //Unexpected input.

    RemoveWhitespace(in,c);
        
    if(c!=':')
        return 1; //Unexpected input.

    c = ReadChar(in);
    RemoveWhitespace(in,c);

    if(ReadLargeInteger(in,L,c)==false)
        return 1; //Unexpected input.

    try
    {
        output = L.ConvertToLong();
    }
    catch(ErrorOverflow& e)
    {
        return 2; //L is too large for 'output'.
    }

    return 0; //Success.
}

//This function is used in parsing the input file. It removes any whitespaces then tries to
//read in the blue coloured vertices of a hypercube. The function assumes g.Dim has been
//set and uses g.Dim to verify the vertices are valid. The vertex set should take the form
//"{<vertex>,<vertex>,...,<vertex>}", where <vertex> is an integer representing the vertex.
//For example "{0, 2, 7}" or the empty set "{}". The integer representing the vertex should
//be between 0 and 2^g.Dim-1 (inclusive). Whitespaces either side of <vertex> are ignored. 
//The function returns true if it succeeds, and false otherwise. 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 ReadHypercube(istream& in, CGraph& g, char& c)
{
    long i, j;

    //Initialize the vertex colours.
    for(i=0;i<32;i++)
        g.C[i] = false;
    
    //The first non-whitespace character should be '{'.
    RemoveWhitespace(in,c);

    if(c!='{')
        return false;

    c = ReadChar(in);
    RemoveWhitespace(in,c);

    //Check for empty edge set
    if(c=='}')
    {
        c = ReadChar(in);
        return true;
    }

    while(true)
    {
        //Read vertex.
        i = -1;
        while(true)
        {
            if('0'<=c && c<='9')
            {
                if(i==-1)
                    i = 0;
                
                i = 10*i+c-'0';
                
                c = ReadChar(in);
            }
            else
                break;
        }
        
        if(i==-1)
            return false;
                
        //Check i is a valid vertex.
        if(i>=(1<<g.Dim))
            return false;
        
        //Set the vertex.
        g.C[i] = true;

        //The next non-whitespace character should be '}' or ','.
        RemoveWhitespace(in,c);
        
        if(c=='}')
        {
            c = ReadChar(in);
            return true;
        }

        if(c!=',')
            return false;

        c = ReadChar(in);
        RemoveWhitespace(in,c);
    }
}

//This function is used in parsing the input file. It removes any whitespaces then tries to
//read in a linear combination of flags and store it in Term.Basis[basisIndex][]. The
//function assumes memory has already been allocated for Term.Basis[][]. It also assumes
//Term.noOfFlags has been set and uses this to verify the flags read in are valid. The linear
//combination of flags should take the form
//"<Fraction>(<Flag>)<+ or -><Fraction>(<Flag>)<+ or ->...".
//For example "1/2(F1)-4/3(F2)+2/3(F5)+1/7(F6)". An ErrorOverflow object may be thrown if an
//integer that is read is too large to be stored.
//The function returns 0 if it succeeds,
//                     1 if the input is incorrectly formatted,
//                     2 if a flag index does not lie in [1,Term.noOfFlags],
//                     3 if a coefficient causes a division by zero.
//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.
long ReadBasis(istream& in, CTermInfo& Term, long basisIndex, char& c)
{
    long           i;
    bool           bFirstTerm;
    CLargeInt      L;
    CLargeFraction Frac;
    long           index;
    
    //Zero out basis.
    for(i=0;i<Term.noOfFlags;i++)
        Term.Basis[basisIndex][i] = 0;

    //bFirstTerm keeps track of whether we are reading in the first term of the linear
    //combination.
    bFirstTerm = true;
    while(true)
    {
        RemoveWhitespace(in,c);

        //Check if the linear combination has ended.
        if(bFirstTerm==true)
        {
            if(c!='+' && c!='-' && c!='(' && (c<'0' || '9'<c))
                break; //Linear combination is "empty".
        }
        else
        {
            //bFirstTerm = false, hence we have read something that looks like "...(F2)". The
            //next character should be + or - if there is another term we have yet to read.
            if(c!='+' && c!='-')
                break;
        }

        //Assume the next term has a coefficient of 1.
        Frac = 1;

        if(bFirstTerm==false)
        {
            //Read the + or - that separates the terms.
            if(c=='-')
                Frac = -Frac;

            c = ReadChar(in);
            RemoveWhitespace(in,c);
        }

        //The term starts with +, -, (, or 0-9.

        if(c=='+' || c=='-')
        {
            if(c=='-')
                Frac = -Frac;

            c = ReadChar(in);
            RemoveWhitespace(in,c);

            //The next character could be ( or 0-9. E.g. "-(F1)...", or "-3(F1)...".
            if(c!='(' && (c<'0' || '9'<c))
                return 1; //Incorrect format.
        }

        if('0'<=c && c<='9')
        {
            //Read the integer coefficient.
            if(ReadLargeInteger(in,L,c)==false)
                return 1; //Incorrect format.

            Frac = Frac*L;

            RemoveWhitespace(in,c);

            //The next character could be / or (. E.g. "3(F1)..." or "3/2(F1)...".
            if(c=='/')
            {
                c = ReadChar(in);
                RemoveWhitespace(in,c);

                if(ReadLargeInteger(in,L,c)==false)
                    return 1; //Incorrect format.

                if(L==0)
                    return 3; //Division by zero.

                Frac = Frac/L;

                RemoveWhitespace(in,c);
            }

            //The next character should be (.
        }

        if(c!='(')
            return 1; //Incorrect format.

        c = ReadChar(in);
        RemoveWhitespace(in,c);

        //Read the flag.
        if(CheckString(in,"F",c)==false)
            return 1; //Incorrect format.

        if(c=='+')
            return 1; //Incorrect format.

        if(ReadLargeInteger(in,L,c)==false)
            return 1; //Incorrect format.

        if(L<1 || Term.noOfFlags<L)
            return 2; //Flag index was out of bounds.

        index = L.ConvertToLong();

        Term.Basis[basisIndex][index-1] = Term.Basis[basisIndex][index-1]+Frac;

        RemoveWhitespace(in,c);

        if(c!=')')
            return 1; //Incorrect format.

        c = ReadChar(in);

        bFirstTerm = false;
    }

    return 0; //Finished successfully.
}

//This function parses the input file given by filename. It initializes, and where
//appropriate allocates memory, for the global variables dimOfH, noOfInducedGraphs,
//InducedGraph[], noOfTerms, and Term[]. It returns true if it succeeds and false if an error
//occurred. It makes use of the functions NumberOfDigits(), ReadChar(), IsWhitespace(),
//RemoveWhitespace(), ReadRepeatingChar(), CheckString(), ReadBool(), ReadLargeInteger(),
//ReadLong(), ReadHypercube(), and ReadBasis().
bool ReadDataFile(const char* filename)
{
    long      i, n, r;
    long      row, col;
    long      nTerm, nFlag, nBasis;
    ifstream  FileIn;
    char      c;
    CLargeInt L;

    cout << "Reading the data file....." << endl;
    
    FileIn.open(filename);

    if(FileIn.good()==false)
    {
        cout << endl;
        cout << "*** ERROR: Cannot open the file. ***" << endl;
        cout << endl;
        return false;
    }
    
    try
    {
        c = ReadChar(FileIn);

        //Check Program.
        {
            RemoveWhitespace(FileIn,c);

            if(CheckString(FileIn,"Program",c)==false)
                throw "0:Expected to read \"Program : VertexInducedHypercubeDensityChecker\".";


            RemoveWhitespace(FileIn,c);
        
            if(c!=':')
                throw "0:Expected to read \"Program : VertexInducedHypercubeDensityChecker\".";

            c = ReadChar(FileIn);
            RemoveWhitespace(FileIn,c);
        
            if(CheckString(FileIn,"VertexInducedHypercubeDensityChecker",c)==false)
                throw "0:Expected to read \"Program : VertexInducedHypercubeDensityChecker\".";

            cout << endl;
            cout << "Program                      : VertexInducedHypercubeDensityChecker" << endl;
        }

        //Read dimOfH.
        {
            r = ReadLong(FileIn,"Dimension of subcubes H",dimOfH,c);
            
            if(r!=0 || dimOfH<0 || 5<dimOfH)
                throw "0:Expected to read \"Dimension of subcubes H : "
                      "<An integer in [0,5]>\".";
            
            cout << "Dimension of subcubes H      : " << dimOfH << endl;
        }

        //Read noOfInducedGraphs.
        {
            r = ReadLong(FileIn,"Number of induced hypercubes",noOfInducedGraphs,c);
            
            if(r==2)
                throw "1:";

            if(r!=0 || noOfInducedGraphs<1)
                throw "0:Expected to read \"Number of induced hypercubes : "
                      "<A positive integer>\".";

            cout << "Number of induced hypercubes : " << noOfInducedGraphs << endl;
            
            InducedGraph = new CGraph[noOfInducedGraphs];
        }

        //Read InducedGraph[].
        {
            //Read "Induced hypercubes :".
            RemoveWhitespace(FileIn,c);

            if(CheckString(FileIn,"Induced hypercubes",c)==false)
                throw "0:Expected to read \"Induced hypercubes : "
                      "<A list of hypercubes>\".";

            RemoveWhitespace(FileIn,c);
        
            if(c!=':')
                throw "0:Expected to read \"Induced hypercubes : "
                      "<A list of hypercubes>\".";

            c = ReadChar(FileIn);

            cout << "Induced hypercubes           :" << endl;

            //Read the graphs.
            for(i=0;i<noOfInducedGraphs;i++)
            {
                RemoveWhitespace(FileIn,c);

                //Ignore the character if it is ','.
                if(i!=0 && c==',')
                    c = ReadChar(FileIn);

                RemoveWhitespace(FileIn,c);

                //Read the dimension.
                if(ReadLargeInteger(FileIn,L,c)==false)
                    throw "2:";

                try
                {
                    InducedGraph[i].Dim = L.ConvertToLong();
                }
                catch(ErrorOverflow& e)
                {
                    throw "2:";
                }

                if(InducedGraph[i].Dim>5 || InducedGraph[i].Dim<0)
                    throw "2:";

                if(InducedGraph[i].Dim>dimOfH)
                    throw "0:Dimension of the induced hypercube is too high.";

                RemoveWhitespace(FileIn,c);
        
                if(c!=':')
                    throw "2:";

                //Read the edge set.
                c = ReadChar(FileIn);

                if(ReadHypercube(FileIn,InducedGraph[i],c)==false)
                    throw "2:";
                
                cout << InducedGraph[i] << endl;
            }
        }

        //Read noOfTerms.
        {
            r = ReadLong(FileIn,"Number of terms",noOfTerms,c);
            
            if((r==0 && noOfTerms<0) || r==1)
                throw "0:Expected to read \"Number of terms : <A non-negative integer>\".";

            if(r==2)
                throw "1:";

            cout << "Number of terms              : " << noOfTerms << endl;
            
            if(noOfTerms==0)
            {
                cout << endl;
                cout << "Finished reading the data file." << endl;
                cout << endl;

                FileIn.close();
                return true;
            }

            //Allocate memory for the term data.
            Term = new CTermInfo[noOfTerms];
        }

        cout << endl;

        try
        {
            //Read the Term[] data.
            for(nTerm=1;nTerm<=noOfTerms;nTerm++)
            {
                //Read the header e.g. "--- Term 1 ---".
                {
                    cout << "\rReading term " << nTerm << " : ";
                    cout << "Reading the header." << flush;

                    RemoveWhitespace(FileIn,c);

                    if(ReadRepeatingChar(FileIn,'-',c)==false)
                        throw "3:";
            
                    RemoveWhitespace(FileIn,c);

                    if(CheckString(FileIn,"Term",c)==false)
                        throw "3:";

                    RemoveWhitespace(FileIn,c);

                    if(ReadLargeInteger(FileIn,L,c)==false)
                        throw "3:";

                    if(L!=nTerm)
                        throw "3:";
                
                    RemoveWhitespace(FileIn,c);

                    if(ReadRepeatingChar(FileIn,'-',c)==false)
                        throw "3:";
                }

                //Read Term[nTerm-1].Type.Dim.
                {
                    cout << "\rReading term " << nTerm << " : ";
                    cout << "Reading the dimension of the type." << flush;

                    r = ReadLong(FileIn,"Dimension of type",Term[nTerm-1].Type.Dim,c);
                
                    if(r!=0 || Term[nTerm-1].Type.Dim<0 || dimOfH<Term[nTerm-1].Type.Dim)
                        throw "0:Expected to read \"Dimension of type : "
                              "<An integer in [0,<Dimension of subcubes H>]>\".";

                    //Initialize the vertex colours of the type.
                    for(i=0;i<32;i++)
                        Term[nTerm-1].Type.C[i] = false;
                }
                
                //Read the dimension of the flags.
                {
                    cout << "\rReading term " << nTerm << " : ";
                    cout << "Reading the dimension of the flags." << flush;

                    //Temporarily store the dimension in 'n'.
                    r = ReadLong(FileIn,"Dimension of flags",n,c);

                    if(r!=0 || 4<n || n<Term[nTerm-1].Type.Dim
                    || (dimOfH+Term[nTerm-1].Type.Dim)/2<n)
                        throw "0:Expected to read \"Dimension of flags : <An integer less\n"
                              "than 5 and in [<Dimension of type>, <Dimension of type +\n"
                              "Dimension of subcubes H>/2]>\".";
                }

                //Read Term[nTerm-1].noOfFlags.
                {
                    cout << "\rReading term " << nTerm << " : ";
                    cout << "Reading the number of flags.       " << flush;

                    r = ReadLong(FileIn,"Number of flags",Term[nTerm-1].noOfFlags,c);
                
                    if(r==2)
                        throw "1:";

                    if(r!=0 || Term[nTerm-1].noOfFlags<0)
                        throw "0:Expected to read \"Number of flags : "
                              "<A non-negative integer>\".";
                }

                //Allocate memory for Term[nTerm-1].FlagList[].
                {
                    cout << "\rReading term " << nTerm << " : ";
                    cout << "Allocating memory for the flags." << flush;

                    if(0<Term[nTerm-1].noOfFlags)
                        Term[nTerm-1].FlagList = new CGraph[Term[nTerm-1].noOfFlags];
                }

                //Read Term[nTerm-1].MatrixDim.
                {
                    cout << "\rReading term " << nTerm << " : ";
                    cout << "Reading the dimension of the matrix." << flush;

                    r = ReadLong(FileIn,"Matrix dimension",Term[nTerm-1].MatrixDim,c);

                    if(r==2)
                        throw "1:";

                    if(r!=0 || Term[nTerm-1].MatrixDim<0)
                        throw "0:Expected to read \"Matrix dimension : "
                              "<A non-negative integer>\".";

                    if(Term[nTerm-1].noOfFlags==0 && Term[nTerm-1].MatrixDim>0)
                        throw "0:Expected the dimension to be 0, "
                              "as the number of flags is 0.";
                }

                //Allocate memory for Term[nTerm-1].Matrix[][] and Term[nTerm-1].Basis[][].
                if(Term[nTerm-1].MatrixDim>0)
                {
                    cout << "\rReading term " << nTerm << " : ";
                    cout << "Allocating memory for the matrix.   " << flush;

                    Term[nTerm-1].Matrix = new CLargeFraction*[Term[nTerm-1].MatrixDim];

                    for(i=0;i<Term[nTerm-1].MatrixDim;i++)
                        Term[nTerm-1].Matrix[i] = NULL;

                    for(i=0;i<Term[nTerm-1].MatrixDim;i++)
                        Term[nTerm-1].Matrix[i] =
                            new CLargeFraction[Term[nTerm-1].MatrixDim];

                    cout << "\rReading term " << nTerm << " : ";
                    cout << "Allocating memory for the basis. " << flush;

                    Term[nTerm-1].Basis = new CLargeFraction*[Term[nTerm-1].MatrixDim];

                    for(i=0;i<Term[nTerm-1].MatrixDim;i++)
                        Term[nTerm-1].Basis[i] = NULL;

                    for(i=0;i<Term[nTerm-1].MatrixDim;i++)
                        Term[nTerm-1].Basis[i] = new CLargeFraction[Term[nTerm-1].noOfFlags];
                }

                //Read Term[nTerm-1].FlagList[].
                {
                    cout << "\rReading term " << nTerm << " : ";
                    cout << "Reading the flags.              " << flush;
                
                    RemoveWhitespace(FileIn,c);
                
                    if(CheckString(FileIn,"Flags",c)==false)
                        throw "0:Expected to read \"Flags :\".";

                    RemoveWhitespace(FileIn,c);

                    if(c!=':')
                        throw "0:Expected to read \"Flags :\".";

                    c = ReadChar(FileIn);

                    cout << "\rReading term " << nTerm << " : ";
                    cout << "                  ";

                    for(nFlag=1;nFlag<=Term[nTerm-1].noOfFlags;nFlag++)
                    {
                        
                        cout << "\rReading term " << nTerm << " : ";
                        cout << "Reading flag " << nFlag << "." << flush;

                        RemoveWhitespace(FileIn,c);

                        if(CheckString(FileIn,"F",c)==false)
                            throw "4:";

                        if(c=='+')
                            throw "4:";

                        if(ReadLargeInteger(FileIn,L,c)==false)
                            throw "4:";

                        if(L!=nFlag)
                            throw "4:";

                        RemoveWhitespace(FileIn,c);

                        if(c!=':')
                            throw "4:";

                        c = ReadChar(FileIn);

                        Term[nTerm-1].FlagList[nFlag-1].Dim = n;
                        if(ReadHypercube(FileIn,Term[nTerm-1].FlagList[nFlag-1],c)==false)
                            throw "4:";

                        //Set the type's colours set using the first flag's colours.
                        if(nFlag==1)
                        {
                            for(i=0;i<(1<<Term[nTerm-1].Type.Dim);i++)
                                Term[nTerm-1].Type.C[i] = Term[nTerm-1].FlagList[0].C[i];
                        }

                        //Check the flag has the same type and dimension as the other flags.
                        if(Term[nTerm-1].FlagList[nFlag-1].Dim!=
                            Term[nTerm-1].FlagList[0].Dim)
                            throw "0:The flag has a different dimension to that of \"F1\".";

                        for(i=0;i<(1<<Term[nTerm-1].Type.Dim);i++)
                        {
                            if(Term[nTerm-1].FlagList[nFlag-1].C[i]!=Term[nTerm-1].Type.C[i])
                                throw "0:The flag's induced \"type\" is not the same\n"
                                      "as those of the other flags. The labelled\n"
                                      "vertices that form the type are vertices\n"
                                      "0 to 2^<Dimension of type>-1.";
                        }
                    }

                    cout << "\rReading term " << nTerm << " : ";
                    cout << "              ";
                    for(i=0;i<NumberOfDigits(Term[nTerm-1].noOfFlags);i++)
                        cout << " ";
                }
                
                //Read Term[nTerm-1].Basis[][].
                {
                    cout << "\rReading term " << nTerm << " : ";
                    cout << "Reading the basis." << flush;
                
                    RemoveWhitespace(FileIn,c);

                    if(CheckString(FileIn,"Basis",c)==false)
                        throw "0:Expected to read \"Basis :\".";

                    RemoveWhitespace(FileIn,c);

                    if(c!=':')
                        throw "0:Expected to read \"Basis :\".";

                    c = ReadChar(FileIn);

                    for(nBasis=1;nBasis<=Term[nTerm-1].MatrixDim;nBasis++)
                    {
                        cout << "\rReading term " << nTerm << " : ";
                        cout << "Reading basis member " << nBasis << "." << flush;
            
                        RemoveWhitespace(FileIn,c);

                        if(CheckString(FileIn,"B",c)==false)
                            throw "5:";

                        if(c=='+')
                            throw "5:";

                        if(ReadLargeInteger(FileIn,L,c)==false)
                            throw "5:";

                        if(L!=nBasis)
                            throw "5:";

                        RemoveWhitespace(FileIn,c);

                        if(c!=':')
                            throw "5:";

                        c = ReadChar(FileIn);

                        r = ReadBasis(FileIn,Term[nTerm-1],nBasis-1,c);

                        if(r==1)
                            throw "5:"; //Incorrectly formatted input.

                        if(r==2)
                            throw "0:Expected to read \"F<An integer in "
                                  "[1,<Number of flags>]>\".";

                        if(r==3)
                            throw "0:A denominator of one of the fractions is zero.";
                    }

                    cout << "\rReading term " << nTerm << " : ";
                    cout << "                      ";
                    for(i=0;i<NumberOfDigits(Term[nTerm-1].MatrixDim);i++)
                        cout << " ";
                }

                //Read Term[nTerm-1].Matrix[][].
                {
                    cout << "\rReading term " << nTerm << " : ";
                    cout << "Reading the matrix." << flush;

                    RemoveWhitespace(FileIn,c);

                    if(CheckString(FileIn,"Matrix",c)==false)
                        throw "0:Expected to read \"Matrix :\".";

                    RemoveWhitespace(FileIn,c);

                    if(c!=':')
                        throw "0:Expected to read \"Matrix :\".";

                    c = ReadChar(FileIn);
                    
                    for(row=1;row<=Term[nTerm-1].MatrixDim;row++)
                    for(col=1;col<=Term[nTerm-1].MatrixDim;col++)
                    {
                        cout << "\rReading term " << nTerm << " : ";
                        cout << "Reading matrix entry (" << row << ", " << col << ").";
                        if(col==1 && row!=1)
                        {
                            n = NumberOfDigits(row-1)+NumberOfDigits(Term[nTerm-1].MatrixDim)
                                -NumberOfDigits(row)-1;

                            for(i=0;i<n;i++)
                                cout << " ";
                        }

                        cout << flush;

                        RemoveWhitespace(FileIn,c);

                        if(ReadLargeInteger(FileIn,L,c)==false)
                            throw "0:Expected to read a fraction.";

                        Term[nTerm-1].Matrix[row-1][col-1] = L;

                        RemoveWhitespace(FileIn,c);

                        //Check if it is a fraction.
                        if(c=='/')
                        {
                            c = ReadChar(FileIn);
                            RemoveWhitespace(FileIn,c);

                            if(ReadLargeInteger(FileIn,L,c)==false)
                                throw "0:Expected to read a fraction.";

                            if(L==0)
                                throw "0:The denominator of the fraction is zero.";

                            Term[nTerm-1].Matrix[row-1][col-1] =
                                Term[nTerm-1].Matrix[row-1][col-1]/L;
                        }
                    }

                    cout << "\rReading term " << nTerm << " : ";
                    cout << "                          ";
                    for(i=0;i<2*NumberOfDigits(Term[nTerm-1].MatrixDim);i++)
                        cout << " ";
                }
            }
        }
        catch(...)
        {
            cout << endl;
            throw;
        }

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

        cout << "\rFinished reading the data file." << endl;
        cout << endl;
        
        FileIn.close();
        return true;
    }
    catch(const char* str)
    {
        cout << endl;
        
        switch(str[0])
        {
        case '0':
            cout << "*** ERROR: Unexpected input. ***" << endl;
            cout << str+2 << endl;
            break;

        case '1':
            cout << "*** ERROR: Overflow. ***" << endl;
            cout << "Integer read was too large to store in a \"long int\"." << endl; 
            break;

        case '2':
            cout << "*** ERROR: Unexpected input. ***" << endl;
            cout << "Expected to read a vertex-coloured hypercube in the form" << endl;
            cout << "\"<The dimension - an integer in [0,5]> : <The set of blue vertices>\".\n";
            cout << "For example \"3 : {0, 2, 7}\"." << endl;
            break;

        case '3':
            cout << "*** ERROR: Unexpected input. ***" << endl;
            cout << "Expected to read \"--- Term " << nTerm << " ---\"." << endl;
            break;

        case '4':
            cout << "*** ERROR: Unexpected input. ***" << endl;
            cout << "Expected to read \"F" << nFlag << " : ";
            cout << "<The flag's blue vertex set>\"."<< endl;
            cout << "An example of a blue vertex set for a dimension 3" << endl;
            cout << "hypercube is \"{0, 2, 7}\"." << endl; 
            break;

        case '5':
            cout << "*** ERROR: Unexpected input. ***" << endl;
            cout << "Expected to read \"B" << nBasis;
            cout << " : <A linear combination of flags>\"." << endl;
            cout << "The linear combination of flags should be in the form" << endl;
            cout << "\"<Fraction>(<Flag>)<+ or -><Fraction>(<Flag>)<+ or ->...\"." << endl;
            cout << "For example \"1/2(F1)-4/3(F2)+2/3(F5)+1/7(F6)\"." << endl;
            break;

        default:
            cout << "*** ERROR: Unknown. ***" << endl;
            cout << "Thrown data : \"" << str << "\"" << endl; 
        }
    }
    catch(ErrorOverflow& e)
    {
        cout << endl;
        cout << "*** ERROR: Overflow. ***" << endl;
        cout << "Integer read was too large to store." << endl;
        cout << "Try increasing LI_MAX in the source code." << endl;        
    }
    catch(bad_alloc& e)
    {
        cout << endl;
        cout << "*** ERROR: Bad allocation. ***" << endl;
        cout << "Probably caused by lack of memory." << endl;
        cout << "Try decreasing LI_MAX in the source code." << endl;
    }
    catch(...)
    {
        cout << endl;
        cout << "*** ERROR: Unknown. ***" << endl;
    }
    
    cout << endl;

    FileIn.close();
    return false;
}

//Returns true if all matrices in Term[] are positive semidefinite (psd). The function
//assumes Term[] and noOfTerms have been initialized. A bad_alloc or ErrorOverflow object may
//be thrown.
bool IsPSD()
{
    long             i, j, n;
    long             nTerm;
    CLargeFraction** M;
    long             dim;
    CLargeFraction   F;
    CTermInfo        Temp;

    for(nTerm=1;nTerm<=noOfTerms;nTerm++)
    {
        cout << "\rChecking matrix " << nTerm << " of " << noOfTerms << " : ";
        cout << "Freeing memory.   " << flush;

        Temp.Free();

        cout << "\rChecking matrix " << nTerm << " of " << noOfTerms << " : ";
        cout << "Allocating memory." << flush;

        if(Term[nTerm-1].MatrixDim<=0)
            continue;

        //Copy Term[nTerm-1].Matrix[][] to Temp.Matrix[][]. Use dim and M to simplify the
        //expressions.
        dim            = Term[nTerm-1].MatrixDim;
        Temp.MatrixDim = dim;

        M           = new CLargeFraction*[dim];
        Temp.Matrix = M;

        for(i=0;i<dim;i++)
            M[i] = NULL;

        for(i=0;i<dim;i++)
            M[i] = new CLargeFraction[dim];

        cout << "\rChecking matrix " << nTerm << " of " << noOfTerms << " : ";
        cout << "Copying matrix.   " << flush;

        for(i=0;i<dim;i++)
        for(j=0;j<dim;j++)
            M[i][j] = Term[nTerm-1].Matrix[i][j];

        //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)
                return false; //Not positive semidefinite.

            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[][].
                        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].
            F = M[n][n];
            for(i=n;i<dim;i++)
                M[n][i] = M[n][i]/F;

            //Subtract row n from other rows.
            for(i=n+1;i<dim;i++)
            {
                F = M[i][n];
                for(j=n;j<dim;j++)
                    M[i][j] = M[i][j]-F*M[n][j];
            }

            //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 << "\rAll matrices are positive semidefinite." << endl;
    cout << endl;

    return true;
}

//Returns the number of blue vertices in the hypercube g.
long NumberOfVertices(CGraph& g)
{
    long i;
    long vertices;

    vertices = 0;
    for(i=0;i<(1<<g.Dim);i++)
        if(g.C[i]==true)
            vertices++;

    return vertices;
}

//Returns true if g1 and g2 are of the same dimension and have identically coloured vertices.
bool IsCopy(CGraph& g1, CGraph& g2)
{
    long i;

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

    for(i=0;i<(1<<g1.Dim);i++)
        if(g1.C[i]!=g2.C[i])
            return false;

    return true;
}

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

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

    if(gSub.Dim<0)
        return true;

    //Create all isomorphic copies of gSuper. We will do this by first choosing 'root', a
    //vertex that will be 0 in our new hypercube. Then we will permute through the gSuper.Dim
    //dimensions of the hypercube.
    for(root=0;root<(1<<gSuper.Dim);root++)
    for(p=0;p<Factorial[gSuper.Dim];p++)
    {
        //Due to the way we constructed Iso[][],
        //{Iso[p][0], Iso[p][1], ..., Iso[p][gSuper.Dim-1]} is a permutation of 
        //{0, 1, ..., gSuper.Dim-1}.

        //Construct the mapping between gSuper and the isomorphic hypercube.
        for(i=0;i<(1<<gSuper.Dim);i++)
        {
            //Initialize the mapping of vertex i to root.
            mapping[i] = root;

            //Check the component of dimension j in i.
            for(j=0;j<gSuper.Dim;j++)
            {
                //The j-th dimension is mapped to the Iso[p][j]-th dimension.

                if(((i>>j)&1)==1)
                {
                    //The j-th bit of i is 1, so to create our mapping we flip the
                    //Iso[p][j]-th bit of mapping[i].
                    mapping[i] ^= 1<<Iso[p][j];
                }

                //If the j-th bit of i is 0 we do not flip the Iso[p][j]-th bit.
            }
        }

        //Assume gSub is a (induced) subcube. If at some point we've shown it isn't, then try
        //a different mapping.
        bIsSub = true;
        for(i=0;i<(1<<gSub.Dim);i++)
        {
            if(bInduced==true)
            {
                //Check if the colours are identical.
                if(gSub.C[i]!=gSuper.C[mapping[i]])
                {
                    //Try a different mapping.
                    bIsSub = false;
                    break;
                }
            }
            else
            {
                //Check if gSub is a subcube.
                if(gSub.C[i]==true && gSuper.C[mapping[i]]==false)
                {
                    //Try a different mapping.
                    bIsSub = false;
                    break;
                }
            }
        }

        if(bIsSub==true)
            return true;
    }

    //Tested all mappings with no success.
    return false;
}

//Each blue coloured vertex set can be assigned a unique integer. FindMinIso() creates gOut a
//graph with the property that the flag (gOut, dimOfType) is isomorphic to the flag 
//(gIn, dimOfType) and has the smallest integer associated with its blue coloured vertex set
//out of all such graphs. This function is used by GenerateH() to avoid storing isomorphic
//copies of graphs, and by CalculateDensity() to speed up searching through lists of flags.
//In order for GenerateH() to function correctly the way FindMinIso() constructs the integer
//representing the blue vertex set is important. The function assumes Iso[][] has been
//initialized.
void FindMinIso(CGraph gIn, CGraph& gOut, long dimOfType)
{
    long   i, j, p, np;
    long   mapping[32];
    long   root;
    CGraph gTemp;

    //gOut holds the graph with the smallest integer associated with C[] so far.
    gOut = gIn;

    if(gIn.Dim<=0)
        return;

    //We wish to test all isomorphic flags of (gIn, dimOfType). We will do this by creating
    //all isomorphic copies of gIn. First we will choose 'root', a vertex that will be 0 in
    //our new hypercube. Then we will permute through the gIn.Dim dimensions of the
    //hypercube.
    for(root=0;root<(1<<gIn.Dim);root++)
    {
        //When dimOfType>=0 we wish to preserve the positions of vertices
        //0,...,2^dimOfType-1.
        if(dimOfType>=0 && root!=0)
            break;

        //Set 'np' the number of dimensions we wish to permute.
        if(dimOfType<0)
            np = gIn.Dim;
        else
            np = gIn.Dim-dimOfType;

        for(p=0;p<Factorial[np];p++)
        {
            //Construct the mapping between gIn and the isomorphic hypercube.
            for(i=0;i<(1<<gIn.Dim);i++)
            {
                //Initialize the mapping of vertex i to root.
                mapping[i] = root;

                //Check the component of dimension j in i.
                for(j=0;j<gIn.Dim;j++)
                {
                    //The j-th dimension is mapped to the (gIn.Dim-1-Iso[p][gIn.Dim-1-j])-th
                    //dimension. Due to the way we constructed Iso[][], this will cause the
                    //j-th dimension to be mapped to the j-th dimension when j<dimOfType,
                    //and permute the higher dimensions.

                    if(((i>>j)&1)==1)
                    {
                        //The j-th bit of i is 1, so to create our mapping we flip the
                        //(gIn.Dim-1-Iso[p][gIn.Dim-1-j])-th bit of mapping[i].
                        mapping[i] ^= 1<<(gIn.Dim-1-Iso[p][gIn.Dim-1-j]);
                    }

                    //If the j-th bit of i is 0 we do not flip the
                    //(gIn.Dim-1-Iso[p][gIn.Dim-1-j])-th bit.
                }
            }

            //Construct gTemp the isomorphic flag described by mapping[].
            gTemp.Dim = gIn.Dim;
            for(i=0;i<(1<<gIn.Dim);i++)
                gTemp.C[i] = gIn.C[mapping[i]];

            //We construct the unique integer associated with C[] in the following way.
            //Construct an array B[], such that B[i] = 1 if C[i]=true, and 0 if C[i]=false.
            //The unique integer in binary has bits B[0]B[1]...B[2^gIn.Dim-1] (where B[0] is
            //the most significant bit). We will not explicitly construct the integer, it
            //will be enough to compare values of C[].

            for(i=0;i<(1<<gIn.Dim);i++)
            {
                if(gOut.C[i]==false && gTemp.C[i]==true)
                    break; //gOut is smaller.

                if(gOut.C[i]==true && gTemp.C[i]==false)
                {
                    //gTemp is smaller.
                    gOut = gTemp;
                    break;
                }

                //gOut.C[i]==gTemp.C[i], try the next bit.
            }

            //Try the next mapping.
        }
    }

    return;
}

//Generates 'H' a linked list representing $\mathcal{H}$ the family of admissible hypercubes
//of dimension dimOfH. The function assumes dimOfH has been initialized. A bad_alloc object
//may be thrown.
void GenerateH()
{
    long        i, n;
    CLinkedList OldList;
    CLinkedList NewList;
    CListNode*  pNode;
    CGraph      g, gMin;
    long        b, p;
    long        graphIndex, total;
    bool        bAddZeroCube[2];

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

    //Generate the dimension zero hypercubes.
    cout << "\rGenerating dimension 0 hypercubes :" << flush;

    //bAddZeroCube[0] and bAddZeroCube[1] hold whether we should add the hypercubes which
    //consist solely of a single vertex coloured red and blue respectively.
    bAddZeroCube[0] = true;
    bAddZeroCube[1] = true;

    //for(i=0;i<noOfForbiddenGraphs;i++)
    //    if(ForbiddenGraph[i].Dim==0)
    //    {
    //        if(ForbiddenGraph[i].C[0]==true)
    //            bAddZeroCube[1] = false; //Forbid blue vertices.
    //        else
    //        {
    //            bAddZeroCube[0] = false; //Forbid red vertices.
    //
    //            if(bForbidOnlyInducedSubcubes==false)
    //                bAddZeroCube[1] = false; //Forbid blue vertices.
    //        }
    //    }
    
    if(bAddZeroCube[0]==true)
    {
        g.Dim  = 0;
        g.C[0] = false;
        OldList.Add(g);
    }

    if(bAddZeroCube[1]==true)
    {
        g.Dim  = 0;
        g.C[0] = true;
        OldList.Add(g);
    }

    for(n=1;n<=dimOfH;n++)
    {
        //Generate NewList, a list of hypercubes of dimension n, from OldList a list of
        //hypercubes of dimension n-1.

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

        pNode = OldList.pHead;
        while(pNode!=NULL)
        {
            graphIndex++;
            cout << "\rGenerating dimension " << n << " hypercubes : ";
            cout << "Extending hypercube " << graphIndex << " of " << total << "." << flush;

            //Create a new graph which is pNode->g with an extra dimension.
            g = pNode->g;
            g.Dim++;

            //There are 2^(g.Dim-1) new vertices which we have to colour. This means there
            //are 2^(2^(g.Dim-1)) possible ways to colour these vertices. We let b encode
            //the colouring via the value of its bits.
            for(b=0;b<(1<<(1<<(g.Dim-1)));b++)
            {
                for(p=0;p<(1<<(g.Dim-1));p++)
                {
                    if(((b>>p)&1)==1)
                        g.C[(1<<(g.Dim-1))+p] = true;
                    else
                        g.C[(1<<(g.Dim-1))+p] = false;
                }

                //Check if the graph is admissible.
                //for(i=0;i<noOfForbiddenGraphs;i++)
                //    if(IsSubcube(ForbiddenGraph[i],g,bForbidOnlyInducedSubcubes)==true)
                //        break;
                //
                //if(i!=noOfForbiddenGraphs)
                //    continue; //A forbidden graph was found.

                //We are guaranteed to generate all hypercubes which are minimally
                //represented due to the way FindMinIso() chooses a minimal representive of
                //the isomorphic copies of a hypercube and the way we are inductively
                //building up our list of hypercubes from smaller dimension hypercubes which
                //are minimally represented. Hence if the hypercube is not already the
                //minimal representation of all its isomorphic copies we can safely ignore
                //it.
                FindMinIso(g,gMin,-1);

                if(IsCopy(g,gMin)==false)
                    continue;
    
                NewList.Add(g);
            }

            //Move to the next hypercube in OldList.
            pNode = pNode->pNext;
        }

        //Make NewList the new OldList.
        OldList.Free();
        OldList.pHead = NewList.pHead;
        NewList.pHead = NULL;

        cout << "\rGenerating dimension " << n << " hypercubes :                          ";

        for(i=0;i<2*NumberOfDigits(total);i++)
            cout << " ";
    }

    //Move OldList to H.
    H.Free();
    H.pHead       = OldList.pHead;
    OldList.pHead = NULL;
    
    //Display how large H is.
    n = H.Size();

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

    cout << endl;

    return;
}

//Displays 'Density' as both a fraction and a decimal. This function is used by
//CalculateDensity() to display an upper bound of the induced density. It may throw an
//ErrorOverflow object.
void DisplayDensity(CLargeFraction Density)
{
    long      i;
    bool      bNegative;
    bool      bRounded;
    CLargeInt Digits, PowerOf10;
    
    Density.Simplify();

    cout << "\rUpper bound of the induced 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;
}

//Calculates then displays the upper bound of the induced density. The function assumes
//Iso[][], Term[], noOfTerms, noOfInducedGraphs, InducedGraph[], H, and dimOfH have been
//initialized. To save memory the function modifies Term[]. A bad_alloc or ErrorOverflow
//object may be thrown.
void CalculateDensity()
{
    long           i, j, p, q, r;
    long           nTerm;
    long           root;
    long           dimOfFlag;
    long           dimOfType;
    long           VertexMapping[2][32];
    long           DimensionMapping[2][5];
    long           index[2];
    CGraph         f[2];
    CListNode*     pNode;
    CLargeFraction Density;
    CTermInfo      TempMatrix;

    if(H.pHead==NULL)
    {
        DisplayDensity(0);
        return;
    }

    //Calculate the densities of the graphs in H.
    {
        cout << "\rCalculating the densities of the subcubes." << flush;

        pNode = H.pHead;
        while(pNode!=NULL)
        {
            pNode->coeff = 0;

            //Look at each permutation of the cube to check for copies of the induced graphs. 
            for(root=0;root<(1<<dimOfH);root++)
            for(p=0;p<Factorial[dimOfH];p++)
            {
                //Create a mapping from the vertices of a subcube of H to the vertex set of two
                //flags f[0], f[1]. Vertex v in flag f[r] will be the vertex
                //'VertexMapping[r][v]' in the subcube of H.
                
                //Create the vertex mappings.
                for(i=0;i<(1<<dimOfH);i++)
                {
                    //Initialize the mapping of vertex i to root.
                    VertexMapping[0][i] = root;

                    //Check the component of dimension j in i.
                    for(j=0;j<dimOfH;j++)
                    {
                        if(((i>>j)&1)==1)
                        {
                            //The j-th bit of i is 1, so to create our mapping we flip the
                            //Iso[p][j]-th bit of VertexMapping[0][i].
                            VertexMapping[0][i] ^= 1<<Iso[p][j];
                        }
                    }
                }

                //Create the new graph.
                f[0].Dim = dimOfH;
                for(i=0;i<(1<<dimOfH);i++)
                    f[0].C[i] = pNode->g.C[VertexMapping[0][i]];

                for(i=0;i<noOfInducedGraphs;i++)
                {
                    f[0].Dim = InducedGraph[i].Dim;
                    if(IsSubcube(InducedGraph[i],f[0],true)==true)
                        pNode->coeff = pNode->coeff+1;
                }
            }
            
            //Scale to convert the count into a density.
            pNode->coeff = pNode->coeff/(Factorial[dimOfH]<<dimOfH); 
            
            pNode = pNode->pNext;
        }

        cout << "\r                                          ";
    }

    //Calculate the contribution of the terms.
    for(nTerm=1;nTerm<=noOfTerms;nTerm++)
    {
        CTermInfo& T = Term[nTerm-1];

        cout << "\rEvaluating term " << nTerm << " of " << noOfTerms << " : " << flush;

        if(T.MatrixDim<=0 || T.noOfFlags<=0)
        {
            T.Free();
            continue;
        }

        //Replace the flags with their "minimal" isomporhic copy. 
        {
            cout << "\rEvaluating term " << nTerm << " of " << noOfTerms << " : ";
            cout << "Initializing flags." << flush;

            dimOfType = T.Type.Dim;
            dimOfFlag = T.FlagList[0].Dim;

            for(i=0;i<T.noOfFlags;i++)
                FindMinIso(T.FlagList[i],T.FlagList[i],dimOfType);
        }

        //Calculate the effect of the basis on the matrix.
        {
            cout << "\rEvaluating term " << nTerm << " of " << noOfTerms << " : ";
            cout << "Initializing matrix." << flush;

            //Allocate memory for the new matrix.
            TempMatrix.MatrixDim = T.noOfFlags;
            TempMatrix.Matrix    = new CLargeFraction*[TempMatrix.MatrixDim];

            for(i=0;i<TempMatrix.MatrixDim;i++)
                TempMatrix.Matrix[i] = NULL;

            for(i=0;i<TempMatrix.MatrixDim;i++)
                TempMatrix.Matrix[i] = new CLargeFraction[TempMatrix.MatrixDim];

            //Set TempMatrix.Matrix = T.Basis^T T.Matrix T.Basis.
            for(i=0;i<TempMatrix.MatrixDim;i++)
            for(j=0;j<TempMatrix.MatrixDim;j++)
            {
                TempMatrix.Matrix[i][j] = 0;

                for(p=0;p<T.MatrixDim;p++)
                for(q=0;q<T.MatrixDim;q++)
                {
                    //To avoid costly CLargeFraction operations check if we are adding zero.
                    if(T.Basis[p][i]==0 || T.Basis[q][j]==0 || T.Matrix[p][q]==0)
                        continue;

                    TempMatrix.Matrix[i][j] = TempMatrix.Matrix[i][j]
                                              +T.Basis[p][i]*T.Matrix[p][q]*T.Basis[q][j];
                }
            }

            //Free T.Matrix[][] and T.Basis[][] (as they are no longer required).
            for(i=0;i<T.MatrixDim;i++)
                delete[] T.Matrix[i];

            delete[] T.Matrix;
            T.Matrix = NULL;
    
            for(i=0;i<T.MatrixDim;i++)
                delete[] T.Basis[i];

            delete[] T.Basis;
            T.Basis = NULL;

            //Store TempMatrix.Matrix[][] in T.Matrix[][].
            T.MatrixDim          = TempMatrix.MatrixDim;
            T.Matrix             = TempMatrix.Matrix;
            TempMatrix.MatrixDim = 0; 
            TempMatrix.Matrix    = NULL;

            //To calculate the contribution of the term we will multiply elements of the
            //matrix by probabilities. These probabilities will always be of the form
            //<an integer>/(dimOfH!*2^dimOfH). To avoid costly CLargeFraction operations we
            //will pre-divide the matrix by (dimOfH!*2^dimOfH).
            for(i=0;i<T.MatrixDim;i++)
            for(j=0;j<T.MatrixDim;j++)
                T.Matrix[i][j] = T.Matrix[i][j]/(Factorial[dimOfH]*(1<<dimOfH));
        }

        //Count how many times FlagList[a] intersected with FlagList[b] appears in a subcube
        //H. Using this count we can calculate $E[p(F_a,F_b,\theta;H)]$, and hence the
        //contribution Term[nTerm-1] makes to the upper bound.

        //Generate all isomorphic copies of H by first choosing 'root', a vertex that will
        //be 0 in our new hypercube, then choosing a permutation of the 'dimOfH' dimensions.
        for(root=0;root<(1<<dimOfH);root++)
        for(p=0;p<Factorial[dimOfH];p++)
        {
            cout << "\rEvaluating term " << nTerm << " of " << noOfTerms << " : ";
            cout << "Checking mapping " << root*Factorial[dimOfH]+p+1 << " of ";
            cout << Factorial[dimOfH]*(1<<dimOfH) << "." << flush;

            //Create a mapping from the vertices of a subcube of H to the vertex set of two
            //flags f[0], f[1]. Vertex v in flag f[r] will be the vertex
            //'VertexMapping[r][v]' in the subcube of H.

            //Set the dimensions which will form the type.
            for(i=0;i<dimOfType;i++)
            {
                DimensionMapping[0][i] = Iso[p][i];
                DimensionMapping[1][i] = Iso[p][i];
            }

            //Map the remaining dimensions.
            for(i=dimOfType;i<dimOfFlag;i++)
            {
                DimensionMapping[0][i] = Iso[p][i];
                DimensionMapping[1][i] = Iso[p][i+dimOfFlag-dimOfType];
            }

            //Create the vertex mappings.
            for(i=0;i<(1<<dimOfFlag);i++)
            {
                //Initialize the mapping of vertex i to root.
                VertexMapping[0][i] = root;
                VertexMapping[1][i] = root;

                //Check the component of dimension j in i.
                for(j=0;j<dimOfFlag;j++)
                {
                    if(((i>>j)&1)==1)
                    {
                        //The j-th bit of i is 1, so to create our mapping we flip the
                        //VertexMapping[][j]-th bit of VertexMapping[][i].
                        VertexMapping[0][i] ^= 1<<DimensionMapping[0][j];
                        VertexMapping[1][i] ^= 1<<DimensionMapping[1][j];
                    }
                }
            }

            //Apply VertexMapping[][] to each subcube in the list H to get two flags f[0],
            //f[1].
            pNode = H.pHead;
            while(pNode!=NULL)
            {
                //Create the two flags.
                for(r=0;r<2;r++)
                {
                    f[r].Dim = dimOfFlag;

                    for(i=0;i<(1<<dimOfFlag);i++)
                        f[r].C[i] = pNode->g.C[VertexMapping[r][i]];
                }

                //Check types match to T.Type, by temporarily reducing the dimension of f[0].
                f[0].Dim = dimOfType;
                if(IsCopy(f[0],T.Type)==false)
                {
                    pNode = pNode->pNext;
                    continue;
                }

                f[0].Dim = dimOfFlag;

                //Find the isomorphic copy of each flag in T.FlagList[], and store the
                //position at which it appears in index[].
                for(r=0;r<2;r++)
                    FindMinIso(f[r],f[r],dimOfType);

                for(index[0]=0;index[0]<T.noOfFlags;index[0]++)
                    if(IsCopy(f[0],T.FlagList[index[0]])==true)
                    {
                        for(index[1]=0;index[1]<T.noOfFlags;index[1]++)
                            if(IsCopy(f[1],T.FlagList[index[1]])==true)
                            {
                                //Rather than keep track of how many times f[0] intersected
                                //with f[1] occurs then use this to update pNode->coeff, we
                                //will update pNode->coeff directly.
                                pNode->coeff = pNode->coeff+T.Matrix[index[0]][index[1]];
                            }
                    }

                pNode = pNode->pNext;
            }
        }

        cout << "\rEvaluating term " << nTerm << " of " << noOfTerms << " : ";
        cout << "                      ";
        for(i=0;i<2*NumberOfDigits(Factorial[dimOfH]*(1<<dimOfH));i++)
            cout << " ";

        T.Free();
    }

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

    //Find the upper bound of the Turan density.
    {
        cout << "\rFinding the maximum coefficient." << flush;

        Density = 0;
        pNode = H.pHead;
        while(pNode!=NULL)
        {
            if(Density<pNode->coeff)
                Density = pNode->coeff;
        
            pNode = pNode->pNext;
        }

        DisplayDensity(Density);
    }

    return;
}

//The main function. Expects argc to be 2, argv[0] to be the name of the executable and
//argv[1] to be the name of the input file.
int main(int argc, char* argv[])
{
    long i;
    bool bHasSpaces;
    
    try
    {
        cout << endl;

        //Check that a file was passed to the program via the command line.
        if(argc!=2)
        {
            //Check if argv[0][] has spaces.
            bHasSpaces = false;
            i          = 0;
            while(argv[0][i]!='\0')
            {
                if(argv[0][i]==' ')
                {
                    bHasSpaces = true;
                    break;
                }

                i++;
            }

            cout << "Usage : " << endl;
            if(bHasSpaces==true)
                cout << "\"" << argv[0] << "\"";
            else
                cout << argv[0];

            cout << " <Input Data File>" << endl;
            cout << endl;

            return 0;
        }

        //Read the input file.
        if(ReadDataFile(argv[1])==false)
        {
            Cleanup();
            return 0;
        }

        //Initialize Iso[][].
        cout << "\rAllocating memory for the isomorphism mapping." << flush;
        
        InitIso();
        
        cout << "\rChecking the matrices are positive semidefinite....." << endl;
        
        if(IsPSD()==false)
        {
            cout << endl;
            cout << endl;
            cout << "*** ERROR: Matrix is not positive semidefinite. ***" << endl;
            cout << "For the calculation to return an upper bound of the" << endl;
            cout << "induced density all matrices must be positive semidefinite." << endl;
            cout << endl;

            Cleanup();
            return 0;
        }

        cout << "\rGenerating admissible subcubes of dimension " << dimOfH;
        cout << "....." << endl;

        GenerateH();

        cout << "\rCalculating the upper bound of the induced density....." << endl;

        CalculateDensity();

        Cleanup();
        return 0; //Calculation completed successfully.
    }
    catch(bad_alloc& e)
    {
        cout << endl;
        cout << endl;
        cout << "*** ERROR: Bad allocation. ***" << endl;
        cout << "Probably caused by lack of memory." << endl;
        cout << "Try decreasing LI_MAX in the source code." << endl;
    }
    catch(ErrorOverflow& e)
    {
        cout << endl;
        cout << endl;
        cout << "*** ERROR: Overflow. ***" << endl;
        cout << "The calculation encountered an integer that was too large to store.\n";
        cout << "Try increasing LI_MAX in the source code." << endl;
    }
    catch(ErrorDivisionByZero& e)
    {
        cout << endl;
        cout << endl;
        cout << "*** ERROR: Division by zero. ***" << endl;
    }
    catch(...)
    {
        cout << endl;
        cout << endl;
        cout << "*** ERROR: Unknown. ***" << endl;
    }

    cout << endl;

    Cleanup();
    return 0;
}

