/*
Copyright (C) 2021 Alessandro Languasco
*/

/*
Theorems 2 and 4 numerical verifications;
finding the q_0(r) values in Remark 17 using Lemma 8
*/

\\ function D
{D(q)= 
     3.125 * min( 2*Pi , 0.5*log(q) )
};

\\ function W(y)
{W(y)= local(R);
       R = 9.645908801;
       sqrt( log(y) / R )
};

\\ function y(q)
{y(q)=local(R);
       R = 9.645908801;
	   exp( 1.44 * R * (log(q))^2 )
};

\\ numerator of the first quotient
{num1(y,D,q)= local(c1);
		 c1=1.015;
		 c1 * (log(q))^2 * y^( - D/ (sqrt(q) * (log(q))^2 ) )
};

\\ denominator of the first quotient
{den1(D,q) = 
			D*sqrt(q) };		

\\  first quotient			
{f1(q)= 
      num1(y(q),D(q),q) / den1(D(q),q) };

\\ numerator of the second quotient      
{num2(W) =local(R);
        R = 9.645908801; 
        (8/9) * ( 2*R*W^2 + (4*R+1)*W + 4*R ) };     

\\ denominator of the second quotient      
{den2(W) = 
		exp(W) };       

\\  second quotient			
{f2(q)= 
      num2(W(y(q))) / den2(W(y(q))) };


         
\\  NEW general estimates for S(m,q)      
{newboundS(C,m,q)= local(alpha, alpha1, res, logtwo, h, M); 
		res=0;
		logtwo = log(2);
		h=(q-1)/m;
	    \\M=factor(h); if (M[1,1] == 2, v = M[1,2]);  e = v-1; 
		\\e = floor(log(h) / log(2)) - 1; \\ upper bound for nu_2(h)-1 
		\\print(m);
		\\alpha =  C * ( log(q) +  log (m));
		alpha =  C *  log(q) ;
		alpha = max(alpha,3);
		alpha = min(alpha, q-1);
		alpha1 = min(log(alpha)/(2*logtwo) , (log(q-1)/logtwo -1) /2  ); \\ might use e/2 as second term, but it oscillates
		res = alpha1 * log(q-1) / ((q-1)^(1/m) -1 ) 
        			+ logtwo/2 * q^2 / ( 2^(alpha/2) - 1 )
        	  		+ (8*log(m *alpha^2)) / (7*q^(1/m)) 
        	  		+ 1.053/(q-(2+0.0001));
        if (m > 1, res = res + log(q-1)/ ((q-1)^(2/m) -1));
        return(res)
     };
             
{checknewLemma8(C,m,q)= local(alpha, alpha1, res, logtwo, h, M); 
		res=0;
		logtwo = log(2);
		h=(q-1)/m;
	    M=factor(h); if (M[1,1] == 2, v = M[1,2]);  e = v-1; 
		e = floor(log(h) / logtwo) - 1; \\ upper bound for nu_2(h)-1 
		\\print(e);
		\\alpha =  C * ( log(q) +  log (m));
		alpha =  C *  log(q) ;
		alpha = max(alpha,3);
		alpha = min(alpha, q-1);	
		\\print(alpha);	
		alpha1 = min(log(alpha)/(2*logtwo) , e /2  ); \\ 
		\\print(alpha1)
		res = alpha1 * log(q-1) / ((q-1)^(1/m) -1 ) 
        			+ logtwo/2 * q^2 / ( 2^(alpha/2) - 1 )
        	  		+ (8*log(m/2*alpha^2)) / (7*q^(1/m)) 
        	  		+ 1.053/(q-(2+0.0001));
        if (m > 1, res = res + log(q-1)/ ((q-1)^(2/m) -1));
        return(res)
     }; 
      

\\ NEW general case formula; right hand side
{newF(C,m,q) = 
        Euler  - m * ( log(y(q)) / (q-1) + f1(q) + f2(q) ) - newboundS(C,m,q)
        };

        
\\ using the bounds on gammaK2 and gammaK1 and the general estimates for S(m,q)

\\{boundgammas(q) = ( 2*1.4262639307  - 0.3145570795 ) * log(q)
\\ };
 
 {boundgammas(q) = ( 2*1.4263  - 0.3145 ) * log(q)
 };
 

\\ general case formula; with gammaK2 and gammaK1 bounds;  right hand side
{newG1(C,q) = 
		Euler - log(q) / (q-1)^2 - boundgammas(q)/(q-1) - newboundS(C,1,q)
		}; 

 
		
		
        
/************
RESULTS
--------
ChecknewLemma8

? q=3;r=1;C=10
%19 = 10
? checknewLemma8(C,r,q)
%20 = 4.4363236918355480679422592649845056847
? q=5;r=1;C=10
%21 = 10
? checknewLemma8(C,r,q)
%22 = 3.9454749367223008192684229223923909408
? q=5;r=2;C=10
%23 = 10
? checknewLemma8(C,r,q)
%24 = 5.1182966403141144874073341447318486223

--------
// general case estimates but r=1; q q>= 10000 

r=1 
C=10, directly
  
? for(C=10,10,forprime(q=10000,40000,if(newF(C,r,q)>0.5,print("C = ", C, " q = ", q, " newF(", C,",", r,",q) = ",newF(C, r, q)); break )))
C = 10 q = 28537 newF(10,1,q) = 0.50001016504330466014731973552983831978
   
? for(C=10,10,forprime(q=3,30000,if(newG1(C,q)>0.5,print("C = ", C, " q = ", q, " newG1(C,q) = ",newG1(C,q)); break )))
C = 10 q = 599 newG1(C,q) = 0.50027116639395622778298528489337643909

--------
r = 2   
 C=10, directly
 
 for(C=10,10,forprimestep(q=10000,10^8,Mod(1,2*r), if(newF(C,r,q)>0.5,print("q = ", q, " newF(", C,",", r,",q) = ",newF(C, r, q)); break )))
q = 1160893 newF(10,2,q) = 0.50000084502053473734588587699425680547


--------
r = 3 
C=10, directly
 
? for(C=10,10,forprimestep(q=1160893,10^8,Mod(1,2*r), if(newF(C,r,q)>0.5,print("q = ", q, " newF(", C,",", r,",q) = ",newF(C, r, q)); break )))
? for(C=10,10,forprimestep(q=10^8,10^9,Mod(1,2*r), if(newF(C,r,q)>0.5,print("q = ", q, " newF(", C,",", r,",q) = ",newF(C, r, q)); break )))
? for(C=10,10,forprimestep(q=10^9,10^(10),Mod(1,2*r), if(newF(C,r,q)>0.5,print("q = ", q, " newF(", C,",", r,",q) = ",newF(C, r, q)); break )))
q = 2089575931 newF(10,3,q) = 0.50000000004587560010937129102625012157 

--------
r = 4
C=10, directly
 
? for(C=10,10,forprimestep(q=2089575931,10^(10),Mod(1,2*r), if(newF(C,r,q)>0.5,print("q = ", q, " newF(", C,",", r,",q) = ",newF(C, r, q)); break )))

NO OUTPUT for q up to 10^(10) so there are no such q-s such that newF(4,q) > 0.5

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