Home | History | Annotate | Download | only in libtommath 
      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 (at) gmail.com, http://math.libtomcrypt.com
     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, 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.3 $ */
    170 /* $Date: 2006/03/31 14:18:44 $ */
    171