Botan  1.11.15
src/lib/math/numbertheory/numthry.cpp
Go to the documentation of this file.
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 }