The La Plata Mountains as seen from above the author’s
            home.


Durango Bill’s

“C” Program to generate Cabtaxi Numbers
(Source Code)



Return to the main Ramanujan page

Web page generated via Sea Monkey's Composer HTML editor
within  a Linux Cinnamon Mint 18 operating system.
(Goodbye Microsoft)


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