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:
lisawbill@mydurango.net)
*/
/* (Remove the "w" for a valid
E-mail)
*/
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
KompoZer