/* 
 * gemelos2a.php
 *
 * 09-10-2016 JAZ
 *
 * Estudio de los primos gemelos
 * 
 * Nueva versión para reducir el tiempo y el consumo de memoria.  SOLO SE PUEDE
 * UTILIZAR PARA VALORES DE pi<=2147483647 (si es PHP de 32 bits o 
 * pi<=9223372036854775807 si es PHP de 64 bits).
 * 
 * Para valores mayores de pi hay que utilizar la versión anterior.
 *
 * ==========================================================================
 */

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

unsigned int numParesGemelosIntervalo(mpz_t n1, mpz_t n2);
long double calcF2Inicial(mpz_t p);

int main(int argc, char **argv)
{
	unsigned int num= 100;	/* Número de primos consecutivos a calcular */
	unsigned int pInicial= 3;	/* primo mayor o igual a 3 */

	char s1[1024];
	char s2[1024];
	char s3[1024];

	mpz_t	p,t;
	long double	f2;
	unsigned int p0;
	mpz_t	n1,n2;
	unsigned int	numPares;
	int	flag;

	if(2==argc)
		num = atoi(argv[1]);
	else if(3==argc)
	{
		num = atoi(argv[1]);
		pInicial = atoi(argv[2]);
	}

/*
Para cada primo pi calcula el número de pares de gemelos en [pi^2, pi^2+2*pi]
 */

	fprintf(stdout, "pi;f1(pi);f4(pi);\n");
	mpz_init(p);
	mpz_init(t);
	mpz_init(n1);
	mpz_init(n2);

	mpz_set_ui(p,pInicial);

	mpz_sub_ui(t,p,1);
	mpz_nextprime(p,t);	// Garantizar que empieza en primo
	
	f2 = calcF2Inicial(p);
	while( 0 < num-- )
	{
		mpz_mul(n1,p,p);
		mpz_mul_ui(n2,p,2);
		mpz_add(n2,n1,n2);
		
		numPares = numParesGemelosIntervalo(n1,n2);

		fprintf(stdout, "%s;%d;%.13Lf;\n",
			mpz_get_str(s1, 10, p),
			numPares,
			(long double) 0.79 * f2 );
			
		p0 = mpz_get_ui(p);
		mpz_nextprime(p,p);
		f2 *= (long double)(mpz_get_ui(p)-2)/(long double) p0;
	}
}


long double calcF2Inicial(mpz_t p)
{
	mpz_t pAct;
	unsigned int p0;
	long double f2= 1;

	mpz_init(pAct);
	mpz_set_ui(pAct,3);

	while(0<mpz_cmp(p,pAct))
	{
		p0 = mpz_get_ui(pAct); 
		mpz_nextprime(pAct,pAct);
		f2 *= (long double)(mpz_get_ui(pAct)-2)/(long double) p0;
	}

	return f2;
}
	
/*==========================================================================*/

unsigned int numParesGemelosIntervalo(mpz_t n1, mpz_t n2)
{
	/* Devuelve el número de primos gemelos en el intervalo ($n1,$n2) */
	mpz_t p3,p4,p5;
	unsigned int num=0;
	mpz_init (p3);
	mpz_init (p4);
	mpz_init (p5);

	mpz_nextprime(p3,n1);
	while(0<mpz_cmp(n2,p3))
	{
		mpz_nextprime(p4,p3);
		mpz_add_ui(p5, p3, 2);
		if(0==mpz_cmp(p4,p5))
			num++;

		mpz_set(p3, p4);
	}
	
	return num;
}