lib-src/numeric/rabin.cc

   1 //////////////////////////////////////////////////////////////////////////////
   2 // NOTICE:
   3 //
   4 // ADLib, Prop and their related set of tools and documentation are in the
   5 // public domain.   The author(s) of this software reserve no copyrights on
   6 // the source code and any code generated using the tools.  You are encouraged
   7 // to use ADLib and Prop to develop software, in both academic and commercial
   8 // settings, and are free to incorporate any part of ADLib and Prop into
   9 // your programs.
  10 //
  11 // Although you are under no obligation to do so, we strongly recommend that
  12 // you give away all software developed using our tools.
  13 //
  14 // We also ask that credit be given to us when ADLib and/or Prop are used in
  15 // your programs, and that this notice be preserved intact in all the source
  16 // code.
  17 //
  18 // This software is still under development and we welcome any suggestions
  19 // and help from the users.
  20 //
  21 // Allen Leung
  22 // 1994
  23 //////////////////////////////////////////////////////////////////////////////
  24
  25 #include <AD/numeric/primes.h>
  26
  27 ////////////////////////////////////////////////////////////////////////////
  28 // Rabin primality checking\cite{RSA, random-algorithms, number-theory}
  29 ////////////////////////////////////////////////////////////////////////////
  30
  31 Bool Rabin_isPrime(const BigInt& n, int trials)
  32 {
  33    BigInt n_minus_1 = n - 1;
  34
  35    //////////////////////////////////////////////////////////////////////
  36    // First, try some small primes to prune out blatant composites.
  37    //////////////////////////////////////////////////////////////////////
  38
  39    if (pow_mod(2 * 3 * 5 * 7 * 11 * 13 * 17 * 23 * 31,n_minus_1,n) != 1)
  40       return false;
  41
  42
  43    for ( ; trials > 0; trials--) {
  44       BigInt x = make_random(n_minus_1) + 1; // get a random number from 1 to n-1
  45
  46       //
  47       // Try testing for composite using Fermat's(Euler's) theorem:
  48       //    x^(phi(n)) = 1 (mod n) iff n is relatively prime to x.
  49       //  where phi(n) is the number of elements relatively prime to
  50       //  n from 1 to n.
  51       //
  52
  53       if (pow_mod(x,n_minus_1,n) != 1) return false;
  54
  55       //
  56       // Then try Lehmer's heuristics:
  57       //  n is prime iff
  58       //        (1) x^(n-1)     =  1 (mod n) and
  59       //        (2) x^((n-1)/p) <> 1 (mod n) for each prime factor p of n-1
  60       //
  61    }
  62
  63    /////////////////////////////////////////////////////////////////////
  64    //  N is a prime with probability >= 1 - 1/(2^trials)
  65    /////////////////////////////////////////////////////////////////////
  66
  67    return true;
  68 }