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-2024 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 /* Work around <https://gcc.gnu.org/bugzilla/show_bug.cgi?id=109635>.  */
 114 #if 13 <= __GNUC__
 115 # pragma GCC diagnostic ignored "-Wanalyzer-use-of-uninitialized-value"
 116 #endif
 117
 118 static void
 119 output_primes (const struct prime *primes, unsigned nprimes)
 120 {
 121   unsigned i;
 122   unsigned p;
 123   int is_prime;
 124
 125   /* Compute wide_uint_bits by repeated shifting, rather than by
 126      multiplying sizeof by CHAR_BIT, as this works even if the
 127      wide_uint representation has holes.  */
 128   unsigned wide_uint_bits = 0;
 129   wide_uint mask = -1;
 130   for (wide_uint_bits = 0; mask; wide_uint_bits++)
 131     mask >>= 1;
 132
 133   puts ("/* Generated file -- DO NOT EDIT */\n");
 134   printf ("#define WIDE_UINT_BITS %u\n", wide_uint_bits);
 135
 136   for (i = 0, p = 2; i < nprimes; i++)
 137     {
 138       unsigned int d8 = i + 8 < nprimes ? primes[i + 8].p - primes[i].p : 0xff;
 139       if (255 < d8) /* this happens at 668221 */
 140         abort ();
 141       printf ("P (%u, %u,\n   (", primes[i].p - p, d8);
 142       print_wide_uint (primes[i].pinv, 0, wide_uint_bits);
 143       printf ("),\n   UINTMAX_MAX / %u)\n", primes[i].p);
 144       p = primes[i].p;
 145     }
 146
 147   printf ("\n#undef FIRST_OMITTED_PRIME\n");
 148
 149   /* Find next prime */
 150   do
 151     {
 152       p += 2;
 153       for (i = 0, is_prime = 1; is_prime; i++)
 154         {
 155           if (primes[i].p * primes[i].p > p)
 156             break;
 157           if (p * primes[i].pinv <= primes[i].lim)
 158             {
 159               is_prime = 0;
 160               break;
 161             }
 162         }
 163     }
 164   while (!is_prime);
 165
 166   printf ("#define FIRST_OMITTED_PRIME %u\n", p);
 167 }
 168
 169 ATTRIBUTE_MALLOC
 170 static void *
 171 xalloc (size_t s)
 172 {
 173   void *p = malloc (s);
 174   if (p)
 175     return p;
 176
 177   fprintf (stderr, "Virtual memory exhausted.\n");
 178   exit (EXIT_FAILURE);
 179 }
 180
 181 int
 182 main (int argc, char **argv)
 183 {
 184   int limit;
 185
 186   char *sieve;
 187   size_t size, i;
 188
 189   struct prime *prime_list;
 190   unsigned nprimes;
 191
 192   if (argc != 2)
 193     {
 194       fprintf (stderr, "Usage: %s LIMIT\n"
 195                "Produces a list of odd primes <= LIMIT\n", argv[0]);
 196       return EXIT_FAILURE;
 197     }
 198   limit = atoi (argv[1]);
 199   if (limit < 3)
 200     return EXIT_SUCCESS;
 201
 202   /* Make limit odd */
 203   if ( !(limit & 1))
 204     limit--;
 205
 206   size = (limit - 1) / 2;
 207   /* sieve[i] represents 3+2*i */
 208   sieve = xalloc (size);
 209   memset (sieve, 1, size);
 210
 211   prime_list = xalloc (size * sizeof (*prime_list));
 212   nprimes = 0;
 213
 214   for (i = 0; i < size;)
 215     {
 216       unsigned p = 3 + 2 * i;
 217       unsigned j;
 218
 219       process_prime (&prime_list[nprimes++], p);
 220
 221       for (j = (p * p - 3) / 2; j < size; j += p)
 222         sieve[j] = 0;
 223
 224       while (++i < size && sieve[i] == 0)
 225         ;
 226     }
 227
 228   output_primes (prime_list, nprimes);
 229
 230   free (sieve);
 231   free (prime_list);
 232
 233   if (ferror (stdout) + fclose (stdout))
 234     {
 235       fprintf (stderr, "write error: %s\n", strerror (errno));
 236       return EXIT_FAILURE;
 237     }
 238
 239   return EXIT_SUCCESS;
 240 }