Durango Bill’s
“C” Program to generate Cabtaxi Numbers
(Source Code)



Return to the main Ramanujan page

Web page generated via KompoZer


//*****************************************************************************
//
//                                Cabtaxi.c
//                                    by
//                                Bill Butler
//
//          May 5, 2008 version (as opposed to multiple earlier versions)
//
//    This program finds CabTaxi numbers where the number can be formed in at
//    least 7 different ways by either the sum of two cubes and/or the difference
//    of two cubes. The tables below illustrate low orders of CabTaxi numbers.
//    (Specifically CabTaxi(2) and CabTaxi(3)).
//
//    For the sum of two cubes, imagine that the following table is an Excel
//    speadsheet. The first row across the top counts upward for J = 1, 2, 3, 4,
//    etc. The second row gives the cube for each of these numbers: 1, 8, 27, 64,
//    etc. Similarly, the first column counts upward for I = 0, 1, 2, 3, 4, etc.
//    and the 2nd column gives the cube for each of these numbers: 0, 1, 8, 27,
//    64, etc. The rest of the cells in the spreadsheet are the sum of two cubes
//    using the respective rows and columns.
//
//           J    1    2    3    4    5    6    7    8    9   10   11   12
//           J^3  1    8   27   64  125  216  343  512  729 1000 1331 1728
//  I  I^3   -------------------------------------------------------------
//  0    0   |    1    8   27   64  125  216  343  512  729 1000 1331 1728
//  1    1   |    2    9   28   65  126  217  344  513  730 1001 1332 1729
//  2    8   |    9   16   35   72  133  224  351  520  737 1008 1339 1736
//  3   27   |   28   35   54   91  152  243  370  539  756 1027 1358 1755
//  4   64   |   65   72   91  128  189  280  407  576  793 1064 1395 1792
//  5  125   |  126  133  152  189  250  341  468  637  854 1125 1456 1853
//  6  216   |  217  224  243  280  341  432  559  728  945 1216 1547 1944
//  7  343   |  344  351  370  407  468  559  686  855 1072 1343 1674 2071
//  8  512   |  513  520  539  576  637  728  855 1024 1241 1512 1843 2240
//  9  729   |  730  737  756  793  854  945 1072 1241 1458 1729 2060 2457
// 10 1000   | 1001 1008 1027 1064 1125 1216 1343 1512 1729 2000 2331 2728
// 11 1331   | 1332 1339 1358 1395 1456 1547 1674 1674 1843 2060 2331 2662
// 12 1728   | 1729 1736 1755 1792 1853 1944 2071 2240 2457 2728 3059 3456
//
//    We note that except for the "0" row, the upper right triangle is a mirror
//    image of the lower left triangle. Thus, from here on, we will ignore the
//    lower left triangle.
//
//    This next table is similar to the table above except the values are the
//    differences of two cubes. Specifically, the values in the table equal
//    I^3 - J^3. If the upper right triangle were included, it would just give
//    the negative values of those seen in the lower left triangle.
//
//           J    0    1    2    3    4    5    6    7    8    9   10   11   12
//           J^3  0    1    8   27   64  125  216  343  512  729 1000 1331 1728
//  I  I^3   ------------------------------------------------------------------
//  0    0   |    0
//  1    1   |    1    0
//  2    8   |    8    7    0
//  3   27   |   27   26   19    0
//  4   64   |   64   63   56   37    0
//  5  125   |  125  124  117   98   61    0
//  6  216   |  216  215  208  189  152   91    0
//  7  343   |  343  342  335  316  279  218  127    0
//  8  512   |  512  511  504  485  448  387  296  169    0
//  9  729   |  729  728  721  702  665  604  513  386  217    0
// 10 1000   | 1000  999  992  973  936  875  784  657  488  271    0
// 11 1331   | 1331 1330 1323 1304 1267 1206 1115  988  819  602  331    0
// 12 1728   | 1728 1727 1720 1701 1664 1603 1512 1385 1216  999  728  397   0
//
//    Note: Additional diagrams in the NextGroup() routine show how these tables
//    are processed.
//
//    CabTaxi(2) is defined as the lowest number that can be formed by the sum
//    and/or difference of 2 cubes in two different ways. (Ignore the zeros in
//    table 2.)
//
//    If we search the upper right portion of the first table and also search the
//    lower table, we find that the lowest number that appears exactly twice is
//    91. It appears at row 3 column 4 in the first table and row 6 column 5 in
//    the second table.
//
//    Thus CabTaxi(2) = 91
//    = 3^3 + 4^3
//    = 6^3 - 5^3
//
//    Similarly the two tables can be searched for the lowest number that appears
//    three times - which turns out to be 728.
//
//    Thus CabTaxi(3) = 728
//    = 6^3 + 8^3            Row 6 column 8 in the first table
//    = 9^3 - 1^3            Row 9 column 1 in the second table
//    = 12^3 - 10^3        Row 12 column 10 in the second table
//
//    The process can be extended to higher orders for CabTaxi(N). Of course the
//    tables above would have to be extended to much larger sizes.
//
//    If we are to extend the process of finding additional CabTaxi numbers to
//    higher order matches (e.g. CabTaxi(10)), three major problems become
//    apparent.
//    1) The tables will become very large.
//    2) The search time will become prohibitive.
//    3) The size of the numbers in the cells will run into precision problems.
//    (Too many digits for standard computer hardware.)
//
//    We note that numbers in the two tables are confined to specific areas. In
//    the first table, numbers from 1 to 1,000 appear only in the upper left
//    corner. Numbers between 1,000 and 2,000 occur in a band that extends from
//    upper right to lower left. Numbers between 2,000 and 3,000 exist only in
//    another band that begins off the right side of the table and extends down
//    to the left. If we are looking for possible number matches, we only have
//    to look at one band at a time.
//
//    Similar sloping bands exist in the second table, but these slope from upper
//    left to lower right.. We also note that bands for the difference of two
//    cubes will tend to parallel the main diagonal. Thus processing within the
//    program for the difference of two cubes will sequentially search down the
//    diagonals whereas the search for the sum of two cubes portion will
//    sequentially search down the columns.
//
//    The program sequentially generates these bands. The new bands (one from the
//    I^3 + J^3 table and another from the I^3 - J^3 table) are used for the
//    current search process. Then these old bands are discarded and new bands
//    are generated for the next batch. The problem then is to search all of the
//    numbers within the bands to see if duplicates exist. The example above uses
//    a band width of 1,000 for the search area. The program initially uses
//    ~600 billion for the band width (see IncLkUp[] table), but since the
//    density of both the I^3 + J^3 and I^3 - J^3 pairs decreases with larger
//    numbers, the magnitude of the numbers that define a band width is allowed
//    to expand with time.
//
//    The "NextGroup()" routine generates about 10 to 40 million (depends on
//    MaxHashData) I^3 and J^3 sums / differences for each new band/group.
//    (Stabilizes at this rate after a few iterations.) As these numbers are
//    generated they are inserted into a hash table. (The hash index uses bits
//    16 to 39 (see below) (16 to 40 on the Skulltrail) of the sum / differences
//    of the two cubes.) Actual storage is via link lists that begin at each hash
//    index.
//
//    After a group has been generated, the "CheckGroup()" routine checks the
//    number of entries that have been stored in each of the link lists that
//    begin at each hash index. If there are at least "N" entries in any link
//    list ("N" was entered by the user), then the CheckGroup() routine also
//    checks the link list to see how many of these also match the rest of the
//    bits for the sum-of-cubes / difference-of-cubes. If there are at least "N"
//    of these, then a solution has been found.
//
//    Components within the I^3 + J^3 sums and I^3 - J^3 differences become much
//    larger than any standard computer hardware can handle. Thus these large
//    cubes are separated into sections that can be handled by ordinary 16, 32
//    and 64 bit integer arithmetic. Large numbers are represented as follows:
//
//
//                       <---- Bit ID's within variables ---->
//
//        63                        24 23          0    15           0
//        ------------------------------------------    --------------
//       |            high  bits      |     Hash    |  | Low 16 bits  |
//        ------------------------------------------    --------------
//        79                        40 39         16    15           0
//
//                        <---- Bit ID's within a number ---->
//
//    The above variables are unsigned 64 and 16 bit integers. The Skulltrail
//    uses a larger hash index of 0 to 24 instead of the 0 to 23 shown above.
//
//    The combined 64 + 16 = 80 bits can represent numbers up to 2^80 ~= 1.2E+24
//
//    Note: Only the low 16 bits, the Hash bits, and the next 32 bits are saved
//    for the result so the output results can only go to 2^(16 + 24 + 32) = 2^72
//    = 4.72E21.
//
//    However, the program will crash if/when the search goes past ArrayLimit^3.
//    (ArrayLimit can be increased as needed.)
//
//    The search can begin at any location above 1.0E+14
//
//
//    Programmer's notes:
//
//    The program will run as is if it is compiled and run via
//    the lcc win32 "C" compiler.
//    http://www.cs.virginia.edu/~lcc-win32/
//    http://www.q-software-solutions.de/downloaders/get_name
//
//    The compiled version will require over 400 MB of RAM. If you have more RAM
//    memory, the program will run faster if you increase the numbers in the
//    "#define" statements as suggested below. No other changes to the code are
//    needed.
//
//    The program uses "Global Variables" as they are an efficient memory storage
//    system if the program is a stand-alone system. Global variables are not
//    recommended for larger projects.
//
//    CabTaxi(7) = 11,302,198,488 ~= 1.13E+10
//    CabTaxi(8) = 137,513,849,003,496 ~= 1.37E+14
//    CabTaxi(9) = 424,910,390,480,793,000 ~= 4.249E+17
//    CabTaxi(10) <= 933,528,127,886,302,221,000 <= 9.336E+20
//    CabTaxi(11) <= 261,858,398,098,545,372,249,216 <= 2.62E+23
//
//***************************************************************************



#include <stdheaders.h>            //    The usual stdio.h, stdlib.h, etc

            //    "ArrayLimit" controls how far you want to search. The program
            //    will search for CabTaxis up to ArrayLimit^3. (and then crash
            //    without warning.) Both ArrayLimit and MaxHashData can be
            //    increased as memory allows. (Optional increases in IncLkUp[].)
            //    No other modifications to the code are needed


                                    //    Skulltrail is using 9800000
#define ArrayLimit 6000000          //    This is the Pentium 4 version that
                                    //    only searches to 2.16E20
                                    //    Program will crash if/when

                                    //    search passes ArrayLimit^3
                                    //    Must stay below (2^80)^(1/3)
                                    //    = 106,528,681

                                    //    Increase this number if you have
                                    //    the available RAM.
#define MaxHashData 13000000        //    Skulltrail is using 40000000
                                    //    If possible/convenient, "MaxHashData"
                                    //    will be optimal if it is a little
                                    //    larger than "HashSize"

                                    //    Increase this number if you have
                                    //    the available RAM.
#define HashSize 16777216           //    Skulltrail is using 33554432 (2^25)
#define HashMask (HashSize-1)
#define NbrHashBits 24              //    Programmer must make sure that these
                                    //    are coordinated powers of 2.


                                    //    Procedure declarations
void InitSys(void);                 //    Initialize the program
void NextGroup(void);               //    Generate next band of numbers
void CheckGroup(void);              //    Search for duplicate entries
void Solution(int, int);            //    Called by CheckGroup() to output a
                                    //    solution..
void Bin2Dbl(unsigned long long, unsigned short);    //    Converts 2-piece binary
                                    //    to 2-piece "double"
void UpdtUpperLim(void);            //    Update Upper Limit and Increment
void MakeCube(int);                 //    Cubes the integer and places the
                                    //    results in global variables
                                    //    "High64bits", and "Low16bits"
void CalcCubeDiff(unsigned long long, int);
                                    //    Two integers are passed to the routine.
                                    //    "X" and "dx". The routine calculates
                                    //    the difference of (X + dx)^3  -  X^3
                                    //    and places the result in High64bits and
                                    //    Low16bits.
void CalcIandJ(int);                //    Reconstructs the I and J values that
                                    //    were part of the original I^3 +/- J^3
                                    //    sums/differences
long long GCD(long long, long long);
                                    //    Finds the Greatest Common Divisor of
                                    //    the two parameters and returns this GCD

void PauseMsg(void);                //    Pauses the program at key points



                                    //    NcubedHigh[i] = Bits 16 to 79 of I^3
unsigned long long NcubedHigh[ArrayLimit+1];
unsigned short NcubedLow[ArrayLimit+1];    //    NcubedLow[i] = 16 rightmost bits of
                                    //    I^3. The above arrays are used to look
                                    //    up the cubes of numbers when forming
                                    //    the I^3 + J^3 sums.

int NextI[ArrayLimit+1];            //    Keeps track of the next "I" to try when
                                    //    forming trial I^3 + J^3 sums.
int LowestJ;                        //    Lowest J to use for each call to
                                    //    NextGroup(). "J" is an index into
                                    //    NextI[]

long long NextJforDiag[ArrayLimit+1];    //    Note: The rows and columns in
                                    //    the I^3 - J^3 table will extend beyond
                                    //    2^32. Thus can't use 32 bit integers.
                                    //    Keeps track of where to start for each
                                    //    diagonal when generating the I^3 - J^3
                                    //    differences. Processing usage is
                                    //    similar to that for NextI[] except I^3
                                    //    - J^3 goes down the diagonals.
                                    //    If "I" - "J" = 1, then diagonal = 1
                                    //    If "I" - "J" = 2, then diagonal = 2
                                    //    etc.
unsigned long long NextDiagI3J3[ArrayLimit+1];        //    When I^3-J^3 becomes too
                                    //    big for current Diagonal, save it for
                                    //    the next call to NextGroup()
unsigned short NextDiagLow16[ArrayLimit+1];    //    Same for the low 16 bits

                                    //    The size of the next 2 arrays matches
                                    //    the hash index/mask size. Link lists
                                    //    for the match tests (below) start here.
unsigned HashHd[HashSize];          //    Link heads for hash table.
unsigned short ChainLen[HashSize];  //    Chain length.
                                    //    As trial I^3 + J^3 and I^3 - J^3 numbers
                                    //    are formed for each new "band"/group,
                                    //    pertinent data is stored in these link
                                    //    lists.
                                    //    The following arrays can contain up to
                                    //    13,000,000 elements on the Pentium 4,
                                    //    and 40,000,000 for the Skulltrail
int IorDiag[MaxHashData];           //    For I^3 + J^3 sums, this contains the
                                    //    value for "I". For I^3 - J^3, the value
                                    //    of the diagonal is stored as a negative
                                    //    number. Specifically this will equal
                                    //    J - I.
int Link[MaxHashData];              //    Links for link lists.
unsigned HighBits[MaxHashData];     //    Bits above the hash index are
                                    //    stored here.
unsigned short Lowest16bits[MaxHashData];    //    Low16bits are bits 0 to 15 from
                                    //    the I^3 + J^3 and I^3 - J^3 results.


unsigned long long High64bits;      //    The MakeCube()and CalcCubeDiff()
unsigned short Low16bits;           //    routines store their results here.

long long Ival;                     //    The CalcIandJ() routine reconstructs
long long Jval;                     //    the I and J that were used for I^3 +/-
                                    //    J^3. The results are stored here.

double DblHigh;                     //    2 piece binary numbers are converted to
double DblLow;                      //    decimal numbers such that the true
                                    //    number equals DblHigh * 1.0E9 + DblLow.
                                    //    DblLow will contain the decimal digits
                                    //    for units, thousands, and millions,
                                    //    while DblHigh will contain the decimal
                                    //    digits for "billions" and higher.

long long PairAtI[100];             //    Solutions are stored here
long long PairAtJ[100];

char Databuff[100];                 //    Dummy array for user input.

unsigned NbrPairs;                  //    User input for minimum number of pairs
                                    //    required for output display. (>= 7)

unsigned long long UpperLim;        //    For each new group, find cube sums/
                                    //    differences up to this limit.
                                    //    UpperLim increases by
unsigned long long Increment;       //    Increment for each new group. In turn,
                                    //    "Increment" intermittently increases as
                                    //    array storage allows. "Increment" is
                                    //    the same as the "band width" phrase
                                    //    used in the Excel Spreadsheet example.
                                    //    Note, the lowest 16 bits of the I^3,
                                    //    J^3 calcs are truncated for these
                                    //    calculations.


                        //    The Increment Look Up array is used to give
                        //    "Increment" a good running head start for any
                        //    arbitrary starting point.
                        //    Note: If you have to decrease "MaxHashData" to
                        //    something below 13000000 due to memory
                        //    restrictions, you will also have to decrease the
                        //    numbers in this look-up table
                        //    Optional increase if you increase "MaxHashData"
                        //    Skulltrail is about double the following
double IncLkUp[] = {    //    Multiply these values by 65,536 to get
                        //    the true size for "Increment".
    1.0E7,              //    Start = 1.0E14 to 1.0E15
    2.0E7,              //    Start = 1.0E15 to 1,0E16
    4.5E7,              //    Start = 1.0E16 to 1.0E17
    1.0E8,              //    Start = 1.0E17 to 1,0E18
    2.0E8,              //    Start = 1.0E18 to 1.0E19
    4.5E8,              //    Start = 1.0E19 to 1,0E20
    1.0E9,              //    Start = 1.0E20 to 1.0E21
    2.0E9,              //    Start = 1.0E21 to 1,0E22 (Not used)
    4.0E9,              //    Start = 1.0E22 to 1.0E23 (Not used)
    8.0E9               //    Start = 1.0E23 to 1.0E24 (Not used)
    };

int NbrStored;                      //    Nbr of items in hash table/I3J3 arrays

int StatusFlag;                     //    Set by user for display/no-display of
                                    //    status data while program runs.

int MaxMultiple;                    //    Ask user for the largest multiple
                                    //    to output in the solutions.
                                    //    e.g. Enter "1" for primitives only
                                    //    or 10,000 for all solutions.

int SolNbr;                         //    This is output info that counts the
                                    //    number of "N-way" solutions for the
                                    //    current program run.

time_t ltime;                       //    Used to get the time

double OldTime, NewTime;            //    Used to calculate the search rate
double OldDblHigh, NewDblHigh;





int main(void) {

  InitSys();                        //    Initialize system.

  while(1) {                        //    Do until user stops the program
    NextGroup();                    //    Form next group of cube sums & diff.
    CheckGroup();                   //    Look for matches.
                                    //    Process up to 13,000,000
    UpdtUpperLim();                 //    trial I^3 + J^3 plus I^3 - J^3
                                    //    each iter. (40,000,000 on Skulltrail)
                                    //    Increment increases by 1/128 as
                                    //    array space allows.
    if (StatusFlag) {               //    Optional status check
      Bin2Dbl(UpperLim, 0);
      printf("Upper Limit is now %15.7g\n", DblHigh * 1.0E9);
      Bin2Dbl(Increment, 0);
      printf("Increment is now %15.7g\n", DblHigh * 1.0E9);
      PauseMsg();                   //    Can optionally take this out, but be
    }                               //    ready to use the "Pause" key
  }
  return 0;
}



//****************************************************************************
//
//                                    InitSys
//
//    This routine initializes the system. The Ncubed[] arrays are filled with
//    cubes such that NcubedHigh[i]*(2^16) + NcubedLow[i] = i^3. There is no
//    overflow / precision problem up to 2^80 ~= 1.2089E+24. Alternately, this
//    could allow "ArrayLimit" to have as many as 106,528,681 elements.
//
//    Note: The program allows you to begin your search at any arbitrary
//    starting point. Thus, 2 or more computers/cores could each be running the
//    program to simultaneously search different zones in the number field.
//
//    When the I^3 + J^3 sums and differences are formed, they are split into two
//    pieces. The lowest 16 bits are in an unsigned short variable and the
//    remaining bits are in a 64 bit unsigned long long variable. The lowest 24
//    bits of this long long variable are used as a hash index. (Skulltrail
//    version uses 25 bits.)
//
//    The hash index system greatly speeds up the process of finding matching
//    pairs of I^3 + J^3 = K^3 + L^3, = M^3 - N^3, etc.
//
//    The NextGroup() routine will generate a large group of trial I^3 + J^3
//    sums and I^3 - J^3 differences.
//
//    For the I^3 + J^3 portion of the search, for each column, the search starts
//    at either
//
//    1) row 0
//    or
//    2) where it stopped on the last iteration
//
//    and sequences down until either
//
//    1) the I^3 + J^3 sum exceeds the Upper Limit
//    or
//    2) the downward search reaches the main diagonal.
//
//    This routine initializes the NextI[] array such that the search will
//    "start" at "where it stopped on the last iteration".
//
//    The I^3 - J^3 portion of the search searches sequentially down each
//    diagonal. The first diagonal searched is to the left of the main diagonal,
//    the next one is to the left of the previous, etc. Within each diagonal,
//    each call to NextGroup() starts where it left off last time (remembered in
//    NextJforDiag[]) and continues down the diagonal until the I^3 - J^3
//    difference reaches the current Upper Limit. This routine also initializes
//    the NextJforDiag[] array using the user entered start point.
//
//    The search can start at any number >= 1.0E+14. A good value for "Increment"
//    is found via the IncLkUp[] table for any arbitrary starting point.
//    Subsequently, "UpperLim" will be increased by "Increment" to include larger
//    trial values of I^3 +/- J^3. The density of I^3 +/- J^3 numbers becomes
//    sparser as the number field expands. Therefore, "Increment" is allowed to
//    become larger with time. Thus, the number of trial I^3 + J^3 numbers that
//    will be looked at will cluster just under 13 million per iteration. (Just
//    under 40,000,000 on the Skulltrail.)
//
//**************************************************************************

void InitSys(void) {

  int j, diag;
  double a, b, c, result;
  double start, limit, Jdbl, Icubed;
  long long NextJ;


  puts("              cabtaxi.c");
  puts("                 by");
  puts("             Bill Butler\n");

  puts("This program finds CabTaxi numbers.");
  puts("It will run until you manually stop the program.\n");

  puts("The number that you enter below will determine how many pairs");
  puts("are required to display the associated CabTaxi number.\n");

  puts("Enter 7 for output of 7 or more pairs");
  puts("Enter 8 for output of 8 or more pairs");
  puts("Etc.");

  gets(Databuff);
  NbrPairs = atoi(Databuff);


  do {
    puts("\nEnter a number for the largest multiple of a");
    puts("primitive solution that you want to output.");
    puts("For example: Enter \"1\" for primitives only");
    puts("or Enter \"10000\" to output everything.");
    puts("Hint: Use to control the rate that solutions are displayed");
    gets(Databuff);
    MaxMultiple = atoi(Databuff);
  } while (MaxMultiple < 1);


  limit = ArrayLimit;                       //    Convert to double for calc
  limit *= 0.9999 * limit * limit;

  do {
    printf("\nWhere do you want to start the search?  1.0E+14 <= start <= %g\n",
        limit);
    puts("(Suggest a very small percent below true start point.)");
    gets(Databuff);
    start = atof(Databuff);
  } while ((start < 1.0E+14) || (start > limit));

  limit = 1.0e15;                           //    Look up value for "Increment"
  j = 0;                                    //    given this starting point.
  while (start >= limit) {
    j++;
    limit *= 10.0;
  }
  Increment = IncLkUp[j];                   //    Both convert to long long
  UpperLim = IncLkUp[j] + floor(start/65536.0);
                                            //    Multiply both "Increment" and
                                            //    UpperLim by 2^16 to get true
                                            //    size.

  puts("\nDo you want to display status results?   1 for yes  0 for no");
  puts("(If yes, don't forget that real data may scroll off the screen.)");
  gets(Databuff);
  StatusFlag = atoi(Databuff);

  puts("\nInitializing the cubes table\n");

  for (j = 1; j <= ArrayLimit; j++) {
    MakeCube(j);                            //    Calculates j^3
    NcubedHigh[j] = High64bits;             //    Leftmost 64 bits
    NcubedLow[j] = Low16bits;               //    Rightmost 16 bits
  }
                                            //    Initialize the NextI[] array

  LowestJ = ceil(cbrt(start/2.0));          //    LowestJ marks where
                                            //    to start search on each call to
                                            //    NextGroup().

    //    Overwrite the default "0's" in NextI[] with the calculated next "I's"
    //    until the upward sloping curved line intersects the Row 0 line in the
    //    "Excel" diagram. The search will start at column J = LowestJ, and move
    //    right. The curved line traces: I^3 + J^3 = start. Thus for each "J"
    //    that will be used, we have to calculate the starting "I's" until the
    //    calculated "I" drops below 0.

  j = LowestJ;                              //    This is the index into NextI[]
  Jdbl = j;                                 //    Floating point version for calcs
  while (1) {
    Icubed = start - Jdbl * Jdbl * Jdbl;    //    From: I^3 + J^3 = "start"
    if    (Icubed < 0.0)                    //    Not valid if result will be < 0
      break;
    NextI[j] = ceil(cbrt(Icubed));
    j++;
    Jdbl += 1.0;
  }
                //    Initialize arrays for diagonals
    //    To initialize the search for the I^3 - J^3 table, we have to find the
    //    location of the first element in each diagonal that is greater than the
    //    user entered "start" value. Each entry in the I^3 - J^3 array is
    //    calculated by I^3 - J^3 where the diagonal is I - J = "D".
    //
    //    If we use an algegraic "X" for "J", then for each diagonal, we have to
    //    solve (X + D)^3 - X^3 - "start" = 0 for "X"
    //    and then use the next highest integer
    //
    //    (X + D)^3 - X^3 - start = 0 expands to:
    //    (3)(X^2)(D) + (3)(X)(D^2) + D^3 - start = 0
    //
    //    This can be solved for "X" via the "Quadratic Formula" where:
    //    A = 3D
    //    B = 3(D^2)
    //    C = D^3 - start
    //
    //    Finally we find the integer ceiling for this result.

  puts("Initializing the diagonal arrays\n");

  for (diag = 1; ; diag++) {            //    Sequentially move left until "J"
    a = 3.0 * diag;                     //    reaches column 0.
    b = a * diag;
    c = diag;                           //    Start by converting to a "double"
    c *= c * c;
    c -= start;
    result = (-b + sqrt(b*b - 4*a*c)) / (2*a);    //    Quadratic formula
    if (result < 0.0)                   //    (You should have paid attention)
      break;                            //    (during your high school algebra)
    NextJ = ceil(result);               //    (class)
    NextJforDiag[diag] = NextJ;
    CalcCubeDiff(NextJ, diag);          //    Also set up for first round of
    NextDiagI3J3[diag] = High64bits;    //    I^3 - J^3 differences.
    NextDiagLow16[diag] = Low16bits;
  }
  for ( ;diag <= ArrayLimit; diag++) {  //    Initialize rest of I^3 - J^3 arrays
    NextDiagI3J3[diag] = NcubedHigh[diag];
    NextDiagLow16[diag] = NcubedLow[diag];
  }

  NewTime = (double) clock();           //    Initialize for search rate
  NewTime = NewTime / (double) CLOCKS_PER_SEC;
  NewDblHigh = start / 1.0E9;

  puts("\nThe search started at:");
  time(&ltime);
  printf("%s\n", ctime(&ltime));
}



//*************************************************************************
//
//                                NextGroup
//
//    This routine generates the next group of I^3 + J^3 sums and I^3 - J^3
//    differences, and places these values in the hash arrays. The magnitude of
//    the generated values will be larger than those generated in the last group,
//    but will be limited by the current Upper Limit. The quantity of these trial
//    "I^3 + J^3" sums and I^3 - J^3 differences will average just under 13
//    million for each call to the routine. (Skulltrail version is just under 40
//    million.)
//
//    The diagrams below illustrate how the tables are processed to generate the
//    current group. Both diagrams assume that the last round ended with an old
//    Upper Limit of 501. The "Increment" is assumed to be 500 which produces a
//    new Upper Limit of 1001.
//
//           J    1    2    3    4    5    6    7    8    9   10   11   12
//           J^3  1    8   27   64  125  216  343  512  729 1000 1331 1728
//  I  I^3   -------------------------------------------------------------
//  0    0   |    1    8   27   64  125  216  343   |    |   \|/ 1331 1728
//  1    1   |    2    9   28   65  126  217  344   |    |  1001 1332 1729
//  2    8   |    9   16   35   72  133  224  351   |    |  1008 1339 1736
//  3   27   |   28   35   54   91  152  243  370   |    |  1027 1358 1755
//  4   64   |   65   72   91  128  189  280  407   |    |  1064 1395 1792
//  5  125   |  126  133  152  189  250  341  468   |    |  1125 1456 1853
//  6  216   |  217  224  243  280  341  432   |    |   \|/ 1216 1547 1944
//  7  343   |  344  351  370  407  468  559  \|/  \|/ 1072 1343 1674 2071
//  8  512   |  513  520  539  576  637  728  855 1024 1241 1512 1843 2240
//  9  729   |  730  737  756  793  854  945 1072 1241 1458 1729 2060 2457
// 10 1000   | 1001 1008 1027 1064 1125 1216 1343 1512 1729 2000 2331 2728
// 11 1331   | 1332 1339 1358 1395 1456 1547 1674 1674 1843 2060 2331 2662
// 12 1728   | 1729 1736 1755 1792 1853 1944 2071 2240 2457 2728 3059 3456
//
//    The routine processes the columns in left to right order. The first element
//    in each column that is processed is denoted by the tail of each arrow. If
//    the old Upper Limit was 501, then all values less than 501 were processed
//    in earlier rounds. In the current call to the routine, each column is
//    processed downward until the new Upper Limit of 1001 is reached or until
//    the main diagonal is reached. When processing reaches this limit, the
//    number for the next lower row is saved in NextI[] to show where processing
//    should begin on the subsequent call to the routine.
//
//    The second table is similar to the table above except the values are the
//    differences of two cubes. Specifically, the values in the table represent
//    I^3 - J^3. Here, processing is along each diagonal from upper left to lower
//    right. The rightmost diagonal is processed first followed by the other
//    diagonals in right to left order. Diagonals are identified by a small
//    integer (the ID number of the diagonal - which is equal I - J) instead of
//    the arrows shown in the first diagram. Again, the elements that are
//    processed are those that are >= the old Upper Limit of 501 and less than
//    the new Upper Limit of 1001.
//
//           J    0    1    2    3    4    5    6    7    8    9   10   11   12
//           J^3  0    1    8   27   64  125  216  343  512  729 1000 1331 1728
//  I  I^3   ------------------------------------------------------------------
//  0    0   |    0
//  1    1   |    1    0
//  2    8   |    8    7    0
//  3   27   |   27   26   19    0
//  4   64   |   64   63   56   37    0
//  5  125   |  125  124  117   98   61    0
//  6  216   |  216  215  208  189  152   91    0
//  7  343   |  343  342  335  316  279  218  127    0
//  8  512   |   8    7    6   485  448  387  296  169    0
//  9  729   |   9    8    7    6    5    4    3   386  217    0
// 10 1000   |  10    9    8    7    6    5    4    3   488  271    0
// 11 1331   | 1331 1330 1323 1304 1267 1206 1115   4    3    2   331    0
// 12 1728   | 1728 1727 1720 1701 1664 1603 1512 1385 1216   3    2   397   0
//                                                        (Long tail forms here)
//    Note: The long tail will extend past 2^32. Thus long long integers are used.
//
//    Each result that is used for the current round is stored in the hash arrays.
//    The I^3 + J^3 or I^3 - J^3 result is split up as follows. Bits 0 to 15 are
//    saved in the Lowest16bits[] array for later matching. Bits 16 to 39 are
//    used for the hash index. (The Skulltrail version uses bits 16 to 40 for the
//    hash index.) Bits 40 to 71 are saved in HighBits[]. For the I^3 + J^3
//    calculations, "I" is saved as the first part of an I/J pair. For the
//    I^3 - J^3 portion, the diagonal (equals I - J) is saved as a negative
//    value.
//
//*************************************************************************

void NextGroup(void) {

  int Count, i, j, diag;
  long long LLj;
  unsigned long long i3j3high;
  unsigned work32;


  Count = 0;    //    Count nbr. of I^3, J^3 results in the group and where to
                //    put data. Note: Usage of the "i" & "j" vars. in this loop
                //    matches the "I" & "J" row/col labels in the "Excel
                //    Spreadsheet" examples.

    //    The routine processes I^3 + J^3 numbers first and then I^3 - J^3
    //    numbers

    //    The main loop for I^3 + J^3 portion has been split into two parts. In
    //    the first portion, the iterative step down each column in the "Excel
    //    spreadsheet" will run into the "main diagonal" which will terminate
    //    operations for this column. After the first few columns, the only exit
    //    condition will be the "if (i3j3high < UpperLim)" test. Thus the
    //    second portion of the code only tests this condition (slight reduction
    //    in code size speeds up the search).

  j = LowestJ;                            //    Start at this J index

  do {      //    Start "J" loop. Each "j" increment moves 1 col to right in the
            //    Excel table. Will process all "i" for this col. until either
            //    i > j (reach the diagonal) or reach UpperLim. This first
            //    portion of the main loop will normally execute just a few times.

    for (i = NextI[j]; i <= j; i++) {
      work32 = NcubedLow[i] + NcubedLow[j]; //    Start sum by adding low bits
      i3j3high = (work32 >> 16);            //    Forms "Carry"
      i3j3high += NcubedHigh[i] + NcubedHigh[j];    //    Then add high cubes
                                            //    If within upper limit then
      if (i3j3high < UpperLim) {            //    include it in the group
        Lowest16bits[++Count] = work32;     //    Will use these for matches
        HighBits[Count] = (i3j3high >> NbrHashBits);
        IorDiag[Count] = i;                 //    Save the I that was used
        work32 = i3j3high;                  //    Construct hash index
        work32 &= HashMask;                 //    Just use low 24 bits
        Link[Count] = HashHd[work32];       //    Update link list
        HashHd[work32] = Count;
        ChainLen[work32]++;                 //    Increment count for link list
      }
      else                                  //    Repeat until too big
        break;
    }
        //    If loop ended by "i" test ("i" increased until it dropped below the
        //    diagonal in the "spreadsheet"), then "i" will equal "j" + 1. Else
        //    exit from loop was via "break" (sum of cubes exceeded upper limit),
        //    and "i" will be <= "j".

    NextI[j] = i;                       //    Save for next group.

    j++;                                //    Move right 1 col in spreadsheet.
    if (i >= j)                         //    If won't use this "j" again.
      LowestJ = j;

  } while ((i + 2) > j);                //    Loop exits after the first couple
                                        //    of columns have been processed.
                                        //    "i + 2" is needed to prevent
                                        //    overshoot problem.

        //    When processing gets to here, then "i" did not reach the spreadsheet
        //    diagonal. Thus from here, the only test required will be that I^3
        //    plus J^3 reaches the upper limit.

  do {                                  //    Continue with the other columns
                                        //    until I^3 + J^3 is > UpperLimit
                                        //    even with "i" = 0. This portion of
                                        //    the loop will execute millions of
                                        //    times.
    i = NextI[j];
    work32 = NcubedLow[i] + NcubedLow[j];     //    Start sum by adding low bits
    i3j3high = (work32 >> 16);                //    Forms "Carry"
    i3j3high += NcubedHigh[i] + NcubedHigh[j];//    Then add high cubes

                                        //    Process all I^3 + J^3 in this
                                        //    column until too big, then move
                                        //    1 column to the right.
    while (i3j3high < UpperLim) {
      Lowest16bits[++Count] = work32;       //    Will use these for matches
      HighBits[Count] = (i3j3high >> NbrHashBits);
      IorDiag[Count] = i;                   //    Save the I that was used.
      work32 = i3j3high;                    //    Construct hash index
      work32 &= HashMask;                   //    Just use low 24 bits
      Link[Count] = HashHd[work32];         //    Update link list
      HashHd[work32] = Count;
      ChainLen[work32]++;                   //    Increment count for link list
                                            //    Form I^3 + J^3 for next "I"
      i++;
      work32 = NcubedLow[i] + NcubedLow[j];  //    Start sum by adding low bits
      i3j3high = (work32 >> 16);             //    Forms "Carry"
      i3j3high += NcubedHigh[i] + NcubedHigh[j];    //    Then add cubes
    }                                        //    Repeat until too big


    NextI[j] = i;                       //    Save for next group.

    j++;                                //    Move right 1 col in spreadheet.

  } while (i);                          //    Loop exits when right edge of
                                        //    band/group intercepts
                                        //    Row = 0 in the spreadsheet.



  //    Now process the I^3 - J^3 portion

  diag = 1;                                 //    Init for first diagonal
  do {                                      //    For all active diagonals
    LLj = NextJforDiag[diag];               //    Get next "J" for this diagonal
    i3j3high = NextDiagI3J3[diag];          //    Get the i3j3 & Low16bits that
    Low16bits = NextDiagLow16[diag];        //    were too big on the last round
                                            //    Move down diagonal until I3J3
    while (i3j3high < UpperLim) {           //    is too big
      Lowest16bits[++Count] = Low16bits;    //    Will use these for matches
      HighBits[Count] = (i3j3high >> NbrHashBits);
      IorDiag[Count] = -diag;               //    Save the diagonal as a negative
      work32 = i3j3high;                    //    Construct hash index
      work32 &= HashMask;
      Link[Count] = HashHd[work32];         //    Update link list
      HashHd[work32] = Count;
      ChainLen[work32]++;                   //    Increment count for link list
                                            //    Calculate next trial i3j3
      LLj++;                                //    Move down diagonal
      CalcCubeDiff(LLj, diag);              //    Find difference of cubes
      i3j3high = High64bits;                //    Get high 64 bits
    }                                       //    Loop continues down diagonal
                                            //    until i3j3 is larger than U. L.
    NextJforDiag[diag] = LLj;               //    When too big, save data for
    NextDiagI3J3[diag] = i3j3high;          //    next call to routine
    NextDiagLow16[diag] = Low16bits;
    diag++;                                 //    Move left to next diagonal
  } while (LLj);                            //    Repeat until at left edge

  NbrStored = Count;
                                        //    Optional status check.
                                        //    Averages just under 13,000,000 per
                                        //    crack. (40,000,000 for Skulltrail)
  if (StatusFlag) {
    printf("This iter. stored %'d trial I^3 +/- J^3 sums\n", NbrStored);
    printf("Must stay < %'d\n", MaxHashData);
  }
}



//************************************************************************
//
//                                CheckGroup
//
//    This routine checks the hash arrays to see if any matches exist. If the
//    quantity of numbers stored in any hash link list is >= the user defined
//    display number, an actual CabTaxi solution may exist. A count is then
//    made to see how many exact matches exist. If the number of exact matches
//    is >= "NbrPairs" (entered by the user) then a solution exists and the
//    Solution() routine is called to display it.
//
//    Statistician's note: The number of elements in the link list chains
//    (= ChainLen[HashIndex]) will closely approximate "Poisson Distribution"
//    where the mean will equal "NbrStored" / "HashSize". In turn "NbrStored"
//    will average 3 to 4 percent under "MaxHashData".
//    (The vast majority of the time, the portion of the code inside the
//    "if (ChainLen[HashIndex] >= LocalPairs)" test will not execute.)
//
//************************************************************************

void CheckGroup(void) {

  int Link1, Link2, Count, StopCount;
  int register HashIndex, LocalPairs;

  LocalPairs = NbrPairs;                    //    Setting this up as a local
                                            //    variable may speed up program

                                            //    For all hash heads
  for (HashIndex = HashMask; HashIndex >= 0; HashIndex--) {

    if (ChainLen[HashIndex] >= LocalPairs) {//    If chain is long enough,
                                            //    check for a possible solution.
//      if (ChainLen[HashIndex] > 50)       //    See if hashing is OK
//        printf("ChainLen[%'d] = %d\n", HashIndex, ChainLen[HashIndex]);

                                            //    Check the list for this hash.
                                            //    StopCount controls how many times
                                            //    that "Link1" moves forward
      StopCount = ChainLen[HashIndex] - LocalPairs;

      for (Link1 = HashHd[HashIndex]; StopCount >= 0;
                StopCount--, Link1 = Link[Link1]) {

        Count = 1;                          //    Will count number of high and
                                            //    low matches
        for (Link2 = Link[Link1]; Link2; Link2 = Link[Link2]) {
          if (Lowest16bits[Link1] != Lowest16bits[Link2])
            continue;
          if (HighBits[Link1] == HighBits[Link2])
            Count++;                        //    Increment count if complete match
        }

        if (Count >= LocalPairs)            //    If true, then solution exists
          Solution(HashIndex, Link1);
      }                                     //    End of search link list loop
    }                                       //    End of check ChainLen[] test

    HashHd[HashIndex] = 0;                  //    Clear garbage for next round
    ChainLen[HashIndex] = 0;
  }
}




//*************************************************************************
//
//                                    Solution
//
//    This routine is called when a solution exists starting at the data
//    referenced by the passed parameters.
//
//*************************************************************************

void Solution(int HashIndex, int Link1) {

  int i, Count, Link2, GCDivInt;
  long long GCDiv, work;


                                            //    Calculate the Cabtaxi number
  work = HighBits[Link1];                   //    Reconstruct all the high bits
  work <<= NbrHashBits;
  work += HashIndex;
  Bin2Dbl(work, Lowest16bits[Link1]);       //    Convert binary to "English"

  Count = 1;                                //    Count number of matches.

  CalcIandJ(IorDiag[Link1]);                //    Reconstruct the original I & J
  PairAtI[Count] = Ival;                    //    I for the first pair
  PairAtJ[Count] = Jval;                    //    J for the first pair

  for (Link2 = Link[Link1]; Link2; Link2 = Link[Link2]) {
                                            //    Must match both sections
    if (Lowest16bits[Link2] != Lowest16bits[Link1])
      continue;
    if (HighBits[Link2] == HighBits[Link1]) {
      Count++;
      CalcIandJ(IorDiag[Link2]);            //    Reconstruct the original I and
      PairAtI[Count] = Ival;                //    J for the other pairs
      PairAtJ[Count] = Jval;
    }
  }
                                            //    Get Greatest Common Divisor
  GCDiv = GCD(PairAtI[Count], PairAtJ[Count]);
  for (i = Count - 1; i; i--) {
    GCDiv = GCD(GCDiv, PairAtI[i]);
    GCDiv = GCD(GCDiv, PairAtJ[i]);
  }
  GCDivInt = GCDiv;                         //    Convert back to simple integer
                                            //    If Greatest Common Divisor is
                                            //    <= the user-entered multiple
  if (GCDivInt <= MaxMultiple) {            //    size, then output the result.
    if (Count >= 9)                         //    "Wake up" something significant
      puts("\n***************   AT LEAST A 9-WAY SOLUTION   ****************");
    SolNbr++;
    printf("\nSolution Number %d\n", SolNbr);

    printf("%d pairs exist at: %'.0lf * E9 + %'.0lf\n", Count, DblHigh, DblLow);
    for (i = Count; i ; i--)                //    Display all pairs
      printf("     Pair at I = %16'lld     J = %16'lld\n", PairAtI[i], PairAtJ[i]);
    printf("The Greatest Common Divisor for this solution is: %'8d\n", GCDivInt);

    time(&ltime);
    printf("%s", ctime(&ltime));

    OldTime = NewTime;                      //    Calculate and output
    NewTime = (double) clock();             //    the search rate
    NewTime = NewTime / (double) CLOCKS_PER_SEC;
    OldDblHigh = NewDblHigh;
    NewDblHigh = DblHigh;
    if (NewTime > (OldTime + 1.0))
      printf("Search rate = %g per second\n\n",
             1.0E9 * (NewDblHigh - OldDblHigh) / (NewTime - OldTime));

  }
}



//*************************************************************************
//
//                                    MakeCube
//
//    This routine cubes the passed integer number and places the lowest 16
//    binary bits of the result in the global integer variable "Low16bits"
//    and the remaing bits into High64bits (long long). The InitSys() routine
//    loads these results into the Ncubed arrays.
//
//    The multiplication process is carried out via 32 bit and 64 bit binary
//    arithmetic with a final split of the lowest 16 bits into "Low16bits",
//    and the remainder into "High64bits".
//
//    The passed number "AnInteger" is known to be < 2^32. This is converted
//    to a "long long" number which can be safely squared without any precision
//    problems. This result is then split into two seperate 32-bit sections that
//    are again converted to "long long" binary numbers. Each of these can be
//    safely multiplied by the original "AnInteger" again without any precision
//    error.
//
//    Finally, the lowest 16 bits are placed in the global variable "Low16bits",
//    and everything else is placed in "High64bits".
//
//**************************************************************************

void MakeCube(int AnInteger) {

  unsigned long long UpperHalf, LowerHalf;
  unsigned work;
                                        //    "AnInteger" will be < 2^32. Thus
                                        //    there is enough precision in an
                                        //    ordinary "long long" variable for
                                        //    the first multiplication.
  UpperHalf = AnInteger;                //    Convert to "long long"
  UpperHalf *= UpperHalf;               //    Square it
                                        //    Split it into two pieces
  work = UpperHalf;                     //    Get lowest 32 bits
  LowerHalf = work;                     //    Lower Half has the lower 32 bits
  UpperHalf >>= 32;                     //    Only the upper 32 bits remain
  UpperHalf *= AnInteger;               //    Form cube of AnInteger
  LowerHalf *= AnInteger;
  Low16bits = LowerHalf;                //    Only uses last 16 bits
  LowerHalf >>= 16;                     //    Get rid of lowest 16 binary bits
  UpperHalf <<= 16;                     //    Align upper and lower bits
  UpperHalf += LowerHalf;               //    Add the two sections
  High64bits = UpperHalf;               //    Store final result

}




//**************************************************************************
//
//                                CalcCubeDiff
//
//    This routine calculates the difference between two cubes and places the
//    results in global variables High64bits and Low16bits. A long long integer
//    "X" and a regular integer "dx" are passed to the routine. The routine
//    calculates (X + dx)^3 - X^3 and places the 64 high bits of the result in
//    High64bits and the low 16 bits in Low16bits.
//
//    For precision purposes, it is known that the result will be less than
//    1.0E24 ~= 2^80. Thus all component arithmetic will be OK as long as the
//    max size stays within 80 bits.
//
//    Instead of trying to calculate (X + dx)^3 - X^3 directly, the routine uses
//    the following algebraic expansion:
//
//    (X + dx)^3 - X^3
//    = 3 * X^2 * dx  +  3 * X * dx^2 + dx^3
//    = (3 * X * dx) * (X + dx) + dx^3
//
//**************************************************************************

void CalcCubeDiff(unsigned long long X, int dx) {

  unsigned long long SumHigh, SumLow;
  unsigned long long Work;

                                        //    First, calculate (3 * X * dx)
  SumHigh = X;                          //    Maximum of 35 bits
  SumHigh *= (3 * dx);                  //    Maximum of 35 + 2 + 24 = 61 bits
  SumLow = SumHigh & 65535;             //    Get low 16 bits
  SumHigh >>= 16;                       //    Upper bits. SumHigh can never
                                        //    overflow since it is known that the
                                        //    result will be <= 80 bits, and
                                        //    SumHigh doesn't have the low 16 bits
  Work = X + dx;                        //    Next term (Maximum of 35 bits)
  SumHigh *= Work;                      //    Multiply both by (X + dx)
  SumLow *= Work;                       //    Maximum of 16 + 35 bits
  SumHigh += NcubedHigh[dx];            //    Add both parts of dx^3
  SumLow += NcubedLow[dx];
  Low16bits = SumLow;                   //    Get the low 16 bits
  SumLow >>= 16;                        //    Get rid of low 16 bits
  SumHigh += SumLow;                    //    Add a big carry
  High64bits = SumHigh;                 //    Save rest of I^3 - J^3

}



//***************************************************************************
//
//                                CalcIandJ
//
//       In the NextGroup() routine, instead of saving both the I and J values
//    (and both would have to be long long variables), the routine only stored
//    I or the diagonal involved as normal integers. (If the diagonal was saved,
//    it was saved as a negative value.) This saves a large amount of array space.
//    However, this requires that the original "I" and "J" values must be
//    reconstructed on the rare occasions when a solution is found.
//
//    If the value in IorDiag[] was a positive value, then the original "I" was
//    saved. In this case, the value for I will be less than 2^24. Thus to
//    reconstruct the value for J, we use the equation I^3 + J^3 = solution
//    value. "I" as stated above, was saved in IorDiag[] and the solution value
//    has already been calculated as floating point numbers which have been
//    stored in global variables "DblHigh" and "DblLow".
//
//    If the value in IorDiag[] is negative, then the diagonal was saved as a
//    negative value. Thus the original values for I and J in
//    I^3 - J^3 = solution-value have to be calculated.
//
//    In this case the value of the diagonal (as a positive value) equals I - J.
//
//    The original equation
//    I^3 - J^3 = solution value
//    expands to:
//    (J + diagonal)^3 - J^3 - solution value = 0
//    which becomes:
//    (3)(J^2)(diagonal) + (3)(J)(diagonal^2) + diagonal^3 - solution value = 0
//    which can be solved for "J" via the quadratic formula.
//    a = 3*diagonal
//    b = 3*diagonal^2
//    c = diagonal^3 - solution value
//    J = (-b+sqrt(b*b - 4*a*c))/(2*a)
//    and "I" will equal J + diagonal
//
//    In both of the above cases, the resulting values for I and J are stored in
//    Ival and Jval where they can be accessed by the solution routine.
//
//****************************************************************************

void CalcIandJ(int IorDiagonal) {

  double Icubed, Jcubed, DblIorDiag;
  double DblJ, a, b, c;

  DblIorDiag = IorDiagonal;                 //    Convert to a double

                                            //    If IorDiagonal is positive,
  if (DblIorDiag >= 0.0) {                  //    then get I and J for I^3 + J^3
    Icubed = DblIorDiag * DblIorDiag * DblIorDiag;
                                            //    Then use: solution = I^3 + J^3
    Jcubed = DblHigh*1.0E9 + DblLow - Icubed;
    Jval = round(cbrt(Jcubed));             //    Saves result and converts to LL
    Ival = IorDiagonal;                     //    Also converts to long long
  }
  else {                                    //    Else the diagonal was saved
    DblIorDiag = -DblIorDiag;               //    Convert to a positive double
    a = 3.0 * DblIorDiag;
    b = a * DblIorDiag;
    c = DblIorDiag * DblIorDiag * DblIorDiag - 1.0E9 * DblHigh - DblLow;
    DblJ = (-b + sqrt(b*b - 4.0*a*c)) / (2.0*a);
    Jval = -round(DblJ);                    //    Save as a negative value
    Ival = round(DblJ + DblIorDiag);
  }
}




//************************************************************************
//
//                            Bin2Dbl
//
//    This routine converts the 2-piece binary number that is passed to the
//    routine to a 2-piece floating point number (each result piece is double).
//    The two unsigned integers passed to the routine are components that
//    would be a complete binary number represented by the following:
//    Binary number = BigInteger * (2^16) + SmallInt.
//    It is known that BigInteger is a 64 bit unsigned long long and
//    SmallInt is an unsigned short. The whole purpose is to end up with
//    something recognizable as a decimal number.
//
//************************************************************************

void Bin2Dbl(unsigned long long BigInteger, unsigned short SmallInt) {

  double BilResult, UnitResult;
  double BilIncr, UnitIncr;
  unsigned long long work;
                                //    First get the decimal value of the
                                //    small integer
  BilResult = 0.0;
  UnitResult = SmallInt;        //    SmallInt is < 2^16.
                                //    Initialize increments for decimal numbers
  BilIncr = 0.0;                //    These two are the decimal equivalent
  UnitIncr = 65536.0;           //    of 2^16.
  work = BigInteger;

  while (work) {

    if (work & 1) {                //    If rightmost binary digit in work is a 1
      BilResult += BilIncr;        //    then add in the current increment
      UnitResult += UnitIncr;
      if (UnitResult >= 1.0E9) {   //    Check for carry
        UnitResult -= 1.0E9;
        BilResult += 1.0;
      }
    }
    UnitIncr += UnitIncr;          //    Double the increment values for the next
    BilIncr += BilIncr;            //    round.
    if (UnitIncr >= 1.0E9) {
      UnitIncr -= 1.0E9;
      BilIncr += 1.0;
    }
    work >>= 1;                    //    Shift next bit into position.
  }
  DblHigh = BilResult;             //    Put result in global variables
  DblLow = UnitResult;
}




//**********************************************************************
//
//                                UpdtUpperLim
//
//    After each group of cube pairs has been generated and checked for
//    duplicates, the number field is expanded so that a new group of paired
//    cubes may be processed. The new group of candidates will exist in the
//    number range bounded by "Old Upper Limit" to "New Upper Limit" where
//    "New Upper Limit" = "Old Upper Limit" plus the current "Increment". The
//    current Increment is allowed to increase with time so that the number of
//    trial sums that are generated will stabilize at just under 13,000,000
//    new candidates per iteration. (Just under 40,000,000 for Skulltrail.)
//
//    Note: If UpperLimit and Increment were treated as full 80 bit binary
//    numbers, the lowest 16 bits would have to be included. Since the program
//    only deals with large numbers, the lowest 16 bits are truncated for the
//    Upper Limit tests. As long as "Increment" is >= 2^7 (and it is initialized
//    to much bigger than this), the incrementing process will proceed normally.
//
//**********************************************************************

void UpdtUpperLim(void) {

  unsigned long long WorkBits;

                            //    If the number of trial cube pairs stored on
                            //    the last iteration was significantly less than
                            //    the maximum that could be used without a
                            //    MaxHashData overflow, then the increment can
                            //    be increased by about 1/128 so that the next
                            //    iteration will generate a larger group size.
  if (NbrStored < (MaxHashData - MaxHashData/25)) {
    WorkBits = Increment;
    WorkBits >>= 7;         //    Get 1/128 of the former Increment
    Increment += WorkBits;
  }
                            //    Update UpperLimit with increment
  UpperLim += Increment;
}



//******************************************************************
//
//                    Misc. Routines
//
//******************************************************************

    /*  This routine finds the greatest common divisor.         */

long long GCD(long long IntOne, long long IntTwo) {

  long long larger, smaller, Remainder;

  if (IntOne < 0)
    IntOne = -IntOne;

  if (IntTwo < 0)
    IntTwo = -IntTwo;

  if (IntOne > IntTwo) {
    larger = IntOne;
    smaller = IntTwo;
  }
  else {
    smaller = IntOne;
    larger = IntTwo;
  }

  if (!smaller)
    return(larger);

  do {
    Remainder = larger - smaller * (larger/smaller);
    larger = smaller;
    smaller = Remainder;
  } while (smaller);
  return(larger);
}


void PauseMsg(void) {               //    Pause before displaying
                                    //    CRT data.
  puts("\nPress RETURN to continue");
  gets(Databuff);
}