## Emirps

### November 2, 2010

We give three solutions. The first uses a function `emirp?`

to identify emirps and enumerates them by testing each *n* in range:

`(define (rev n) (undigits (reverse (digits n))))`

(define (palin? n) (= n (rev n)))

(define (emirp? n) (and (prime? n) (prime? (rev n)) (not (palin? n))))

(define (emirps1 n) (filter emirp? (range 2 n)))

There are 36 emirps less than a thousand:

`> (emirps1 1000)`

(13 17 31 37 71 73 79 97 107 113 149 157 167 179 199 311 337

347 359 389 701 709 733 739 743 751 761 769 907 937 941 953

967 971 983 991)

`Emirps1`

operates in linear time. But `prime?`

is fairly slow, and sieving is fast, so it makes sense to identify the emirps using a sieve:

`(define (emirps2 n)`

(let ((n10 (next-power-of-ten n)))

(let loop ((ps (filter (complement palin?) (primes n10))) (es '()))

(cond ((null? ps) (take-while (right-section < n) (sort < es)))

((member (rev (car ps)) (cdr ps))

(loop (cdr ps) (cons (car ps) (cons (rev (car ps)) es))))

(else (loop (cdr ps) es))))))

It is tempting, but wrong, to sieve up to *n* instead of 10^{⌈log10n⌉}. (You can probably guess that I succumbed to the temptation in my first version of the function.) The problem is that primes whose reversal is bigger than *n* will mistakenly disappear from the result. The calculation `next-power-of-ten`

was originally written as `(inexact->exact (expt 10 (ceiling (log10 n))))))`

, but was changed due to a nagging insecurity that somewhere, someday, `(log10 1000)`

might be 3.0000000000000004, leading to an incorrect answer:

`(define (next-power-of-ten n)`

(let loop ((i 1) (ten^i 10))

(if (<= n ten^i) (expt 10 i)

(loop (+ i 1) (* ten^i 10)))))

Testing shows that `emirps2`

is faster than `emirps1`

for *n*=1000:

`> (time (length (emirps1 1000)))`

(time (length (emirps1 1000)))

no collections

14 ms elapsed cpu time

22 ms elapsed real time

76344 bytes allocated

36

> (time (length (emirps2 1000)))

(time (length (emirps2 1000)))

no collections

1 ms elapsed cpu time

1 ms elapsed real time

29032 bytes allocated

36

But consider what happens when we increase *n* from a thousand to a million:

`> (time (length (emirps1 1000000)))`

(time (length (emirps1 1000000)))

111 collections

16760 ms elapsed cpu time, including 289 ms collecting

17245 ms elapsed real time, including 312 ms collecting

465986352 bytes allocated, including 469509920 bytes reclaimed

11184

> (time (length (emirps2 1000000)))

(time (length (emirps2 1000000)))

5 collections

24402 ms elapsed cpu time, including 33 ms collecting

26126 ms elapsed real time, including 34 ms collecting

22172848 bytes allocated, including 19376512 bytes reclaimed

11184

Oops! The problem is that `member`

performs linear search, making `emirps2`

quadratic rather than linear. But the idea of sieving for primes instead of testing for them is good. Here is a third version of `emirps`

:

`(define (emirps3 n)`

(let* ((n10 (next-power-of-ten n))

(ps (filter (complement palin?) (primes n10)))

(qs (sort < (append ps (map rev ps)))))

(let loop ((qs qs) (es '()))

(cond ((or (< n (car qs)) (null? (cdr qs))) (reverse es))

((= (car qs) (cadr qs))

(loop (cddr qs) (cons (car qs) es)))

(else (loop (cdr qs) es))))))

Here `ps`

is the non-palindromic primes, and `qs`

combines those primes with their reverses, in order; for `n`

=1000, the beginning of the `qs`

list is (13 13 14 16 17 17 19 23 …). Then a simple scan through `qs`

, stopping at *n*, identifies twins, which indicate that both a prime and its reverse are in the list. This is a linear algorithm, like `emirps1`

, but much faster, since it uses a sieve rather than testing each number for primality, and since the loop is only over the primes rather than the complete range. Here are some timings, in milliseconds, for `(length (emirps `

*n*`))`

:

`n emirps1 emirps2 emirps3 length`

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

1000 14 0 36

10000 163 15 6 240

50000 1020 724 60 980

150000 2529 23734 664 2954

500000 8427 24415 671 6700

1000000 16760 24402 690 11184

By comparison, `(length (primes 1000000))`

takes 123ms and finds 78498 primes.

We used `take-while`

, `filter`

, `range`

, `sort`

, `digits`

, `undigits`

, `log10`

, `complement`

and `right-section`

from the Standard Prelude, and `primes`

and `prime?`

from two previous exercises. You can run the program at http://programmingpraxis.codepad.org/WCYfLGjL.

Pages: 1 2

[…] Praxis – Emirps By Remco Niemeijer In today’s Programming Praxis, our task is to enumerate all the non-palindrome prime numbers that are still […]

My Haskell solution (see http://bonsaicode.wordpress.com/2010/11/02/programming-praxis-emirps/ for a version with comments):

Here is my Haskell solution: http://www.gleocadie.net/?p=178&lang=en

Comments, advices are welcomed.

I’m loving this website! Thanks a lot for the problems.

I’m loving this website, thanks a lot!

Haskell:

Love the site.

Sorry, line 14 should read

if prime(int(str(n)[::-1])) and str(n) != str(n)[::-1]:

@Joe: your prime() function isn’t quite correct. For instance, prime(10) returns True.

Wrote it in python (though for speed, a different language may be better). In the interest of speed, I wrote a function called reverse(n) that writes n backwards using only numerical operations (avoiding strings). Interestingly, trying to memoize my is_prime() function made everything slower (perhaps due to memory usage?). Running ./emirps.py 1000000 took just under 15 seconds on my aging laptop. For more speed I used pypy (Python with a just in time compiler); it finished in just under 3 seconds.

Changing line 36 in main(n) from

to

yields a decent speed up; normal execution finishes in just under 11 seconds, while pypy finishes it in under 2.

Probably if you avoid using math.log10(n), you’ll end up with a faster execution.

My C Implementation

http://codepad.org/xVdzeVs6

This one was easy one. ;)

My C Implementation

http://codepad.org/xVdzeVs6

This one was easy one. ;)

We can build on the sieve of Erastosthenes exercise to create a list of primes < 1,000,000. Since the output array of the sieve is sorted, we can use binary search on the table for O(log n) search of the sieve, when checking for whether each reversal is prime. The Factor code for sieve is already posted on this blog and won't be reproduced here.

USING: kernel sequences vectors math math.parser locals

binary-search sieve ;

IN: emirp

: emirp? ( n vec — ? )

swap

number>string dup reverse 2dup =

[ 3drop f ]

[ nip string>number swap sorted-member? ] if ;

:: emirp-filter ( primes — semirp )

V{ } clone

primes

[ dup primes emirp?

[ suffix ]

[ drop ] if

] each ;

: Semirp ( n — vec )

primes emirp-filter ;

Factor session:

( scratchpad ) 100000 Semirp

— Data stack:

V{ 13 17 31 37 71 73 79 97 107 113 149 157 167 179 199…

( scratchpad ) length .

1646

Sorry, somehow the code block didn’t work on last post, probably got a tag wrong somewhere…

We can build on the sieve of Erastosthenes exercise to create a list of primes < 1,000,000. Since the output array of the sieve is sorted, we can use binary search on the table for O(log n) search of the sieve, when checking for whether each reversal is prime. The Factor code for sieve is already posted on this blog and won't be reproduced here.

USING: kernel sequences vectors math math.parser locals

binary-search sieve ;

IN: emirp

: emirp? ( n vec — ? )

swap

number>string dup reverse 2dup =

[ 3drop f ]

[ nip string>number swap sorted-member? ] if ;

:: emirp-filter ( primes — semirp )

V{ } clone

primes

[ dup primes emirp?

[ suffix ]

[ drop ] if

] each ;

: Semirp ( n — vec )

primes emirp-filter ;

Factor session:

( scratchpad ) 100000 Semirp

— Data stack:

V{ 13 17 31 37 71 73 79 97 107 113 149 157 167 179 199…

( scratchpad ) length .

1646

ruby solution : http://codepad.org/0ZlE1j4F