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

   1 /*
   2   Name:     iprime.c
   3   Purpose:  Pseudoprimality testing routines
   4   Author:   M. J. Fromberger <http://spinning-yarns.org/michael/>
   5   Info:     $Id: iprime.c 635 2008-01-08 18:19:40Z sting $
   6
   7   Copyright (C) 2002-2008 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
  51 #ifdef IMATH_LARGE_PRIME_TABLE
  52     , 1009, 1013, 1019, 1021, 1031, 1033,
  53     1039, 1049, 1051, 1061, 1063, 1069, 1087, 1091,
  54     1093, 1097, 1103, 1109, 1117, 1123, 1129, 1151,
  55     1153, 1163, 1171, 1181, 1187, 1193, 1201, 1213,
  56     1217, 1223, 1229, 1231, 1237, 1249, 1259, 1277,
  57     1279, 1283, 1289, 1291, 1297, 1301, 1303, 1307,
  58     1319, 1321, 1327, 1361, 1367, 1373, 1381, 1399,
  59     1409, 1423, 1427, 1429, 1433, 1439, 1447, 1451,
  60     1453, 1459, 1471, 1481, 1483, 1487, 1489, 1493,
  61     1499, 1511, 1523, 1531, 1543, 1549, 1553, 1559,
  62     1567, 1571, 1579, 1583, 1597, 1601, 1607, 1609,
  63     1613, 1619, 1621, 1627, 1637, 1657, 1663, 1667,
  64     1669, 1693, 1697, 1699, 1709, 1721, 1723, 1733,
  65     1741, 1747, 1753, 1759, 1777, 1783, 1787, 1789,
  66     1801, 1811, 1823, 1831, 1847, 1861, 1867, 1871,
  67     1873, 1877, 1879, 1889, 1901, 1907, 1913, 1931,
  68     1933, 1949, 1951, 1973, 1979, 1987, 1993, 1997,
  69     1999, 2003, 2011, 2017, 2027, 2029, 2039, 2053,
  70     2063, 2069, 2081, 2083, 2087, 2089, 2099, 2111,
  71     2113, 2129, 2131, 2137, 2141, 2143, 2153, 2161,
  72     2179, 2203, 2207, 2213, 2221, 2237, 2239, 2243,
  73     2251, 2267, 2269, 2273, 2281, 2287, 2293, 2297,
  74     2309, 2311, 2333, 2339, 2341, 2347, 2351, 2357,
  75     2371, 2377, 2381, 2383, 2389, 2393, 2399, 2411,
  76     2417, 2423, 2437, 2441, 2447, 2459, 2467, 2473,
  77     2477, 2503, 2521, 2531, 2539, 2543, 2549, 2551,
  78     2557, 2579, 2591, 2593, 2609, 2617, 2621, 2633,
  79     2647, 2657, 2659, 2663, 2671, 2677, 2683, 2687,
  80     2689, 2693, 2699, 2707, 2711, 2713, 2719, 2729,
  81     2731, 2741, 2749, 2753, 2767, 2777, 2789, 2791,
  82     2797, 2801, 2803, 2819, 2833, 2837, 2843, 2851,
  83     2857, 2861, 2879, 2887, 2897, 2903, 2909, 2917,
  84     2927, 2939, 2953, 2957, 2963, 2969, 2971, 2999,
  85     3001, 3011, 3019, 3023, 3037, 3041, 3049, 3061,
  86     3067, 3079, 3083, 3089, 3109, 3119, 3121, 3137,
  87     3163, 3167, 3169, 3181, 3187, 3191, 3203, 3209,
  88     3217, 3221, 3229, 3251, 3253, 3257, 3259, 3271,
  89     3299, 3301, 3307, 3313, 3319, 3323, 3329, 3331,
  90     3343, 3347, 3359, 3361, 3371, 3373, 3389, 3391,
  91     3407, 3413, 3433, 3449, 3457, 3461, 3463, 3467,
  92     3469, 3491, 3499, 3511, 3517, 3527, 3529, 3533,
  93     3539, 3541, 3547, 3557, 3559, 3571, 3581, 3583,
  94     3593, 3607, 3613, 3617, 3623, 3631, 3637, 3643,
  95     3659, 3671, 3673, 3677, 3691, 3697, 3701, 3709,
  96     3719, 3727, 3733, 3739, 3761, 3767, 3769, 3779,
  97     3793, 3797, 3803, 3821, 3823, 3833, 3847, 3851,
  98     3853, 3863, 3877, 3881, 3889, 3907, 3911, 3917,
  99     3919, 3923, 3929, 3931, 3943, 3947, 3967, 3989,
 100     4001, 4003, 4007, 4013, 4019, 4021, 4027, 4049,
 101     4051, 4057, 4073, 4079, 4091, 4093, 4099, 4111,
 102     4127, 4129, 4133, 4139, 4153, 4157, 4159, 4177,
 103     4201, 4211, 4217, 4219, 4229, 4231, 4241, 4243,
 104     4253, 4259, 4261, 4271, 4273, 4283, 4289, 4297,
 105     4327, 4337, 4339, 4349, 4357, 4363, 4373, 4391,
 106     4397, 4409, 4421, 4423, 4441, 4447, 4451, 4457,
 107     4463, 4481, 4483, 4493, 4507, 4513, 4517, 4519,
 108     4523, 4547, 4549, 4561, 4567, 4583, 4591, 4597,
 109     4603, 4621, 4637, 4639, 4643, 4649, 4651, 4657,
 110     4663, 4673, 4679, 4691, 4703, 4721, 4723, 4729,
 111     4733, 4751, 4759, 4783, 4787, 4789, 4793, 4799,
 112     4801, 4813, 4817, 4831, 4861, 4871, 4877, 4889,
 113     4903, 4909, 4919, 4931, 4933, 4937, 4943, 4951,
 114     4957, 4967, 4969, 4973, 4987, 4993, 4999
 115 #endif
 116 };
 117 static const int s_ptab_size = sizeof(s_ptab)/sizeof(s_ptab[0]);
 118
 119 /* {{{ mp_int_is_prime(z) */
 120
 121 /* Test whether z is likely to be prime:
 122    MP_TRUE  means it is probably prime
 123    MP_FALSE means it is definitely composite
 124  */
 125 mp_result mp_int_is_prime(mp_int z)
 126 {
 127   int       i;
 128   mp_small  rem;
 129   mp_result res;
 130
 131   /* First check for divisibility by small primes; this eliminates a
 132      large number of composite candidates quickly
 133    */
 134   for(i = 0; i < s_ptab_size; ++i) {
 135     if((res = mp_int_div_value(z, s_ptab[i], NULL, &rem)) != MP_OK)
 136       return res;
 137
 138     if(rem == 0)
 139       return MP_FALSE;
 140   }
 141
 142   /* Now try Fermat's test for several prime witnesses (since we now
 143      know from the above that z is not a multiple of any of them)
 144    */
 145   {
 146     mpz_t  tmp;
 147
 148     if((res = mp_int_init(&tmp)) != MP_OK) return res;
 149
 150     for(i = 0; i < 10 && i < s_ptab_size; ++i) {
 151       if((res = mp_int_exptmod_bvalue(s_ptab[i], z, z, &tmp)) != MP_OK)
 152         return res;
 153
 154       if(mp_int_compare_value(&tmp, s_ptab[i]) != 0) {
 155         mp_int_clear(&tmp);
 156         return MP_FALSE;
 157       }
 158     }
 159
 160     mp_int_clear(&tmp);
 161   }
 162
 163   return MP_TRUE;
 164 }
 165
 166 /* }}} */
 167
 168 /* {{{ mp_int_find_prime(z) */
 169
 170 /* Find the first apparent prime in ascending order from z */
 171 mp_result mp_int_find_prime(mp_int z)
 172 {
 173   mp_result  res;
 174
 175   if(mp_int_is_even(z) && ((res = mp_int_add_value(z, 1, z)) != MP_OK))
 176     return res;
 177
 178   while((res = mp_int_is_prime(z)) == MP_FALSE) {
 179     if((res = mp_int_add_value(z, 2, z)) != MP_OK)
 180       break;
 181
 182   }
 183
 184   return res;
 185 }
 186
 187 /* }}} */
 188
 189 /* Here there be dragons */