lib-src/numeric/solovay.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 //  Primality testing by Solovay and Strassen\cite{random-algorithms}
  29 //
  30
  31 Bool Solovay_Strassen_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 the Euclidean algorithm
  48       //
  49       if (gcd(x,n) != 1) return false;
  50
  51       //
  52       // Then try testing for quadratic
  53       // residues\cite{number-theory} (chapter 3)
  54       //
  55    }
  56
  57    //
  58    //  N is a prime with probability >= 1 - 1/(2^trials)
  59    //
  60
  61    return true;
  62 }