Bugs item #1050304, was opened at 20041019 22:41
Message generated for change (Comment added) made by andrejv
You can respond by visiting:
https://sourceforge.net/tracker/?func=detail&atid=104933&aid=1050304&group_id=4933
Please note that this message will contain a full copy of the comment thread,
including the initial issue submission, for this request,
not just the latest update.
Category: None
Group: None
>Status: Closed
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^n2 for being prime from n=1 up to
100. In MuPAD the results just need seconds, but in
maxima5.9.1 (with clisp or cmucl) it takes really long
for some n. For example n=37.

>Comment By: Andrej Vodopivec (andrejv)
Date: 20060321 04:39
Message:
Logged In: YES
user_id=1179910
The new version of primep uses MillerRabin test and has no
problems with these computations.
Andrej

Comment By: HG Graebe (hggraebe)
Date: 20050110 10: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 RabinMillerTest. 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:m1, t:0, isComposite:FALSE,
while evenp(q) do (q:q/2, t:t+1),
/* print(q,t, m1=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 (t1) 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^n2,n,1,10);
map(RabinMillerTest, SmallTestExamples);
map(primep, SmallTestExamples);
BigTestExamples:makelist(3^n2,n,40,50);
map(RabinMillerTest, BigTestExamples);
map(primep, BigTestExamples);
Compare with MuPAD's result:
BigTestExamples:=[3^n2 $ 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
