src/make-prime-list.c

   1 /* Factoring of uintmax_t numbers. Generation of needed tables.
   2
   3    Contributed to the GNU project by Torbjörn Granlund and Niels Möller
   4    Contains code from GNU MP.
   5
   6 Copyright 2012 Free Software Foundation, Inc.
   7
   8 This program is free software; you can redistribute it and/or modify it under
   9 the terms of the GNU General Public License as published by the Free Software
  10 Foundation; either version 3 of the License, or (at your option) any later
  11 version.
  12
  13 This program is distributed in the hope that it will be useful, but WITHOUT ANY
  14 WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
  15 PARTICULAR PURPOSE.  See the GNU General Public License for more details.
  16
  17 You should have received a copy of the GNU General Public License along with
  18 this program.  If not, see http://www.gnu.org/licenses/.  */
  19
  20 #include <config.h>
  21
  22 #include <stdint.h>
  23 #include <stdio.h>
  24 #include <string.h>
  25 #include <stdlib.h>
  26 #include <errno.h>
  27
  28 /* Deactivate config.h's "rpl_"-prefixed definitions of these symbols.  */
  29 #undef fclose
  30 #undef strerror
  31
  32 struct prime
  33 {
  34   unsigned p;
  35   uintmax_t pinv; /* Inverse mod b = 2^{bitsize of uintmax_t} */
  36   uintmax_t lim; /* floor(UINTMAX_MAX / p) */
  37 };
  38
  39 static uintmax_t _GL_ATTRIBUTE_CONST
  40 binvert (uintmax_t a)
  41 {
  42   uintmax_t x = 0xf5397db1 >> (4*((a/2) & 0x7));
  43   for (;;)
  44     {
  45       uintmax_t y = 2*x - x*x*a;
  46       if (y == x)
  47         return x;
  48       x = y;
  49     }
  50 }
  51
  52 static void
  53 process_prime (struct prime *info, unsigned p)
  54 {
  55   info->p = p;
  56   info->pinv = binvert (p);
  57   info->lim = UINTMAX_MAX / (uintmax_t) p;
  58 }
  59
  60 static void
  61 output_primes (const struct prime *primes, unsigned nprimes)
  62 {
  63   unsigned i;
  64   unsigned p;
  65   int is_prime;
  66   const char *suffix;
  67
  68   puts ("/* Generated file -- DO NOT EDIT */\n");
  69
  70   if (sizeof (uintmax_t) <= sizeof (unsigned long))
  71     suffix = "UL";
  72   else if (sizeof (uintmax_t) <= sizeof (unsigned long long))
  73     suffix = "ULL";
  74   else
  75     {
  76       fprintf (stderr, "Don't know how to write uintmax_t constants.\n");
  77       exit (EXIT_FAILURE);
  78     }
  79
  80 #define SZ (int)(2*sizeof (uintmax_t))
  81
  82   for (i = 0, p = 2; i < nprimes; i++)
  83     {
  84       unsigned int d8 = i + 8 < nprimes ? primes[i + 8].p - primes[i].p : 0xff;
  85       if (255 < d8) /* this happens at 668221 */
  86         abort ();
  87       printf ("P (%2u, %3u, 0x%0*jx%s, 0x%0*jx%s) /* %d */\n",
  88               primes[i].p - p, d8,
  89               SZ, primes[i].pinv, suffix,
  90               SZ, primes[i].lim, suffix, primes[i].p);
  91       p = primes[i].p;
  92     }
  93
  94   printf ("\n#undef FIRST_OMITTED_PRIME\n");
  95
  96   /* Find next prime */
  97   do
  98     {
  99       p += 2;
 100       for (i = 0, is_prime = 1; is_prime; i++)
 101         {
 102           if (primes[i].p * primes[i].p > p)
 103             break;
 104           if ( (uintmax_t) p * primes[i].pinv <= primes[i].lim)
 105             {
 106               is_prime = 0;
 107               break;
 108             }
 109         }
 110     }
 111   while (!is_prime);
 112
 113   printf ("#define FIRST_OMITTED_PRIME %d\n", p);
 114 }
 115
 116 static void *
 117 xalloc (size_t s)
 118 {
 119   void *p = malloc (s);
 120   if (p)
 121     return p;
 122
 123   fprintf (stderr, "Virtual memory exhausted.\n");
 124   exit (EXIT_FAILURE);
 125 }
 126
 127 int
 128 main (int argc, char **argv)
 129 {
 130   int limit;
 131
 132   char *sieve;
 133   size_t size, i;
 134
 135   struct prime *prime_list;
 136   unsigned nprimes;
 137
 138   if (argc != 2)
 139     {
 140       fprintf (stderr, "Usage: %s LIMIT\n"
 141                "Produces a list of odd primes <= LIMIT\n", argv[0]);
 142       return EXIT_FAILURE;
 143     }
 144   limit = atoi (argv[1]);
 145   if (limit < 3)
 146     exit (EXIT_SUCCESS);
 147
 148   /* Make limit odd */
 149   if ( !(limit & 1))
 150     limit--;
 151
 152   size = (limit-1)/2;
 153   /* sieve[i] represents 3+2*i */
 154   sieve = xalloc (size);
 155   memset (sieve, 1, size);
 156
 157   prime_list = xalloc (size * sizeof (*prime_list));
 158   nprimes = 0;
 159
 160   for (i = 0; i < size;)
 161     {
 162       unsigned p = 3+2*i;
 163       unsigned j;
 164
 165       process_prime (&prime_list[nprimes++], p);
 166
 167       for (j = (p*p - 3)/2; j < size; j+= p)
 168         sieve[j] = 0;
 169
 170       while (i < size && sieve[++i] == 0)
 171         ;
 172     }
 173
 174   output_primes (prime_list, nprimes);
 175
 176   if (ferror (stdout) + fclose (stdout))
 177     {
 178       fprintf (stderr, "write error: %s\n", strerror (errno));
 179       return EXIT_FAILURE;
 180     }
 181
 182   return EXIT_SUCCESS;
 183 }