Botan
1.11.15
|
00001 /* 00002 * Number Theory Functions 00003 * (C) 1999-2011 Jack Lloyd 00004 * 00005 * Botan is released under the Simplified BSD License (see license.txt) 00006 */ 00007 00008 #include <botan/numthry.h> 00009 #include <botan/reducer.h> 00010 #include <botan/internal/bit_ops.h> 00011 #include <botan/internal/mp_core.h> 00012 #include <algorithm> 00013 00014 namespace Botan { 00015 00016 /* 00017 * Return the number of 0 bits at the end of n 00018 */ 00019 size_t low_zero_bits(const BigInt& n) 00020 { 00021 size_t low_zero = 0; 00022 00023 if(n.is_positive() && n.is_nonzero()) 00024 { 00025 for(size_t i = 0; i != n.size(); ++i) 00026 { 00027 const word x = n.word_at(i); 00028 00029 if(x) 00030 { 00031 low_zero += ctz(x); 00032 break; 00033 } 00034 else 00035 low_zero += BOTAN_MP_WORD_BITS; 00036 } 00037 } 00038 00039 return low_zero; 00040 } 00041 00042 /* 00043 * Calculate the GCD 00044 */ 00045 BigInt gcd(const BigInt& a, const BigInt& b) 00046 { 00047 if(a.is_zero() || b.is_zero()) return 0; 00048 if(a == 1 || b == 1) return 1; 00049 00050 BigInt x = a, y = b; 00051 x.set_sign(BigInt::Positive); 00052 y.set_sign(BigInt::Positive); 00053 size_t shift = std::min(low_zero_bits(x), low_zero_bits(y)); 00054 00055 x >>= shift; 00056 y >>= shift; 00057 00058 while(x.is_nonzero()) 00059 { 00060 x >>= low_zero_bits(x); 00061 y >>= low_zero_bits(y); 00062 if(x >= y) { x -= y; x >>= 1; } 00063 else { y -= x; y >>= 1; } 00064 } 00065 00066 return (y << shift); 00067 } 00068 00069 /* 00070 * Calculate the LCM 00071 */ 00072 BigInt lcm(const BigInt& a, const BigInt& b) 00073 { 00074 return ((a * b) / gcd(a, b)); 00075 } 00076 00077 namespace { 00078 00079 /* 00080 * If the modulus is odd, then we can avoid computing A and C. This is 00081 * a critical path algorithm in some instances and an odd modulus is 00082 * the common case for crypto, so worth special casing. See note 14.64 00083 * in Handbook of Applied Cryptography for more details. 00084 */ 00085 BigInt inverse_mod_odd_modulus(const BigInt& n, const BigInt& mod) 00086 { 00087 BigInt u = mod, v = n; 00088 BigInt B = 0, D = 1; 00089 00090 while(u.is_nonzero()) 00091 { 00092 const size_t u_zero_bits = low_zero_bits(u); 00093 u >>= u_zero_bits; 00094 for(size_t i = 0; i != u_zero_bits; ++i) 00095 { 00096 if(B.is_odd()) 00097 { B -= mod; } 00098 B >>= 1; 00099 } 00100 00101 const size_t v_zero_bits = low_zero_bits(v); 00102 v >>= v_zero_bits; 00103 for(size_t i = 0; i != v_zero_bits; ++i) 00104 { 00105 if(D.is_odd()) 00106 { D -= mod; } 00107 D >>= 1; 00108 } 00109 00110 if(u >= v) { u -= v; B -= D; } 00111 else { v -= u; D -= B; } 00112 } 00113 00114 if(v != 1) 00115 return 0; // no modular inverse 00116 00117 while(D.is_negative()) D += mod; 00118 while(D >= mod) D -= mod; 00119 00120 return D; 00121 } 00122 00123 } 00124 00125 /* 00126 * Find the Modular Inverse 00127 */ 00128 BigInt inverse_mod(const BigInt& n, const BigInt& mod) 00129 { 00130 if(mod.is_zero()) 00131 throw BigInt::DivideByZero(); 00132 if(mod.is_negative() || n.is_negative()) 00133 throw Invalid_Argument("inverse_mod: arguments must be non-negative"); 00134 00135 if(n.is_zero() || (n.is_even() && mod.is_even())) 00136 return 0; // fast fail checks 00137 00138 if(mod.is_odd()) 00139 return inverse_mod_odd_modulus(n, mod); 00140 00141 BigInt u = mod, v = n; 00142 BigInt A = 1, B = 0, C = 0, D = 1; 00143 00144 while(u.is_nonzero()) 00145 { 00146 const size_t u_zero_bits = low_zero_bits(u); 00147 u >>= u_zero_bits; 00148 for(size_t i = 0; i != u_zero_bits; ++i) 00149 { 00150 if(A.is_odd() || B.is_odd()) 00151 { A += n; B -= mod; } 00152 A >>= 1; B >>= 1; 00153 } 00154 00155 const size_t v_zero_bits = low_zero_bits(v); 00156 v >>= v_zero_bits; 00157 for(size_t i = 0; i != v_zero_bits; ++i) 00158 { 00159 if(C.is_odd() || D.is_odd()) 00160 { C += n; D -= mod; } 00161 C >>= 1; D >>= 1; 00162 } 00163 00164 if(u >= v) { u -= v; A -= C; B -= D; } 00165 else { v -= u; C -= A; D -= B; } 00166 } 00167 00168 if(v != 1) 00169 return 0; // no modular inverse 00170 00171 while(D.is_negative()) D += mod; 00172 while(D >= mod) D -= mod; 00173 00174 return D; 00175 } 00176 00177 word monty_inverse(word input) 00178 { 00179 word b = input; 00180 word x2 = 1, x1 = 0, y2 = 0, y1 = 1; 00181 00182 // First iteration, a = n+1 00183 word q = bigint_divop(1, 0, b); 00184 word r = (MP_WORD_MAX - q*b) + 1; 00185 word x = x2 - q*x1; 00186 word y = y2 - q*y1; 00187 00188 word a = b; 00189 b = r; 00190 x2 = x1; 00191 x1 = x; 00192 y2 = y1; 00193 y1 = y; 00194 00195 while(b > 0) 00196 { 00197 q = a / b; 00198 r = a - q*b; 00199 x = x2 - q*x1; 00200 y = y2 - q*y1; 00201 00202 a = b; 00203 b = r; 00204 x2 = x1; 00205 x1 = x; 00206 y2 = y1; 00207 y1 = y; 00208 } 00209 00210 // Now invert in addition space 00211 y2 = (MP_WORD_MAX - y2) + 1; 00212 00213 return y2; 00214 } 00215 00216 /* 00217 * Modular Exponentiation 00218 */ 00219 BigInt power_mod(const BigInt& base, const BigInt& exp, const BigInt& mod) 00220 { 00221 Power_Mod pow_mod(mod); 00222 00223 /* 00224 * Calling set_base before set_exponent means we end up using a 00225 * minimal window. This makes sense given that here we know that any 00226 * precomputation is wasted. 00227 */ 00228 pow_mod.set_base(base); 00229 pow_mod.set_exponent(exp); 00230 return pow_mod.execute(); 00231 } 00232 00233 namespace { 00234 00235 bool mr_witness(BigInt&& y, 00236 const Modular_Reducer& reducer_n, 00237 const BigInt& n_minus_1, size_t s) 00238 { 00239 if(y == 1 || y == n_minus_1) 00240 return false; 00241 00242 for(size_t i = 1; i != s; ++i) 00243 { 00244 y = reducer_n.square(y); 00245 00246 if(y == 1) // found a non-trivial square root 00247 return true; 00248 00249 if(y == n_minus_1) // -1, trivial square root, so give up 00250 return false; 00251 } 00252 00253 return true; // fails Fermat test 00254 } 00255 00256 size_t mr_test_iterations(size_t n_bits, size_t prob, bool random) 00257 { 00258 const size_t base = (prob + 2) / 2; // worst case 4^-t error rate 00259 00260 /* 00261 * For randomly chosen numbers we can use the estimates from 00262 * http://www.math.dartmouth.edu/~carlp/PDF/paper88.pdf 00263 * 00264 * These values are derived from the inequality for p(k,t) given on 00265 * the second page. 00266 */ 00267 if(random && prob <= 80) 00268 { 00269 if(n_bits >= 1536) 00270 return 2; // < 2^-89 00271 if(n_bits >= 1024) 00272 return 4; // < 2^-89 00273 if(n_bits >= 512) 00274 return 5; // < 2^-80 00275 if(n_bits >= 256) 00276 return 11; // < 2^-80 00277 } 00278 00279 return base; 00280 } 00281 00282 } 00283 00284 /* 00285 * Test for primaility using Miller-Rabin 00286 */ 00287 bool is_prime(const BigInt& n, RandomNumberGenerator& rng, 00288 size_t prob, bool is_random) 00289 { 00290 if(n == 2) 00291 return true; 00292 if(n <= 1 || n.is_even()) 00293 return false; 00294 00295 // Fast path testing for small numbers (<= 65521) 00296 if(n <= PRIMES[PRIME_TABLE_SIZE-1]) 00297 { 00298 const u16bit num = n.word_at(0); 00299 00300 return std::binary_search(PRIMES, PRIMES + PRIME_TABLE_SIZE, num); 00301 } 00302 00303 const size_t test_iterations = mr_test_iterations(n.bits(), prob, is_random); 00304 00305 const BigInt n_minus_1 = n - 1; 00306 const size_t s = low_zero_bits(n_minus_1); 00307 00308 Fixed_Exponent_Power_Mod pow_mod(n_minus_1 >> s, n); 00309 Modular_Reducer reducer(n); 00310 00311 for(size_t i = 0; i != test_iterations; ++i) 00312 { 00313 const BigInt a = BigInt::random_integer(rng, 2, n_minus_1); 00314 00315 BigInt y = pow_mod(a); 00316 00317 if(mr_witness(std::move(y), reducer, n_minus_1, s)) 00318 return false; 00319 } 00320 00321 return true; 00322 } 00323 00324 }