[Maxima-bugs] [ maxima-Bugs-1050304 ] primep() sometimes takes for ever

[SourceForge.net <no...@so...>]

Bugs item #1050304, was opened at 2004-10-19 20:41
Message generated for change (Comment added) made by hg-graebe
You can respond by visiting: 
https://sourceforge.net/tracker/?func=detail&atid=104933&aid=1050304&group_id=4933

Category: None
Group: None
Status: Open
Resolution: None
Priority: 5
Submitted By: Frank Thieme (lefloyd)
Assigned to: Nobody/Anonymous (nobody)
Summary: primep() sometimes takes for ever

Initial Comment:
we had to test 3^n-2 for being prime from n=1 up to
100. In MuPAD the results just need seconds, but in
maxima-5.9.1 (with clisp or cmucl) it takes really long
for some n. For example n=37.

----------------------------------------------------------------------

Comment By: HG Graebe (hg-graebe)
Date: 2005-01-10 09:41

Message:
Logged In: YES 
user_id=1193894

Frank is right, and I'm wondering about that. Here is as a
first hack
an implementation of the Rabin-Miller-Test.  May be it is
better to
have it in the lisp part, but I did not (yet) realize how
and where.

div(a,b):=first(divide(a,b));

modpower(a,n,m):=
  if n<0 then error("modpower only for n>0")
  else if n=0 then 1
  else if n=1 then a
  else if oddp(n) then mod(a*modpower(mod(a*a,m),div(n,2),m),m)
  else mod(modpower(mod(a*a,m),div(n,2),m),m);

RabinMillerTest(m):=block([q,t,j,isComposite],
  if m<10 then return(primep(m)),
  /* avoids some trivialities, a good prime test makes anyway 
     a smallPrimesTest for a list of precompiled small
primes */ 
  q:m-1, t:0, isComposite:FALSE,
  while evenp(q) do (q:q/2, t:t+1),
  /* print(q,t, m-1=2^t*q), */
  thru 20 do ( 
    block([a,b],
      /* print("start a try"), */
      a:random(m), b:modpower(a,q,m),
      /* print(a,b), */
      if b=1 or b=-1 then return(0),
      for j:1 thru (t-1) do (
        b:mod(b*b,m),
        /* print("loop with j=",j,", b=",b), */
        if b=1 then (isComposite:TRUE,return(0)),
        if b=-1 then j:t
      ),
      /* print("loop finished with b=",b,isComposite), */
      if isComposite then return(0),
      if b#-1 then (isComposite:TRUE,return(0))
      /* , print(a,"no information in the loop") */
    ),
    if isComposite then return(0)
  ),
  not(isComposite));

Programming first time a Maxima piece was a hard task. May
be I did not
yet understand it well but my opinion was that Maxima has an
ugly old
fashioned syntax, uncommon organization of control flow etc.
Also
I missed a description of the elementary arithmetic
functions in the
Reference Manual and the DESCRIBE part (even a remark how to
include a
comment), so I defined my own. 

And now test it:
          
RabinMillerTest(79);
       
CarmichaelZahlen:[1105,1729,294409, 56052361, 118901521,
172947529,
216821881, 228842209, 1299963601];

map(RabinMillerTest, CarmichaelZahlen);

SmallTestExamples:makelist(3^n-2,n,1,10);
map(RabinMillerTest, SmallTestExamples);
map(primep, SmallTestExamples);

BigTestExamples:makelist(3^n-2,n,40,50);
map(RabinMillerTest, BigTestExamples);
map(primep, BigTestExamples);

Compare with MuPAD's result: 

BigTestExamples:=[3^n-2 $ n=40..50];
map(BigTestExamples, isprime);

[FALSE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE,
FALSE, FALSE, FALSE]


----------------------------------------------------------------------

You can respond by visiting: 
https://sourceforge.net/tracker/?func=detail&atid=104933&aid=1050304&group_id=4933