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@mywdurango.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