Project Euler's "Harshad Numbers" Project Euler Problem 387 states:


  A Harshad or Niven number is a number that is divisible by the sum of its digits. 
  201 is a Harshad number because it is divisible by 3 (the sum of its digits.) 
  When we truncate the last digit from 201, we get 20, which is a Harshad number. 
  When we truncate the last digit from 20, we get 2, which is also a Harshad number. 
  Let's call a Harshad number that, while recursively truncating the last digit, always results in a Harshad number a right truncatable Harshad number.
  
  Also: 
  201/3=67 which is prime. 
  Let's call a Harshad number that, when divided by the sum of its digits, results in a prime a strong Harshad number.
  
  Now take the number 2011 which is prime. 
  When we truncate the last digit from it we get 201, a strong Harshad number that is also right truncatable. 
  Let's call such primes strong, right truncatable Harshad primes.


Given a number n, determine the first n strong, right truncatable Harshad primes.

Input: An integer, n, where n&lt;50. (This means your program will need to be somewhat efficient.)

Output: The first n strong, right truncatable Harshad primes, with some form of whitespace in between each value.

Reference Table:
The first 20 strong, right truncatable Harshad primes are:
1-10:181 211 271 277 421 457 631 2011 2017 2099
11-20:2473 2477 4021 4027 4073 4079 4231 4813 4817 6037 8011

Shortest code wins. Good luck!
Python 188 bytes

r=range(10)
h=zip(r,r)
w=input()
p=lambda n:all(n%~i for i in range(1,int(n**.5)))%n
while w:
 n,s=h.pop(1)
 for i in r:
    d=s+i;i+=10*n
    if i%d&lt;1:h+=[(i,d)]
    if p(n/s)and p(i):print i;w-=1


This could be made shorter, but my primary concern was efficiency. This code maintains a list of the 'right truncable harshad numbers', as well as the their digital sum, and then uses this list to extend itself by one digit iteratively, allowing only new truncable numbers. I don't see any way to obtain the efficiency required without maintaining such a list. Runtime for n&lt;=50 is less than one second.

A less efficient version, which calculates everthing as it goes (139 bytes, has an acceptable runtime for n&lt;=31):

s=9
w=input()
while w:
 n=s=s+2;x=2;f=0
 while n&gt;9:n/=10;z=sum(map(int,`n`));x-=n%z;f=f or(n/z&lt;3)+n/z
 if(1&lt;&lt;s)%s==(1&lt;&lt;f)%f==x:print s;w-=1