polly/lib/External/isl/imath/examples/randprime.c

   1 /*
   2   Name:     randprime.c
   3   Purpose:  Generate a probable prime at random.
   4   Author:   M. J. Fromberger
   5
   6   Usage:  randprime [-s] <bits> [<outfile>]
   7
   8   Generate a randomly-chosen probable prime having <bits> significant bits, and
   9   write it to the specified output file or to the standard output.  If the "-s"
  10   option is given, a prime p is chosen such that (p - 1) / 2 is also prime.
  11
  12   A prime is obtained by reading random bits from /dev/random, setting the
  13   low-order bit, and testing for primality.  If the first candidate is not
  14   prime, successive odd candidates are tried until a probable prime is found.
  15
  16   Copyright (C) 2002-2008 Michael J. Fromberger, All Rights Reserved.
  17
  18   Permission is hereby granted, free of charge, to any person obtaining a copy
  19   of this software and associated documentation files (the "Software"), to deal
  20   in the Software without restriction, including without limitation the rights
  21   to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
  22   copies of the Software, and to permit persons to whom the Software is
  23   furnished to do so, subject to the following conditions:
  24
  25   The above copyright notice and this permission notice shall be included in
  26   all copies or substantial portions of the Software.
  27
  28   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  29   IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  30   FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN NO EVENT SHALL THE
  31   AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  32   LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
  33   OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
  34   SOFTWARE.
  35  */
  36
  37 #include <errno.h>
  38 #include <limits.h>
  39 #include <stdio.h>
  40 #include <stdlib.h>
  41 #include <string.h>
  42
  43 #include <getopt.h>
  44 #include <unistd.h>
  45
  46 #include "imath.h"
  47 #include "iprime.h"
  48
  49 /* Load the specified buffer with random bytes */
  50 int randomize(unsigned char *buf, size_t len);
  51
  52 /* Overwrite the specified value with n_bits random bits */
  53 mp_result mp_int_randomize(mp_int a, mp_size n_bits);
  54
  55 /* Find a prime starting from the given odd seed */
  56 mp_result find_prime(mp_int seed, FILE *fb);
  57 mp_result find_strong_prime(mp_int seed, FILE *fb);
  58
  59 typedef mp_result (*find_f)(mp_int, FILE *);
  60
  61 int main(int argc, char *argv[]) {
  62   int opt, modbits;
  63   FILE *ofp = stdout;
  64   mp_result res;
  65   find_f find_func = find_prime;
  66   char tag = 'p';
  67   mpz_t value;
  68
  69   /* Process command-line arguments */
  70   while ((opt = getopt(argc, argv, "s")) != EOF) {
  71     switch (opt) {
  72       case 's':
  73         find_func = find_strong_prime;
  74         tag = 'P';
  75         break;
  76       default:
  77         fprintf(stderr, "Usage: randprime [-s] <bits> [<outfile>]\n");
  78         return 1;
  79     }
  80   }
  81
  82   if (optind >= argc) {
  83     fprintf(stderr,
  84             "Error:  You must specify the number of significant bits.\n");
  85     fprintf(stderr, "Usage: randprime [-s] <bits> [<outfile>]\n");
  86     return 1;
  87   }
  88   modbits = (int)strtol(argv[optind++], NULL, 0);
  89   if (modbits < CHAR_BIT) {
  90     fprintf(stderr, "Error:  Invalid value for number of significant bits.\n");
  91     return 1;
  92   }
  93   if (modbits % 2 == 1) ++modbits;
  94
  95   /* Check if output file is specified */
  96   if (optind < argc) {
  97     if ((ofp = fopen(argv[optind], "wt")) == NULL) {
  98       fprintf(stderr,
  99               "Error:  Unable to open output file for writing.\n"
 100               " - Filename: %s\n"
 101               " - Error:    %s\n",
 102               argv[optind], strerror(errno));
 103       return 1;
 104     }
 105   }
 106
 107   mp_int_init(&value);
 108   if ((res = mp_int_randomize(&value, modbits - 1)) != MP_OK) {
 109     fprintf(stderr,
 110             "Error:  Unable to generate random start value.\n"
 111             " - %s (%d)\n",
 112             mp_error_string(res), res);
 113     goto EXIT;
 114   }
 115   fprintf(stderr, "%c: ", tag);
 116   find_func(&value, stderr);
 117   fputc('\n', stderr);
 118
 119   /* Write the completed value to the specified output file */
 120   {
 121     int len;
 122     char *obuf;
 123
 124     len = mp_int_string_len(&value, 10);
 125     obuf = malloc(len);
 126     mp_int_to_string(&value, 10, obuf, len);
 127     fputs(obuf, ofp);
 128     fputc('\n', ofp);
 129
 130     free(obuf);
 131   }
 132
 133 EXIT:
 134   fclose(ofp);
 135   mp_int_clear(&value);
 136   return 0;
 137 }
 138
 139 int randomize(unsigned char *buf, size_t len) {
 140   FILE *rnd = fopen("/dev/random", "rb");
 141   size_t nr;
 142
 143   if (rnd == NULL) return -1;
 144
 145   nr = fread(buf, sizeof(*buf), len, rnd);
 146   fclose(rnd);
 147
 148   return (int)nr;
 149 }
 150
 151 mp_result mp_int_randomize(mp_int a, mp_size n_bits) {
 152   mp_size n_bytes = (n_bits + CHAR_BIT - 1) / CHAR_BIT;
 153   unsigned char *buf;
 154   mp_result res = MP_OK;
 155
 156   if ((buf = malloc(n_bytes)) == NULL) return MP_MEMORY;
 157
 158   if ((mp_size)randomize(buf, n_bytes) != n_bytes) {
 159     res = MP_TRUNC;
 160     goto CLEANUP;
 161   }
 162
 163   /* Clear bits beyond the number requested */
 164   if (n_bits % CHAR_BIT != 0) {
 165     unsigned char b_mask = (1 << (n_bits % CHAR_BIT)) - 1;
 166     unsigned char t_mask = (1 << (n_bits % CHAR_BIT)) >> 1;
 167
 168     buf[0] &= b_mask;
 169     buf[0] |= t_mask;
 170   }
 171
 172   /* Set low-order bit to insure value is odd */
 173   buf[n_bytes - 1] |= 1;
 174
 175   res = mp_int_read_unsigned(a, buf, n_bytes);
 176
 177 CLEANUP:
 178   memset(buf, 0, n_bytes);
 179   free(buf);
 180
 181   return res;
 182 }
 183
 184 mp_result find_prime(mp_int seed, FILE *fb) {
 185   mp_result res;
 186   int count = 0;
 187
 188   if (mp_int_is_even(seed)) {
 189     if ((res = mp_int_add_value(seed, 1, seed)) != MP_OK) {
 190       return res;
 191     }
 192   }
 193
 194   while ((res = mp_int_is_prime(seed)) == MP_FALSE) {
 195     ++count;
 196
 197     if (fb != NULL && (count % 50) == 0) {
 198       fputc('.', fb);
 199     }
 200     if ((res = mp_int_add_value(seed, 2, seed)) != MP_OK) {
 201       return res;
 202     }
 203   }
 204
 205   if (res == MP_TRUE && fb != NULL) fputc('+', fb);
 206
 207   return res;
 208 }
 209
 210 mp_result find_strong_prime(mp_int seed, FILE *fb) {
 211   mp_result res = MP_OK;
 212   mpz_t t;
 213
 214   mp_int_init(&t);
 215   for (;;) {
 216     if (find_prime(seed, fb) != MP_TRUE) break;
 217     if (mp_int_copy(seed, &t) != MP_OK) break;
 218
 219     if (mp_int_mul_pow2(&t, 1, &t) != MP_OK ||
 220         mp_int_add_value(&t, 1, &t) != MP_OK) {
 221       break;
 222     }
 223
 224     if ((res = mp_int_is_prime(&t)) == MP_TRUE) {
 225       if (fb != NULL) fputc('!', fb);
 226
 227       res = mp_int_copy(&t, seed);
 228       break;
 229     } else if (res != MP_FALSE)
 230       break;
 231
 232     if (fb != NULL) fputc('x', fb);
 233     if (mp_int_add_value(seed, 2, seed) != MP_OK) break;
 234   }
 235
 236   mp_int_clear(&t);
 237   return res;
 238 }
 239
 240 /* Here there be dragons */