#include <stdio.h>
#include <gmp.h>
#include <stdlib.h>

static int const ifactors[17]=
{
    11,41,101,251,271,3541,5051,9091,
    21401,25601,27961,60101,7019801,
    0,0,0,0 /* too large */
};
static char const *factors[17]=
{
    "11","41","101","251","271","3541","5051","9091",
    "21401","25601","27961","60101","7019801",
    "182521213001","14103673319201","78875943472201","1680588011350901"
};
/* remember there are 2 threes too */
static mpz_t gfactors[3][1<<17];
static unsigned int red85085[3][1<<17];

/* Could sort them? Yeah - why not! */
static int indices[3<<17];


int derefcmp(const void* pa, const void* pb)
{
    int i1=*(const int*)pa;
    int i2=*(const int*)pb;
    return mpz_cmp(gfactors[i1&3][i1>>2], gfactors[i2&3][i2>>2]);
}

int main(int argc, char**argv)
{
    int threes;
    int k;
    
    mpz_t gq, gt; /* temporary for q, and the power therefrom */

    /* d=divisors(10^100-1); */
    for(threes=0; threes<=2;++threes)
    {
	static const int threepow[3]={1,3,9};
	int family;
	mpz_init_set_ui(gfactors[threes][0], threepow[threes]);
	indices[threes<<17]=threes;

	for(family=0; family<17; ++family)
	{
	    int familybase=1<<family;
	    int offset;

	    mpz_init_set_str(gfactors[threes][familybase],
			     factors[family],
			     10);
	    if(threes)
	    {
		mpz_mul_ui(gfactors[threes][familybase],
			   gfactors[threes][familybase],
			   threepow[threes]);
	    }

	    indices[(threes<<17)+familybase]=threes+(familybase<<2);

	    for(offset=1; offset<(1<<family); ++offset)
	    {
		mpz_init(gfactors[threes][familybase+offset]);
		mpz_mul(gfactors[threes][familybase+offset],
			gfactors[threes][offset], /* all smaller factors */
			gfactors[0][familybase]); /* only the largest factor */

		indices[(threes<<17)+familybase+offset]=threes+((familybase+offset)<<2);
	    }
	}
    }

    /* Sort the indices */
    qsort(indices, 3<<17, sizeof(int), derefcmp);

    /* create reduced versions */
    for(threes=0; threes<=2; ++threes)
    {
	for(k=0; k<(1<<17); ++k)
	{
	    red85085[threes][k] = mpz_tdiv_ui(gfactors[threes][k], 85085);
	}
    }

    mpz_init(gq);
    mpz_init(gt);
    /* for(k=1,3000,print("Doing k=",k); */
    for(k=argv[1]?atoi(argv[1]):1857; k<(argc>2?atoi(argv[2]):48000); ++k)
    {
	int must3=k%3?0:2; /* if there's a 3, demand 2 threes from d */
	int mustother=0;
	int i;
	for(i=0; i<17; ++i)
	{
	    if(!ifactors[i] || ifactors[i]>k) break;
	    if(!(k%ifactors[i])) mustother|=(1<<i);
	}
	printf("Doing k=%i, requires f=%d:%05x\n",k, must3, mustother);
        fflush(stdout);

	/* for(n=1,length(d),p=d[n];q=2*p*k+1; */
	for(i=0; i<(3<<17); ++i)
	{
	    int index=indices[i];
	    int n3=index&3, n=index>>2;

	    /* throw away ones done at a lower k */
	    if((n3<must3) || (mustother & ~n))
	    {
		continue;
	    }
	    
	    /* get rid of easy comosites using the reduced form */
	    if(  (red85085[n3][n]*(k%5) % 5 == 2) 
	       ||(red85085[n3][n]*(k%7) % 7 == 3) 
	       ||(red85085[n3][n]*(k%11) % 11 == 5) 
	       ||(red85085[n3][n]*(k%13) % 13 == 6) 
               ||(red85085[n3][n]*(k%17) % 17 == 8))
	    {
		/* There's a chance that the first few values are divisible by small primes as they _are_ the small primes */
		if(k>2 || index>4)
		{
		    continue;
		}
		/* don't risk skipping the really small ones */
	    }

	    /* p=gfactors[n3][n] */
	    mpz_mul_ui(gq, gfactors[n3][n], 2*k);
	    mpz_add_ui(gq, gq, 1);

	    if(mpz_probab_prime_p(gq, 1))
	    {
		/* if(lift(Mod(10,q)^p+1)==0&&ispseudoprime(q),print(q)))) */
		mpz_set_ui(gt, 10);
		mpz_powm(gt, gt, gfactors[n3][n], gq);
		mpz_add_ui(gt, gt, 1);
	    
		if(mpz_cmp(gt, gq)==0) /* power+1 == q */
		{
		    mpz_out_str(stdout, 10, gq);
		    printf("\n");
		}
	    }
	}
    }
    return 0;
}