home *** CD-ROM | disk | FTP | other *** search
/ NetNews Usenet Archive 1992 #26 / NN_1992_26.iso / spool / comp / programm / 3126 < prev    next >
Encoding:
Text File  |  1992-11-11  |  5.6 KB  |  157 lines
  1. Newsgroups: comp.programming
  2. Path: sparky!uunet!mcsun!Germany.EU.net!Urmel.Informatik.RWTH-Aachen.DE!martin
  3. From: martin@math.rwth-aachen.de (  Martin Schoenert)
  4. Subject: Re: Probabilistic Prime Numbers
  5. Message-ID: <martin.721557289@bert>
  6. Sender: news@Urmel.Informatik.RWTH-Aachen.DE (Newsfiles Owner)
  7. Nntp-Posting-Host: bert.math.rwth-aachen.de
  8. Organization: Rechnerbetrieb Informatik  /  RWTH Aachen
  9. References: <2185@snap>
  10. Date: 12 Nov 92 08:34:49 GMT
  11. Lines: 144
  12.  
  13. paj@uk.co.gec-mrc (Paul Johnson) writes:
  14.  
  15.     I am trying to write a small prime number generator using Algorithm P
  16.     in Knuth.
  17.  
  18.     Note for those without a copy of Knuth: this algorithm uses a randomly
  19.     generated number to perform a primality test.  It can return one of
  20.     two results: "Not prime" or "Probably prime".  If the latter, then the
  21.     chance of error is < 0.25.  Repeated tests with different random
  22.     numbers can reduce this probability very quickly.
  23.  
  24.     I cannot understand the discussion of this test in Knuth, nor can I
  25.     understand what his variables are.  Is `n' the number under test, or
  26.     is it `q'?  If the former, why the statement "Let n = 1 + (2^k)q" just
  27.     before the statement of Algorithm P.
  28.  
  29. n is the number to test.  The statement "Let  n = 1  + (2^k)q" means that
  30. you  should  find  an integer k and an odd  integer q such that the above
  31. equation holds.  For example, if n = 273, then k = 4 and q = 17.
  32.  
  33. He coninues:
  34.  
  35.     What I would really like is some C or Pascal code which implements
  36.     Algorithm P.  Failing that, a restatement with a more careful
  37.     declaration of variables would be gratefully received.  Sedgwick and
  38.     Numerical Recipies do not deal with primes.  Can anyone help?
  39.  
  40. I can't  provide  C or Pascal code,  but  it  shouldn't be  difficult  to
  41. translate the  following into C or  Pascal. (Note  that you will  need  a
  42. large integer  arithmetic.  If  you don't  have that you  are limited  to
  43. numbers less  than 2^32 and are better of just trial dividing  by all odd
  44. numbers up to 2^16).
  45.  
  46. #############################################################################
  47. ##
  48. #V  Primes[]  . . . . . . . . . . . . . . . . . . . . . . . .  list of primes
  49. ##
  50. ##  'Primes' is a set, i.e., sorted list, of the 168 primes less than 1000.
  51. ##
  52. ##  This is used in 'IsPrimeInt' and 'FactorsInt' to cast  out  small  primes
  53. ##  quickly.
  54. ##
  55. Primes := [   2,  3,  5,  7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53,
  56.      59, 61, 67, 71, 73, 79, 83, 89, 97,101,103,107,109,113,127,131,137,139,
  57.     149,151,157,163,167,173,179,181,191,193,197,199,211,223,227,229,233,239,
  58.     241,251,257,263,269,271,277,281,283,293,307,311,313,317,331,337,347,349,
  59.     353,359,367,373,379,383,389,397,401,409,419,421,431,433,439,443,449,457,
  60.     461,463,467,479,487,491,499,503,509,521,523,541,547,557,563,569,571,577,
  61.     587,593,599,601,607,613,617,619,631,641,643,647,653,659,661,673,677,683,
  62.     691,701,709,719,727,733,739,743,751,757,761,769,773,787,797,809,811,821,
  63.     823,827,829,839,853,857,859,863,877,881,883,887,907,911,919,929,937,941,
  64.     947,953,967,971,977,983,991,997 ];
  65.  
  66.  
  67. #############################################################################
  68. ##
  69. #F  PowerModInt(<r>,<e>,<m>)  . . . . . . power of one integer modulo another
  70. ##
  71. PowerModInt := function ( r, e, m )
  72.     local   pow, f;
  73.  
  74.     # reduce r initially
  75.     r := r mod m;
  76.  
  77.     # handle special case
  78.     if e = 0  then
  79.         return 1;
  80.     fi;
  81.  
  82.     # now use the repeated squaring method (right-to-left)
  83.     pow := 1;
  84.     f := 2 ^ (LogInt( e, 2 ) + 1);
  85.     while 1 < f  do
  86.         pow := (pow * pow) mod m;
  87.         f := QuoInt( f, 2 );
  88.         if f <= e  then
  89.             pow := (pow * r) mod m;
  90.             e := e - f;
  91.         fi;
  92.     od;
  93.  
  94.     # return the power
  95.     return pow;
  96. end;
  97.  
  98.  
  99. #############################################################################
  100. ##
  101. #F  IsPrimeInt( <n> ) . . . . . . . . . . . . . . . . . . .  test for a prime
  102. ##
  103. ##  'IsPrimeInt' returns 'true' if the  integer <n>  is  a prime  and 'false'
  104. ##  otherwise.  By convention 'IsPrimeInt( 0 ) = IsPrimeInt( 1 ) = false' and
  105. ##  we define 'IsPrimeInt( -<n> ) = IsPrimeInt( <n> )'.
  106. ##
  107. IsPrimeInt := function ( n )
  108.     local  p, e, o, x, i, d;
  109.  
  110.     # make $n$ positive and handle trivial cases
  111.     if n < 0         then n := -n;       fi;
  112.     if n in Primes   then return true;   fi;
  113.     if n <= 1000     then return false;  fi;
  114.  
  115.     # do trial divisions by the primes less than 1000
  116.     # faster than anything fancier because $n$ mod <small int> is very fast
  117.     for p  in Primes  do
  118.         if n mod p = 0  then return false;  fi;
  119.         if n < (p+1)^2  then return true;   fi;
  120.     od;
  121.  
  122.     # find $e$ and $o$ odd such that $n-1 = 2^e * o$
  123.     e := 0;  o := n-1;   while o mod 2 = 0  do e := e+1;  o := o/2;  od;
  124.  
  125.     # try bases [2,3,5,7]
  126.     for d  in [2,3,5,7]  do
  127.  
  128.         # look at the seq $2^o, 2^{2 o}, 2^{4 o}, .., 2^{2^e o}=2^{n-1}$
  129.         x := PowerModInt( d, o, n );
  130.         i := 0;
  131.         while i < e  and x <> 1  and x <> n-1  do
  132.             x := x * x mod n;
  133.             i := i + 1;
  134.         od;
  135.  
  136.         # if it is not of the form $.., -1, 1, 1, ..$ then $n$ is composite
  137.         if not (x = n-1 or (i = 0 and x = 1))  then
  138.             return false;
  139.         fi;
  140.  
  141.         # there are no strong pseudo-primes to base 2 smaller than 2047
  142.         if d = 2  and n < 2047  then
  143.             return true;
  144.         fi;
  145.  
  146.     od;
  147.  
  148.     # after those tests $n$ is probably a prime
  149.     return "n maybe a prime";
  150. end;
  151.  
  152. Martin.
  153.  
  154. -- .- .-. - .. -.  .-.. --- ...- . ...  .- -. -. .. -.- .-
  155. Martin Sch"onert,   Martin.Schoenert@Math.RWTH-Aachen.DE,  +49 241 804551
  156. Lehrstuhl D f"ur Mathematik, Templergraben 64, RWTH, D 51 Aachen, Germany
  157.