/* 
  Copyright (C) 2023 Alessandro Languasco
*/ 


{genseq_fixedpq(boundprime, boundexp = 5, N1 = 1 ) = local(num, denom , value, m , 
M, element, millisec, minutes, seconds, C1, C2, logratio, length_trivial,
setseq, L, rightelement, delta, aux, indexmin, indexmax, numprimes, ppow, qpow,
seqlength, name, namefile, seqfile, logp, logq, a, pn, qn, MINVALUE, MAXVALUE,
u,v,index, u1,v1, u2,v2, x, uright, vright, pmin, qmin, pmax, qmax, prec,
accuracy, eps, umin, vmin, umax, vmax, string, mat, boundseq, jj, n, Y, N, Sums, avg,
ustart, vstart, m_trivial, M_trivial, m_contfrac, M_contfrac, logtwo,
starttimeconstcomp_trivial, elaptimeconstcomp_trivial, totcases,
starttimeconstcomp_contfrac, elaptimeconstcomp_contfrac, time,
starttimeconstcompcomp_global, elaptimeconstcomp_global);

accuracy = 15;
prec = 19;
default(realprecision,prec);   
eps = 1.*10^(-accuracy);

elaptimeconstcomp_trivial = 0;
elaptimeconstcomp_contfrac = 0;

starttimeconstcompcomp_global = getwalltime(); 

\\default(parisizemax,"32G");

print("---------");		   
\\print("Trivial part: generating and sorting the sequence n = p^u*q^v; n >= N; u,v >= 0");
print("N = max ( N1,  exp((log p) * (log q); N1 = ", N1 );
print("Accuracy = ", accuracy, " decimal digits" );
print("---------");

\\ bound for the exponents to be used in the trivial part and as starting point
\\boundexp = 10; 

\\ for the continued fractions part

numprimes = 0;
forprime(p=2, boundprime,
	forprime(q=p+1,boundprime,
	numprimes +=1
	));

MAXVALUE = 10^1000000;
MINVALUE = 0;

boundseq = ceil((boundexp+1)*(boundexp+1)*numprimes)+1; 

name="sequence";
namefile = Strprintf("%s-%s-%s.csv", name, boundprime, boundexp); 
seqfile = fileopen(namefile, "w");	
string = "p;q;C1pq;C2pq;mupq"; 
filewrite(seqfile, string);

logtwo = log(2);		
Y = [0];		
forprime(p=2, boundprime, 
	logp = log(p);
	
	forprime(q=p+1,boundprime, 
		ppow = 1;
		qpow = 1;
		n = ppow * qpow;

		starttimeconstcomp_trivial = getwalltime(); 			
		setseq = vector(boundseq);
		jj = 1;
		N = max(N1, exp(log(p)*log(q)) );
		for(u=0,boundexp,  
			for(v=0,boundexp, 
					if (n >= N,
						setseq[jj] = n;
						jj += 1);
				qpow *= q; 		
				n = ppow * qpow;
				);
			qpow = 1;
			ppow *= p; 	
			n = ppow;
			);
		\\ sort the S-units	
		setseq = Set(setseq);
		\\ get rid of the extra zero element
		setseq = setdelta(setseq, Y);				
		\\print(setseq);
		\\ number of S-units generated
		L = #setseq;
		
		M = MINVALUE;
		m = MAXVALUE;
		
		Sums = 0;
		totcases = 0;
		for(n = 1, L-1, 
			element = setseq[n];
			rightelement = setseq[n+1]; 
			delta = rightelement - element;
			if( delta < element, 
			
				totcases += 1;
				aux = log(element);
				num = aux - log( delta );
				denom = log( aux ); 
				value = num / denom ;
				Sums += value;

				if (value < m, m = value; indexmin = n );
				if (value > M, M = value; indexmax = n );
				);
			);		
		
		time = getwalltime() - starttimeconstcomp_trivial;
		elaptimeconstcomp_trivial += time;
		
		starttimeconstcomp_contfrac = getwalltime(); 
			
		m_trivial = m ;
		M_trivial = M ;
		\\print("length of the trivial part = ", totcases);
		length_trivial = totcases;

		setseq = []; \\ cleaning memory
		ustart = boundexp;
		vstart = boundexp;
    
		M = MINVALUE;
		m = MAXVALUE;

		seqlength = 0; 
		logq = log(q);

		a = contfrac(logp/logq);
		\\ generate all the convergents		
		mat = contfracpnqn(a,#a); \\ matrix of all the convergents
		pn = mat[1,]; \\ numerators of the convergents
		qn = mat[2,]; \\ denominators of the convergents
		L = #qn; \\ number of convergents
		
		u = ustart;
		v = vstart;			
		\\rightelement = p^u * q^v;

		seqlength +=1 ;	
		for ( jj = 1, L,
			\\element = rightelement;
			\\ looking for a denominator <= u
			index = setsearch(qn, u);
			\\ if u is not in the list of denominators, setsearch (with two parameters) returns 0
			if (index == 0, 
					\\ u was not in the list of denominators; setsearch (with three parameters) 
					\\ returns the position in which u should be inserted in such a list
					\\ so this value minus 1 is the largest denominator < u
					index = setsearch(qn, u , 1 ) - 1 );
				
				\\ upper approximants; p_n/q_n ; n odd; stored with an even index by pari/gp
			if ( index <= 1, 						
							\\filewrite1(seqfile, p);
							\\filewrite1(seqfile,";");
							\\filewrite1(seqfile, q);
							\\filewrite1(seqfile,";");
							\\filewrite(seqfile, "no qn <= u;no qn <= u"); 
							\\fileflush(seqfile);
							next(1));
			if ( index % 2  == 0 ,
							\\ index is even; it stores an odd convergent
							v1 = qn[index]; u1 = pn[index],
							\\ index is odd; it stores an even convergent
							\\ so the largest odd convergent is in position (index - 1)								
							v1 = qn[index - 1]; u1 = pn[index - 1];
					);


			\\ looking for a numerator <= v							
			index = setsearch(pn, v);
			\\ if v is not in the list of numerators, setsearch (with two parameters) returns 0				
			if (index == 0, 
					\\ v was not in the list of numerators; setsearch (with three parameters) 
					\\ returns the position in which v should be inserted in such a list
					\\ so this value minus 1 is the largest numerator < v
					index = setsearch(pn, v , 1 ) - 1 );
 
			\\ lower approximants; p_n/q_n ; n even; stored with an odd index by pari/gp					
			if ( index <= 0, 
							\\filewrite1(seqfile,";");
							\\filewrite1(seqfile, q);
							\\filewrite1(seqfile,";");
							\\filewrite(seqfile, "no pn <= v;no pn <= v"); 
							\\fileflush(seqfile);
							next(1));
			if ( index % 2  == 1 ,
							\\ index is odd; it stores an even convergent
							u2 = pn[index]; v2 = qn[index],
							\\ index is even; it stores an odd convergent;
							\\ so the largest even convergent is in position (index - 1)
							u2 = pn[index - 1]; v2 = qn[index - 1];
					);		
														
			x = abs (v1*logp - u1*logq) - abs(v2*logp - u2*logq);
			if ( x  < 0 ,
					uright = u - v1; vright = v + u1;
					logratio = u1 * logq - v1 * logp;	
					\\ratio = q^u1/p^v1;
					\\rightelement  = element * ratio
					\\else
					,
					uright = u + v2; vright = v - u2;	
					logratio = v2 * logp - u2 * logq;	
					\\ratio = p^v2/q^u2;
					\\rightelement  = element * ratio;
					);
						
			if (uright < 0 || vright < 0, print("break - uv");  
							next(1));
				
			\\ updating counter
			seqlength +=1 ;	
			\\print("seqlength = " ,seqlength);
			\\ delta = rightelement - element;		
			\\ delta = element * (ratio - 1 );

								
			\\ if( delta >= element, print("delta > element")); 
			\\ if( delta < element , 
			if( logratio < logtwo , 
				\\aux = log(element);
				aux = u * logp + v * logq ;				
				\\ delta = element * (ratio -1); log( delta ) = log(element) + log(ratio-1) = aux + log(ratio-1) 
				\\num = aux - log( delta );
				if (logratio < eps,
					num = - log(logratio);
					, \\else
					num =  - log( exp(logratio) - 1 );
					);
				denom = log( aux ); 
				value = num / denom ;
				Sums += value;
				
				if (value < m, m = value; indexmin = seqlength; pmin = p; qmin = q; umin=u; vmin =v );
				if (value > M, M = value; indexmax = seqlength; pmax = p; qmax = q; umax=u; vmax =v );
				);
				
			\\ right neighbor exponents
			u = uright ;
			v = vright ;				

			);	
		m_contfrac = m ;
		M_contfrac = M ;	
		time = getwalltime()-starttimeconstcomp_contfrac;
		elaptimeconstcomp_contfrac += time;
		print("---------");
		
		print("C1(p,q) - C2(p,q) - mu((p,q)");

		totcases = seqlength + length_trivial;
		avg = Sums/totcases;

		print("p = ", p, " q = ", q," n = p^u*q^v; n >= N; u,v >= 0");
		print("N = max ( N1,  exp((log p) * (log q); N1 = ", N1 );
		print("Trivial part a_0, b_0 = ",  boundexp);
		print("length of the trivial part = ", length_trivial);
		
		print("length of the cont.frac. part = ", seqlength);	
		print("Number of convergents = ", L);
		print("Generating up to ", L, " right neighbors in the sequence");
		
		print ("the generated n_i sequence has ", totcases , " elements"); 
		print ("p runs from 2 to ", boundprime);
		print ("q runs from p+1 to ", boundprime);
		C1 = min(m_trivial, m_contfrac);
		C2 = max(M_trivial, M_contfrac);
		filewrite1(seqfile, p);
		filewrite1(seqfile,";");
		filewrite1(seqfile, q);
		filewrite1(seqfile,";");
		
		if( C1 < MAXVALUE,
			print ("C1(", p,",", q,") = ", C1 );
			print ("C2(", p,",", q,") = ", C2 );
			print ("Sum(", p,",", q,") = ", Sums);
			print ("mu(", p,",", q,") = ", avg);	

			filewrite1(seqfile, C1);
			filewrite1(seqfile,";");	
			filewrite1(seqfile, C2);
			filewrite1(seqfile,";");
			filewrite(seqfile, avg),
			\\else
			filewrite(seqfile, ";; ");
			);
										
	);
	fileflush(seqfile);
); 
fileflush(seqfile);
print("---------");
fileclose(seqfile); 

time = getwalltime() - starttimeconstcompcomp_global; 

elaptimeconstcomp_global = time;
seconds=floor(elaptimeconstcomp_global /1000)%60;
minutes=floor(elaptimeconstcomp_global /60000);
millisec=elaptimeconstcomp_global - minutes*60000 - seconds*1000;
print("Constants computation time [GLOBAL TIME] (output time included): ", minutes, " min, ", seconds, " sec, ", millisec, " millisec");

seconds=floor(elaptimeconstcomp_trivial /1000)%60;
minutes=floor(elaptimeconstcomp_trivial /60000);
millisec=elaptimeconstcomp_trivial - minutes*60000 - seconds*1000;
print("Constants computation time [TRIVIAL PART TIME] (output time included): ", minutes, " min, ", seconds, " sec, ", millisec, " millisec");

seconds=floor(elaptimeconstcomp_contfrac /1000)%60;
minutes=floor(elaptimeconstcomp_contfrac /60000);
millisec=elaptimeconstcomp_contfrac - minutes*60000 - seconds*1000;
print("Constants computation time [CONTFRAC PART TIME] (output time included): ", minutes, " min, ", seconds, " sec, ", millisec, " millisec");


print("******  END PROGRAM   ********");   
}


/******************

gp2c-run -pmy_ -g -W seq_combined_fixed_primes.gp >results-1000.txt <istruzione.txt

******************/