## 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

### 17 Responses to “Emirps”

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

```import Data.Numbers.Primes

emirps :: [Integer]
emirps = [p | p <- primes, let rev = reverse (show p)
, isPrime (read rev), show p /= rev]

main :: IO ()
main = print \$ takeWhile (< 1000000) emirps
```

```import Data.List

isPrime l =
isPrimeHelper l primes

isPrimeHelper a (p:ps)
| p*p > a = True
| a `mod` p == 0 = False
| otherwise = isPrimeHelper a ps

primes = 2 : filter isPrime [3,5..]

permirps = drop 4 (takeWhile (<1000000) primes)

isEmirps x =
let sx = show x in
let rev = reverse sx in
sx /= rev &&  isPrimeHelper (read rev) primes

emirps = filter isEmirps permirps

main = print \$ emirps
```

4. fengshaun said

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

```import Data.Numbers.Primes

emirps :: [Integer]
emirps = filter (\a -> isEmirp a) . filter (\a -> not . isPalindrome \$ a) . takeWhile (< 1000000) \$ primes
where
isPalindrome :: Integer -> Bool
isPalindrome x = (reverse . show \$ x) == (show x)

isEmirp :: Integer -> Bool
isEmirp x = isPrime . read . reverse . show \$ x

main :: IO ()
main = print . show \$ emirps
```
5. fengshaun said

I’m loving this website, thanks a lot!

```import Data.Numbers.Primes

emirps :: [Integer]
emirps = filter (\a -> isEmirp a) . filter (\a -> not . isPalindrome \$ a) . takeWhile (< 1000000) \$ primes
where
isPalindrome :: Integer -> Bool
isPalindrome x = (reverse . show \$ x) == (show x)

isEmirp :: Integer -> Bool
isEmirp x = isPrime . read . reverse . show \$ x

main :: IO ()
main = print . show \$ emirps
```
6. Joe Eisenberg said

Love the site.

```from math import sqrt

def prime(number):
if number == 2:
return True
f = round(sqrt(number))
for n in range(3,int(f)+1,2):
if number % n == 0:
return False
return True

for n in range(1000001,11,-2):
if prime(n):
if prime(int(str(n)[::-1])):
print str(n) + ":" + str(n)[::-1]
```
7. Joe Eisenberg said

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

8. Graham said

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

9. Graham said

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.

```#!/usr/bin/env python2.6

import math

def reverse(n):
"""
Given integer n, returns n written backwards.
Note: uses only numerical operations (not string ones) for speed.
"""
d, r = 0, 0
while n > 0:
l = int(math.log10(n))
r += (10 ** d) * (n // (10 ** (l)))
d += 1
n %= (10**l)
return r

def is_prime(n):
"""
Simple check for primality, testing n mod k for odd k up to 1 + sqrt(n).
"""
if n == 2:
return True
elif n % 2 == 0:
return False
else:
for k in xrange(3, 1 + int(math.sqrt(n)), 2):
if n % k == 0:
return False
return True

def main(n):
"""
Prints out all emirps less than n.
"""
for p in xrange(2, n):
r = reverse(p)
if (is_prime(p)) and (r != p) and (is_prime(r)):
print p

if __name__ == '__main__':
import sys
main(int(sys.argv[1]))
```
10. Graham said

Changing line 36 in main(n) from

```for p in xrange(2, n):
```

to

```for p in xrange(13, n, 2):
```

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

11. turuthok said

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

12. My C Implementation
This one was easy one. ;)

13. My C Implementation
This one was easy one. ;)

14. David said

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:

— Data stack:
V{ 13 17 31 37 71 73 79 97 107 113 149 157 167 179 199…
1646

15. David said

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:

— Data stack:
V{ 13 17 31 37 71 73 79 97 107 113 149 157 167 179 199…
1646

16. sealfin said

November 2nd, 2010.c:

```#ifndef LEONARDO
#error "This program requires Leonardo IDE to run."
#endif

#include "seal_bool.h" /* <http://GitHub.com/sealfin/C-and-C-Plus-Plus/blob/master/seal_bool.h> */
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <time.h>

#define N 1000000L

#if N <= 0
#error "N ≤ 0."
#endif

long f_ReverseDigits_A( long p )
{
long result = 0;
while( p != 0 )
{
result *= 10;
result += ( p % 10 );
p /= 10;
}
return result;
}

size_t f_NumberOfDigits( const long p )
{
long i = 10;
size_t number_of_digits = 1;
while( p / i != 0 )
{
i *= 10;
number_of_digits ++;
}
return number_of_digits;
}

char *g_s;

long f_ReverseDigits_B( const long p )
{
static char *s = NULL;
size_t i = 0, k;

if( s == NULL )
{
s = ( char* )malloc( sizeof( char ) * ( f_NumberOfDigits( N ) + 1 ));
g_s = s;
}
sprintf( s, "%ld", p );
k = strlen( s ) - 1;
while( i < k )
{
const char c = s[ i ];
s[ i ] = s[ k ];
s[ k ] = c;
i ++;
k --;
}
return atol( s );
}

bool g_isPrime[ N ];
long g_numberOfPrimes = 0, *g_primes;

bool f_IsPrime( const long p )
{
if( p % 2 == 0 )
return p == 2;
else
return g_isPrime[ p ];
}

void p_EnumerateEmirps( long ( *p_reverseDigits )( long ), const char * const p_fileName, long *p_secondsTaken )
{
const time_t beginning_time = time( NULL );
long i = 0;
FILE *f = fopen( p_fileName, "w" );

for( ; i < g_numberOfPrimes; i ++ )
{
const long prime = g_primes[ i ];
const long reversed_prime = p_reverseDigits( prime );
if( f_IsPrime( reversed_prime ) && ( reversed_prime != prime ))
fprintf( f, "%ld\n", prime );
}
fclose( f );
*p_secondsTaken = time( NULL ) - beginning_time;
}

char *f_AllocateMemoryForAndSetMessage( const long p_secondsTaken, const char p_variant, const char p_terminator )
{
size_t message_length = strlen( "• " );
char *message;
message_length += f_NumberOfDigits( p_secondsTaken );
message_length += strlen( " second" );
if( p_secondsTaken != 1 )
message_length ++;
message_length += strlen( " using the f_ReverseDigits_  function " );
message = ( char* )malloc( sizeof( char ) * ( message_length + 1 ));
sprintf( message, "• %ld second%s using the f_ReverseDigits_%c function%c", p_secondsTaken, ( p_secondsTaken != 1 )?"s":"", p_variant, p_terminator );
return message;
}

void main( void )
{
long i = 0, k, seconds_taken_A, seconds_taken_B;

/* Firstly, let's determine the prime numbers in the range [ 0, N ) so that later we can test if a reversed prime number is still a prime number. */
for( ; i < N; i ++ )
g_isPrime[ i ] = true;
if( N >= 1 )
g_isPrime[ 0 ] = false;
if( N >= 2 )
g_isPrime[ 1 ] = false;
if( N >= 3 )
g_numberOfPrimes ++;
for( i = 3; i < N; i += 2 )
if( g_isPrime[ i ] )
{
g_numberOfPrimes ++;
for( k = i + i; k < N; k += i )
g_isPrime[ k ] = false;
}

/* Secondly, let's create a list of the prime numbers in the range [ 0, N ) so that later we can iterate through that list of prime numbers, testing if each prime number is also an emirp. */
g_primes = ( long* )malloc( sizeof( long ) * g_numberOfPrimes );
k = 0;
if( N >= 3 )
g_primes[ k ++ ] = 2;
for( i = 3; k < g_numberOfPrimes; i += 2 )
if( f_IsPrime( i ))
g_primes[ k ++ ] = i;

p_EnumerateEmirps( f_ReverseDigits_A, "November 2nd, 2010 (A).out", &seconds_taken_A );
p_EnumerateEmirps( f_ReverseDigits_B, "November 2nd, 2010 (B).out", &seconds_taken_B );

free( g_primes );
free( g_s );

{
float height_A, height_B;
char *message = ( char* )malloc( sizeof( char ) * ( strlen( "To enumerate the emirps less than " ) + f_NumberOfDigits( N ) + strlen( " took:" ) + 1 )), *message_A, *message_B;

sprintf( message, "To enumerate the emirps less than %ld took:", N );

if( seconds_taken_A > seconds_taken_B )
{
height_A = 100;
height_B = 100 * (( float )seconds_taken_B / seconds_taken_A );
}
else
{
height_A = 100 * (( float )seconds_taken_A / seconds_taken_B );
height_B = 100;
}

// Set-up the view.
/**
View( Out 0 );
ViewOrigin( Out 100, Out 112, 0 );

SmallText( Out 0, Out 0, Out 11, Out String, Out L, Out Left, 0 )
Assign L = message;
**/

// Set-up the bar chart for the f_ReverseDigits_A function.
message_A = f_AllocateMemoryForAndSetMessage( seconds_taken_A, 'A', ';' );
/**
RectangleFrameColor( 1, Out Cyan, 0 );
Rectangle( Out 1, Out -36, Out Y, Out 30, Out H, 0 )
Assign Y = -1 * height_A H = height_A;

LineColor( 1, Out Cyan, 0 );

Line( Out 1, Out -30, Out Y1, Out 0, Out Y2, 0 )
Assign Y1, Y2 = -1 * height_A - 6;
Line( Out 1, Out 0, Out Y, Out 0, Out -6, 0 )
Assign Y = -1 * height_A - 6;

Line( Out 1, Out -36, Out Y1, Out -30, Out Y2, 0 )
Assign Y1 = -1 * height_A Y2 = -1 * height_A - 6;
Line( Out 1, Out -6, Out Y1, Out 0, Out Y2, 0 )
Assign Y1 = -1 * height_A Y2 = -1 * height_A - 6;
Line( Out 1, Out -6, Out -1, Out 0, Out -6, 0 );

SmallTextColor( 1, Out Cyan, 0 );
SmallText( Out 1, Out 0, Out 22, Out String, Out L, Out Left, 0 )
Assign L = message_A;
**/

// Set-up the bar chart for the f_ReverseDigits_B function.
message_B = f_AllocateMemoryForAndSetMessage( seconds_taken_B, 'B', '.' );
/**
RectangleFrameColor( 2, Out Magenta, 0 );
Rectangle( Out 2, Out 6, Out Y, Out 30, Out H, 0 )
Assign Y = -1 * height_B H = height_B;

LineColor( 2, Out Magenta, 0 );

Line( Out 2, Out 12, Out Y1, Out 42, Out Y2, 0 )
Assign Y1, Y2 = -1 * height_B - 6;
Line( Out 2, Out 42, Out Y, Out 42, Out -6, 0 )
Assign Y = -1 * height_B - 6;

Line( Out 2, Out 6, Out Y1, Out 12, Out Y2, 0 )
Assign Y1 = -1 * height_B Y2 = -1 * height_B - 6;
Line( Out 2, Out 36, Out Y1, Out 42, Out Y2, 0 )
Assign Y1 = -1 * height_B Y2 = -1 * height_B - 6;
Line( Out 2, Out 36, Out -1, Out 42, Out -6, 0 );

SmallTextColor( 2, Out Magenta, 0 );
SmallText( Out 2, Out 0, Out 33, Out String, Out L, Out Left, 0 )
Assign L = message_B;
**/

free( message );
free( message_A );
free( message_B );
}
}```

The solution interpreted using Leonardo IDE 3.4.1 on an Apple Power Mac G4 (AGP Graphics) (450MHz processor, 1GB memory) running Mac OS 9.2.2 (International English).

“November 2nd, 2010 (A).out” & “November 2nd, 2010 (B).out”:

```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
[11,112 lines omitted]
998281
998311
998513
998539
998617
998629
998633
998651
998743
998819
998831
998857
998861
998909
998941
998947
998951
999029
999043
999101
999133
999149
999199
999239
999331
999377
999431
999623
999631
999653
999667
999769
999773
999853
999931
999983```

I interpreted the instruction that I “should strive for maximum speed” as an invitation to compare two functions for reversing the digits of an integer. The results were surprising: the function `f_ReverseDigits_B`, which reverses the digits of an integer by converting that integer into a string, reversing the characters of that string, and converting that string back into an integer, was as fast – and on occasion faster and on no occasion slower – than the function `f_ReverseDigits_A`, which reverses the digits of an integer arithmetically.

(I’m just trying to solve the problems posed by this ‘site whilst I try to get a job; I’m well aware that my solutions are far from the best – but, in my defence, I don’t have any traditional qualifications in computer science :/ )