/*
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        : DensityChecker.
//
//Author      : Rahil Baber.
//
//Date        : 10th November 2010.
//
//Compilation : This program has been successfully compiled under Visual C++ 2010, and
//              gcc 4.4.1. E.g. in Linux type "g++ -O3 DensityChecker.cpp -o dcheck"
//
//Usage       : <executable> <input file>
//              E.g. in Linux "./dcheck 2-9-prf.txt"
//
//Description : Using Razborov's method it calculates a bound on the Turan density for a
//              given family of 3-graphs. 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 3-graphs of order at most 9. Users may
//              wish/need to change the hard coded values 'LI_MAX' and 'noOfDecimalPlaces'
//              specified below. Note that while some of our proofs can be verified in under
//              10 minutes others may take several hours to verify.

#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 (increase) 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 3-graphs that have at most nine vertices. Because C++ indexes
//arrays starting from 0, we take our vertex set to be {0,1,...,nV-1}, where nV is the number
//of vertices in our graph. We can represent a flag as (g, n) a CGraph g and an integer n, by
//taking the labelled vertices to be 0,1,...,n-1 and to avoid confusion we label them
//0,1,...,n-1 respectively (instead of 1,2,...,n).
class CGraph
{
public:
    //Number of vertices.
    long nV;

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

//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 graphs
    //of order "orderOfH"). 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 graphs of order "orderOfH"). 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 order of the graphs $H_i$. Denoted by the letter $l$ in the thesis. Due to the way the
//function GenerateH() and ReadGraphEdgeSet() has been implemented orderOfH can be at most 9.
//In order for the density calculation to be valid orderOfH must be at least 3. It is
//initialized by ReadDataFile().
long orderOfH;

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

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

//If bForbidOnlyInducedSubgraphs is set to false, the admissible graphs are those that do not
//contain a subgraph from our forbidden family $\mathcal{F}$. If set to true then the
//admissible graphs are those that do not contain an *induced* subgraph from $\mathcal{F}$.
//It is initialized by ReadDataFile().
bool bForbidOnlyInducedSubgraphs;

//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[10] = {1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880};

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

//*** Function declarations. ***

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

//Allocates memory for Iso[][] and initializes it. The function assumes orderOfH has been set
//to a value between 3 and 9. Iso[][] is initialized to an orderOfH! by orderOfH array
//representing all permutations of {0,1,...,orderOfH-1}. Iso[p][n] will be the image of n
//under permutation p. Iso[][] will be ordered in such a way that {Iso[p][] : p<=r!} is the
//set of permutations in which only {0,1,..,r-1} are permuted and {r,...,orderOfH-1} are left
//in place. 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 edge set of a graph. The function assumes g.nV has been set and uses this to
//verify the edges are valid. The edge set should take the form "{<edge>,<edge>,...,<edge>}",
//where <edge> is a three digit integer representing the edge. For example "{123, 325, 167}"
//or the empty set "{}". The three integers representing the edge should be distinct and be
//between 1 and g.nV (inclusive). Whitespaces either side of the <edge> sequences 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 ReadGraphEdgeSet(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 orderOfH, noOfForbiddenGraphs,
//ForbiddenGraph[], bForbidOnlyInducedSubgraphs, 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(), ReadGraphEdgeSet(), 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 edges in the graph g.
long NumberOfEdges(CGraph& g);

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

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

//Each edge set of a graph can be assigned a unique integer. FindMinIso() creates gOut a
//graph with the property that the flag (gOut, orderOfType) is isomorphic to the flag 
//(gIn, orderOfType) and has the smallest integer associated with its edge set out of all
//such graphs. This function is used by 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 edge set is important. The function assumes Iso[][] has been initialized.
void FindMinIso(CGraph gIn, CGraph& gOut, long orderOfType);

//Generates 'H' a linked list representing $\mathcal{H}$ the family of admissible graphs of
//order orderOfH. The function assumes orderOfH, ForbiddenGraph[], noOfForbiddenGraphs, and
//bForbidOnlyInducedSubgraphs, have 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 Turan density. It may throw an
//ErrorOverflow object.
void DisplayDensity(CLargeFraction Density);

//Calculates then displays the upper bound of the Turan density. The function assumes
//Iso[][], Term[], noOfTerms, H, and orderOfH 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  = (1L<<(8*sizeof(long)-2))-1;
        max += 1L<<(8*sizeof(long)-2);
        min  = -max;
        min -= 1;

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

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

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

    Size      = 0;
    bPositive = true;

    if(input==0)
        return;

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

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

    return;
}

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

    //A and B have the same sign.

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

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

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

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

    return 0; //A = B.
}

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

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

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

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

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

    bitN = i;

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

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

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

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

        //Move to next bit.
        rBit--;

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

        if(rByte==-1)
            break;

        nBit--;

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

    //outputR has been initialized.

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

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

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

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

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

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

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

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

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

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

        //Calculate the next bit of outputQ.
    }

    //outputR holds the magnitude of the remainder.

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

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

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

        outputQ.Size--;
    }

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

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

        outputR.Size--;
    }

    return;
}

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

    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, j, k;
    bool bFirst;

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

    //Display the edges.
    out << "Edges = {";

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

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

    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, ForbiddenGraph[], Term[], and Iso[][]. To deallocate
//the memory associated with Iso[][] it requires orderOfH to be set.
void Cleanup()
{
    long i;

    H.Free();

    delete[] ForbiddenGraph;
    ForbiddenGraph = NULL;

    delete[] Term;
    Term = NULL;

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

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

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

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

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

    //Initialize Iso[][].

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

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

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

        p++;
    }
    
    return;
}

//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 edge set of a graph. The function assumes g.nV has been set and uses this to
//verify the edges are valid. The edge set should take the form "{<edge>,<edge>,...,<edge>}",
//where <edge> is a three digit integer representing the edge. For example "{123, 325, 167}"
//or the empty set "{}". The three integers representing the edge should be distinct and be
//between 1 and g.nV (inclusive). Whitespaces either side of the <edge> sequences 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 ReadGraphEdgeSet(istream& in, CGraph& g, char& c)
{
    long i, j, k;

    //Initialize the adjacency matrix.
    for(i=0;i<9;i++)
    for(j=0;j<9;j++)
    for(k=0;k<9;k++)
        g.M[i][j][k] = 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 edge.
        if('1'<=c && c<='9')
            i = c-'1';
        else
            return false;

        c = ReadChar(in);

        if('1'<=c && c<='9')
            j = c-'1';
        else
            return false;

        c = ReadChar(in);

        if('1'<=c && c<='9')
            k = c-'1';
        else
            return false;

        c = ReadChar(in);
        
        //Check ijk is a valid edge.
        if(i==j || i==k || j==k)
            return false;

        if(i>=g.nV || j>=g.nV || k>=g.nV)
            return false;

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

        //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 orderOfH, noOfForbiddenGraphs,
//ForbiddenGraph[], bForbidOnlyInducedSubgraphs, 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(), ReadGraphEdgeSet(), and ReadBasis().
bool ReadDataFile(const char* filename)
{
    long      i, j, k, 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 : DensityChecker\".";

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

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

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

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

        //Read noOfForbiddenGraphs.
        {
            r = ReadLong(FileIn,"Number of forbidden graphs",noOfForbiddenGraphs,c);
            
            if(r==2)
                throw "1:";

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

            cout << "Number of forbidden graphs : " << noOfForbiddenGraphs << endl;
            
            ForbiddenGraph = new CGraph[noOfForbiddenGraphs];
        }

        //Read ForbiddenGraph[].
        {
            //Read "Forbidden graphs :".
            RemoveWhitespace(FileIn,c);

            if(CheckString(FileIn,"Forbidden graphs",c)==false)
                throw "0:Expected to read \"Forbidden graphs : <List of graphs>\".";

            RemoveWhitespace(FileIn,c);
        
            if(c!=':')
                throw "0:Expected to read \"Forbidden graphs : <List of graphs>\".";

            c = ReadChar(FileIn);

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

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

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

                RemoveWhitespace(FileIn,c);

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

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

                if(ForbiddenGraph[i].nV>9 || ForbiddenGraph[i].nV<0)
                    throw "2:";

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

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

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

        //Read bForbidOnlyInducedSubgraphs.
        {            
            RemoveWhitespace(FileIn,c);

            if(CheckString(FileIn,"Forbid only induced subgraphs",c)==false)
                throw "0:Expected to read \"Forbid only induced subgraphs : <yes or no>\".";

            RemoveWhitespace(FileIn,c);
        
            if(c!=':')
                throw "0:Expected to read \"Forbid only induced subgraphs : <yes or no>\".";

            c = ReadChar(FileIn);
            
            if(ReadBool(FileIn,bForbidOnlyInducedSubgraphs,c)==false)
                throw "0:Expected to read \"Forbid only induced subgraphs : <yes or no>\".";
                    
            cout << "Forbid induced subgraphs   : ";
            cout << (bForbidOnlyInducedSubgraphs?"yes":"no") << 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.nV.
                {
                    cout << "\rReading term " << nTerm << " : ";
                    cout << "Reading the order of the type." << flush;

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

                    //Initialize the edge set of the type.
                    for(i=0;i<9;i++)
                    for(j=0;j<9;j++)
                    for(k=0;k<9;k++)
                        Term[nTerm-1].Type.M[i][j][k] = false;
                }

                //Read the order of the flags.
                {
                    cout << "\rReading term " << nTerm << " : ";
                    cout << "Reading the order of the flags." << flush;

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

                    if(r!=0 || n<Term[nTerm-1].Type.nV
                    || (orderOfH+Term[nTerm-1].Type.nV)/2<n)
                        throw "0:Expected to read \"Order of flags : <An integer in [<Order "
                              "of type>,\n<Order of type + Order of subgraphs 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];

                        //Set the order of the flags.
                        for(i=0;i<Term[nTerm-1].noOfFlags;i++)
                            Term[nTerm-1].FlagList[i].nV = n;
                    }
                }

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

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

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

                        //Check the flag has the same induced type as the other flags.
                        for(i=0;i<Term[nTerm-1].Type.nV;i++)
                        for(j=0;j<Term[nTerm-1].Type.nV;j++)
                        for(k=0;k<Term[nTerm-1].Type.nV;k++)
                            if(Term[nTerm-1].FlagList[nFlag-1].M[i][j][k]!=
                               Term[nTerm-1].Type.M[i][j][k])
                                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"
                                      "1 to <Order of type>.";
                    }

                    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 graph in the form" << endl;
            cout << "\"<The order - an integer in [0,9]> : <The edge set>\"." << endl;
            cout << "For example \"4 : {123, 124, 234}\"." << 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 << " : <The flag's edge set>\"." << endl;
            cout << "An example of an edge set for an order 4 graph is" << endl;
            cout << "\"{123, 124, 234}\"." << 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()
{
    //Overview of the algorithm:
    //
    //This function checks if the matrices stored in Term[] are indeed positive semidefinite.
    //
    //It takes the matrix M and essentially tries to Cholesky decompose it, in particular to
    //avoid square root operations we choose the variant form of U^T D U where D is diagonal
    //(and U is upper triangle), however we don't really care about U or D just that the
    //process is possible.
    //
    //The algorithm works by doing symmetric row and column operations on M. If C is the
    //operation of adding lambda times column i to column j then, C^T M C is the operation of
    //adding lambda times column i to column j, and lambda times row i to row j. Note that
    //C^T M C is symmetric and the property of whether it is PSD or not is preserved
    //(provided i does not equal j). Using these symmetric row and column operations we can
    //try to diagonalise M, and in the process prove whether it is PSD or not.
    //
    //If M is a 1 by 1 matrix then it is easy to check if it is PSD or not. If it isn't a
    //1 by 1 matrix then do the following: 
    //
    //If M[0][0]<0 then the matrix is not PSD.
    //
    //If M[0][0]=0 then either there exists an i such that M[0][i] is not 0 in which case the
    //matrix is not PSD, or all M[0][i]=M[i][0]=0, in which case it is equivalent to remove
    //row 0 and column 0 to create the matrix M' and to check whether M' is PSD or not by
    //restarting the process with M=M'.
    //
    //If M[0][0]>0 then subtract rows and columns to make M[i][0]=M[0][i]=0 for all i (except
    //i=0). If we throw away row 0 and column 0 the resulting matrix M' is PSD if and only if
    //M is. We can check M' by restarting the process with M=M'.
    //
    //It is costly to create M' each time we deal with row 0 and column 0 so instead we will
    //continue to modify M and keep track of how many rows and columns we have discarded
    //through the use of the variable n.  
    //
    //Note that if we keep track of the operations we would get the U^T D U Cholesky
    //decomposition, however the entries are likely ro be extremely large, which is why the
    //data files do not store the matrices in this form.

    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
        //the property of being 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 edges in the graph g.
long NumberOfEdges(CGraph& g)
{
    long i, j, k;
    long edges;

    edges = 0;

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

    return edges;
}

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

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

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

    return true;
}

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

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

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

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

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

                if(bIsSub==false)
                    break;
            }

            if(bIsSub==false)
                break;
        }

        if(bIsSub==true)
            return true;
    }

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

//Each edge set of a graph can be assigned a unique integer. FindMinIso() creates gOut a
//graph with the property that the flag (gOut, orderOfType) is isomorphic to the flag 
//(gIn, orderOfType) and has the smallest integer associated with its edge set out of all
//such graphs. This function is used by 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 edge set is important. The function assumes Iso[][] has been initialized.
void FindMinIso(CGraph gIn, CGraph& gOut, long orderOfType)
{
    long i, j, k, i2, j2, k2, p;
    bool bFoundSmaller, bBreak;

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

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

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

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

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

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

                    //Compare bits

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

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

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

                if(bBreak==true)
                    break;
            }

            if(bBreak==true)
                break;
        }

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

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

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

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

    return;
}

//Generates 'H' a linked list representing $\mathcal{H}$ the family of admissible graphs of
//order orderOfH. The function assumes orderOfH, ForbiddenGraph[], noOfForbiddenGraphs, and
//bForbidOnlyInducedSubgraphs, have been initialized. A bad_alloc object may be thrown.
void GenerateH()
{
    long        i, j, k, n;
    CLinkedList OldList;
    CLinkedList NewList;
    CListNode*  pNode;
    CGraph      g, gMin;
    long        noOfEdgesToSet;
    long        b, p;
    bool        e;
    long        graphIndex, total; 

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

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

    //Check if order 0 graphs are forbidden.
    for(i=0;i<noOfForbiddenGraphs;i++)
        if(ForbiddenGraph[i].nV==0)
        {
            cout << "\rFound 0 graphs.                       " << endl;
            cout << endl;
            return;
        }

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

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

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

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

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

            //Create a new graph which is pNode->g with an extra vertex v.
            g = pNode->g;
            g.nV++;

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

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

                    //Update the adjacency matrix.
                    g.M[i][j][k] = e;
                    g.M[i][k][j] = e;
                    g.M[j][i][k] = e;
                    g.M[j][k][i] = e;
                    g.M[k][i][j] = e;
                    g.M[k][j][i] = e;
            
                    p++;
                }

                //Check if the graph is admissible.
                for(i=0;i<noOfForbiddenGraphs;i++)
                    if(IsSubgraph(ForbiddenGraph[i],g,bForbidOnlyInducedSubgraphs)==true)
                        break;

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

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

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

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

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

        cout << "\rGenerating order " << n << " graphs :                      ";
        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 graph." << endl;
    else
        cout << "\rFound " << n << " graphs."<< 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 Turan 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 Turan density" << endl;
    
    if(Density.D==1)
    {
        cout << " = " << Density.N << endl;
        cout << endl;
        return;
    }

    cout << " = " << Density;
    
    //Display Density as a decimal to 'noOfDecimalPlaces' decimal places.

    //Check if Density is negative.
    if(Density<0)
    {
        Density   = -Density;
        bNegative = true;
    }
    else
        bNegative = false;

    //Calculate the digits that will make up the decimal.
    
    //Calculate 10^noOfDecimalPlaces.
    PowerOf10 = 1;
    for(i=0;i<noOfDecimalPlaces;i++)
        PowerOf10 = PowerOf10*10;

    Density = Density*PowerOf10;

    //'Digits' holds Density*10^noOfDecimalPlaces rounded down. We wish to round up.
    Digits = Density.N/Density.D;
    
    //Check if rounding up is necessary.
    if(Density.N%Density.D==0)
        bRounded = false;
    else
    {
        //Rounding is necessary, currently the magnitude has been rounded down. If the result
        //is negative then this is fine. If the result is positive then we add one so that
        //the magnitude is rounded up.
        bRounded = true;

        if(bNegative==false)
            Digits = Digits+1;
    }

    //Start displaying the digits.
    cout << endl;
    cout << " = ";
    
    if(bNegative==true)
        cout << "-";

    //Display the digits that lie to the left of the decimal point.
    cout << Digits/PowerOf10 << ".";

    //Display the digits that lie to the right of the decimal point, one at a time.
    Digits = Digits%PowerOf10;
    for(i=0;i<noOfDecimalPlaces;i++)
    {
        //Display the digit (10*Digits)/PowerOf10 rounded down.
        Digits = Digits*10;

        cout << Digits/PowerOf10;

        Digits = Digits%PowerOf10;
    }

    if(bRounded==true)
        cout << " (rounded up)";
    
    cout << endl;
    cout << endl;

    return;
}

//Calculates then displays the upper bound of the Turan density. The function assumes
//Iso[][], Term[], noOfTerms, H, and orderOfH have been initialized. To save memory the
//function modifies Term[]. A bad_alloc or ErrorOverflow object may be thrown.
void CalculateDensity()
{
    long           i, j, k, p, q, r;
    long           nTerm;
    long           orderOfFlag;
    long           orderOfType;
    long           mapping[2][9];
    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 subgraphs." << flush;

        pNode = H.pHead;
        while(pNode!=NULL)
        {
            //Density is (number of edges)/(orderOfH choose 3).
            pNode->coeff = CLargeFraction(NumberOfEdges(pNode->g)*6)
                           /(orderOfH*(orderOfH-1)*(orderOfH-2)); 
            
            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;

            orderOfType = T.Type.nV;
            orderOfFlag = T.FlagList[0].nV;

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

        //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>/(orderOfH!). To avoid costly CLargeFraction operations we will
            //pre-divide the matrix by (orderOfH!).
            for(i=0;i<T.MatrixDim;i++)
            for(j=0;j<T.MatrixDim;j++)
                T.Matrix[i][j] = T.Matrix[i][j]/Factorial[orderOfH];
        }

        //Count how many times FlagList[a] intersected with FlagList[b] appears in a subgraph
        //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.

        //We will use Iso[p][] to enumerate through the possible ways the flags could appear
        //in a subgraph H.
        for(p=0;p<Factorial[orderOfH];p++)
        {
            cout << "\rEvaluating term " << nTerm << " of " << noOfTerms << " : ";
            cout << "Checking mapping " << p+1 << " of " << Factorial[orderOfH] << ".";
            cout << flush;

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

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

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

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

            //Apply mapping[][] to each graph 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].nV = orderOfFlag;

                    for(i=0;i<orderOfFlag;i++)
                    for(j=0;j<orderOfFlag;j++)
                    for(k=0;k<orderOfFlag;k++)
                        f[r].M[i][j][k] =
                            pNode->g.M[mapping[r][i]][mapping[r][j]][mapping[r][k]];
                }

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

                f[0].nV = orderOfFlag;

                //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],orderOfType);

                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[orderOfH]);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 << "Turan density all matrices must be positive semidefinite." << endl;
            cout << endl;

            Cleanup();
            return 0;
        }

        cout << "\rGenerating admissible subgraphs of order " << orderOfH << "....." << endl;

        GenerateH();

        cout << "\rCalculating the upper bound of the Turan 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;
}