source4/heimdal/lib/hcrypto/libtommath/bn_mp_prime_next_prime.c

   1 #include <tommath.h>
   2 #ifdef BN_MP_PRIME_NEXT_PRIME_C
   3 /* LibTomMath, multiple-precision integer library -- Tom St Denis
   4  *
   5  * LibTomMath is a library that provides multiple-precision
   6  * integer arithmetic as well as number theoretic functionality.
   7  *
   8  * The library was designed directly after the MPI library by
   9  * Michael Fromberger but has been written from scratch with
  10  * additional optimizations in place.
  11  *
  12  * The library is free for all purposes without any express
  13  * guarantee it works.
  14  *
  15  * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
  16  */
  17
  18 /* finds the next prime after the number "a" using "t" trials
  19  * of Miller-Rabin.
  20  *
  21  * bbs_style = 1 means the prime must be congruent to 3 mod 4
  22  */
  23 int mp_prime_next_prime(mp_int *a, int t, int bbs_style)
  24 {
  25    int      err, res = MP_NO, x, y;
  26    mp_digit res_tab[PRIME_SIZE], step, kstep;
  27    mp_int   b;
  28
  29    /* ensure t is valid */
  30    if (t <= 0 || t > PRIME_SIZE) {
  31       return MP_VAL;
  32    }
  33
  34    /* force positive */
  35    a->sign = MP_ZPOS;
  36
  37    /* simple algo if a is less than the largest prime in the table */
  38    if (mp_cmp_d(a, ltm_prime_tab[PRIME_SIZE-1]) == MP_LT) {
  39       /* find which prime it is bigger than */
  40       for (x = PRIME_SIZE - 2; x >= 0; x--) {
  41           if (mp_cmp_d(a, ltm_prime_tab[x]) != MP_LT) {
  42              if (bbs_style == 1) {
  43                 /* ok we found a prime smaller or
  44                  * equal [so the next is larger]
  45                  *
  46                  * however, the prime must be
  47                  * congruent to 3 mod 4
  48                  */
  49                 if ((ltm_prime_tab[x + 1] & 3) != 3) {
  50                    /* scan upwards for a prime congruent to 3 mod 4 */
  51                    for (y = x + 1; y < PRIME_SIZE; y++) {
  52                        if ((ltm_prime_tab[y] & 3) == 3) {
  53                           mp_set(a, ltm_prime_tab[y]);
  54                           return MP_OKAY;
  55                        }
  56                    }
  57                 }
  58              } else {
  59                 mp_set(a, ltm_prime_tab[x + 1]);
  60                 return MP_OKAY;
  61              }
  62           }
  63       }
  64       /* at this point a maybe 1 */
  65       if (mp_cmp_d(a, 1) == MP_EQ) {
  66          mp_set(a, 2);
  67          return MP_OKAY;
  68       }
  69       /* fall through to the sieve */
  70    }
  71
  72    /* generate a prime congruent to 3 mod 4 or 1/3 mod 4? */
  73    if (bbs_style == 1) {
  74       kstep   = 4;
  75    } else {
  76       kstep   = 2;
  77    }
  78
  79    /* at this point we will use a combination of a sieve and Miller-Rabin */
  80
  81    if (bbs_style == 1) {
  82       /* if a mod 4 != 3 subtract the correct value to make it so */
  83       if ((a->dp[0] & 3) != 3) {
  84          if ((err = mp_sub_d(a, (a->dp[0] & 3) + 1, a)) != MP_OKAY) { return err; };
  85       }
  86    } else {
  87       if (mp_iseven(a) == 1) {
  88          /* force odd */
  89          if ((err = mp_sub_d(a, 1, a)) != MP_OKAY) {
  90             return err;
  91          }
  92       }
  93    }
  94
  95    /* generate the restable */
  96    for (x = 1; x < PRIME_SIZE; x++) {
  97       if ((err = mp_mod_d(a, ltm_prime_tab[x], res_tab + x)) != MP_OKAY) {
  98          return err;
  99       }
 100    }
 101
 102    /* init temp used for Miller-Rabin Testing */
 103    if ((err = mp_init(&b)) != MP_OKAY) {
 104       return err;
 105    }
 106
 107    for (;;) {
 108       /* skip to the next non-trivially divisible candidate */
 109       step = 0;
 110       do {
 111          /* y == 1 if any residue was zero [e.g. cannot be prime] */
 112          y     =  0;
 113
 114          /* increase step to next candidate */
 115          step += kstep;
 116
 117          /* compute the new residue without using division */
 118          for (x = 1; x < PRIME_SIZE; x++) {
 119              /* add the step to each residue */
 120              res_tab[x] += kstep;
 121
 122              /* subtract the modulus [instead of using division] */
 123              if (res_tab[x] >= ltm_prime_tab[x]) {
 124                 res_tab[x]  -= ltm_prime_tab[x];
 125              }
 126
 127              /* set flag if zero */
 128              if (res_tab[x] == 0) {
 129                 y = 1;
 130              }
 131          }
 132       } while (y == 1 && step < ((((mp_digit)1)<<DIGIT_BIT) - kstep));
 133
 134       /* add the step */
 135       if ((err = mp_add_d(a, step, a)) != MP_OKAY) {
 136          goto LBL_ERR;
 137       }
 138
 139       /* if didn't pass sieve and step == MAX then skip test */
 140       if (y == 1 && step >= ((((mp_digit)1)<<DIGIT_BIT) - kstep)) {
 141          continue;
 142       }
 143
 144       /* is this prime? */
 145       for (x = 0; x < t; x++) {
 146           mp_set(&b, ltm_prime_tab[t]);
 147           if ((err = mp_prime_miller_rabin(a, &b, &res)) != MP_OKAY) {
 148              goto LBL_ERR;
 149           }
 150           if (res == MP_NO) {
 151              break;
 152           }
 153       }
 154
 155       if (res == MP_YES) {
 156          break;
 157       }
 158    }
 159
 160    err = MP_OKAY;
 161 LBL_ERR:
 162    mp_clear(&b);
 163    return err;
 164 }
 165
 166 #endif
 167
 168 /* $Source: /cvs/libtom/libtommath/bn_mp_prime_next_prime.c,v $ */
 169 /* $Revision: 1.4 $ */
 170 /* $Date: 2006/12/28 01:25:13 $ */