#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#define uchar unsigned char
#define uint unsigned int
#ifndef _unlikely
#define _unlikely(x) (x)
#endif
#ifndef ifunlike
#define ifunlike(x) if(_unlikely(x))
#endif

int main(int argc, char** argv){
	if(argc<2){
		puts("Program usage: quadratfrei p N [nmin]\n\n"
			"Let a number which is divisible by any of 4, 9, ... p^2 be called Q-number,\n"
			"and a number which is not, be called F-number. The program computes\n"
			"which numbers in [1, N] are Q-numbers and compares the total number\n"
			"of F-numbers in each interval [1, n] with the expected value, defined as\n"
			"n*(1-1/4)*...*(1-1/p^2). Let the difference of those two quantities be called\n"
			"the F-excess at n. The program outputs two files: one listing the successive\n"
			"maxima of the F-excesses and another one listing the succesive maxima of\n"
			"Q-excesses.\n"
			"    If the value of N is not provided the program will employ a default value\n"
			"for it (currently, 900). The maximum allowed value for p is 97.\n"
			"    If the extra argument nmin is provided, the maxima are kept starting\n"
			"from nmin."
			"    The output files include in each row of text `n´, the number of F- or\n"
			"Q-numbers in [1, n], and its the deviation from the expected vaue. The\n"
			"file is closed by a last line in which the values of N and the number of\n"
			"F- or Q-numbers in [1, N] is written");
		return 0;
	}

	static const uint Primes[]={2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,-1U};
	static const uint Squares[]={4,9,25,49,121,169,289,361,529,841,961,1369,1681,1849,2209,2809,3481,3721,4489,5041,5329,6241,6889,7921,9409,-1U};

	char *pstr;
	unsigned int p, n1, N;
	FILE *fmaxQ, *fmaxF;
	char filename[40];
	double r; //1-(1-1/4)*(1-1/9)*...*(1-1/p^2)
	uint Q; //Number of Q-numbers so far found.
	double devM, devm; //max. and min. deviations of M from the expected value.
	uint Nextsq[]={4,9,25,49,121,169,289,361,529,841,961,1369,1681,1849,2209,2809,3481,3721,4489,5041,5329,6241,6889,7921,9409,-1U};
	uint *pfin; //Pointer to within the array nextsq signalling the end of the squares to be checked

	p=atoi(argv[1]);
	pstr=argv[1]; while(*pstr>='0' && *pstr<='9') pstr++; *pstr='\0'; //Leave the string giving p '\0'-terminated,
	pstr=argv[1];														//so that it can be copied to other place.
	if(argc>=3) N=atoi(argv[2]);
	else N=900;
	if(argc>=4) n1=atoi(argv[3]);
	else n1=1;

	if(p<2 || p>97){	puts("Error: value of p out of range"); return 1;}
	if(N<2){puts("Error: value of N too small"); return 1;}
	if(N<n1){puts("Error: The value specified for n1 is larger than N"); return 1;}

	sprintf(filename,"maxdevQ_%s-%u.txt",pstr,n1);
	fmaxQ=fopen(filename,"w");
	if(fmaxQ==NULL){printf("The file %s could not be created for writing. Aborting program.",filename); return 1;}
	//
	sprintf(filename,"maxdevF_%s-%u.txt",pstr,n1);
	fmaxF=fopen(filename,"w");
	if(fmaxF==NULL){printf("The file %s could not be created for writing. Aborting program.",filename); return 1;}

	//Initial values
	p*=p;  //p will not be used anymore
	for(pfin=Nextsq;*pfin<=p;pfin++);
	r=1; for(const uint *pt=Nextsq; pt<pfin; pt++) r*=1.0-1.0/ *pt;
	r=1.0-r;
	//
	devM=-1.0E-12; //We want to avoid negative maxima. We want the first to be listed that in which the number of Q-values is >= the expected value
	devm=0; //The first minima will appear at n= 1,2,3,7... (if p=2, only 1,2 and 3).

	//Up to n1 just increment Q, not keeping track of maximum and minimum
	Q=0;
	uint n;
	for(n=1; n<n1; n++){
		//Check if n equals any of the "next-square" and update the latter where needed.
		int b=0;
		for(uint *pt=Nextsq; pt<pfin; pt++){
			if(*pt!=n) continue;
			b=1;
			*pt+=Squares[(uchar)(pt-Nextsq)];
		}
		Q+=b;
	}

	//r-=0.3828125; //=49/128. See the code using r in the loop for an explanation for this
	for(;n<=N; n++){
		//Check if n equals any of the "next-square" and update the latter where needed.
		int b=0;
		for(uint *pt=Nextsq; pt<pfin; pt++){
			if(*pt!=n) continue;
			b=1;
			*pt+=Squares[(uchar)(pt-Nextsq)];
		}
		/* not worth doing.
		double dev=Q-0.375*n; //No floating point error expected here
		dev-=0.0078125*n;  //Now dev is smaller than Q. We gain a couple of decimal places right
		dev-=r*n; //another one here*/
		//Instead of expe+=r at each iteration of the loop, in order to avoid the accumulation of rounding error.
		double expe=r*n;
		if(b){
			ifunlike(++Q-expe>devM){
				devM=Q-expe;
				fprintf(fmaxQ,"%u %u %.12f\n",n,Q,devM);
			}
		}else{
			ifunlike(Q-expe<devm){
				devm=Q-expe;
				fprintf(fmaxF,"%u %u %.12f\n",n,n-Q,-devm);
			}
		}
	}
	fprintf(fmaxQ,"%u %u\n",N,Q); fclose(fmaxQ);
	fprintf(fmaxF,"%u %u\n",N,N-Q); fclose(fmaxF);

	return 0;
}
