//
// Original code by Phil Carmody: http://fatphil.org/maths/googolplex/gp+10s.c
// Some comments and quadratic residue check added by Bertram Felgenhauer.
//
// Assume that q | b^e + 1 for a prime q. Then the order of b modulo q
// must be even and divides 2e, so we may assume that order equals 2d
// where d | e. The order must also divide phi(q) = q-1, so that
// q = 2k d + 1 where k = (q-1)/(2d). We turn k into a parameter.
//
// Thus q is a prime factor of b^e + 1 if
//   d | e
//   q = 2k d + 1 is prime
//   b^d + 1 = 0 (mod q)
// Note that if d is odd, then b must be a square modulo q iff k is even.
//
// We take
//   e = 10^100-1
//   b = 10
//
// Optimizations
//   3 | k        ==> 9 | d  -- avoids overlaps with smaller k
//   p | k, p | e ==> p | d  -- ditto
//   k d != 2 (mod 5)        -- otherwise 5 | q
//   k d != 5 (mod 11)       -- otherwise 11 | q
//   k d != 6 (mod 13)       -- otherwise 13 | q
//   k d != 8 (mod 17)       -- otherwise 17 | q
//
// + 10 must be a square modulo q depending on k's parity (d is always odd).
//
#include <stdio.h>
#include <gmp.h>
#include <stdlib.h>

//
// e = 10^100-1 = 3^2
//   11 41 101 251 271 3541 5051 9091 21401 25601 27961 60101 7019801
//   182521213001 14103673319201 78875943472201 1680588011350901
//
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];
//
// 85085 = 5 7 11 13 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 */
            }

            {
                /* check quadraticity of 10 */
                /* (5|2kd+1) (2|2kd+1) (-1)^(kd) = (2kd+1|5) (-1)^(kd(kd-1)/2) */
                static const int qr[5] = {0, 1, 42, 1, 0};
                if (qr[red85085[n3][n]*k%5] ==
                    ((mpz_get_ui(gfactors[n3][n])*(unsigned)k & 3) >> 1))
                    continue;
            }

            /* 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;
}