C++ 389 bytes
#include <iostream>
#include <boost/multiprecision/cpp_int.hpp>
#include <boost/multiprecision/miller_rabin.hpp>
using namespace boost::random;typedef boost::multiprecision::cpp_int Z;int main(int,char**v){mt19937 m(clock());independent_bits_engine<mt11213b,256,Z>g(m);Z n{v[1]},p;while(p++<=n)if(miller_rabin_test(p,25,g)&&p.convert_to<std::string>().find( "666" )!=-1)std::cout<<p<<" ";}
This is a full program!
You'll need Boost to compile it. (Or copy and paste into your favorite online C++ shell.)
Run it from the command-line giving n as argument.
Ungolfed:
#include <iostream>
#include <boost/multiprecision/cpp_int.hpp>
#include <boost/multiprecision/miller_rabin.hpp>
using namespace boost::random;
typedef boost::multiprecision::cpp_int integer;
int main( int argc, char** argv )
{
mt19937 mt( clock() );
independent_bits_engine <mt11213b, 256, integer> rng( mt );
integer input {argv[ 1 ]};
integer possible;
while (possible++ <= input)
if (
// is_prime( possible )
miller_rabin_test( possible, 25, rng )
&&
// possible has "666" in it
(possible.convert_to <std::string> ().find( "666" ) != std::string::npos))
std::cout << possible << " ";
}
Shortcuts were made in terms of random number testing. The original code started testing possible primes at 6661 and incremented by two. You'll also get a compiler warning because of that (-1) there instead of npos.
Still, this runs pretty quickly. It only took about 40 seconds to find all 214 satan primes under 1,000,000 on my old AMD Sempron 130.
:^D
output the nth satan primechallenge... \$\endgroup\$