prime numbers - fix



A few general observations and suggestions:

1.  Add checks for invalid input;  possibly each argument to $strongprime 
must be
an positive integer? 

As a general project, I  suggest that we build better tools for invalid 
input.  The 
function require-set in  nset might  be a starting point. Until we have 
standard tools, write your own.

2.  Once you are sure that n and b are both positive integers, you can use
Common Lisp arithmetic operators on them  -- for example, replace  (cexpt 
b s) with (expt b s).

3.  Write r-test style testing code.  If your testing code is 
test-prime.mac, you can run
your testing code through the undocumented  feature of load:

(c1)  load("test-prime.mac", 'test);

Be sure to test your code for all kinds of crazy inputs -- CRE 
expressions, bigfloats, values on the 
boundary, etc.

4.  Is code fast if nobody uses it?    Write user-documenation in texi 
format; if it isn't documented,
few will ever use it.

Barton




Andrei Zorine <zoav1@uic.nnov.ru>
Sent by: maxima-admin@math.utexas.edu
12/15/2003 05:29 AM

 
        To:     Maxima <Maxima@www.ma.utexas.edu>
        cc: 
        Subject:        [Maxima] prime numbers - fix


The last line of $strongprime should be nil) and not '$FALSE)). This way 
is correctly identifies Mersenne numbers at least, and does fast. More 
benchmarks are needed (but I believe now it's stable enough).
--
A.Z.

> (DEFMFUN $STRONGPRIME(N B)
>   (LET ((S (SUB1 N)) (NM1 (SUB1 N)) (R 0) X)
>        (DO () ((MODDP S) S)
>                   (SETQ R (+ 1 R) S (DIV S 2)))
>        (LET ((MODULUS N))
>                    (SETQ X (CEXPT B S)))
>        (IF (OR (= X 1) (= X -1)) 
>                   (RETURN-FROM $STRONGPRIME T))
>        (SETQ R (SUB1 R))
>        (DO ((J 0 (+ 1 J)))
>                   ((> J R) NIL)
>                   (SETQ X (MOD (* X X) N))
>                   (COND
>                    ((LIKE X NM1) (RETURN-FROM $STRONGPRIME T))
>                    ((LIKE X 1) (RETURN-FROM $STRONGPRIME NIL))))
>        '$FALSE))
          ^^^^^^^^^
Should be: NIL))

_______________________________________________
Maxima mailing list
Maxima@www.math.utexas.edu
http://www.math.utexas.edu/mailman/listinfo/maxima