######################################
#######################################
## Computation about GA1 numbers
######################################
######################################

Digits:=20;

read "lisSA300.m";

read "lisCA300.m";

expga:=evalf(exp(gamma));

Gronw:=proc(n);
# computes the Gronwall function G(n)
numtheory[sigma](n)/n/log(log(evalf(n)));
end;


############################################
############################################
## Tests for deciding whether an integer is GA1
##################################################
############################################

critGA1:=proc(n) local lis,lisp,noli,i,maxi,bou,p,pm;
# This is a crude test for proper GA1 numbers
# i.e. for GA1 numberS N satisfying OMEGA(N) > 2.
if isprime(n) then return [n,false,[[n,1]]]; fi;
if irem(n,2)=0 and isprime(n/2) then return [n,false,[[2,1],[n/2,1]]]; fi;
lis:=ifactors(n)[2];
lisp:=map(x->x[1],lis);
maxi:=Gronw(n); bou:=true;
noli:=nops(lisp);
for i from noli by -1 to 1 while bou do
   p:=lisp[i];
   if Gronw(n/p) >= maxi then bou:=false; pm:=p;fi;
od;
# n is the number to be tested, bou=true or false, lis is the factorization
# of n and, if bou=false, pm is the largest prime factor of n 
# for which G(n/p) >= G(n).
if bou then [n,bou,lis] else [n,bou,lis,pm]; fi;
end;


fastGA1:=proc(n,eps) local lis,lisso,lisOm,nf,logn,Om,omeg,p,pf,r,pr,Q,Qp,i,sig,lln,epsp1;
# fast test for proper GA1 numbers; it uses the set S(n).
# more exactly, it tests whether the Gronwall's quotient Q(n) 
# is smaller than 1+eps
epsp1:=1.+eps;
lis:=ifactors(n)[2]; omeg:=nops(lis);
lisso:=sort(lis,proc(a,b) evalb(a[1] <= b[1]) end);
lisOm:=map(x-> x[2],lisso);
Om:=add(i,i=lisOm); # Om is equal to OMEGA(n)
if Om <= 2 then return [n,false,`OMEGA(n)=`,Om]; fi;
nf:=evalf(n); logn:=log(nf); # here we have OMEGA(n) >= 3
lln:=log(logn);
p:=lisso[omeg][1]; # p is the largest prime factor of n
pf:=evalf(p); r:=lisOm[omeg]; # r is equal to v_p(n)
if r=1 then 
   sig:=evalf(p/(p+1));
else   
   pr:=p^(r+1); sig:=evalf((pr-p)/(pr-1));
fi;
Q:=sig*lln/log(logn-log(pf));
if Q > epsp1 then return [n,false,`Q(n) >=`,Q,p]; fi;
for i from omeg -1 by -1 to 1 do
   if lisOm[i] > r then 
      p:=lisso[i][1]; pf:=evalf(p); r:=lisOm[i];  
      pr:=p^(r+1); sig:=evalf((pr-p)/(pr-1));
      Qp:=sig*lln/log(logn-log(pf));
      if Qp > epsp1 then return [n,false,`Q(n) >= `,Qp,p]; fi;
      Q:=max(Q,Qp);
   fi;
od;
[n,true,`Gronwall quotient=`,Q]; # pay attention, if Q > 1 then n is not GA1.
end;


gronwallquotient:=proc(n) 
local lis,lisso,lisOm,nf,logn,Om,omeg,p,pf,r,pr,Q,Qp,i,sig,lln,resu,pQ;
# computes the Gronwall quotient Q(n) = max for p|n of G(n/p)/G(n)
# defines the set S(n) and uses the formula Q(n) = max for p in S(n) of G(n/p)/G(n)
# This procedure is very similar to "fastGA1"
lis:=ifactors(n)[2]; omeg:=nops(lis);
lisso:=sort(lis,proc(a,b) evalb(a[1] <= b[1]) end);
lisOm:=map(x-> x[2],lisso);
Om:=add(i,i=lisOm);
if Om <= 1 then return [n,"the Gronwall quotient is not defined",op(lis)]; fi;
if Om = 2 then
   if nops(lisso) = 1 then # here, n is the square of a prime
      p:=lisso[omeg][1]; Q:=Gronw(p)/Gronw(n);
   else # here, n is the product of two distinct primes
      p:=lisso[omeg][1]; Q:=max(Gronw(p)/Gronw(n),Gronw(n/p)/Gronw(n));
   fi;
   return [Q,"Om=2",op(lis)];
fi;
nf:=evalf(n); logn:=log(nf); # here we have OMEGA(n) >= 3
lln:=log(logn);
p:=lisso[omeg][1]; pf:=evalf(p); r:=lisOm[omeg];
if r=1 then 
   sig:=evalf(p/(p+1));
else   
   pr:=p^(r+1); sig:=evalf((pr-p)/(pr-1));
fi;
Q:=sig*lln/log(logn-log(pf)); resu:=[[p,r,Q]]; pQ:=p;
for i from omeg-1 by -1 to 1 do
   if lisOm[i] > r then 
      p:=lisso[i][1]; pf:=evalf(p); r:=lisOm[i];  
      pr:=p^(r+1); sig:=evalf((pr-p)/(pr-1));
      Qp:=sig*lln/log(logn-log(pf));
      resu:=[[Qp,p,r],op(resu)];
      if Qp > Q then Q:=Qp; pQ:=p; fi;
   fi;
od;
[Q,Om,pQ,op(resu)];
end;

############################################
#######################################################
##   Some technical procedures
#######################################################
############################################

cherche:=proc(liste,m) local a,b,c,ifi;
# "cherche" means look for; "liste" is a list of numbers in ascending order
# ifi is the number of elements of "liste"
# If  m > liste[ifi], returns  ifi+1; if m=liste[ifi] then returns ifi;
# if not,  returns c such that ini <= c < ifi and liste[c] <=m < liste[ifi].
if nops(liste)=0 or liste[1] >= m  then return 0; fi;
ifi:=nops(liste); 
if m > liste[ifi] then return ifi +1; fi;
if m = liste[ifi] then return ifi; fi;
a:=1; b:=ifi; # here, liste[1] < m < liste[ifi] 
while b-a > 1 do
   c:=iquo(a+b,2);
   if liste[c] <= m then a:=c; else b:=c; fi; # liste <= m < liste[b] holds
od;
a;
end;

shortenleft:=proc(liste,m) local a,b,c,ifi;
# "liste" is a list of integers in ascending order
# we want to take off the terms of this list smaller than "m"
ifi:=nops(liste);
if ifi = 0 or m <= liste[1] then 
   return liste;
elif m > liste[ifi] then 
   return [];
fi;
a:=1; b:=ifi; # here, liste[a] < m and liste[b] >= m hold
while b-a > 1 do
   c:=iquo(a+b,2);
   if liste[c] < m then 
      a:=c; 
   else 
      b:=c; # here liste[a] < m and liste[b] >= m hold
   fi;
od;
# here we have 1 <= a < a+1 = b <= nops(liste) and liste[a] < m and liste[b] >= m
liste[b..ifi];
end;

shortenright:=proc(liste,m) local a,b,c,ifi;
# "liste" is a list of integers in ascending order
# we want to take off the terms of this list greater than "m"
ifi:=nops(liste);
if ifi = 0 or m >= liste[ifi] then 
   return liste;
elif m < liste[1] then 
   return [];
fi;
a:=1; b:=ifi; # here, liste[a] <= m and liste[b] > m hold
while b-a > 1 do
   c:=iquo(a+b,2);
   if liste[c] <= m then 
      a:=c; 
   else 
      b:=c; # here liste[a] <= m and liste[b] > m hold
   fi;
od;
# here we have 1 <= a < a+1 = b <= nops(liste) and liste[a] <= m and liste[b] > m
liste[1..a];
end;

facto:=proc(n,pmin,pmax,P0,P) local a,b,d,m,mot,inc,ins,motN;
global lisCA300,lisSA300;
# prints the factorization of the number n provided that
# all prime factors of n are <= pmax
# all prime factors p >= pmin satisfies of v_p(n) <= 1
# pmax is the largest prime dividing n 
# P0 is the product of the primes < pmin
# P is the product of the primes <= pmax
inc:=cherche(lisCA300,n); # Here, CA and SA numbers are labelled
if n=lisCA300[inc] then
   motN:=`**N=`; # these numbers are CA and thus SA
else
  ins:=cherche(lisSA300,n);
   if n=lisSA300[ins] then
      motN:=`*N=`; # these numbers are SA
   else
      motN:=`N=`;  # these numbers are not SA
   fi;
fi;       
a:=P/P0;
d:=igcd(a,n);
b:=a/d;
m:=n/d;
if b=1 then mot:=`no`; else mot:=ifactor(b); fi;
[motN,ifactor(m),[pmin,`...`,pmax],`abs=`,mot];
end;

GA1:=proc(x,lisGA1,imp) local com,n,nf,pmin,pminm1,pmax,p,Pro,Pr,xf,lis,fa,
kmaxf,kmax,k,nu,S,tabp,lisxy,ix,iy;
# print the factorisation of the numbers of "lisGA1"
# and computes the number of GA1-numbers with OMEGA(N)=k
# and the numpers og GA1-numbers N such that P(N) = p
print();
print("PROPER GA1 NUMBERS N up to x=",evalf(x,2));
print();
print("*SA numbers are starred, **CA numbers are double starred.");
print("In the square bracket [pmin...pmax], pmax is the largest prime factor of N");
print(" and pmin is the smallest prime s.t. for all p's satisfying pmin <= p <= pmax, p^2 does not divide N.");
print(" The absent 'abs' are those primes p satisfying pmin <= p <= pmax and not dividing N");
print(" After the absent primes, the Gronwall quotient Q(N) is given");
print(   "followed by the prime q s.t. Q(N)=G(N/q)/G(N).");
print();print();print();
# iy:=cherche(lisGA1,y); if iy=0 then iy:=1; fi;
ix:=cherche(lisGA1,x); if ix > nops(lisGA1) then ix:=nops(lisGA1); fi;
lisxy:=lisGA1[1..ix];
xf:=evalf(x)*1.000001;
kmax:=log(xf)/log(2.); kmaxf:=log(xf)/log(log(xf)); 
for k from 1 to kmax do nu[k]:=0; od;
pmax:=prevprime(ceil(log(xf))+1);
Pr:=1; p:=2;
while p<= pmax do # Pro[p] is the product of the primes up to p
   Pr:=Pr*p; Pro[p]:=Pr; tabp[p]:=0;
   p:=nextprime(p);
od;
p:=2;
com:=0; # numbers the GA1 numbers  
for n in lisxy do
   com:=com+1; nf:=evalf(n);
   lis:=gronwallquotient(n);
   k:=lis[2];
   if k <= kmax then
      nu[k]:=nu[k]+1;
   fi;
   p:=pmax;
   while irem(n,p) > 0  do
      p:=prevprime(p);
   od; # here, p is the largest prime factor of n
   tabp[p]:=tabp[p]+1;
   pmin:=nextprime(ceil(sqrt(2.*log(nf)))+1);
   while irem(n,pmin^2) > 0 do 
      pmin:=prevprime(pmin);
   od;
   pmin:=nextprime(pmin);
   # here, pmin is the smallest prime satisfying v_pmin <= 1
   fa:=facto(n,pmin,p,Pro[pmin]/pmin,Pro[p]);
   if imp >= 6 then 
      print(com,op(fa),evalf(lis[1],9),lis[3]);
   fi;
od;
print();print();
print("x=",evalf(x,2),"log x/log 2=", evalf(kmax,5), "log x /log log x =",evalf(kmaxf,5));
print();print();print();
S:=0;
for k from 3 to kmax while S < com do 
   S:=S+nu[k]; 
   print(`k=`,k,"# of proper GA1 numbers N<=x s.t. OMEGA(N)= k (resp. <= k)",nu[k],S); 
od;
print(); print();
p:=2;
while p <= pmax do
   print(`p=`,p,"# of proper GA1 numbers N<=x s.t. P(N)=p:",tabp[p]);
   p:=nextprime(p); 
od;
end;


############################################
################################################
## Computation of GA1 numbers up to x
#################################################
############################################

algo1GA1:=proc(x,imp) local n,lis,lisdiv,M,p,r,pr,logx,ti,com,lisGA1;
# crude code to find proper GA1 numbers up to x
# it tests all divisors of M
# it takes 1/2 hour on my ebookfor x=10^20
p:=2; M:=1;ti:=time();
logx:=log(evalf(x))+0.001;
while p <= logx do # computation of M=M(x)
   pr:=p; r:=2;
   while p^r/r <= logx do r:=r+1; od;
   r:=r-1; M:=M*p^r;
   p:=nextprime(p);
od;
lisdiv:=numtheory[divisors](M);
print(`PROPER GA1 <= x=`,evalf(x,5),`M=`,M,ifactor(M),
      `nomber of divisors=`,nops(lisdiv),`time=`,time()-ti);
com:=0; lisGA1:=[];
for n in lisdiv do
   com:=com+1;
   if irem(com,50000)=0 then
      print(`com=`,com,time()-ti);
   fi;
   if n<= x and n > 4 then
      lis:=fastGA1(n,0.000001);  
      if imp > 2 then print(`n=`,n,lis); fi;
      if lis[2] and lis[4] <= 1 then 
         lisGA1:=[n,op(lisGA1)];
         print(`n=`,n, `is GA1`);
      fi;
   fi;
od;
print(`total time=`,time()-ti);
sort(lisGA1);
end;

# We do not give a code for the second algorithm for computing
# proper GA1 numbers, described in 7.3.

algo3GA1:=proc(y,x,eps,imp) 
local ti,epsp1,logx,M,M1,p0,p1,p,pf,ps,ppr,pri,logpri,i1,i2,P0,P,logpf,
thet,theta,pmth,di,mindi,Delta,listep,listeprac,lisdiv,lisdivrac,
lisM1,iM1,pe,e,ie,prev,iep,iey,
r,n,logn,lln,bou,restn,ri,qun,lis,sig,priprp1,Q,comGA1,lisGA1,lisGA1f,j,i,
epsQ,
epsmin1,epsmin2,Nepsmin1,Nepsmin2,comptep;
# computes all proper GA1 numbers N satisfying y < N <= x 
# by using as few as possible divisors of M
print();
print("GA1 NUMBERS N s.t.N > y=",evalf(y,5),"and <= x=",evalf(x,5));
print();print();
###########################################
# Computation of M and M1
###########################################
M:=1; M1:=1; epsp1:=1.+eps; epsmin2:=infinity; epsmin1:=-infinity; ti:=time();
p:=2; P0:=1; i1:=0; i2:=0; thet:=0.; 
logx:=ceil(log(evalf(x))+eps);
while p <= logx do # computation of M = M(x)
   ppr:=p^2; r:=2; pf:=evalf(p); i2:=i2+1;
   while ppr <= r*logx do r:=r+1; ppr:=ppr*p; od;
   r:=r-1; ppr:=ppr/p; M:=M*ppr;
   if r > 1 then 
      p0:=p; M1:=M1*ppr; P0:=P0*p; # P0 is the product of the primes up to p0
      i1:=i1+1; pri[i1]:=p; logpri[i1]:=log(pf); # pri[i] is the i-th prime
   fi;     
   thet:=thet+log(pf); pmth[p]:=pf-thet;
   theta[p]:=thet; # theta is the Chebichev's function in p
   p:=nextprime(p);
od;   
# here, p0 = pri[i1] is the largest prime factor of M1
# and the number of prime factors of M is i2
ps:=prevprime(p); # "ps" is the largest prime factor of M
print("M1=", evalf(M1,5),ifactor(M1),"# of divisors of M1=",numtheory[tau](M1));
print("M=",evalf(M,5),"prime factors of M/M1=",[nextprime(p0),`...`,ps]);
print("number of prime factors of M/M1=",i2-i1);
print();
####################################################
# Computation of Delta[p]
#####################################################
mindi:=pmth[ps]; di[ps]:=infinity; Delta[ps]:=infinity; pf:=evalf(ps);
if imp >= 4 then
   print(`ps=`,ps,`ps-th=`,evalf(pmth[ps],5),`di[ps]=`,di[ps],
         `P=`,evalf(exp(theta[ps]),6),`Delta=`,Delta[ps],`e^p=`,evalf(exp(pf),6));
fi;
p:=prevprime(ps); pf:=evalf(p); di[p]:=mindi; 
Delta[p]:=exp(theta[p]+mindi)/epsp1;
if imp >= 4 then
   print(`p=`,p,`th(p)=`,evalf(theta[p],5),`p-th=`,evalf(pmth[p],5),
         `di[p]=`,evalf(di[p],5),`Delta=`,evalf(Delta[p],5));
fi;
while p > 2 do
   mindi:=min(pmth[p],mindi);
   p:=prevprime(p); pf:=evalf(p);
   di[p]:=mindi; Delta[p]:=exp(theta[p]+mindi)/epsp1;
   if imp >= 4 then
      print(`p=`,p,`th[p]=`,evalf(theta[p],5),`p-th=`,evalf(pmth[p],5),
            `di[p]=`,evalf(di[p],5),`Delta=`,evalf(Delta[p],5)); 
   fi;      
od;
###################################################
# Study of the divisors of M1
###################################################
comGA1:=0; lisGA1:=[]; 
comptep:=0; # will count the numbers of GA1 numbers N with P(N)=p
lisM1:=ifactors(M1)[2];
lisM1:=sort(lisM1,proc(a,b) evalb(a[1] <= b[1]) end);
lisdiv:=[1];
for iM1 from 1 to nops(lisM1) do
   p:=lisM1[iM1][1]; e:=lisM1[iM1][2]; pf:=evalf(p);
   prev:=lisdiv;
   pe:=1;
   for ie from 1 to e do
      pe:=p*pe;
      lisdiv:=[op(lisdiv),op(map(a->a*pe,prev))];
   od;
   lisdiv:=sort(lisdiv);
   iep:=cherche(lisdiv,exp(pf)/epsp1);
   iey:=cherche(lisdiv,y);
   print(`p=`,p,`#lisdiv=`,nops(lisdiv),`iep=`,evalf(iep,6),`Delta=`,evalf(Delta[p],5));
   for j from max(iep,iey+1) to nops(lisdiv) do
      n:=lisdiv[j]; # GA1 test of the numbers of "lisdiv"
      if fastGA1(n,eps)[2] then
         lisGA1:=[n,op(lisGA1)];
         lis:=gronwallquotient(n);
         epsQ:=lis[1]-1.;
         if epsQ >= 0 then
            if epsQ < epsmin2 then 
               epsmin2:=epsQ; Nepsmin2:=ifactor(n);
            fi;
         else
            if epsQ > epsmin1 then 
               epsmin1:=epsQ; Nepsmin1:=ifactor(n);
            fi;
         fi;
         if   lis[1] >= 1. then 
            if imp >= 3 then 
               print(`FALSE GA1*********`,n,ifactor(n),lis[1]);
            fi;
         else
            comGA1:=comGA1+1;
            if imp >= 3 then print(comGA1,n,ifactor(n),lis[1]); fi;
         fi;
      fi;
   od;
   print(`p=`,p,`comptep=`,comGA1-comptep,`temps=`,time()-ti); 
   comptep:=comGA1;
od;
##################################################
# The divisors of M not dividing M1
##################################################
lisdiv:=shortenleft(lisdiv,Delta[p0]); # take off the divisors of M1 smaller than Delta[p0]
p1:=nextprime(p0); comptep:=0; P:=P0; # P is the product of primes <= p
p:=p1;
while p <= logx do # loop on large p's
   pf:=evalf(p); logpf:=log(pf); P:=P*p;
   listep:=map(t->t*p,lisdiv); 
   # listep contains the divisors of M with largest prime factor p
   # which can be mutiplied by primes > p to get a GA1 number
   listep:=shortenright(listep,x); # here, we take off the divisors greater than x
   iep:=cherche(listep,exp(pf)/epsp1);
   iey:=cherche(listep,y);
   print(`p=`,p,`#listep=`,nops(listep),`iep=`,evalf(iep,6),`P=`,evalf(P,6),
         `Delta=`,evalf(Delta[p],6),`exp=`,evalf(exp(pf),6),time()-ti);
   for j from max(iep,iey+1) to nops(listep) do  # GA1 test of the numbers of "listep"
      n:=listep[j];
      logn:=log(evalf(n)); lln:=log(logn);
      sig:=pf/(pf+1.);
      Q:=sig*lln/log(logn-logpf);
      bou:=evalb(Q < epsp1);
      r:=1;
      for i from i1 by -1 to 1 while bou do
         restn:=irem(n,pri[i],'qun');
         if restn=0 then
            ri:=1; # calculation of the exponent "ri" of "pri[i]" in "n"
            while irem(qun,pri[i],'qun')=0 do ri:=ri+1; od;
            if ri > r then # here, pri[i] belongs to the set S(n)
               r:=ri;
               priprp1:=pri[i]^(r+1); sig:=evalf((priprp1-pri[i])/(priprp1-1));
               Q:=sig*lln/log(logn-logpri[i]);
               if Q > epsp1 then bou:=false; fi;
            fi;
         fi;
      od;      
      if bou then 
         lisGA1:=[n,op(lisGA1)];     
         lis:=gronwallquotient(n);
         epsQ:=lis[1]-1.;
         if epsQ >= 0 then
            if epsQ < epsmin2 then 
               epsmin2:=epsQ; Nepsmin2:=op(facto(n,p1,p,P0,P));
            fi;
         else
            if epsQ > epsmin1 then 
               epsmin1:=epsQ; Nepsmin1:=op(facto(n,p1,p,P0,P));
            fi;
         fi;
         if   lis[1] >= 1. then 
            if imp >= 3 then 
               print(`FALSE GA1*******`,op(facto(n,p1,p,P0,P)),lis[1]);
            fi;
         else
            comGA1:=comGA1+1;
            if imp >= 3 then print(comGA1,op(facto(n,p1,p,P0,P)),lis[1]); fi;
         fi;
      fi;
   od;
   ############################################
   # Elimination  of bad  divisors
   ############################################
   lisdivrac:=shortenleft(lisdiv,Delta[p]);
   listeprac:=shortenleft(listep,Delta[p]);
   lisdiv:=[op(listeprac),op(lisdivrac)];
   lisdiv:=sort(lisdiv);
   lisdiv:=shortenright(lisdiv,x);
   # here, lisdiv is the list of divisors <= x of M with largest prime factor <= p
   # and not smaller than Delta[p]. Some of them can be > x
   print(`p=`,p,"comptep=",comGA1-comptep,"temps=",time()-ti);
   print("#lisdivrac=",nops(lisdivrac),"#listeprac=",nops(listeprac),
         "#lisdiv=",nops(lisdiv));
   p:=nextprime(p); comptep:=comGA1;
od;
#########################################
# Conclusion
###########################################
print("epsmin2=",epsmin2,"Nepsmin2=",Nepsmin2);
# 1+epsmin2 is the minimal value of the Gronwall quotient Q(n) 
# for all tested n s.t. Q(n) > 1.
print("epsmin1=",epsmin1,"Nepsmin1=",Nepsmin1);
# 1-epsmin1 is the maximal value of the Gronwall quotient Q(n) 
# for all tested n s.t. Q(n) < 1.
Digits:=40;
lisGA1f:=[]; # Here, we remove the selected n's with Q(n) > 1
for n in lisGA1 do
   lis:=fastGA1(n,0);
   if lis[2] then lisGA1f:=[n,op(lisGA1f)]; fi;
od;
Digits:=20;
lisGA1:=sort(lisGA1f);
save lisGA1,"lisGA1.m";
print(`total number of GA1=`,nops(lisGA1),`total time=`,time()-ti);
lisGA1;
end;