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(<ime);
printf("%s\n", ctime(<ime));
}
//*************************************************************************
//
// 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(<ime);
printf("%s", ctime(<ime));
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);
}