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-2013 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 <limits.h>
  23 #include <stdint.h>
  24 #include <inttypes.h>
  25 #include <stdio.h>
  26 #include <string.h>
  27 #include <stdlib.h>
  28 #include <errno.h>
  29
  30 /* Deactivate config.h's "rpl_"-prefixed definitions of these symbols.  */
  31 #undef fclose
  32 #undef malloc
  33 #undef strerror
  34
  35 /* An unsigned type that is no narrower than 32 bits and no narrower
  36    than unsigned int.  It's best to make it as wide as possible.
  37    For GCC 4.6 and later, use a heuristic to guess whether unsigned
  38    __int128 works on your platform.  If this heuristic does not work
  39    for you, please report a bug; in the meantime compile with, e.g.,
  40    -Dwide_uint='unsigned __int128' to override the heuristic.  */
  41 #ifndef wide_uint
  42 # if 4 < __GNUC__ + (6 <= __GNUC_MINOR__) && ULONG_MAX >> 31 >> 31 >> 1 != 0
  43 typedef unsigned __int128 wide_uint;
  44 # else
  45 typedef uintmax_t wide_uint;
  46 # endif
  47 #endif
  48
  49 struct prime
  50 {
  51   unsigned p;
  52   wide_uint pinv; /* Inverse mod b = 2^{bitsize of wide_uint} */
  53   wide_uint lim; /* floor ((wide_uint) -1 / p) */
  54 };
  55
  56 static wide_uint _GL_ATTRIBUTE_CONST
  57 binvert (wide_uint a)
  58 {
  59   wide_uint x = 0xf5397db1 >> (4*((a/2) & 0x7));
  60   for (;;)
  61     {
  62       wide_uint y = 2*x - x*x*a;
  63       if (y == x)
  64         return x;
  65       x = y;
  66     }
  67 }
  68
  69 static void
  70 process_prime (struct prime *info, unsigned p)
  71 {
  72   wide_uint max = -1;
  73   info->p = p;
  74   info->pinv = binvert (p);
  75   info->lim = max / p;
  76 }
  77
  78 static void
  79 print_wide_uint (wide_uint n, int nesting, unsigned wide_uint_bits)
  80 {
  81   /* Number of bits per integer literal.  8 is too many, because
  82      uintmax_t is 32 bits on some machines so we cannot shift by 32 bits.
  83      So use 7.  */
  84   int hex_digits_per_literal = 7;
  85   int bits_per_literal = hex_digits_per_literal * 4;
  86
  87   unsigned remainder = n & ((1 << bits_per_literal) - 1);
  88
  89   if (n != remainder)
  90     {
  91       int needs_parentheses = n >> bits_per_literal >> bits_per_literal != 0;
  92       if (needs_parentheses)
  93         printf ("(");
  94       print_wide_uint (n >> bits_per_literal, nesting + 1, wide_uint_bits);
  95       if (needs_parentheses)
  96         printf (")\n%*s", nesting + 3, "");
  97       printf (" << %d | ", bits_per_literal);
  98     }
  99   else if (nesting)
 100     {
 101       printf ("(uintmax_t) ");
 102       hex_digits_per_literal
 103         = ((wide_uint_bits - 1) % bits_per_literal) % 4 + 1;
 104     }
 105
 106   printf ("0x%0*xU", hex_digits_per_literal, remainder);
 107 }
 108
 109 static void
 110 output_primes (const struct prime *primes, unsigned nprimes)
 111 {
 112   unsigned i;
 113   unsigned p;
 114   int is_prime;
 115
 116   /* Compute wide_uint_bits by repeated shifting, rather than by
 117      multiplying sizeof by CHAR_BIT, as this works even if the
 118      wide_uint representation has holes.  */
 119   unsigned wide_uint_bits = 0;
 120   wide_uint mask = -1;
 121   for (wide_uint_bits = 0; mask; wide_uint_bits++)
 122     mask >>= 1;
 123
 124   puts ("/* Generated file -- DO NOT EDIT */\n");
 125   printf ("#define WIDE_UINT_BITS %u\n", wide_uint_bits);
 126
 127   for (i = 0, p = 2; i < nprimes; i++)
 128     {
 129       unsigned int d8 = i + 8 < nprimes ? primes[i + 8].p - primes[i].p : 0xff;
 130       if (255 < d8) /* this happens at 668221 */
 131         abort ();
 132       printf ("P (%u, %u,\n   (", primes[i].p - p, d8);
 133       print_wide_uint (primes[i].pinv, 0, wide_uint_bits);
 134       printf ("),\n   UINTMAX_MAX / %d)\n", primes[i].p);
 135       p = primes[i].p;
 136     }
 137
 138   printf ("\n#undef FIRST_OMITTED_PRIME\n");
 139
 140   /* Find next prime */
 141   do
 142     {
 143       p += 2;
 144       for (i = 0, is_prime = 1; is_prime; i++)
 145         {
 146           if (primes[i].p * primes[i].p > p)
 147             break;
 148           if (p * primes[i].pinv <= primes[i].lim)
 149             {
 150               is_prime = 0;
 151               break;
 152             }
 153         }
 154     }
 155   while (!is_prime);
 156
 157   printf ("#define FIRST_OMITTED_PRIME %u\n", p);
 158 }
 159
 160 static void *
 161 xalloc (size_t s)
 162 {
 163   void *p = malloc (s);
 164   if (p)
 165     return p;
 166
 167   fprintf (stderr, "Virtual memory exhausted.\n");
 168   exit (EXIT_FAILURE);
 169 }
 170
 171 int
 172 main (int argc, char **argv)
 173 {
 174   int limit;
 175
 176   char *sieve;
 177   size_t size, i;
 178
 179   struct prime *prime_list;
 180   unsigned nprimes;
 181
 182   if (argc != 2)
 183     {
 184       fprintf (stderr, "Usage: %s LIMIT\n"
 185                "Produces a list of odd primes <= LIMIT\n", argv[0]);
 186       return EXIT_FAILURE;
 187     }
 188   limit = atoi (argv[1]);
 189   if (limit < 3)
 190     exit (EXIT_SUCCESS);
 191
 192   /* Make limit odd */
 193   if ( !(limit & 1))
 194     limit--;
 195
 196   size = (limit-1)/2;
 197   /* sieve[i] represents 3+2*i */
 198   sieve = xalloc (size);
 199   memset (sieve, 1, size);
 200
 201   prime_list = xalloc (size * sizeof (*prime_list));
 202   nprimes = 0;
 203
 204   for (i = 0; i < size;)
 205     {
 206       unsigned p = 3+2*i;
 207       unsigned j;
 208
 209       process_prime (&prime_list[nprimes++], p);
 210
 211       for (j = (p*p - 3)/2; j < size; j+= p)
 212         sieve[j] = 0;
 213
 214       while (i < size && sieve[++i] == 0)
 215         ;
 216     }
 217
 218   output_primes (prime_list, nprimes);
 219
 220   if (ferror (stdout) + fclose (stdout))
 221     {
 222       fprintf (stderr, "write error: %s\n", strerror (errno));
 223       return EXIT_FAILURE;
 224     }
 225
 226   return EXIT_SUCCESS;
 227 }