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 a band 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