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