/* 
  Copyright (C) 2023 Alessandro Languasco
*/ 
 

{npq_alpha(boundprime, alpha = 2) = local(millisec, minutes, seconds, 
npqalpha, ok, iterations, name, namefile, seqfile, logp, logq, a, pn, qn, 
log_gratio, u, v, index, u1,v1, u2,v2, x, uleft, vleft, string, mat,
ustart, vstart, coeff, golden_ratio, loglogalpha, time, starttime, elaptime,
logoneoveralpha, logratio, uend, vend
);

starttime = getwalltime(); 
  
print("---------");
if (alpha <= 1, error("alpha must be greater than 1"));

print("Looking for n_{p,q} (alpha); alpha = ", alpha );
print("with the continued fractions algorithm; decreasing order" );
print("---------");


name="npq_data";
namefile = Strprintf("%s-%s-(%1.3f).csv", name, boundprime, alpha); 
seqfile = fileopen(namefile, "w");	
string = "p;q;alpha;u_start;v_start;u_end;v_end;iterations;npqalpha <= "; 
filewrite(seqfile, string);
		
golden_ratio = (1+sqrt(5))/2;		
log_gratio = log(golden_ratio);
loglogalpha = log(log(alpha));
logoneoveralpha = - log(alpha);

forprime(p=2, boundprime, 
	logp = log(p);
	
	forprime(q=p+1,boundprime, 	
		
		\\ starting point
			    
		iterations = 0; 
		logq = log(q);
		
		\\  generate all the convergents;
		a = contfrac(logp/logq);
		mat = contfracpnqn(a, #a);  
		\\ mat = matrix of all the convergents
		
		pn = mat[1,]; \\ numerators of the convergents
		qn = mat[2,]; \\ denominators of the convergents

		\\L = #qn; \\ number of convergents
		\\print("number of the convergents = ", L);
		
		\\ coeff is M in the paper
		coeff = ceil( (log(logq) - loglogalpha )/ log_gratio );
		
		\\ pari/gp stores q_0 in position 1; so there's a shift of one position in storing
		\\ the convergents; s_M is the M-th denominator of the convergent p_M/q_M
		\\ and it is stored in position M+1 in the qn array
		 
		ustart = qn[coeff + 1] + qn[coeff + 2];
		vstart = 0;
		
		\\  n_start is p^ustart * q^0
		u = ustart;
		v = vstart;			 
	
		npqalpha = 0; \\ initialization

		iterations += 1;	\\ updating the iterations counter
		
		ok = 0 ; \\ when ok = 1 we got npq(alpha)
		
		while ( ok == 0, \\ perform until I have a success
		
			\\ GOING IN THE DECREASING ORDER IN THE n_i SEQUENCE	
			
			x = u * logp - v * logq ;
			
			if ( x  < 0 ,
						
				\\ SEARCHING FOR AN UPPER APPROXIMANT 
				\\ 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 index of the largest numerator < v
					index = setsearch(pn, v , 1 ) - 1 );
				
					\\ upper approximants; p_n/q_n ; n odd; stored with an even index by pari/gp
				if ( index <= 1,  
							print("index - error; upper approximant");
							next(1)); \\error; exiting
				if ( index % 2  == 0 ,
							\\ index is even; it stores an odd convergent; good
							u1 = qn[index]; v1 = pn[index],
							\\ else
							\\ index is odd; it stores an even convergent
							\\ so the largest odd convergent is in position (index - 1)								
							u1 = qn[index - 1]; v1 = pn[index - 1]; 
					);
					
				\\exponents for the left neighbor
				uleft = u + u1; 
				vleft = v - v1; 
				\\ ratio  = leftelement/element  
				logratio = u1 * logp - v1 * logq;
				\\ left neighbour
				\\ leftelement = p^uleft * q^vleft;	 	 
				,
 
 				\\ ELSE (so x>0)
 			
 				\\ SEARCHING FOR AN LOWER APPROXIMANT 
				\\ 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 index of the largest denominator < u
					index = setsearch(qn, u , 1 ) - 1 );
 
				\\ lower approximants; p_n/q_n ; n even; stored with an odd index by pari/gp					
				if ( index <= 0, 
							print("index - error; lower approximant");
							next(1)); \\error; exiting 
				if ( index % 2  == 1 ,
							\\ index is odd; it stores an even convergent; good
							u2 = qn[index]; v2 = pn[index],
							\\ index is even; it stores an odd convergent;
							\\ so the largest even convergent is in position (index - 1)
							u2 = qn[index - 1]; v2 = pn[index - 1];
					);		
				
				\\exponents for the left neighbor
				uleft = u - u2; 
				vleft = v + v2;
				\\ ratio = leftelement/element  
				logratio  = v2 * logq - u2 * logp;				
				\\ left neighbour
				\\ leftelement = p^uleft * q^vleft;	  
	
			);
						
			if (uleft < 0 || vleft < 0, print("break - uv");  
							next(1));
			
			\\ updating the iterations counter
			iterations +=1 ;	
			\\print("iterations = " ,iterations);
			 
			
			if ( logratio <= logoneoveralpha , 
								\\ in this case npqalpha is approximated by p^u * q^v 
								\\ npqalpha = p^u * q^v; 
								uend = u;
								vend = v;
								ok = 1,
							\\else
								\\ we can generate another left element
								\\ setting data for the next iteration
								\\ ok = 0;
								u = uleft;
								v = vleft; 
								);		
	
			);	
		print("---------");
  
		print("p = ", p, "; q = ", q,"; n = p^u*q^v; u,v >= 0");
		print("Starting with [p,q]-exponents: ustart = ", ustart, ", vstart = ", vstart);		
		\\print("number of convergents = ", L);		
		print("performed iterations = ", iterations);	
		filewrite1(seqfile, p);
		filewrite1(seqfile,";");
		filewrite1(seqfile, q);
		filewrite1(seqfile,";");
		filewrite1(seqfile, alpha);
		filewrite1(seqfile,";");
		filewrite1(seqfile, ustart);
		filewrite1(seqfile,";");
		filewrite1(seqfile, vstart);
		filewrite1(seqfile,";");
		filewrite1(seqfile, uend);
		filewrite1(seqfile,";");
		filewrite1(seqfile, vend);
		filewrite1(seqfile,";");
		filewrite1(seqfile, iterations);
		filewrite1(seqfile,";");

		
		npqalpha = p^uend * q^vend;
		print ("n_{", p,",", q,"}(", alpha, ") = ", npqalpha);
		print ("with [p,q]-exponents: u = " , uend,"; v = ", vend);
		filewrite(seqfile, npqalpha);
		fileflush(seqfile)								
		
	);
); 
print("---------");
fileclose(seqfile); 

time = getwalltime() - starttime; 

elaptime = time;
seconds = floor(elaptime /1000)%60;
minutes = floor(elaptime /60000);
millisec = elaptime - minutes*60000 - seconds*1000;
print("npqalpha computation time  (output time included): ", minutes, " min, ", seconds, " sec, ", millisec, " millisec");
 

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


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

 gp2c-run -pmy_ -g -W  npq_alpha_contfrac.gp
Reading GPRC: /Users/languasc/.gprc
GPRC Done.

                  GP/PARI CALCULATOR Version 2.15.4 (released)
         arm64 running darwin (aarch64/GMP-6.2.1 kernel) 64-bit version
    compiled: Jul 10 2023, Apple clang version 14.0.3 (clang-1403.0.22.14.1)
                           threading engine: pthread
                 (readline v8.1 enabled, extended help enabled)

                     Copyright (C) 2000-2022 The PARI Group

PARI/GP is free software, covered by the GNU General Public License, and comes 
WITHOUT ANY WARRANTY WHATSOEVER.

Type ? for help, \q to quit.
Type ?18 for how to get moral (and possibly technical) support.

parisizemax = 2048000000, primelimit = 500000, nbthreads = 8

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