Topic: -={ Challenge }=- is Prime

[Lee Mac]

The Challenge:

To construct a function to determine whether a natural number (positive integer) is prime.

Pour construire une fonction pour dιterminer si un nombre naturel (entier positif) est premier.

Для построения функции определить, является ли натуральное число (целое положительное число) является простым.


Definition of a prime number:

A number only evenly divisible by itself and 1

Notes:

0 and 1 are not prime


To kick off:

Code: [Select]

(defun LM:Primep ( n )
  (and (/= n 1)
    (or (= n 2)
      (not
        (vl-some
          (function
            (lambda ( d ) (zerop (rem n d)))
          )
          (
            (lambda ( i r / l )
              (while (<= (setq i (+ 2 i)) r) (setq l (cons i l)))
              (cons 2 l)
            )
            1 (sqrt n)
          )
        )
      )
    )
  )
)

Enjoy!

[VovKa]

no Ukrainian translation, can't participate, sorry

[roy_043]

No new mathematical insights just a slightly different approach. On a number like 100000001 (can be divided by 17) my code is faster than Lee's. But on big numbers that are primes, like 1508138953, performance is similar.
Of course on even numbers my code wins hands down.

Code: [Select]

(defun kg:PrimeP (n / maxI i)
  (and
    (/= n 1)
    (or
      (= n 2)
      (and
        (/= (rem n 2) 0) ; check for even number
        (/= (setq maxI (sqrt n)) (fix maxI)) ; rare occurrence but you never know...
        (progn
          (setq i 1)
          (while (and (/= (rem n (setq i (+ i 2))) 0) (<= i maxI)))
          (/= (rem n i) 0)
        )
      )
    )
  )
)

[Lee Mac]

no Ukrainian translation, can't participate, sorry

Probably best that you don't - we'd have no chance 
 

[Lee Mac]

No new mathematical insights just a slightly different approach. On a number like 100000001 (can be divided by 17) my code is faster than Lee's. But on big numbers that are primes, like 1508138953, performance is similar.
Of course on even numbers my code wins hands down.

Nice coding Roy 

[Lee Mac]

Using theory that all integers can be represented in the form 6k+i for integer k, i=−1, 0, 1, 2, 3, 4 with 2 dividing (6k + 0), (6k + 2), (6k + 4) and 3 dividing (6k + 3)

Code: [Select]

(defun LM:Primep ( n / m u )
  (and (/= n 1)
    (or (= n 2) (= n 7)
      (not
        (or
          (zerop (rem n 2))
          (zerop (rem n 3))
          (zerop (rem (setq u (sqrt n)) 1))
          (progn
            (setq i 0.0)
            (while
              (and (<= (setq m (+ (* 6.0 (setq i (1+ i))) 1.0)) u)
                   (/= 0 (rem n m))
                   (/= 0 (rem n (- m 2.0)))
              )
            )
            (or (zerop (rem n m))
                (zerop (rem n (- m 2.0)))
            )
          )
        )
      )
    )
  )
)


[gile]

Hi,

Here's something I wrote sometimes ago, trying to make a function which 'learns'.
Rather than restart calculations from the first prime numbers, all 'knowed' prime numbers are stored in a list within the function.

Here's a 'demonstrative version'

Code: [Select]

(defun-q
  prime
  (num / lst den loop)
  (setq lst '(2 3))
  (if (and (= (type num) 'INT) (< 0 num))
    (if (or (member num lst) (= 1 num))
      (alert (strcat "I already know: " (itoa num) " is prime."))
      (if (< num (setq den (last lst)))
        (alert (strcat "I already know: " (itoa num) " is not prime."))
        (if (vl-some '(lambda (x) (zerop (rem num x))) lst)
          (alert (strcat "I calculated: " (itoa num) " is not prime."))
          (progn
            (setq loop   T
                  den (+ den 2)
            )
            (while (and loop (<= den num))
              (if (vl-some '(lambda (x) (zerop (rem den x))) lst)
                 (setq den (+ den 2))
                 (progn
                   (setq lst (append lst (list den)))
                   (alert (strcat "I learned a new prime number: "
                                  (itoa den)
                          )
                   )
                   (if (zerop (rem num den))
                     (if (= num den)
                       (alert (strcat "Now, I know "
                                      (itoa num)
                                      " is prime."
                              )
                       )
                       (progn
                         (alert (strcat "I calculated: "
                                        (itoa num)
                                        " is not prime."
                                )
                         )
                         (setq loop nil)
                       )
                     )
                     (setq den (+ den 2))
                   )
                 )
              )
            )
          )
        )
      )
    )
    (alert "Needs a strictly positive number")
  )
  (setq prime
         (cons (car prime)
               (cons (list 'setq 'lst (list 'quote lst))
                     (cddr prime)
               )
         )
  )
  (princ)
)
The same one with a boolean return.
(eval (cadr isPrime)) returns the list of the 'already knowed' prime numbers

Code: [Select]

(defun-q
  isPrime
  (num / lst den loop result)
  (setq lst '(2 3))
  (if (and (= (type num) 'INT) (< 0 num))
    (if (or (member num lst) (= 1 num))
      (setq result T)
      (if (<= (setq den (last lst)) num)
        (if (not (vl-some '(lambda (x) (zerop (rem num x))) lst))
          (progn
            (setq loop   T
                  den (+ den 2)
            )
            (while (and loop (<= den num))
              (if (vl-some '(lambda (x) (zerop (rem den x))) lst)
                 (setq den (+ den 2))
                 (progn
                   (setq lst (append lst (list den)))
                   (if (zerop (rem num den))
                     (if (= num den)
                       (setq result T)
                       (setq loop nil)
                     )
                     (setq den (+ den 2))
                   )
                 )
              )
            )
          )
        )
      )
    )
  )
  (setq isPrime
         (cons (car isPrime)
               (cons (list 'setq 'lst (list 'quote lst))
                     (cddr isPrime)
               )
         )
  )
  result
)

[VovKa]

Code: [Select]

(kg:PrimeP 3)
(LM:Primep 3)
(LM:Primep 5)


[qjchen]

Code: [Select]

;;;by qjchen@gmail.com
(defun is-prime(n)
 (defun find-2 (n i)
  (cond ((> (* i i) n) n)
        ((= (rem n i) 0) i)
        (T (find-2 n (1+ i)))
  )
 )
 (cond ((and (= (type n) 'INT) (> n 1)) (= (find-2 n 2) n))
       (T nil))
)

recently I also join a challenge here
http://bbs.mjtd.com/dispbbs.asp?boardid=37&Id=83019

there is also a is-prime question.

[Lee Mac]

Code: [Select]

(kg:PrimeP 3)
(LM:Primep 3)
(LM:Primep 5)


Ahhh!    How did i miss that 

[Lee Mac]

Same theory as above, but with modulo 12

Code: [Select]

(defun LM:Primep ( n / u m )
  (and (/= n 1)
    (or (member n '(2 3 5 7 11))
      (not
        (or
          (zerop (rem n 2))
          (zerop (rem n 3))
          (zerop (rem (setq u (sqrt n)) 1))
          (progn
            (setq i 0.0)
            (while
              (and
                (<= (setq m (+ (* 12.0 (setq i (1+ i))) 7.0)) u)
                (/= (rem n m) 0)
                (/= (rem n (- m 2)) 0)
                (/= (rem n (- m 6)) 0)
                (/= (rem n (- m 8)) 0)
              )
            )
            (vl-some (function (lambda ( x ) (zerop (rem n (- m x))))) '(0 2 6 8))
          )
        )
      )
    )
  )
)

[Lee Mac]

Hi,

Here's something I wrote sometimes ago, trying to make a function which 'learns'.
Rather than restart calculations from the first prime numbers, all 'knowed' prime numbers are stored in a list within the function.

Nice solution Gile - that is some great lateral thinking 

[roy_043]

Code: [Select]

(kg:PrimeP 3)
(LM:Primep 3)
(LM:Primep 5)


Code: [Select]

(kg:PrimeP 3) => nil
(LM:Primep 3) => T (first version of LM:Primep)
(LM:Primep 5) => T (first version of LM:Primep)
Well spotted VovKa!

[roy_043]

Code: [Select]

;;;by qjchen@gmail.com
(defun is-prime(n)
....


It's always nice to see (and try to understand ) a recursive function. But for very big numbers recursion is a problem:

Code: [Select]

(is-prime 1508138953) => error : LISP recursion stack limit exceeded

[LE3]

To construct a function to determine whether a natural number (positive integer) is prime.

Pour construire une fonction pour dιterminer si un nombre naturel (entier positif) est premier.

Для построения функции определить, является ли натуральное число (целое положительное число) является простым.

Y Ahora? porque nada mas para los que hablen o entiendan esos idiomas....  

Buena onda eh!  

De todas maneras aqui esta una opcion:

Code: [Select]

static bool isPrime(int num)
{
if (num<=0 || num==1) return false;
for (int p=2; p<=num/2; p++) { if (num%p==0) return false; }
return true;
}

[Lee Mac]

Luis, my C++ is muy rusty, but would this work?

Code: [Select]

#include <math.h>

static bool isPrime(int num)
{
if (num<=0 || num==1 || num%2==0) return false;
         if (num==3) return true;
for (int p=1; p<=sqrt(num); p+=2) { if (num%p==0) return false; }
return true;
}


[LE3]

On a first look and test - nope, it returns different values.

Also on this line:
for (int p=1; p<=sqrt(num); p+=2) { if (num%p==0) return false; }

Needs to be something like:
for (int p=1; p<=(int)sqrt((double)num); p+=2) { if (num%p==0) return false; }

Luis, my C++ is muy rusty, but would this work?

Code: [Select]

#include <math.h>

static bool isPrime(int num)
{
if (num<=0 || num==1 || num%2==0) return false;
         if (num==3) return true;
for (int p=1; p<=sqrt(num); p+=2) { if (num%p==0) return false; }
return true;
}


[Chuck Gabriel]

http://www.theswamp.org/index.php?topic=26363.msg318146#msg318146

[David Bethel]

Not real fancy, but should work:

Code: [Select]

(defun is_prime (i / r n)
  (setq n 3
        r (cond ((/= (type i) 'INT) T)
                ((<= i 1)           T)
                ((= i 2)          nil)
                ((zerop (rem i 2))  T)))
  (while (and (not r)
              (<= n (sqrt i)))
         (if (zerop (rem i n))
             (setq r T)
             (setq n (+ n 2))))
  (not r))

-David

Added 'INT and 2 test

[roy_043]

Code: [Select]

(is_prime 2)

[ElpanovEvgeniy]

my version:

Code: [Select]

(defun EE:Primep (n)
 (defun f (n i)
  (cond ((< i 2))
        ((= (rem n i) 0) nil)
        ((f n (1- i)))
  )
 )
 (and (> n 1) (f n (fix (sqrt n))))
)

[ElpanovEvgeniy]

version iteration

Code: [Select]

(defun EE:Primep_i (n / i)
 (setq i (fix (sqrt n)))
 (while (and (> i 1) (< 0 (rem n i)) (setq i (1- i))))
 (and (< i 2) (> n 1))
)

[Lee Mac]

The recursive function is most probably impractical as the stack limit will be hit quite easily for high primes, so the iterative version is preferable. But most are great examples (as usual )

Lee

[Lee Mac]

Perhaps a way to approach the problem is:

For testing a prime 'p'

** Choose a number 'x' with a large number of distinct prime factors

** Represent all integers modulo x

Hence for x=12 (2² x 3) we have 12k+i for integers k and i=0,1,2,3,4,5,6,7,8,9,10,11

** Check division of 'p' by prime factors of 'x'

Hence for x=12:

prime factor 2 we can remove i=0,2,4,6,8,10
prime factor 3 we can remove i=3,6,9

Hence we are left to check all integers of the form 12k+i <= sqrt(p)  for i=1,5,7,11

Severely reducing the number of division tests we have to make.

[Lee Mac]

Perhaps an example of such a function:

Code: [Select]

(defun LM:isPrime ( p / i k x )
  ;; © Lee Mac 2010
  (setq i (fix (sqrt p)) k 0)

  (if (< 1 p)
    (cond
      ( (vl-position p '(2 3 5 7 11)) T )
      ( (vl-some (function (lambda ( d ) (zerop (rem p d)))) '(2 3 5 7 11)) nil )
      ( (progn
          (while (and (<= (setq x (+ (* 12 (setq k (1+ k))) 11)) i)
                      (not (vl-some (function (lambda ( d ) (zerop (rem p d)))) (list x (- x 4) (- x 6) (- x 10))))))

          (> x i)
        )
      )
    )
  )
)


Code: [Select]

_$ (Benchmark '((LM:isPrime 123457) (EE:Primep_I 123457)))
Benchmarking ...............Elapsed milliseconds / relative speed for 4096 iteration(s):

    (LM:ISPRIME 123457)......1077 / 1.36 <fastest>
    (EE:PRIMEP_I 123457).....1467 / 1.00 <slowest>


[VovKa]

Code: [Select]

(LM:isPrime 1,2,25,35,49....)

[Lee Mac]

Code: [Select]

(LM:isPrime 1,2,25,35,49....)

Tweaked. 

[highflyingbird]

Code: [Select]

(defun HFB:isPrime (n / HD LG LH VA ret)
  (if (= (rem n 2) 0)
    nil
    (progn
      (setq hd 8
    lg 6
    lh 8
    va (1- n)
      )
      (setq ret t)
      (while (>= va hd)
(if (= (rem (- va hd) lg) 0)
  (setq ret nil
va hd
  )
)
(setq lh (+ lh 8)
      lg (+ lg 4)
      hd (+ hd lh)
)
      )
      ret
    )
  )
)
found some bug's in Lee's code.

Code: [Select]

(defun c:test(/ I J RANGE T0)

  (setq range (getreal "\nPlease enter a range(<100000000):"))
  (setq range (1- (fix range)))
  ;;Highflybird's
  (setq t0 (getvar "TDUSRTIMER"))
  (setq i 2)
  (setq j 1)    ;because 2 is a prime
  (repeat range
    (if (HFB:isprime i)
      (progn
        (setq j (1+ j))
        ;;(princ (strcat "\n" (itoa i)))
      )
    )
    (setq i (1+ i))
  )
  (princ "\nIt takes")
  (princ (* (- (getvar "TDUSRTIMER") t0) 86400))
  (princ "Second.")
  (princ (strcat "\nThe count of primes is:" (itoa j)))

  ;;Lee Mac's
  (setq t0 (getvar "TDUSRTIMER"))
  (setq i 2)
  (setq j 0)
  (repeat range
    (if (LM:isPrime i)
      (progn
        (setq j (1+ j))
        ;;(princ (strcat "\n" (itoa i)))
      )
    )
    (setq i (1+ i))
  )

  (princ "\nIt takes")
  (princ (* (- (getvar "TDUSRTIMER") t0) 86400))
  (princ "Second")
  (princ (strcat "\nThe count of primes is:" (itoa j)))
 
  (princ)
)

[highflyingbird]

for example:

Lee's  :
1000--------185        primes
10000------1284        primes

but actually:

1000--------168        primes
10000------1229        primes   

[David Bethel]

I came up with the same number of primes:

Command: iptest

Please enter a range(<100000000):1000

It takes0.016Second
The count of primes is:168

Command:
Command:
IPTEST
Please enter a range(<100000000):10000

It takes0.14Second
The count of primes is:1229

Command:
Command:
IPTEST
Please enter a range(<100000000):100000

It takes2.75Second
The count of primes is:9592

Command:

[Lee Mac]

Looks like I'm going to have to take another look at my code 

[Lee Mac]

Code: [Select]

(HFB:isPrime 2)


Code: [Select]

(HFB:isPrime 1)


[Lee Mac]

Hopefully better 

Code: [Select]

(defun LM:isPrime ( p / i r f )
  ;; © Lee Mac 2010
  (if (or (= 2 p) (= 3 p)) T
    (progn
      (setq i (fix (sqrt p))
            r (or (> 2 p) (zerop (rem p 2)) (zerop (rem p 3))) f -1)

      (while (and (not r) (<= (setq f (+ f 6)) i))
        (if (or (zerop (rem p f)) (zerop (rem p (+ f 2))))
          (setq r T)
        )
      )

      (not r)
    )
  )
)
Liked your loop control David

[Lee Mac]

Provided there aren't bugs (highly possible!)

Please enter a range(<100000000):1000

---- HFB ----
It takes 0.007 Seconds
The count of primes is: 168

---- DB ----
It takes 0.01 Seconds
The count of primes is: 168

---- LM -----
It takes 0.006 Seconds
The count of primes is: 168

Please enter a range(<100000000):10000

---- HFB ----
It takes 0.281 Seconds
The count of primes is: 1229

---- DB ----
It takes 0.34 Seconds
The count of primes is: 1229

---- LM -----
It takes 0.21 Seconds
The count of primes is: 1229

Please enter a range(<100000000):100000

---- HFB ----
It takes 6.783 Seconds
The count of primes is: 9592

---- DB ----
It takes 6.639 Seconds
The count of primes is: 9592

---- LM -----
It takes 3.436 Seconds
The count of primes is: 9592

[David Bethel]

Code: [Select]


(defun is_prime (i / r n s)
  (setq n 3
        s (sqrt i)
        r (cond ((/= (type i) 'INT) T)
                ((<= i 1)           T)
                ((= i 2)          nil)
                ((zerop (rem i 2))  T)))
  (while (and (not r)
              (<= n s))
         (if (zerop (rem i n))
             (setq r T)
             (setq n (+ n 2))))
  (not r))

Surprised me how much faster it is when you don't use the (sqrt) function each iteneration  -David

[Lee Mac]

Certainly improves the time David, the less calculated within the loop the better 

Please enter a range(<100000000):1000

---- HFB ----
It takes 0.007 Seconds
The count of primes is: 168

---- DB ----
It takes 0.016 Seconds
The count of primes is: 168

---- LM -----
It takes 0.013 Seconds
The count of primes is: 168

Please enter a range(<100000000):10000

---- HFB ----
It takes 0.184 Seconds
The count of primes is: 1229

---- DB ----
It takes 0.228 Seconds
The count of primes is: 1229

---- LM -----
It takes 0.147 Seconds
The count of primes is: 1229

Please enter a range(<100000000):100000

---- HFB ----
It takes 6.322 Seconds
The count of primes is: 9592

---- DB ----
It takes 5.37 Seconds
The count of primes is: 9592

---- LM -----
It takes 3.316 Seconds
The count of primes is: 9592

Please enter a range(<100000000):200000

---- HFB ----
It takes 18.25 Seconds
The count of primes is: 17984

---- DB ----
It takes 15.028 Seconds
The count of primes is: 17984

---- LM -----
It takes 8.587 Seconds
The count of primes is: 17984

[pkohut]

A dollar short and a day late. Was in Python so...

Code: [Select]

def PrintIsPrime(x):
"""Determines if x is prime and prints result"""
s1 = "Is prime!"
s2 = "Is not prime!"

if x < 2 or (x > 2 and x % 2 == 0):
s = s2
elif x == 2:
s = s1
else:
i = 3
while x % i != 0:
i += 2
s = s1 if i == x else s2
print(x, s)

def IsPrime(x):
"""Returns True/False result of x == prime number"""
if x < 2 or (x > 2 and x % 2 == 0):
return False
elif x == 2:
return True
else:
i = 3
while x % i != 0:
i += 2
return True if i == x else False


>>> PrintIsPrime(3)
3 Is prime!
>>> PrintIsPrime(101)
101 Is prime!
>>> PrintIsPrime(100000001)
100000001 Is not prime!
>>> IsPrime(1508138953)
1508138953 Is prime!  <<<<-------That one took a minute or so to calc.

>>> IsPrime(3)
True
>>> IsPrime(3)
True
>>> IsPrime(21)
False
>>> IsPrime(31)
True

[highflyingbird]

How about" Miller-Rabin primality test"?
when the number < 3.4X10^14 ,the result of test is correct.

Code: [Select]

static unsigned g_aPrimeList[] = { // 100素数表
        2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41,
43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97
};
unsigned Montgomery(INT64 n,unsigned p,unsigned m)
{ //快速计算(n^e)%m的值
unsigned k=1;
n%=m;
while(p!=1)
{
if(0!=(p&1))
k=(k*n)%m;
n=(n*n)%m;
p>>=1;
}
return(n*k)%m;
}

bool RabbinMillerTest( unsigned n )
{
if (n<2)
{ // 小于2的数即不是合数也不是素数
throw 0;
}

const unsigned nPrimeListSize=sizeof(g_aPrimeList)/sizeof(unsigned);//求素数表元素个数
for(int i=0;i<nPrimeListSize;++i)
{// 按照素数表中的数对当前素数进行判断
if (n/2+1<=g_aPrimeList[i])
{// 如果已经小于当前素数表的数，则一定是素数
return true;
}
if (0==n%g_aPrimeList[i])
{// 余数为0则说明一定不是素数
return false;
}
}
// 找到r和m，使得n = 2^r * m + 1;
int r = 0, m = n - 1; // ( n - 1 ) 一定是合数
while ( 0 == ( m & 1 ) )
{
m >>= 1; // 右移一位
r++; // 统计右移的次数
}
const unsigned nTestCnt = 7; // 表示进行测试的次数
for ( unsigned i = 0; i < nTestCnt; ++i )
{ // 利用随机数进行测试，
int a = g_aPrimeList[ i ];
if ( 1 != Montgomery( a, m, n ) )
{
int j = 0;
int e = 1;
for ( ; j < r; ++j )
{
if ( n - 1 == Montgomery( a, m * e, n ) )
{
break;
}
e <<= 1;
}
if (j == r)
{
return false;
}
}
}
return true;
}

RabbinMillerTest (2234223321);
return 0;

because this way used "_int64",so it's not  applicable for LISP . but for a small number ,it's ok.

Code: [Select]


(defun paw_mod( bs power diver)
  (cond
    ((zerop power)
     1
    )
    ((= power 1)
     bs
    )
    ((= (logand power 1) 0)
     (paw_mod (rem (* bs bs) diver) (lsh power -1) diver)
    )
    (t
      (rem (* (paw_mod (rem (* bs bs) diver) (/ power 2) diver) bs) diver)
    )
  )
)

(defun RabbinMiller (base num / D K N RET TT)
  (setq n (1- num))
  (setq d n)
  (while (= (logand d 1) 0)
    (setq d (lsh d -1))
  )
  (setq k (paw_mod base d num))
  (if (or (= k 1) (= k n))
    t
    (progn
      (setq tt (/ n 2))
      (while (/= d tt)
(setq d (lsh d 1))
(if (= (paw_mod base d num) n)
  (setq ret t
d tt
  )
)
ret
      )
    )
  )
)

(defun check(p)
  (if (> p 13)
    (cond
      ( (or (= (rem p 2) 0)
    (= (rem p 3) 0)
    (= (rem p 5) 0)
    (= (rem p 7) 0)
    (= (rem p 11) 0)
    (= (rem p 13) 0)
)
        nil
      )
      ( (and (RabbinMiller 2 p)
     (RabbinMiller 7 p)
     (RabbinMiller 5 p)
)
        T
      )
      (t nil)
    )
    (and (member p '(2 3 5 7 11 13)))
  )
)


(check 23421);
nil;