#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>       /* clock_t, clock, CLOCKS_PER_SEC */
#include "gmp.h"
//#define LIMIT 4000000000U
#define LIMIT 10000000U
#define PRIMEN 40
const int mask[8]={
0x01,
0x02,
0x04,
0x08,
0x10,
0x20,
0x40,
0x80
};
//typedef unsigned int uint128_t __attribute__((mode(TI)));
#define GETNOTPOSSIBLE(x)   (notpossible[(x)>>3]& mask[(x)&7])
#define SETNOTPOSSIBLE(x)   (notpossible[(x)>>3]|=mask[(x)&7])
#define SETPOSSIBLE(x)      (notpossible[(x)>>3]&=~(mask[(x)&7]))z
//#define mulmod(a,b,m)  (unsigned long long int)(((uint128_t)(a)*(uint128_t)(b)) % ((uint128_t)(m)))
int main(void)
{
    printf("Calculating\n");
    mpz_t x;
    mpz_t legp;
    mpz_init(x);
    mpz_init(legp);
    mpz_t mpztemp;
    mpz_init(mpztemp);
    int i=0,j=0,k=0;
    long long majorcounter=1,minorcounter=0;
    unsigned long long int temp;
    mpz_set_ui(x, LIMIT);
    unsigned long int testprimes[PRIMEN]={0};
    unsigned long int partialresidue[PRIMEN];
    for(i=0;i<PRIMEN;++i){
        mpz_nextprime(x,x);
        partialresidue[i]=1;
        testprimes[j++]=mpz_get_ui(x);
    }
    const unsigned int STEP=LIMIT/1000;
    char* notpossible=(char*)calloc(STEP/8+1,1);
    clock_t t = clock();
    while(majorcounter<LIMIT){
        printf("Going from %llu ",majorcounter);
        printf("to %llu, ",majorcounter+STEP);
        printf("%d-1000\n",++k);
        for(j=0;j<PRIMEN;++j){
            mpz_set_ui(legp,testprimes[j]);
            for(minorcounter=0,temp=partialresidue[j],mpz_set_ui(mpztemp,partialresidue[j]);minorcounter<STEP;++minorcounter){//calcular el residuo parcial para el numero primo actual
                i=minorcounter+majorcounter;//actual numero. Nota: Ya que majorcounter empieza en 1, esto nunca sera cero, por lo tanto no cancelara los datos de la variable temp
                //temp*=i;
                //temp%=testprimes[j];
                mpz_mul_ui(mpztemp,mpztemp,i);
                mpz_mod_ui(mpztemp,mpztemp,testprimes[j]);
                //mulmod(temp,i,testprimes[j]);
                if(GETNOTPOSSIBLE(minorcounter)) continue;//Ya ha sido marcado como un candidato imposible, saltarlo
                mpz_set_ui(x,mpz_get_ui(mpztemp)+1);
                int legendre=mpz_legendre(x,legp);//verificar si n!+1 (mod m) es un residuo cuadratico modulo m
                if(legendre==-1) SETNOTPOSSIBLE(minorcounter);//si no lo es, definitivamente no es un cuadrado, eliminarlo como posiblidad
            }
            //partialresidue[j]=temp;
            partialresidue[j]=mpz_get_ui(mpztemp);
        }
        for(i=0;i<STEP;++i){//verificar la verdad de la conjetura
            if(GETNOTPOSSIBLE(i)==0)
                printf("%d\n",majorcounter+i);
        }
        memset(notpossible,0,STEP/8);//Vaciar el array
        majorcounter+=STEP;
    }
    /* free used memory */
    printf("Clicks taken: %d\n",clock()-t);
    mpz_clear(x);
    mpz_clear(legp);
    free(notpossible);
    printf("Finishing");
    return EXIT_SUCCESS;
}