Durango Bill's
"C" Program for Ramanujan Quadruples
Date: 30-MAY-01 Time: 18:58:03 rama4.c
Line Depth
Nbr. {} /* Text
() */
1 0 0 /******************************************************************/
2 0 0 /* */
3 0 0 /* Rama4.c */
4 0 0 /* by */
5 0 0 /* Bill Butler */
6 0 0 /* */
7 0 0 /* This program is dedicated to finding quadruple Ramanujan */
8 0 0 /* numbers. A quadruple Ramanujan is a number such that: */
9 0 0 /* I^3 + J^3 = K^3 + L^3 = M^3 + N^3 = O^3 + P^3. */
10 0 0 /* */
11 0 0 /* The first quadruple occurs at: */
12 0 0 /* */
13 0 0 /* 13,322^3 + 16,630^3 = */
14 0 0 /* 10,200^3 + 18,072^3 = */
15 0 0 /* 5,436^3 + 18,948^3 = */
16 0 0 /* 2,421^3 + 19,083^3 = 6,963,472,309,248 ( 6963472309248 ) */
17 0 0 /* */
18 0 0 /* This program may be used, copied, modified, etc. without any*/
19 0 0 /* obligation by any person for any purpose. I would appreciate */
20 0 0 /* that any published version (modified or original, or any */
21 0 0 /* published results) include a note crediting me with the */
22 0 0 /* program/algorithm. (e.g. "Original algorithm by Bill Butler.")*/
23 0 0 /* */
24 0 0 /* The algorithm (especially the hash index portion) is very */
25 0 0 /* efficient and appears to be at least 10 times faster than the */
26 0 0 /* "heap" algorithm used by David Wilson who published the first */
27 0 0 /* known Ramanujan quintuple. */
28 0 0 /* */
29 0 0 /* The program is in ANSI "C", and runs "as is" without bugs. */
30 0 0 /* (Bad news - this is under DOS on an old 80486 computer). It */
31 0 0 /* should start displaying Ramanujan quadruples within a few */
32 0 0 /* minutes on any contemporary Athlon or Intel processor. The */
33 0 0 /* program will eventually run into a precision problem before it*/
34 0 0 /* gets to the first Ramanujan quintuple. If you have 64 bit */
35 0 0 /* integer arithmetic available, you might try converting all */
36 0 0 /* "double" (real) numbers to 64 bit integers and then run the */
37 0 0 /* program. If you are successful, please let me know how it */
38 0 0 /* works out. (E-mail: lisabill@centurylink.net) */
39 0 0 /* (Using 64 bit integers, it should find the first Ramanujan */
40 0 0 /* quintuple in less than a week.) */
41 0 0 /* */
42 0 0 /* Note: All "int" variables are 32 bits. */
43 0 0 /* */
44 0 0 /******************************************************************/
45 0 0
46 0 0
47 0 0
48 0 0 #include <stdio.h> /* (May not need these) */
49 0 0 #include <stdlib.h> /* Need for atof() */
50 0 0 #include <ctype.h> /* tolower() */
51 0 0 #include <math.h>
52 0 0 #include <string.h>
53 0 0
54 0 0 #define HASHMAX 360000 /* Keep HASHMAX about 20 pct*/
55 0 0 #define AvgGrpSize 300000 /* larger than AvgGrpSize. */
56 0 0
57 0 0 #define MAXi 210000 /* Maximum trial "I" or "J" */
58 0 0 /* that will be used for */
59 0 0 /* I^3 + J^3. To extend the */
60 0 0 /* search, increase this */
61 0 0 /* number. Suggest 400000. */
62 0 0 /* (All other constants are */
63 0 0 /* OK for any size search.) */
64 0 0
65 0 0 double Ncubed[MAXi+1]; /* Ncubed[i] = i^3 */
66 0 0 unsigned Bits[MAXi+1]; /* The rightmost 32 bits */
67 0 0 /* of the integer versions */
68 0 0 /* of these cubes. */
69 0 0 int NextJ[MAXi+1]; /* Keeps track of the next */
70 0 0 /* "J" to try when forming */
71 0 0 /* trial I^3 + J^3. */
72 0 0
73 0 0 /* The size of the next 2 */
74 0 0 /* arrays matches the 19 bit*/
75 0 0 /* hash size. (2^19) */
76 0 0 unsigned HashHd[524288]; /* Heads for hash table. */
77 0 0 int ChainLen[524288]; /* Chain length. */
78 0 0
79 0 0
80 0 0 double I3J3[HASHMAX+1]; /* I^3 + J^3 */
81 0 0 /* "double" will eventually */
82 0 0 /* run into a precision */
83 0 0 /* problem. If you have 64 */
84 0 0 /* bit integers, use them */
85 0 0 /* instead of real numbers. */
86 0 0 int Ival[HASHMAX+1]; /* The "I" in I^3 + J^3 */
87 0 0 int Jval[HASHMAX+1]; /* The "J" in I^3 + J^3 */
88 0 0 int Link[HASHMAX+1]; /* Link lists. */
89 0 0
90 0 0 char Databuff[100]; /* Dummy array for input. */
91 0 0
92 0 0 double UpperLim; /* For each iter., find cube */
93 0 0 /* sums up to this limit. */
94 0 0 double Increment; /* Increases the upper */
95 0 0 /* limit by this amount on */
96 0 0 /* each iter. Note: */
97 0 0 /* "Increment" increases */
98 0 0 /* with time. */
99 0 0 int NbrStored; /* Nbr. of items in hash */
100 0 0 /* table. */
101 0 0 unsigned Bit19Mask = 524287; /* Mask for right 19 bits */
102 0 0
103 0 0
104 1 0 main() {
105 1 0
106 1 0 int i;
107 1 0
108 1 0 InitSys(); /* Initialize system. */
109 1 0
110 2 0 while(1) { /* Do forever. */
111 2 0 NextGroup(); /* Form next group of */
112 2 0 /* cubes. */
113 2 0 CheckGroup(); /* Look for matches. */
114 3 0 if (NbrStored < AvgGrpSize) { /* Process about 300,000 */
115 3 0 /* trial I^3 + J^3 ea. iter.*/
116 3 0 Increment *= 1.1; /* Increase by 10 % as needed*/
117 2 0 }
118 2 0 UpperLim += Increment; /* Set up for next group */
119 1 0 }
120 0 0 }
121 0 0
122 0 0
123 0 0
124 0 0 /******************************************************************/
125 0 0 /* */
126 0 0 /* InitSys */
127 0 0 /* */
128 0 0 /* This routine initializes the system. The Ncubed[] array is */
129 0 0 /* filled with cubes such that Ncubed[i] = I^3. */
130 0 0 /* */
131 0 0 /* The Bits[] array contains the last 32 bits of the integer */
132 0 0 /* representation of all the I^3's. These will be used to form a */
133 0 0 /* hash index. The hash index system greatly speeds up the */
134 0 0 /* process of finding matching pairs of I^3 + J^3 = K^3 + L^3, */
135 0 0 /* etc. (The actual hash index adds the 32 bit portions of */
136 0 0 /* I^3 + J^3, and uses bits 8 to 26 of the result as the hash */
137 0 0 /* index. Note: the least significant bit is bit "0".) */
138 0 0 /* */
139 0 0 /* The NextGroup() routine will generate a large group of */
140 0 0 /* trial I^3 + J^3 numbers. The "J's" that will be used in this */
141 0 0 /* routine will be >= "I" until the sum of I^3 + J^3 reaches the */
142 0 0 /* upper limit for the current group. (Group size is kept within */
143 0 0 /* bounds by stopping the process at "UpperLim".). This routine */
144 0 0 /* initializes the NextJ[] array for this process. */
145 0 0 /* */
146 0 0 /* "UpperLim" is set to 200000000 for the first iteration. */
147 0 0 /* Subsequently it will be increased by "Increment" to include */
148 0 0 /* larger trial values of I^3 + J^3. The density of I^3 + J^3 */
149 0 0 /* numbers becomes sparser as the number field expands. */
150 0 0 /* "Increment" will thus become larger with time. Thus, the */
151 0 0 /* number of trial I^3 + J^3 numbers that will be looked at will */
152 0 0 /* gradually cluster around "AvgGrpSize" (300,000) per iteration.*/
153 0 0 /* */
154 0 0 /* Finally, when a group of trial I^3 + J^3 numbers cluster at */
155 0 0 /* a given hash index, all numbers in this link list group are */
156 0 0 /* copied to the end of the I3J3[] array for final sorting. The */
157 0 0 /* last position in the array is initialized to -1 to aid the */
158 0 0 /* "insertion sort". */
159 0 0 /* */
160 0 0 /******************************************************************/
161 0 0
162 1 0 InitSys() {
163 1 0
164 1 0
165 1 0 unsigned i, Ibits;
166 1 0 double Ifloat, Icubed;
167 1 0
168 1 0 puts("\nThis program finds Ramanujan quadruples. It will run");
169 1 0 puts("until the user manually breaks the program. However it");
170 1 0 puts("will pause anytime a Ramanujan quintuple is found.");
171 1 0 pausemsg();
172 1 0 puts("Initializing the cubes table");
173 2 0 for (i = 1; i <= MAXi; i++) {
174 2 0 Ifloat = i; /* Convert int to real. */
175 2 0 Icubed = Ifloat * Ifloat * Ifloat;
176 2 0 Ncubed[i] = Icubed;
177 2 0 Ibits = i * i * i; /* Calculate the rightmost */
178 2 0 Bits[i] = Ibits; /* 32 bits. (Note "C" */
179 2 0 /* ignores the overflow.) */
180 2 0 NextJ[i] = i; /* Init the "J's" */
181 1 0 }
182 1 0 puts("Starting search");
183 1 0
184 1 0 UpperLim = 200000000.0; /* Initial upper limit. */
185 1 0 Increment = 200000000.0; /* Initial increment. */
186 1 0 I3J3[HASHMAX} = -1.0; /* Used for sort. */
187 1 0 /* Requires a negative */
188 1 0 /* number. */
189 0 0 }
190 0 0
191 0 0
192 0 0
193 0 0 /******************************************************************/
194 0 0 /* */
195 0 0 /* NextGroup */
196 0 0 /* */
197 0 0 /* This routine generates the next group of I^3 + J^3 and */
198 0 0 /* places it in the hash arrays. The number of these trial */
199 0 0 /* I^3 + J^3 sums that are processed will oscillate near */
200 0 0 /* "AvgGrpSize" (about 300,000). */
201 0 0 /* */
202 0 0 /******************************************************************/
203 0 0
204 1 0 NextGroup() {
205 1 0
206 1 0
207 1 0 int Count, i, j;
208 1 0 unsigned IJhash;
209 1 0 double LimNbr, i3j3sum; /* Forces "LimNbr" to a */
210 1 0 /* "Register" position. */
211 1 0 /* Otherwise "UpperLim" */
212 1 0 /* could be used. */
213 1 0
214 2 0 for (i = 524287; i >= 0; i--) { /* Clear old garbage. */
215 2 0 HashHd[i] = 0;
216 2 0 ChainLen[i] = 0;
217 1 0 }
218 1 0 Count = 0;
219 1 0 i = 1;
220 1 0 LimNbr = UpperLim; /* Sets up reg. number */
221 2 0 do { /* Do for all "i" in group */
222 2 0 j = NextJ[i];
223 2 0 /* Do for all "j" such that */
224 3 0 while(1) { /* i^3 + j^3 is < UpperLim */
225 3 0 i3j3sum = Ncubed[i] + Ncubed[j]; /* I^3 + J^3 */
226 4 0 if (i3j3sum < LimNbr) { /* If within limit, then */
227 4 0 Count++; /* include it in the group. */
228 4 0 IJhash = Bits[i] + Bits[j]; /* Generate the */
229 4 0 IJhash = IJhash >> 8; /* hash index. */
230 4 0 IJhash &= Bit19Mask;
231 4 0 I3J3[Count] = i3j3sum; /* Add to the hash arrays. */
232 4 0 Ival[Count] = i;
233 4 0 Jval[Count] = j;
234 4 0 Link[Count] = HashHd[IJhash]; /* Update link list */
235 4 0 HashHd[IJhash] = Count;
236 4 0 ChainLen[IJhash]++;
237 4 0 j++;
238 3 0 }
239 3 0 else /* Repeat until too big */
240 3 0 break;
241 2 0 }
242 2 0 NextJ[i] = j; /* Set up for next round */
243 2 0 i++;
244 1 0 } while (j > i);
245 1 0
246 1 0 NbrStored = Count;
247 1 0 /* Optional status check. */
248 1 0 /* Will average about */
249 1 0 /* 300,000 per crack. */
250 1 1 /*
251 1 1 printf("This iter. stored %d trial I^3+J^3 numbers\n", NbrStored);
252 1 0 */
253 0 0 }
254 0 0
255 0 0
256 0 0
257 0 0 /******************************************************************/
258 0 0 /* */
259 0 0 /* CheckGroup */
260 0 0 /* */
261 0 0 /* This routine checks the hash arrays to see if any matches */
262 0 0 /* exist. If at least 4 numbers exist at any hash index, the */
263 0 0 /* entire link list is copied to the sort arrays where it is */
264 0 0 /* sorted. Note: The end of all arrays dimensioned by "HASHMAX" */
265 0 0 /* are used for sorting. (Usually the list is very short.) */
266 0 0 /* */
267 0 0 /******************************************************************/
268 0 0
269 1 0 CheckGroup() {
270 1 0
271 1 0 int i, j, k, EndLoc, Next;
272 1 0 double TempDbl;
273 1 0
274 2 0 for (i = 524287; i >= 0; i--) {
275 2 0 if (ChainLen[i] < 4) /* Chain is too */
276 2 0 continue; /* short for quads. */
277 2 0 EndLoc = HASHMAX;
278 3 0 for (Next = HashHd[i]; Next; Next = Link[Next]) {
279 3 0 TempDbl = I3J3[Next];
280 3 0 EndLoc--;
281 3 0 /* Use "Insertion sort" */
282 4 0 for (j = EndLoc, k = j + 1; I3J3[k] > TempDbl; j++, k++) {
283 4 0 I3J3[j] = I3J3[k];
284 4 0 Ival[j] = Ival[k];
285 4 0 Jval[j] = Jval[k];
286 3 0 }
287 3 0 I3J3[j] = TempDbl;
288 3 0 Ival[j] = Ival[Next];
289 3 0 Jval[j] = Jval[Next];
290 2 0 }
291 2 0 /* Check for quadruples. */
292 3 0 for (j = HASHMAX-4, k = HASHMAX-1; j >= EndLoc; j--, k--) {
293 3 0 if (I3J3[j] == I3J3[k])
294 4 0 printf("%7d%7d %7d%7d %7d%7d %7d%7d %.0lf\n",
295 4 0 Ival[j], Jval[j], Ival[j+1], Jval[j+1],
296 4 0 Ival[j+2], Jval[j+2], Ival[k], Jval[k],
297 3 0 I3J3[j]);
298 3 0 /* The following can be deleted */
299 3 0 /* unless the program is modified */
300 3 0 /* for 64 bit integers. */
301 4 0 if (I3J3[j] == I3J3[k+1]) {
302 4 0 puts("Ramanujan Quintuple. Press ENTER to continue.");
303 4 0 pausemsg();
304 3 0 }
305 2 0 }
306 1 0 }
307 0 0 }
308 0 0
309 0 0
310 0 0 /******************************************************************/
311 0 0 /* */
312 0 0 /* Misc. Routines */
313 0 0 /* */
314 0 0 /******************************************************************/
315 0 0
316 1 0 pausemsg() {
317 1 0
318 1 0 puts("\nPress RETURN to continue");
319 1 0 gets(Databuff);
320 0 0 }
Return to Ramanujan Number Page
Web page generated via Sea Monkey's Composer HTML editor
within a Linux Cinnamon Mint 18 operating system.
(Goodbye Microsoft)