Durango Bill’s
“C” Program to generate Ramanujan 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)
//*************************************************************************
//
// Ramanujans.c
// by
// Bill Butler
//
// This program finds Ramanujan numbers. A Ramanujan number is a number
// formed by the sum of two cubes in 2 or more different ways. For example:
//
// 12^3 + 1^3 = 9^3 + 10^3 = 1729
//
// There are an infinite number of other paired cubes that have a common sum.
//
// Also, the number of pairs can be extended to higher orders.
//
// For example, the first quadruple, Taxicab(4), occurs at:
//
// 13,322^3 + 16,630^3 =
// 10,200^3 + 18,072^3 =
// 5,436^3 + 18,948^3 =
// 2,421^3 + 19,083^3 = 6,963,472,309,248 ( 6963472309248 )
// (Takes about 30 seconds using this program on a 3GHz processor)
// (Also nails Taxicab(5) in less than 3 1/4 hrs.)
//
// 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 = 1, 2, 3, 4, etc. and the 2nd column gives the
// cube for each of these numbers: 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 13
// J^3 1 8 27 64 125 216 343 512 729 1000 1331 1728 2197
// I I^3 ------------------------------------------------------------------
// 1 1 | 2 9 28 65 126 217 344 513 730 1001 1332 1729 2198
// 2 8 | 9 16 35 72 133 224 351 520 737 1008 1339 1736 2205
// 3 27 | 28 35 54 91 152 243 370 539 756 1027 1358 1755 2224
// 4 64 | 65 72 91 128 189 280 407 576 793 1064 1395 1792 2261
// 5 125 | 126 133 152 189 250 341 468 637 854 1125 1456 1853 2322
// 6 216 | 217 224 243 280 341 432 559 728 945 1216 1547 1944 2413
// 7 343 | 344 351 370 407 468 559 686 855 1072 1343 1674 2071 2540
// 8 512 | 513 520 539 576 637 728 855 1024 1241 1512 1843 2240 2709
// 9 729 | 730 737 756 793 854 945 1072 1241 1458 1729 2060 2457 2926
// 10 1000 | 1001 1008 1027 1064 1125 1216 1343 1512 1729 2000 2331 2728 3197
// 11 1331 | 1332 1339 1358 1395 1456 1547 1674 1674 1843 2060 2331 2662 3059
// 12 1728 | 1729 1736 1755 1792 1853 1944 2071 2240 2457 2728 3059 3456 3925
// 13 2197 | 2198 2205 2224 2261 2322 2413 2540 2709 2926 3197 3528 3925 4394
//
// The table can be extended for as many rows and columns as needed. The lower
// left triangle of the table is a mirror image of the upper right triangle.
// Thus, the lower left triangle can be ignored.
//
// Within the upper right triangle, we want to look for two identical numbers.
// After some searching, we note that the number 1729 appears twice. It is
// found at row 1, column 12; and again at row 9, column 10. Thus 1^3 + 12^3
// equals 9^3 + 10^3 equals 1729. Hence, 1729 is a Ramanujan number. If we
// infinitely increased the size of the table, we could find an infinite
// number of other paired cube sums. Also, we could find an infinite number of
// "triples". The bad news is that the first of these "Ramanujan triples"
// occurs at:
// (Row 228, Column 423), (Row 167, Column 436), and (Row 255, Column 414)
// for 228^3 + 423^3 = 167^3 + 436^3 = 255^3 + 414^3, which all share a
// common cell entry of 87,539,319.
//
// If we are to extend the process of finding additional Ramanujan triples or
// even higher order matches (e.g. Ramanujan quadruples), three major problems
// become apparent. 1) The table 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 table are confined to specific areas. 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.
//
// The program sequentially generates these bands. Each new band is used for
// the current search process. Then this old band is discarded and a new band
// is generated. 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 10 billion for
// the band width (for 4 & 5 pair combinations), but since the density of the
// I^3 + J^3 sums decreases with larger numbers, this band width is allowed
// to expand with time.
//
// The program generates a little over 2.4 million numbers at a time for each
// new band. (Stabilizes at this rate after a few dozen iterations.) These
// numbers are inserted into a hash table, and then the hash table is searched
// to see if duplicate entries exist. If there are enough entries for any hash
// number, then the associated link list of results for this hash number are
// searched for actual solutions.
//
// The algorithm (especially the hash index portion) is very efficient and
// appears to be at least 10 times faster than the "heap" algorithm used by
// David Wilson who published the first Ramanujan quintuple.
//
// If the search is continued to include the lowest known Ramanujan sextuple,
// (= Taxicab(6)), numbers in the table would have 23 decimal digits. No
// simple data representation (or standard computer hardware) has this much
// precision. Thus numbers are split into two pieces. One part includes the
// rightmost 9 decimal digits and is maintained as "int" variables. The
// remaining leftmost 14 decimal digits are maintained as "double" variables.
// The integer portion does double duty as it is used to construct the hash
// indexes.
//
// The program has been updated from an old ANSI "C" version that ran under
// DOS on a 80486 computer. (And that was a copy of a version that ran on a
// 80386 computer, and that was a copy of the original that I wrote for a
// 32032 32-bit National Semiconductor co-processor board that I added to an
// IBM PC-XT back in the early 80's. (Those were the days when you sold your
// soul for 1 MB of direct RAM addressing - no PC/DOS segmentation. Anybody
// remember what "EDLIN" was?)
//
// This program may be used, copied, modified, etc. without any obligation by
// any person for any not-for-profit purpose. I would appreciate that this or
// any derivative version would include a note crediting me with the original
// program/algorithm. (e.g. "Original algorithm by Bill Butler.")
//
// Final note: The program runs best if you have >= 1 GB of RAM. It runs as is
// under Windows XP when compiled by the lcc-win32 C compiler as a "console"
// program.
//
//***************************************************************************
#include <stdheaders.h> // The usual stdio.h, stdlib.h, etc
// "ArrayLimit" controls how far you want to search. The program
// will search for Ramanujans up to ArrayLimit^3. (No other
// modifications to the code are needed.)
// DO NOT use anything > 30000000.
// At 30,000,000 for ArrayLimit, the author was able to run time
// trials up to 2.15E22, and successfully run a test to see if it
// would find the best known candidate for Taxicab(6). (Started a
// short distance lower than this.) However, when the program was
// running with "ArrayLimit" at 30,000,000 at the same time that
// several other programs were running, the author's computer
// (1 GB RAM) crashed and corrupted Norton's "Go Back" files.
// Norton's "System Works" had to be uninstalled and then
// reinstalled. (System "Commit Charge" was not known, but an
// overload is possible.)
#define ArrayLimit 10000000 // Program will crash if/when search passes
// ArrayLimit^3
// Procedure declarations
void InitSys(void); // Initialize the program
void NextGroup(void); // Generate next band of numbers
void CheckGroup(void); // Search for duplicate entries
void MakeCube(int); // Cubes the integer and places the results
// in global variables CubedHigh and CubedLow
void pausemsg(void); // Pauses the program at key points
// Note: Each 10000000 increase in ArrayLimit
// uses another 160 MB of RAM. ArrayLimit of
// 30000000 is good enough for Taxicab(6).
// (But estimated run time is 7 years.)
double NcubedHigh[ArrayLimit+1];// NcubedHigh[i] = leftmost 14 digits of i^3
int NcubedLow[ArrayLimit+1]; // NcubedLow[i] = 9 rightmost digits of i^3
// These arrays are used to look up the cubes
// of numbers when forming the I^3 + J^3 sums
int NextJ[ArrayLimit+1]; // Keeps track of the next "J" to try when
// forming trial I^3 + J^3.
// The size of the next 2 arrays matches
// the 22 bit hash size. Link lists for
// the data arrays below start here.
unsigned HashHd[4194304]; // Link heads for hash table.
int ChainLen[4194304]; // Chain length.
// As trial I^3 + J^3 numbers are generated
// for each new "band"/group, they are
// stored in these link lists. (Purists might
// want to use a structure array, but time
// trials for both encoding methods showed
// that simple arrays execute faster.)
double I3J3High[3000000]; // Leftmost 14 digits of I^3 + J^3
int I3J3Low[3000000]; // 9 rightmost digits of I^3 + J^3
int Ival[3000000]; // The "I" in I^3 + J^3
int Jval[3000000]; // The "J" in I^3 + J^3
int Link[3000000]; // Link lists.
double CubedHigh; // The MakeCube() routine is passed an integer
int CubedLow; // number. It cubes this number and places the
// lowest (rightmost) 9 decimal digits in
// CubedLow and the remaining digits in
// CubedHigh. These are then copied into the
// Ncubed arrays.
char Databuff[100]; // Dummy array for user input.
int NbrPairs; // User input for minimum number of pairs
// required for output display. (=2,3,4,or 5)
double InitValues[10] = { // Initial upper limit & increment. Will keep
0.0, 0.0, 2000.0, // output display in roughly sorted order.
20000000.0,
10000000000.0,
10000000000.0 };
double UpperLim; // For each new group, find cube sums up to
// this limit. UpperLim increases by
double 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.
int NbrStored; // Nbr. of items in hash table/I3J3 arrays.
unsigned Bit22Mask = 4194303; // Mask for 22 bit index into Hash Table
int StatusFlag; // Set by user for display/no-display of
// status data while program runs.
int main(void) {
InitSys(); // Initialize system.
while(1) { // Do until user stops the program
NextGroup(); // Form next group of cube sums.
CheckGroup(); // Look for matches.
// Process a little over 2,400,000
if (NbrStored < 2400000) // trial I^3 + J^3 each iter.
Increment *= 1.05; // Increases by 5 % as array space allows
UpperLim += Increment; // Set up for next group
if (StatusFlag) { // Optional status check
printf("Upper Limit is now %g\n", UpperLim);
printf("Increment is now %g\n", Increment);
}
}
return 0;
}
//****************************************************************************
//
// InitSys
//
// This routine initializes the system. The Ncubed[] arrays are filled with
// cubes such that NcubedHigh[i]*E9 + NcubedLow[i] = I^3. 23+ digits of
// precision can be handled.
//
// Note: The program allows you to begin your search at any arbitrary
// starting point. Thus, 2 or more computers could each be running the program
// to simultaneously search different zones in the number field.
//
// The NcubedLow[] array contains the integer bit version of the last 9
// decimal digits of all the I^3's. These will be used to form a hash index.
// The hash index system greatly speeds up the process of finding matching
// pairs of I^3 + J^3 = K^3 + L^3, etc. (The actual hash index adds the bit
// portions of I^3 + J^3, and uses bits 8 to 29 of the result as the hash
// index. Note: The least significant bit is bit "0".)
//
// The NextGroup() routine will generate a large group of trial I^3 + J^3
// numbers. The "J's" that will be used in this routine will be >= "I" until
// the sum of I^3 + J^3 reaches the upper limit for the current group. (Group
// size is kept within bounds by stopping the process at "UpperLim".). This
// routine initializes the NextJ[] array for this process.
//
// "UpperLim" and "Increment" are initially set to 10 billion for the first
// iteration. (Valid for 4 & 5 pairs. Lower values for pairs = 2 & 3).
// 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 gradually cluster at a little over 2.4 million per
// iteration.
//
//**************************************************************************
void InitSys(void) {
int i;
double start, Difference, Idbl, Jcalc;
puts(" ramanujans.c");
puts(" by");
puts(" Bill Butler\n");
puts("This program finds Ramanujan 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 Ramanujan number.\n");
puts("Enter 2 for output of 2 or more pairs (Pauses on 3 pairs)");
puts("Enter 3 for output of 3 or more pairs (Pauses on 4 pairs)");
puts("Enter 4 for output of 4 or more pairs (Pauses on 5 pairs)");
puts("Enter 5 for output of 5 or more pairs (Pauses on 6 pairs)");
puts("(For output from 2, 3, 4 use \"Pause\" key as needed)");
gets(Databuff);
NbrPairs = atoi(Databuff);
start = ArrayLimit; // Convert to double for calc
printf("\nWhere do you want to start the search? (0 <= start <= %g)\n",
0.999999 * start * start * start);
puts("(Suggest a very small percent below true start point if > 0.)");
gets(Databuff);
start = atof(Databuff);
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");
puts("(Will take a few seconds.)\n");
for (i = 1; i <= ArrayLimit; i++) {
MakeCube(i); // Calculates i^3
NcubedHigh[i] = CubedHigh; // Leftmost 14 digits
NcubedLow[i] = CubedLow; // Rightmost 9 digits
// Initialize NextJ[]
Idbl = i; // Convert to double
Difference = start - Idbl*Idbl*Idbl; // Derived from nbr = I^3+J^3
if (Difference < 0.0) // If calc will be negative
NextJ[i] = i; // then ordinary start point
else { // Else calc which J to use
Jcalc = pow(Difference, 1.0/3.0) + 1.0; // Calculated value for J.
// Adding 1.0 avoids an array
// overflow problem
if (Jcalc > Idbl) // If this is bigger than "i"
NextJ[i] = Jcalc; // use the calculated value
else // Else
NextJ[i] = i; // use the normal value.
}
/* Optional debug check
printf("At i = %d NcubedHigh = %'.0lf * E9 NcubedLow = %'d NextJ = %d\n",
i, NcubedHigh[i], NcubedLow[i], NextJ[i]);
*/
}
UpperLim = InitValues[NbrPairs] + start; // Initial upper limit.
Increment = InitValues[NbrPairs]; // Initial increment.
puts("Starting search\n");
}
//*************************************************************************
//
// NextGroup
//
// This routine generates the next group of I^3 + J^3 sums and places these
// sums in the hash arrays. The number of these trial "I^3 + J^3" sums that
// are processed will average a little over 2.4 million.
//
//*************************************************************************
void NextGroup(void) {
int Count, i, j;
unsigned i3j3sumLow; // Does double duty as hash index
double LimNbr, i3j3sumHigh; // Forces "LimNbr" to a "Register"
// position. Otherwise "UpperLim"
// could be used.
for (i = 4194303; i >= 0; i--) { // Clear old garbage.
HashHd[i] = 0;
ChainLen[i] = 0;
}
Count = 0; // Count nbr. of I^3+J^3 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" example.
i = 1; // First "I" for I^3 (Start on row 1)
LimNbr = UpperLim; // Sets up reg. variable
do { // Do for all "i" in group (Will move
// down until I^3 + J^3 is too big)
j = NextJ[i]; // First "J" for J^3 (Leftmost column
// for current row in current group)
// Do for all "J" such that
while(1) { // I^3 + J^3 is < UpperLim
i3j3sumHigh = NcubedHigh[i] + NcubedHigh[j]; // Forms I^3 + J^3
i3j3sumLow = NcubedLow[i] + NcubedLow[j];
if (i3j3sumLow >= 1000000000) { // If carry overflow (>= 1 billion),
i3j3sumLow -= 1000000000; // decrease integer portion by 1e9
i3j3sumHigh += 1.0; // and increment the "billions"
}
// If within upper limit then
if ((i3j3sumHigh * 1.0E9 + i3j3sumLow) < LimNbr) {
Count++; // include it in the group.
I3J3High[Count] = i3j3sumHigh; // Add to the hash arrays.
I3J3Low[Count] = i3j3sumLow;
Ival[Count] = i;
Jval[Count] = j;
i3j3sumLow >>= 8; // Generate the
i3j3sumLow &= Bit22Mask; // hash index
Link[Count] = HashHd[i3j3sumLow]; // Update link list
HashHd[i3j3sumLow] = Count;
ChainLen[i3j3sumLow]++;
j++; // Move 1 col to right in
} // "Excel spreadsheet"
else // Repeat until too big
break;
}
NextJ[i] = j; // Set up for next group. Move down
i++; // 1 row in "Excel spreadsheet".
} while (j > i); // Loop exits when right edge of
// band/group intercepts the sloping
// diagonal in the spreadsheet. At
// this point, "i" will equal "j".
NbrStored = Count;
// Optional status check.
// Averages ~2,400,000+ per crack.
if (StatusFlag) {
printf("This iter. stored %'d trial I^3 + J^3 sums\n", NbrStored);
puts("Must stay < 3,000,000\n");
}
}
//************************************************************************
//
// 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 Ramanujan solution may exist. If true, the
// link list chain is checked to see if at least "NbrPairs" entries are
// identical. If true, then the information is output.
//
// The variable "NbrPairs" (see start of program) controls how many pairs
// have to exist for the Ramanujan number before this info is output.
//
//************************************************************************
void CheckGroup(void) {
int HashIndex, Ilink, Jlink;
int Count;
int StopNbr;
// For all hash heads
for (HashIndex = 0; HashIndex <= 4194303; HashIndex++) {
if ((ChainLen[HashIndex] - NbrPairs) < 0) // Chain is too short for output.
continue;
// Optional check to see if
// hash algorithm is efficient.
/*
printf("Link list[%d] has %d members\n", HashIndex, ChainLen[HashIndex]);
*/
// Check the list for this hash.
// StopNbr controls how many times
// that "Ilink" moves forward
StopNbr = ChainLen[HashIndex] - NbrPairs;
for (Ilink = HashHd[HashIndex]; StopNbr >= 0;
StopNbr--, Ilink = Link[Ilink]) {
Count = 1; // Will count number of identical
// I3J3 sums
for (Jlink = Link[Ilink]; Jlink; Jlink = Link[Jlink]) {
if (I3J3Low[Jlink] != I3J3Low[Ilink]) // Not a match if either is
continue; // different.
if (I3J3High[Jlink] != I3J3High[Ilink])
continue;
// If to here, the I3J3 sum at Jlink
// matches the I3J3 sum at Ilink
Count++;
}
if (Count < NbrPairs) // If not enough pair matches,
continue; // then keep looking.
printf("\n%d pairs exist at Ramanujan Number: %'.0lf * E+9 + %'d\n",
Count, I3J3High[Ilink], I3J3Low[Ilink]);
for (Jlink = Ilink; Jlink; Jlink = Link[Jlink]) {
if (I3J3Low[Jlink] != I3J3Low[Ilink]) // Both must match for a
continue; // paired solution
if (I3J3High[Jlink] != I3J3High[Ilink])
continue;
printf(" Pair at I = %'d J = %'d\n", Ival[Jlink], Jval[Jlink]);
}
if (Count > NbrPairs) // Pause if big match
pausemsg();
}
}
}
//*************************************************************************
//
// MakeCube
//
// This routine cubes the passed integer number and places the lowest 9
// decimal digits of the result in the global integer variable "CubedLow"
// and the high 14 digits into CubedHigh (double). Subsequently these
// results are loaded into the Ncubed arrays.
//
// The multiplication process is similar to ordinary base 10 arithmetic
// (and its associated "carry") except base 10,000,000 is used.
//
// The passed number "AnInteger" is known to be <= 3E7. This is converted
// to a "double" number which can be safely squared without any precision
// problems. This result is then split into three seperate sections that
// have the equivalent of 7 decimal digits each. Each of these can be safely
// multiplied by the original "AnInteger" again without any precision error.
//
// Finally, the end result is placed into two variables with the lowest 9
// decimal digits stored in one variable and the high 14 digits in the other.
// These two variables are then placed in the global variables "CubedHigh"
// and "CubedLow". (The latter becomes an integer variable.)
//
//**************************************************************************
void MakeCube(int AnInteger) {
double HighDigits, MidDigits, LowDigits, temp, Original;
// "AnInteger" will be <= 3E7. Thus
// there is enough precision in an
// ordinary "double" variable for
// the first multiplication.
Original = AnInteger; // Convert to "double"
temp = Original * Original; // Initially this will be <= 9E14
// Split it into high, mid, low digits
LowDigits = fmod(temp, 10000000.0); // These will be the low 7 digits
MidDigits = (temp - LowDigits) / 10000000.0; // This is the "carry"
// Now reduce MidDigits to 7 digits
temp = fmod(MidDigits, 10000000.0); // This will be the middle 7 digits
HighDigits = (MidDigits - temp) / 10000000.0; // Another "carry"
MidDigits = temp; // At this point "AnInteger" has been
// been squared. The lowest 7 digits are
// in "LowDigits", the next 7 digits are
// in "MidDigits", and the remaining
// digits are in "HighDigits".
LowDigits *= Original; // Mult. all three by "Original" to get
MidDigits *= Original; // cube. Then will have to sort out the
HighDigits *= Original; // pieces. After this multiplication, the
// true cubed result would be:
// HighDigits*E14+MidDigits*E7+LowDigits
temp = fmod(MidDigits, 100.0); // Get last 2 digits from MidDigits
MidDigits -= temp; // Remove last two digits
MidDigits /= 100.0; // Shift right two decimal digits
// so it measures "Billions"
LowDigits += 10000000.0 * temp; // All 9 low digits are in LowDigits
// (But still has some higher digits)
temp = fmod(LowDigits, 1000000000.0); // These are the final lowest 9 digits
MidDigits += (LowDigits - temp) / 1000000000.0; // Add carry to Mid Digits
MidDigits += 100000.0 * HighDigits; // The 14 high digits are in place
CubedHigh = MidDigits; // Store final result
CubedLow = temp; // Also converts it to an integer
}
//******************************************************************
//
// Misc. Routines
//
//******************************************************************
void pausemsg(void) {
puts("\nLarge number of pairs exist");
puts("Press RETURN to continue");
gets(Databuff);
}