source4/heimdal/lib/des/imath/iprime.c

   1 /*
   2   Name:     iprime.c
   3   Purpose:  Pseudoprimality testing routines
   4   Author:   M. J. Fromberger <http://www.dartmouth.edu/~sting/>
   5   Info:     $Id: iprime.c,v 1.5 2007/01/05 21:01:48 lha Exp $
   6
   7   Copyright (C) 2002 Michael J. Fromberger, All Rights Reserved.
   8
   9   Permission is hereby granted, free of charge, to any person
  10   obtaining a copy of this software and associated documentation files
  11   (the "Software"), to deal in the Software without restriction,
  12   including without limitation the rights to use, copy, modify, merge,
  13   publish, distribute, sublicense, and/or sell copies of the Software,
  14   and to permit persons to whom the Software is furnished to do so,
  15   subject to the following conditions:
  16
  17   The above copyright notice and this permission notice shall be
  18   included in all copies or substantial portions of the Software.
  19
  20   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
  21   EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
  22   MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
  23   NONINFRINGEMENT.  IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
  24   BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
  25   ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
  26   CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
  27   SOFTWARE.
  28  */
  29
  30 #include "iprime.h"
  31 #include <stdlib.h>
  32
  33 static const int s_ptab[] = {
  34     3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43,
  35     47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101,
  36     103, 107, 109, 113, 127, 131, 137, 139, 149, 151,
  37     157, 163, 167, 173, 179, 181, 191, 193, 197, 199,
  38     211, 223, 227, 229, 233, 239, 241, 251, 257, 263,
  39     269, 271, 277, 281, 283, 293, 307, 311, 313, 317,
  40     331, 337, 347, 349, 353, 359, 367, 373, 379, 383,
  41     389, 397, 401, 409, 419, 421, 431, 433, 439, 443,
  42     449, 457, 461, 463, 467, 479, 487, 491, 499, 503,
  43     509, 521, 523, 541, 547, 557, 563, 569, 571, 577,
  44     587, 593, 599, 601, 607, 613, 617, 619, 631, 641,
  45     643, 647, 653, 659, 661, 673, 677, 683, 691, 701,
  46     709, 719, 727, 733, 739, 743, 751, 757, 761, 769,
  47     773, 787, 797, 809, 811, 821, 823, 827, 829, 839,
  48     853, 857, 859, 863, 877, 881, 883, 887, 907, 911,
  49     919, 929, 937, 941, 947, 953, 967, 971, 977, 983,
  50     991, 997, 1009, 1013, 1019, 1021, 1031, 1033,
  51     1039, 1049, 1051, 1061, 1063, 1069, 1087, 1091,
  52     1093, 1097, 1103, 1109, 1117, 1123, 1129, 1151,
  53     1153, 1163, 1171, 1181, 1187, 1193, 1201, 1213,
  54     1217, 1223, 1229, 1231, 1237, 1249, 1259, 1277,
  55     1279, 1283, 1289, 1291, 1297, 1301, 1303, 1307,
  56     1319, 1321, 1327, 1361, 1367, 1373, 1381, 1399,
  57     1409, 1423, 1427, 1429, 1433, 1439, 1447, 1451,
  58     1453, 1459, 1471, 1481, 1483, 1487, 1489, 1493,
  59     1499, 1511, 1523, 1531, 1543, 1549, 1553, 1559,
  60     1567, 1571, 1579, 1583, 1597, 1601, 1607, 1609,
  61     1613, 1619, 1621, 1627, 1637, 1657, 1663, 1667,
  62     1669, 1693, 1697, 1699, 1709, 1721, 1723, 1733,
  63     1741, 1747, 1753, 1759, 1777, 1783, 1787, 1789,
  64     1801, 1811, 1823, 1831, 1847, 1861, 1867, 1871,
  65     1873, 1877, 1879, 1889, 1901, 1907, 1913, 1931,
  66     1933, 1949, 1951, 1973, 1979, 1987, 1993, 1997,
  67     1999, 2003, 2011, 2017, 2027, 2029, 2039, 2053,
  68     2063, 2069, 2081, 2083, 2087, 2089, 2099, 2111,
  69     2113, 2129, 2131, 2137, 2141, 2143, 2153, 2161,
  70     2179, 2203, 2207, 2213, 2221, 2237, 2239, 2243,
  71     2251, 2267, 2269, 2273, 2281, 2287, 2293, 2297,
  72     2309, 2311, 2333, 2339, 2341, 2347, 2351, 2357,
  73     2371, 2377, 2381, 2383, 2389, 2393, 2399, 2411,
  74     2417, 2423, 2437, 2441, 2447, 2459, 2467, 2473,
  75     2477, 2503, 2521, 2531, 2539, 2543, 2549, 2551,
  76     2557, 2579, 2591, 2593, 2609, 2617, 2621, 2633,
  77     2647, 2657, 2659, 2663, 2671, 2677, 2683, 2687,
  78     2689, 2693, 2699, 2707, 2711, 2713, 2719, 2729,
  79     2731, 2741, 2749, 2753, 2767, 2777, 2789, 2791,
  80     2797, 2801, 2803, 2819, 2833, 2837, 2843, 2851,
  81     2857, 2861, 2879, 2887, 2897, 2903, 2909, 2917,
  82     2927, 2939, 2953, 2957, 2963, 2969, 2971, 2999,
  83     3001, 3011, 3019, 3023, 3037, 3041, 3049, 3061,
  84     3067, 3079, 3083, 3089, 3109, 3119, 3121, 3137,
  85     3163, 3167, 3169, 3181, 3187, 3191, 3203, 3209,
  86     3217, 3221, 3229, 3251, 3253, 3257, 3259, 3271,
  87     3299, 3301, 3307, 3313, 3319, 3323, 3329, 3331,
  88     3343, 3347, 3359, 3361, 3371, 3373, 3389, 3391,
  89     3407, 3413, 3433, 3449, 3457, 3461, 3463, 3467,
  90     3469, 3491, 3499, 3511, 3517, 3527, 3529, 3533,
  91     3539, 3541, 3547, 3557, 3559, 3571, 3581, 3583,
  92     3593, 3607, 3613, 3617, 3623, 3631, 3637, 3643,
  93     3659, 3671, 3673, 3677, 3691, 3697, 3701, 3709,
  94     3719, 3727, 3733, 3739, 3761, 3767, 3769, 3779,
  95     3793, 3797, 3803, 3821, 3823, 3833, 3847, 3851,
  96     3853, 3863, 3877, 3881, 3889, 3907, 3911, 3917,
  97     3919, 3923, 3929, 3931, 3943, 3947, 3967, 3989,
  98     4001, 4003, 4007, 4013, 4019, 4021, 4027, 4049,
  99     4051, 4057, 4073, 4079, 4091, 4093, 4099, 4111,
 100     4127, 4129, 4133, 4139, 4153, 4157, 4159, 4177,
 101     4201, 4211, 4217, 4219, 4229, 4231, 4241, 4243,
 102     4253, 4259, 4261, 4271, 4273, 4283, 4289, 4297,
 103     4327, 4337, 4339, 4349, 4357, 4363, 4373, 4391,
 104     4397, 4409, 4421, 4423, 4441, 4447, 4451, 4457,
 105     4463, 4481, 4483, 4493, 4507, 4513, 4517, 4519,
 106     4523, 4547, 4549, 4561, 4567, 4583, 4591, 4597,
 107     4603, 4621, 4637, 4639, 4643, 4649, 4651, 4657,
 108     4663, 4673, 4679, 4691, 4703, 4721, 4723, 4729,
 109     4733, 4751, 4759, 4783, 4787, 4789, 4793, 4799,
 110     4801, 4813, 4817, 4831, 4861, 4871, 4877, 4889,
 111     4903, 4909, 4919, 4931, 4933, 4937, 4943, 4951,
 112     4957, 4967, 4969, 4973, 4987, 4993, 4999
 113 };
 114 static const int s_ptab_size = sizeof(s_ptab)/sizeof(s_ptab[0]);
 115
 116
 117 /* {{{ mp_int_is_prime(z) */
 118
 119 /* Test whether z is likely to be prime:
 120    MP_TRUE  means it is probably prime
 121    MP_FALSE means it is definitely composite
 122  */
 123 mp_result mp_int_is_prime(mp_int z)
 124 {
 125   int       i, rem;
 126   mp_result res;
 127
 128   /* First check for divisibility by small primes; this eliminates a
 129      large number of composite candidates quickly
 130    */
 131   for(i = 0; i < s_ptab_size; ++i) {
 132     if((res = mp_int_div_value(z, s_ptab[i], NULL, &rem)) != MP_OK)
 133       return res;
 134
 135     if(rem == 0)
 136       return MP_FALSE;
 137   }
 138
 139   /* Now try Fermat's test for several prime witnesses (since we now
 140      know from the above that z is not a multiple of any of them)
 141    */
 142   {
 143     mpz_t  tmp;
 144
 145     if((res = mp_int_init(&tmp)) != MP_OK) return res;
 146
 147     for(i = 0; i < 10 && i < s_ptab_size; ++i) {
 148       if((res = mp_int_exptmod_bvalue(s_ptab[i], z, z, &tmp)) != MP_OK)
 149         return res;
 150
 151       if(mp_int_compare_value(&tmp, s_ptab[i]) != 0) {
 152         mp_int_clear(&tmp);
 153         return MP_FALSE;
 154       }
 155     }
 156
 157     mp_int_clear(&tmp);
 158   }
 159
 160   return MP_TRUE;
 161 }
 162
 163 /* }}} */
 164
 165 /* {{{ mp_int_find_prime(z) */
 166
 167 /* Find the first apparent prime in ascending order from z */
 168 mp_result mp_int_find_prime(mp_int z)
 169 {
 170   mp_result  res;
 171
 172   if(mp_int_is_even(z) && ((res = mp_int_add_value(z, 1, z)) != MP_OK))
 173     return res;
 174
 175   while((res = mp_int_is_prime(z)) == MP_FALSE) {
 176     if((res = mp_int_add_value(z, 2, z)) != MP_OK)
 177       break;
 178
 179   }
 180
 181   return res;
 182 }
 183
 184 /* }}} */
 185
 186 /* Here there be dragons */