######################################################################
##SALIKHOV.txt: Save this file as  SALIKHOV.txt                      #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read  SALIKHOV.txt<Enter>                          #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#DoronZeil at gmail dot com                                          #
######################################################################
 
#Created: Nov. 2019
 
print(`Created:  Nov. 2019`):
print(` This is SALIKHOV.txt `):
print(`It is one of the Maple package that accompanies the article `):
print(` Automatic Discovery of Irrationailty Proofs and Independence Measures`):
print(`by Doron Zeilberger and Wadim Zudiln`):

print(`and also available from Zeilberger's website`):
print(``):
print(`Please report bugs to DoronZeil at gmail dot com `):
print(``):
 print(`The most current version of this  package and paper`):
 print(` are  available from`):
 print(`http://sites.math.rutgers.edu/~zeilberg/  .`):



print(`---------------------------------------`):
 print(`For a list of the procedures for general families type ezraG();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):

print(`---------------------------------------`):
 print(`For a list of the Supporting procedures type ezra1();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):

print(`---------------------------------------`):

print(`---------------------------------------`):
 print(`For a list of the MAIN procedures type ezra();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

with(combinat):



ezraG:=proc()

if args=NULL then
 print(` The  procedures for the general cases are: `):
 print(` NuEven, NuOdd , PaperG `):
 print(``):

else
ezra(args):
fi:

end:

ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are: Atoms1, Atoms2, dn, ExtractLogs, GrC, MakeNice, OpeNice, Ra1, RaE1, RaE2, `):
 print(` Sabc1, Sabc1f, Sabc2, Sabc2, SalPairs, SalPairsPC, SInt `):
 print(``):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are: deltE, deltT1, deltT2, Info1, Info2, INFOab, GuessK1, GuessK2, Ope, Paper, Sabc1f, Sabc2f, Shors, Zug1, Zug2 `):
 print(` `):




elif nargs=1 and args[1]=Atoms1 then
print(`Atoms1(a,b,c): The atoms of Sabc1(a,b,c,n) for all n. Try:`):
print(`Atoms1(5,7,35);`):


elif nargs=1 and args[1]=Atoms2 then
print(`Atoms2(a,b,c): The atoms of Sabc2(a,b,c,n) for all n. Try:`):
print(`Atoms2(5,7,35);`):

elif nargs=1 and args[1]=deltE then
print(`deltE(f,X): Given f a linear combination of 1, X[1],.., X[nops(X)] and 1`):
print(`with integer coefficients, finds`):
print(` the empirical delta such that abs(f)=1/max(coeff(f,1),coeff(f,X[1]),coeff(f,X[2]))^delta. Try: `):
print(` deltE(Sabc1f(5,7,35,200),[log(2),log(3)]); `):
print(`deltE(Sabc2f(5,7,35,200),[log(2),log(3)]);`):


elif nargs=1 and args[1]=deltT1 then
print(`deltT1(a,b,c): The theoretical delta for Salikhov's Diophantine approximations inspired by Sabc1(a,b,c,n). Try:`):
print(`deltT1(5,7,35);`):


elif nargs=1 and args[1]=deltT2 then
print(`deltT2(a,b,c): The theoretical delta for Salikhov's Diophantine approximations inspired by Sabc2(a,b,c,n). Try:`):
print(`deltT2(5,7,35);`):

elif nargs=1 and args[1]=dn then
print(`dn(n): lcm(1..n). Try: `):
print(`dn(100);`):

elif nargs=1 and args[1]=ExtractLogs then
print(`ExtractLogs(L): Given an expression that is a linear combination of log(k)'s  and 1, outputs the set of k's that show up.`):
print(`Try: `):
print(`ExtractLogs(5*log(3)-11*log(5));`):



elif nargs=1 and args[1]=Info1 then
print(`Info1(a,b,c): inputs positive integers a,b,c, and outputs a list consisting of`):
print(`(i) the set of integers such that Sabc1(a,b,c,n) are linear combinations (with rational number coefficients) of their logs`):
print(`(ii) The conjectured integer K such that K^n*dn(2*n)*Sabc1(a,b,c,n) have integer coefficients`):
print(`(iii) The theoretical delta such that`):
print(`T(a,b,c,n):=Sabc1(a,b,c,n)*K^n*dn(2*n), has the property, that for some universal constant CONST`):
print(`|T(a,b,c,n)|<=CONST/(max(abs(coeff.s))^delta .`):
print(`(iv): The smallest empirical delta between 100 and 150. Try:`):
print(`Info1(5,7,35); `):


elif nargs=1 and args[1]=Info2 then
print(`Info2(a,b,c): inputs positive integers a,b,c, and outputs a list consisting of`):
print(`(i) the set of integers such that Sabc1(a,b,c,n) are linear combinations (with rational number coefficients) of their logs`):
print(`(ii) The conjectured integer K such that K^n*dn(2*n)*Sabc2(a,b,c,n) have integer coefficients`):
print(`(iii) The theoretical delta such that`):
print(`T(a,b,c,n):=Sabc2(a,b,c,n)*K^n*dn(2*n), has the property, that for some universal constant CONST`):
print(`|T(a,b,c,n)|<=CONST/(max(abs(coeff.s))^delta .`):
print(`(iv): The smallest empirical delta between 100 and 150. Try:`):
print(`Info2(5,7,35); `):


elif nargs=1 and args[1]=INFOab then
print(`INFOab(a,b,x): inputs two integers 1<=a<b and outputs the information about a Salikhov type theorem `):
print(`of simultaneous approximations to log(b/(b+1)) and log(a/(a+1)). `):
print(` [b/(b+1), a/(a+1)], [K1,K2, nu], where K1 is such that `):
print(`Zug1(a,b,n)*lcm(1..2*n)*K1^n and Zug2(a,b,n)*lcm(1..2*n)*K2^n are integers`):
print(`and nu is an upper bound for the linear independence measure of 1, log(b/(b+1)), log(a/(a+1))`):
print(` i.e. if Q=max(|q|,|p1|,|p2|), Q>=Q0 is sufficiently large, then `):
print(` |q+p1*log((b-1)/(b+1))+p2*log((a-1)/(a+1)))|>1/Q^nu `):
print(` It also outputs empirical delta1 and delta2.`):
print(`It also outputs the indicial polynomial in the variable x, and its roots. `):
print(` To get a synopsis of the original Salikhov proof, type: `):
print(`INFOab(2,3,x);`):

elif nargs=1 and args[1]=GrC then
print(`GrC(a,b,c): The growth constants of Sabc1(a,b,c,n) and Sabc2(a,b,c,n)`):
print(`It also gives the locations (or rather their squares). Try:`):
print(`GrC(5,7,35); `):


elif nargs=1 and args[1]=GuessK1 then
print(`GuessK1(a,b,c,N): guesses an integer K and a small integer L such that Sabc1(a,b,c,n)*dn(2*n)*K^n*L is an integer for all n`):
print(`between 20 and N. The output is [K,L].`):
print(`Try: `):
print(`GuessK1(5,7,35,100); `):


elif nargs=1 and args[1]=GuessK2 then
print(`GuessK2(a,b,c,N): guesses an integer K and a small integer L such that Sabc2(a,b,c,n)*dn(2*n)*K^n*L is an integer for all n`):
print(`between 20 and N. The output is [K,L].`):
print(`Try: `):
print(`GuessK2(5,7,35,100); `):

elif nargs=1 and args[1]=MakeNice then
print(`MakeNice(ope,n,N): inputs a recurrence operator ope in n, N, outputs Ope1(n,1/N) such that`):
print(`1-Ope1(n,N) is equivalent fo Ope. Try:`):
print(`MakeNice(N^2-n*N-1,n,N);`):


elif nargs=1 and args[1]=NuEven then
print(`NuEven(a,N): a closed form expression for the exponent nu in the linear independence measure of {1,log(a/(a+1)), log((a+1)/(a+2))}`):
print(`when a is a positive even integer. `):		   
print(`It also returns the cubic in N whose largest and smallest roots are used in the expression. Try:`):
print(`NuEven(a,N):`):


elif nargs=1 and args[1]=NuOdd then
print(`NuOdd(a,N): a closed form expression for the exponent nu in in the linear independence measure of {1,log(a/(a+1)), log((a+1)/(a+2))}`):
print(`when a is a positive odd integer that is a>=3. When a=1 it is  20.01872047...`):
print(`It also returns the cubic in N whose largest and smallest roots are used in the expression. Try:`):
print(`NuOdd(a,N):`):

elif nargs=1 and args[1]=Ope then
print(`Ope(a,b,c,n,N): The third-order linear recurrence operator, where N is the forward shift in n annihilating both`):
print(`Sabc1(a,b,c,n) and Sabc2(a,b,c,n) (q.v.) it was obtained via the Maple package EKHAD, procedure`):
print(`AZd, implementing the Almkvist-Zeilberger algorithm. Try:`):
print(`Ope(a,b,c,n,N);`):
print(`Ope(5,7,35,n,N);`):


elif nargs=1 and args[1]=OpeNice then
print(`OpeNice(a,b,c,n,N): the operator in 1/N, let's call it Ope1,  such that Ope(a,b,c,n,N) is equivalent fo 1-Ope1. Try:`):
print(`OpeNice(5,7,35,n,N);`):

elif nargs=1 and args[1]=Paper then
print(`Paper(M): States  and sketches the proof of simultaneous linear diophantine approximation`):
print(`for log((b-1)/(b+1)) and log((a-1)/(a+1)) for all 3<=a<b<M a and b odd. Try:`):
print(`Paper(5); `):


elif nargs=1 and args[1]=PaperG then
print(`PaperG(a,N):inputs symbols a and N and outputs an article with a proof `):
print(`(modulo divisibility lemmas left the reader) of the linear independence `):
print(`of {1,log(a/(a+1)),log((a+1)/(a+2))}, together with explicit upper bounds`):
print(`for the measure of indepenedce (different for a even and odd). Try:`):
print(`PaperG(a,N);`):

elif nargs=1 and args[1]=Ra1 then
print(`Ra1(a,b,c,N): The rations of the sequence of coefficients of the atoms of Sabc1(a,b,c) with N. Try:`):
print(`Ra1(5,7,35,100); `):

elif nargs=1 and args[1]=RaE1 then
print(`RaE1(a,b,c,N):abs(Sabc1f(a,b,c,N)/Sabc1f(a,b,c,N-1)). Try: `):
print(`RaE1(5,7,35,30);`):


elif nargs=1 and args[1]=RaE2 then
print(`RaE2(a,b,c,N):abs(Sabc2f(a,b,c,N)/Sabc2f(a,b,c,N-1)). Try: `):
print(`RaE2(5,7,35,30);`):

elif nargs=1 and args[1]=Sabc1 then
print(`Sabc1(a,b,c,n1): The generalized Salikhov's integral from 0 to a, straight from the definition. in other words`):
print(`int(x^(2*n1)*(x^2-a^2)^n1*(x^2-b^2)^n1/(x^2-c^2)^(2*n1+1),x=0..a):`):
print(`Pretty slow. You should use the much faster Sabc1f(a,b,c,n1). Try `):
print(`[seq(Sabc1(5,7,35,n1),n1=0..20)];`); 

elif nargs=1 and args[1]=Sabc1f then
print(`Sabc1f(a,b,c,n1): The generalized Salikhov's integral from 0 to a, Fast version. in other words`):
print(`int(x^(2*n1)*(x^2-a^2)^n1*(x^2-b^2)^n1/(x^2-c^2)^(2*n1+1),x=0..a): Try: `):
print(`[seq(Sabc1f(5,7,35,n1),n1=0..20)];`); 



elif nargs=1 and args[1]=Sabc2 then
print(`Sabc2(a,b,c,n1): The generalized Salikhov's integral from 0 to b, straight from the definition. in other words`):
print(`int(x^(2*n1)*(x^2-a^2)^n1*(x^2-b^2)^n1/(x^2-c^2)^(2*n1+1),x=0.ba):`):
print(`Pretty slow. You should use the much faster Sabc2f(a,b,c,n1). Try `):
print(`[seq(Sabc2(5,7,35,n1),n1=0..20)];`); 


elif nargs=1 and args[1]=Sabc2f then
print(`Sabc2f(a,b,c,n1): The generalized Salikhov's integral from 0 to b, straight from the RECURRENCE. in other words`):
print(`int(x^(2*n1)*(x^2-a^2)^n1*(x^2-b^2)^n1/(x^2-c^2)^(2*n1+1),x=0..b), but using the recurrence.`):
print(` Try `):
print(`[seq(Sabc2f(5,7,35,n1),n1=0..20)]; `):


elif nargs=1 and args[1]=SalPairs then
print(`SalPairs(N): all the pairs of primes [p,q] p<q, up to the N-th prime such that Atoms1(p,q,p*q)=Atoms2(p,q,p*q)`):
print(`Try: `):
print(`SalPairs(10);`):


elif nargs=1 and args[1]=SalPairsPC then
print(`SalPairsPC():  a pre-computed version of SalPairs(20); `):
print(`Try: `):
print(`SalPairsPC();`):

elif nargs=1 and args[1]=Shors then
print(`Shors(a,b,c,N): The indicial polynomial and three roots of the limiting indicial equation of Ope(a,b,c,n,N). Try:`):
print(`Shors(5,7,35,N);`):

elif nargs=1 and args[1]=SInt then
print(`SInt(a,b,c,n1,x): the integrand in the general Salikhov integral. Try:`):
print(`SInt(5,7,35,n,x);`):

elif nargs=1 and args[1]=Zug1 then
print(`Zug1(a,b,n1): The  pair of rational number [A,B] such that Sabc1(2*a+1,2*b+1,(2*a+1)*(2*b+1),n1)=A+B*log(b/(b+1))`):
print(`Try: `):
print(`Zug1(2,3,10);`):


elif nargs=1 and args[1]=Zug2 then
print(`Zug2(a,b,n1): The  pair of rational number [A,B] such that Sabc2(2*a+1,2*b+1,(2*a+1)*(2*b+1),n1)=A+B*log(a/(a+1))`):
print(`Try: `):
print(`Zug2(2,3,10);`):

else
print(`There is no ezra for`,args):
fi:
 
end:

Digits:=3000:

#SInt(a,b,c,n1,x): the integrand in the general Salikhov integral. Try:
#SInt(5,7,35,n,x);
SInt:=proc(a,b,c,n1,x)
 x^(2*n1)*(x^2-a^2)^n1*(x^2-b^2)^n1/(x^2-c^2)^(2*n1+1):
end:

#Sabc1(a,b,c,n1): The generalized Salikhov's integral from 0 to a. The original is Aabc(5,7,35,n1);
Sabc1:=proc(a,b,c,n1) local x:
option remember:
int(x^(2*n1)*(x^2-a^2)^n1*(x^2-b^2)^n1/(x^2-c^2)^(2*n1+1),x=0..a):
end:


#Sabc2(a,b,c,n): The generalized Salikhov's integral from 0 to a. The original is Aabc(5,7,35,n1);
Sabc2:=proc(a,b,c,n1) local x:
option remember:
int(x^(2*n1)*(x^2-a^2)^n1*(x^2-b^2)^n1/(x^2-c^2)^(2*n1+1),x=0..b):
end:

#Ope(a,b,c,n,N): The third-order linear recurrence operator, where N is the forward shift in n annihilating both
#Sabc1(a,b,c,n) and Sabc2(a,b,c,n) (q.v.) it was obtained via the Maple package EKHAD, procedure
#AZd, implementing the Almkvist-Zeilberger algorithm. Try:
#Ope(a,b,c,n,N);
#Ope(5,7,35,n,N);

Ope:=proc(a,b,c,n,N):
-4*a^4*b^4*(2*n+1)*(n+2)*(n+1)*(a-b)^2*(a+b)^2*(8*a^6*b^6*n^2-48*a^6*b^4*c^2*n^2+96*a^6*b^2*c^4*n^2-64*a^6*c^6*n^2-48*a^4*b^6*c^2*n^2+48*a^4*b^4*c^4*n^2-48*a^4*b^2*
c^6*n^2+72*a^4*c^8*n^2+96*a^2*b^6*c^4*n^2-48*a^2*b^4*c^6*n^2-72*a^2*b^2*c^8*n^2-64*b^6*c^6*n^2+72*b^4*c^8*n^2+32*a^6*b^6*n-204*a^6*b^4*c^2*n+426*a^6*b^2*c^4*n-292*a
^6*c^6*n-204*a^4*b^6*c^2*n+180*a^4*b^4*c^4*n-192*a^4*b^2*c^6*n+330*a^4*c^8*n+426*a^2*b^6*c^4*n-192*a^2*b^4*c^6*n-348*a^2*b^2*c^8*n-292*b^6*c^6*n+330*b^4*c^8*n+30*a^
6*b^6-198*a^6*b^4*c^2+441*a^6*b^2*c^4-315*a^6*c^6-198*a^4*b^6*c^2+126*a^4*b^4*c^4-159*a^4*b^2*c^6+357*a^4*c^8+441*a^2*b^6*c^4-159*a^2*b^4*c^6-408*a^2*b^2*c^8-315*b^
6*c^6+357*b^4*c^8)-2*(n+2)*(64*a^12*b^10*n^4-256*a^12*b^8*c^2*n^4-128*a^12*b^6*c^4*n^4+1792*a^12*b^4*c^6*n^4-2560*a^12*b^2*c^8*n^4+1024*a^12*c^10*n^4+64*a^10*b^12*n
^4-1280*a^10*b^10*c^2*n^4+3648*a^10*b^8*c^4*n^4-6656*a^10*b^6*c^6*n^4+5440*a^10*b^4*c^8*n^4+384*a^10*b^2*c^10*n^4-1152*a^10*c^12*n^4-256*a^8*b^12*c^2*n^4+3648*a^8*b
^10*c^4*n^4-3072*a^8*b^8*c^6*n^4+3200*a^8*b^6*c^8*n^4-7680*a^8*b^4*c^10*n^4+2880*a^8*b^2*c^12*n^4-128*a^6*b^12*c^4*n^4-6656*a^6*b^10*c^6*n^4+3200*a^6*b^8*c^8*n^4+
6656*a^6*b^6*c^10*n^4-1152*a^6*b^4*c^12*n^4+1792*a^4*b^12*c^6*n^4+5440*a^4*b^10*c^8*n^4-7680*a^4*b^8*c^10*n^4-1152*a^4*b^6*c^12*n^4-2560*a^2*b^12*c^8*n^4+384*a^2*b^
10*c^10*n^4+2880*a^2*b^8*c^12*n^4+1024*b^12*c^10*n^4-1152*b^10*c^12*n^4+384*a^12*b^10*n^3-1632*a^12*b^8*c^2*n^3-624*a^12*b^6*c^4*n^3+11328*a^12*b^4*c^6*n^3-16608*a^
12*b^2*c^8*n^3+6720*a^12*c^10*n^3+384*a^10*b^12*n^3-7872*a^10*b^10*c^2*n^3+22704*a^10*b^8*c^4*n^3-43200*a^10*b^6*c^6*n^3+36480*a^10*b^4*c^8*n^3+2112*a^10*b^2*c^10*n
^3-7584*a^10*c^12*n^3-1632*a^8*b^12*c^2*n^3+22704*a^8*b^10*c^4*n^3-16896*a^8*b^8*c^6*n^3+18528*a^8*b^6*c^8*n^3-50592*a^8*b^4*c^10*n^3+19248*a^8*b^2*c^12*n^3-624*a^6
*b^12*c^4*n^3-43200*a^6*b^10*c^6*n^3+18528*a^6*b^8*c^8*n^3+46272*a^6*b^6*c^10*n^3-8016*a^6*b^4*c^12*n^3+11328*a^4*b^12*c^6*n^3+36480*a^4*b^10*c^8*n^3-50592*a^4*b^8*
c^10*n^3-8016*a^4*b^6*c^12*n^3-16608*a^2*b^12*c^8*n^3+2112*a^2*b^10*c^10*n^3+19248*a^2*b^8*c^12*n^3+6720*b^12*c^10*n^3-7584*b^10*c^12*n^3+800*a^12*b^10*n^2-3568*a^
12*b^8*c^2*n^2-904*a^12*b^6*c^4*n^2+24728*a^12*b^4*c^6*n^2-37408*a^12*b^2*c^8*n^2+15344*a^12*c^10*n^2+800*a^10*b^12*n^2-16608*a^10*b^10*c^2*n^2+48264*a^10*b^8*c^4*n
^2-97216*a^10*b^6*c^6*n^2+85584*a^10*b^4*c^8*n^2+3576*a^10*b^2*c^10*n^2-17352*a^10*c^12*n^2-3568*a^8*b^12*c^2*n^2+48264*a^8*b^10*c^4*n^2-28464*a^8*b^8*c^6*n^2+34544
*a^8*b^6*c^8*n^2-115776*a^8*b^4*c^10*n^2+44880*a^8*b^2*c^12*n^2-904*a^6*b^12*c^4*n^2-97216*a^6*b^10*c^6*n^2+34544*a^6*b^8*c^8*n^2+113392*a^6*b^6*c^10*n^2-19656*a^6*
b^4*c^12*n^2+24728*a^4*b^12*c^6*n^2+85584*a^4*b^10*c^8*n^2-115776*a^4*b^8*c^10*n^2-19656*a^4*b^6*c^12*n^2-37408*a^2*b^12*c^8*n^2+3576*a^2*b^10*c^10*n^2+44880*a^2*b^
8*c^12*n^2+15344*b^12*c^10*n^2-17352*b^10*c^12*n^2+672*a^12*b^10*n-3112*a^12*b^8*c^2*n-420*a^12*b^6*c^4*n+22128*a^12*b^4*c^6*n-34702*a^12*b^2*c^8*n+14460*a^12*c^10*
n+672*a^10*b^12*n-14032*a^10*b^10*c^2*n+40980*a^10*b^8*c^4*n-89424*a^10*b^6*c^6*n+82994*a^10*b^4*c^8*n+1968*a^10*b^2*c^10*n-16374*a^10*c^12*n-3112*a^8*b^12*c^2*n+
40980*a^8*b^10*c^4*n-14208*a^8*b^8*c^6*n+22748*a^8*b^6*c^8*n-108996*a^8*b^4*c^10*n+43278*a^8*b^2*c^12*n-420*a^6*b^12*c^4*n-89424*a^6*b^10*c^6*n+22748*a^6*b^8*c^8*n+
116112*a^6*b^6*c^10*n-20136*a^6*b^4*c^12*n+22128*a^4*b^12*c^6*n+82994*a^4*b^10*c^8*n-108996*a^4*b^8*c^10*n-20136*a^4*b^6*c^12*n-34702*a^2*b^12*c^8*n+1968*a^2*b^10*c
^10*n+43278*a^2*b^8*c^12*n+14460*b^12*c^10*n-16374*b^10*c^12*n+180*a^12*b^10-852*a^12*b^8*c^2-54*a^12*b^6*c^4+6786*a^12*b^4*c^6-11115*a^12*b^2*c^8+4725*a^12*c^10+
180*a^10*b^12-3768*a^10*b^10*c^2+11034*a^10*b^8*c^4-27792*a^10*b^6*c^6+27921*a^10*b^4*c^8+60*a^10*b^2*c^10-5355*a^10*c^12-852*a^8*b^12*c^2+11034*a^8*b^10*c^4+1692*a
^8*b^8*c^6+2454*a^8*b^6*c^8-35313*a^8*b^4*c^10+14535*a^8*b^2*c^12-54*a^6*b^12*c^4-27792*a^6*b^10*c^6+2454*a^6*b^8*c^8+42336*a^6*b^6*c^10-7344*a^6*b^4*c^12+6786*a^4*
b^12*c^6+27921*a^4*b^10*c^8-35313*a^4*b^8*c^10-7344*a^4*b^6*c^12-11115*a^2*b^12*c^8+60*a^2*b^10*c^10+14535*a^2*b^8*c^12+4725*b^12*c^10-5355*b^10*c^12)*N-(2*n+3)*(32
*a^10*b^10*n^4+64*a^10*b^8*c^2*n^4-1408*a^10*b^6*c^4*n^4+4352*a^10*b^4*c^6*n^4-5120*a^10*b^2*c^8*n^4+2048*a^10*c^10*n^4+64*a^8*b^10*c^2*n^4-4352*a^8*b^8*c^4*n^4+
15936*a^8*b^6*c^6*n^4-29408*a^8*b^4*c^8*n^4+29440*a^8*b^2*c^10*n^4-11520*a^8*c^12*n^4-1408*a^6*b^10*c^4*n^4+15936*a^6*b^8*c^6*n^4-33024*a^6*b^6*c^8*n^4+38336*a^6*b^
4*c^10*n^4-37440*a^6*b^2*c^12*n^4+17280*a^6*c^14*n^4+4352*a^4*b^10*c^6*n^4-29408*a^4*b^8*c^8*n^4+38336*a^4*b^6*c^10*n^4-10368*a^4*b^4*c^12*n^4+5184*a^4*b^2*c^14*n^4
-7776*a^4*c^16*n^4-5120*a^2*b^10*c^8*n^4+29440*a^2*b^8*c^10*n^4-37440*a^2*b^6*c^12*n^4+5184*a^2*b^4*c^14*n^4+7776*a^2*b^2*c^16*n^4+2048*b^10*c^10*n^4-11520*b^8*c^12
*n^4+17280*b^6*c^14*n^4-7776*b^4*c^16*n^4+224*a^10*b^10*n^3+400*a^10*b^8*c^2*n^3-10072*a^10*b^6*c^4*n^3+32048*a^10*b^4*c^6*n^3-38336*a^10*b^2*c^8*n^3+15488*a^10*c^
10*n^3+400*a^8*b^10*c^2*n^3-31280*a^8*b^8*c^4*n^3+115104*a^8*b^6*c^6*n^3-215912*a^8*b^4*c^8*n^3+220096*a^8*b^2*c^10*n^3-87168*a^8*c^12*n^3-10072*a^6*b^10*c^4*n^3+
115104*a^6*b^8*c^6*n^3-235344*a^6*b^6*c^8*n^3+274544*a^6*b^4*c^10*n^3-277608*a^6*b^2*c^12*n^3+130896*a^6*c^14*n^3+32048*a^4*b^10*c^6*n^3-215912*a^4*b^8*c^8*n^3+
274544*a^4*b^6*c^10*n^3-62928*a^4*b^4*c^12*n^3+33696*a^4*b^2*c^14*n^3-58968*a^4*c^16*n^3-38336*a^2*b^10*c^8*n^3+220096*a^2*b^8*c^10*n^3-277608*a^2*b^6*c^12*n^3+
33696*a^2*b^4*c^14*n^3+60912*a^2*b^2*c^16*n^3+15488*b^10*c^10*n^3-87168*b^8*c^12*n^3+130896*b^6*c^14*n^3-58968*b^4*c^16*n^3+544*a^10*b^10*n^2+872*a^10*b^8*c^2*n^2-\
25204*a^10*b^6*c^4*n^2+83236*a^10*b^4*c^6*n^2-101824*a^10*b^2*c^8*n^2+41696*a^10*c^10*n^2+872*a^8*b^10*c^2*n^2-78408*a^8*b^8*c^4*n^2+288996*a^8*b^6*c^6*n^2-555340*a
^8*b^4*c^8*n^2+581784*a^8*b^2*c^10*n^2-234504*a^8*c^12*n^2-25204*a^6*b^10*c^4*n^2+288996*a^6*b^8*c^6*n^2-574872*a^6*b^6*c^8*n^2+675008*a^6*b^4*c^10*n^2-723084*a^6*b
^2*c^12*n^2+352356*a^6*c^14*n^2+83236*a^4*b^10*c^6*n^2-555340*a^4*b^8*c^8*n^2+675008*a^4*b^6*c^10*n^2-104808*a^4*b^4*c^12*n^2+67572*a^4*b^2*c^14*n^2-158868*a^4*c^16
*n^2-101824*a^2*b^10*c^8*n^2+581784*a^2*b^8*c^10*n^2-723084*a^2*b^6*c^12*n^2+67572*a^2*b^4*c^14*n^2+172152*a^2*b^2*c^16*n^2+41696*b^10*c^10*n^2-234504*b^8*c^12*n^2+
352356*b^6*c^14*n^2-158868*b^4*c^16*n^2+520*a^10*b^10*n+764*a^10*b^8*c^2*n-25074*a^10*b^6*c^4*n+87616*a^10*b^4*c^6*n-111022*a^10*b^2*c^8*n+46428*a^10*c^10*n+764*a^8
*b^10*c^2*n-77924*a^8*b^8*c^4*n+285516*a^8*b^6*c^6*n-572230*a^8*b^4*c^8*n+628432*a^8*b^2*c^10*n-260718*a^8*c^12*n-25074*a^6*b^10*c^4*n+285516*a^6*b^8*c^6*n-533616*a
^6*b^6*c^8*n+633480*a^6*b^4*c^10*n-759822*a^6*b^2*c^12*n+391836*a^6*c^14*n+87616*a^4*b^10*c^6*n-572230*a^4*b^8*c^8*n+633480*a^4*b^6*c^10*n+3552*a^4*b^4*c^12*n+32040
*a^4*b^2*c^14*n-176778*a^4*c^16*n-111022*a^2*b^10*c^8*n+628432*a^2*b^8*c^10*n-759822*a^2*b^6*c^12*n+32040*a^2*b^4*c^14*n+206532*a^2*b^2*c^16*n+46428*b^10*c^10*n-\
260718*b^8*c^12*n+391836*b^6*c^14*n-176778*b^4*c^16*n+150*a^10*b^10+210*a^10*b^8*c^2-7407*a^10*b^6*c^4+29133*a^10*b^4*c^6-39627*a^10*b^2*c^8+17253*a^10*c^10+210*a^8
*b^10*c^2-22998*a^8*b^8*c^4+81501*a^8*b^6*c^6-180831*a^8*b^4*c^8+220317*a^8*b^2*c^10-96759*a^8*c^12-7407*a^6*b^10*c^4+81501*a^6*b^8*c^6-123444*a^6*b^6*c^8+151380*a^
6*b^4*c^10-250341*a^6*b^2*c^12+145431*a^6*c^14+29133*a^4*b^10*c^6-180831*a^4*b^8*c^8+151380*a^4*b^6*c^10+88140*a^4*b^4*c^12-19305*a^4*b^2*c^14-65637*a^4*c^16-39627*
a^2*b^10*c^8+220317*a^2*b^8*c^10-250341*a^2*b^6*c^12-19305*a^2*b^4*c^14+87516*a^2*b^2*c^16+17253*b^10*c^10-96759*b^8*c^12+145431*b^6*c^14-65637*b^4*c^16)*N^2-8*c^2*
(2*n+5)*(2*n+3)*(n+3)*(b-c)*(b+c)*(a-c)*(a+c)*(8*a^6*b^6*n^2-48*a^6*b^4*c^2*n^2+96*a^6*b^2*c^4*n^2-64*a^6*c^6*n^2-48*a^4*b^6*c^2*n^2+48*a^4*b^4*c^4*n^2-48*a^4*b^2*c
^6*n^2+72*a^4*c^8*n^2+96*a^2*b^6*c^4*n^2-48*a^2*b^4*c^6*n^2-72*a^2*b^2*c^8*n^2-64*b^6*c^6*n^2+72*b^4*c^8*n^2+16*a^6*b^6*n-108*a^6*b^4*c^2*n+234*a^6*b^2*c^4*n-164*a^
6*c^6*n-108*a^4*b^6*c^2*n+84*a^4*b^4*c^4*n-96*a^4*b^2*c^6*n+186*a^4*c^8*n+234*a^2*b^6*c^4*n-96*a^2*b^4*c^6*n-204*a^2*b^2*c^8*n-164*b^6*c^6*n+186*b^4*c^8*n+6*a^6*b^6
-42*a^6*b^4*c^2+111*a^6*b^2*c^4-87*a^6*c^6-42*a^4*b^6*c^2-6*a^4*b^4*c^4-15*a^4*b^2*c^6+99*a^4*c^8+111*a^2*b^6*c^4-15*a^2*b^4*c^6-132*a^2*b^2*c^8-87*b^6*c^6+99*b^4*c
^8)*N^3:

end:




#MakeNice(ope,n,N): inputs a recurrence operator ope in n, N, outputs Ope1(n,1/N) such that
#1-Ope1(n,N) is equivalent fo Ope. Try:
#MakeNice(N^2-n*N-1,n,N);
MakeNice:=proc(ope,n,N) local k,i,ope1:
k:=degree(ope,N): 

ope1:=ope/coeff(ope,N,k):

ope1:=1-expand(subs(n=n-k,ope1)/N^k):
add(factor(coeff(ope1,N,-i))*N^(-i),i=1..3):

end:

#OpeNice(a,b,c,n,N): the operator in 1/N, let's call it Ope1,  such that Ope(a,b,c,n,N) is equivalent fo 1-Ope1. Try:
#OpeNice(5,7,35);
OpeNice:=proc(a,b,c,n,N) option remember: MakeNice(Ope(a,b,c,n,N),n,N):end:



#Sabc1f(a,b,c,n1): A much faster version of Sabc1(a,b,c,n) using the recurrence. Try:
#[seq(Sabc1f(5,7,35,i),i=0..20)];
Sabc1f:=proc(a,b,c,n1) local ope,N,n,i:
option remember:
ope:=OpeNice(a,b,c,n,N):

if 0<=n1 and n1<=2 then
 RETURN(Sabc1(a,b,c,n1)):
else
add(subs(n=n1,coeff(ope,N,-i))*Sabc1f(a,b,c,n1-i),i=1..3):
fi:

end:


#Sabc2f(a,b,c,n1): A much faster version of Sabc2(a,b,c,n) using the recurrence. Try:
#[seq(Sabc2f(5,7,35,i),i=0..20)];
Sabc2f:=proc(a,b,c,n1) local ope,N,n,i:
option remember:
ope:=OpeNice(a,b,c,n,N):

if 0<=n1 and n1<=2 then
 RETURN(Sabc2(a,b,c,n1)):
else
add(subs(n=n1,coeff(ope,N,-i))*Sabc2f(a,b,c,n1-i),i=1..3):
fi:

end:

#deltE(f,X): Given f a linear combination of 1, X[1],.., X[nops(X)] and 1
#with integer coefficients, finds
#the empirical delta such that abs(f)=1/max(coeff(f,1),coeff(f,X[1]),coeff(f,X[2]))^delta. Try:
#deltE(Sabc1f(5,7,35),[log(2),log(3)]);
#deltE(Sabc2f(5,7,35),[log(2),log(3)]);

deltE:=proc(f,X) local f1,i,L,E1,gad,gu:

if not (type(X,list) and nops(X)>0) then
 print(`X must be a non-emptylist`):
 RETURN(FAIL):
fi:

f1:=numer(f):

 gu:=f1:

 gad:=abs(coeff(gu,X[1],1)):

 gu:=coeff(f1,X[1],0):

 for i from 2 to nops(X) do
 gad:=max(abs(coeff(gu,X[i],1)),gad):
  gu:=coeff(gu,X[i],0):
 od:

if not type(gu,numeric) then
 RETURN(FAIL):
fi:

L:=max(abs(gu),gad):

E1:=evalf(f1):

evalf(-log(abs(E1))/log(L)):

end:



#dn(n): lcm(1...n). Try: dn(10);
dn:=proc(n) local i:lcm(seq(i,i=1..n)):end:

#deltT1(a,b,c): The theoretical delta for Salikhov's Diophantine approximations inspired by Sabc1(a,b,c,n). Try:
#deltT1(5,7,35);

deltT1:=proc(a,b,c) local K,gu,N,n,ka,ga,i,ku1,ku2,kak,kak1,kak2,lu:

K:=GuessK1(a,b,c,200):

if K=FAIL then
 RETURN(FAIL):
fi:

K:=K[1]:

gu:=Shors(a,b,c,N)[2]:


kak:=RaE1(a,b,c,100):

kak1:=abs(gu[1]-kak):
kak2:=abs(gu[2]-kak):

if kak1<kak2 then
 lu:=gu[1]:
else
 lu:=gu[2]:
fi:



lu:=-(log(lu)+log(K)+2)/(log(gu[3])+log(K)+2):

evalf(lu,20):

end:


#deltT2(a,b,c): The theoretical delta for Salikhov's Diophantine approximations assuming that multiplying by K^n*dn(2*n) clears denominators
#deltT2(5,7,35);

deltT2:=proc(a,b,c) local K,gu,N,n,ka,ga,i,ku1,ku2,kak,kak1,kak2,lu:

K:=GuessK2(a,b,c,200):

if K=FAIL then
 RETURN(FAIL):
fi:

K:=K[1]:

gu:=Shors(a,b,c,N)[2]:


kak:=RaE2(a,b,c,100):

kak1:=abs(gu[1]-kak):
kak2:=abs(gu[2]-kak):

if kak1<kak2 then
 lu:=gu[1]:
else
 lu:=gu[2]:
fi:



lu:=-(log(lu)+log(K)+2)/(log(gu[3])+log(K)+2):

evalf(lu,20):

end:



#deltT(a,b,c): The theoretical delta for Salikhov's Diophantine approximations inspired by Sabc1(a,b,c,n). Try:
#deltT1(5,7,35);

deltT:=proc(a,b,c) local K,K1,K2,gu,N,n,ka,ga,i,ku1,ku2,kak,kak1,kak2,lu:

K1:=GuessK1(a,b,c,200):

if K1=FAIL then
 RETURN(FAIL):
fi:



K2:=GuessK2(a,b,c,200):

if K2=FAIL then
 RETURN(FAIL):
fi:


K:=lcm(K1[1],K2[1]):

gu:=Shors(a,b,c,N)[2]:


lu:=min(
-(log(gu[1])+log(K)+2)/(log(gu[3])+log(K)+2),
-(log(gu[2])+log(K)+2)/(log(gu[3])+log(K)+2)):

evalf(lu,20):

end:


#GuessK1(a,b,c,N): guesses an integer K and a small integer L such that Sabc1(a,b,c,n)*dn(2*n)*K^n*L is an integer for all n
#between 20 and N. The output is [K,L].
#Try:
#GuessK1(5,7,35,100);
GuessK1:=proc(a,b,c,N) local mu,gu,i,K,n1:

if N<=20 then
 print(N, `is too small`):
 RETURN(FAIL):
fi:

gu:=ifactor(denom(Sabc1f(a,b,c,50)*dn(100))):




if type(gu, `*`) then
mu:=1:
for i from 1 to nops(gu) do
if type(op(i,gu),`^`) then

  if op(2,op(i,gu))>=5 then
   mu:=mu*op(1,op(1,op(i,gu)))^round(op(2,op(i,gu))/50):
  fi:
fi:
od:

elif type(gu,`^`) then
   mu:=op(1,op(1,gu))^round(op(2,gu)/50):
else
mu:=1:
fi:



K:=mu:


gu:=lcm(seq(denom(Sabc1f(a,b,c,n1)*dn(2*n1)*K^n1),n1=20..N)):


if {seq(denom(gu*Sabc1f(a,b,c,n1)*dn(2*n1)*K^n1),n1=20..N)}<>{1} then
 RETURN(FAIL):
fi:

[K,gu]:

end:



#GuessK2(a,b,c,N): guesses an integer K and a small integer L such that Sabc2(a,b,c,n)*dn(2*n)*K^n*L is an integer for all n
#between 20 and N. The output is [K,L].
#Try:
#GuessK2(5,7,35,100);
GuessK2:=proc(a,b,c,N) local mu,gu,i,K,n1:

if N<=20 then
 print(N, `is too small`):
 RETURN(FAIL):
fi:

gu:=ifactor(denom(Sabc2f(a,b,c,50)*dn(100))):




if type(gu, `*`) then

mu:=1:

for i from 1 to nops(gu) do
if type(op(i,gu),`^`) then

  if op(2,op(i,gu))>=5 then
   mu:=mu*op(1,op(1,op(i,gu)))^round(op(2,op(i,gu))/50):
  fi:
fi:
od:

elif type(gu,`^`) then
   mu:=op(1,op(1,gu))^round(op(2,gu)/50):
else
mu:=1:
fi:



K:=mu:


gu:=lcm(seq(denom(Sabc2f(a,b,c,n1)*dn(2*n1)*K^n1),n1=20..N)):


if {seq(denom(gu*Sabc2f(a,b,c,n1)*dn(2*n1)*K^n1),n1=20..N)}<>{1} then
 RETURN(FAIL):
fi:

[K,gu]:

end:


#ExtractLogs1(L): Given an expression of the form log(k) or c*log(k) where k is an integer, outputs k.
#Tty:
#ExtractLog1(log(11));  ExtractLog1(7*log(11)); 
ExtractLogs1:=proc(L) local i:

if not ( type(L, function) or type(L,`*`)) then
 RETURN(FAIL):
fi:

if type(L,function) then
 if not op(0,L)=ln then
    RETURN(FAIL):
 else
  RETURN(op(1,L)):
 fi:
fi:


for i from 1 to nops(L) do
 if type(op(i,L),function) then
   RETURN(ExtractLogs1(op(i,L))):
 fi:
od:

FAIL:

end:


#ExtractLogs(L): Given an expression that is a linear combination of log(k)'s  and 1, outputs the set of k's that show up.
#Try:
#ExtractLogs(5*log(3)-11*log(5));
ExtractLogs:=proc(L) local i,gu,mu:

if not type(L,`+`) then
 mu:=ExtractLogs1(L):
  if mu<>FAIL then
   gu:={mu}:
  fi:


else
 gu:={}:

 for i from 1 to nops(L) do
  if not type(op(i,L),numeric) then
   mu:=ExtractLogs1(op(i,L)):
  if mu=FAIL then
    RETURN(FAIL):
  else
  gu:=gu union {mu}:
 fi:
 fi:
od:
fi:

gu:

mu:=L-add(coeff(L,log(gu[i]),1)*log(gu[i]),i=1..nops(gu)):

if not (type(mu,fraction) or type(mu,integer)) then
 RETURN(FAIL):
fi:
gu:


end:

 
#Atoms1(a,b,c): The atoms of Sabc1(a,b,c,n) for all n. Try:
#Atoms1(5,7,35);
Atoms1:=proc(a,b,c) local gu,mu,i:
option remember:
 gu:=ExtractLogs(Sabc1(a,b,c,0)):

 if gu=FAIL then
   RETURN(FAIL):
 fi:

for i from 1 to 2 do
  mu:=ExtractLogs(Sabc1(a,b,c,i)):

 
 if mu=FAIL then
   RETURN(FAIL):
 else
 gu:=gu union mu:
 fi:
od:
gu:
end:


 
#Atoms2(a,b,c): The atoms of Sabc2(a,b,c,n) for all n. Try:
#Atoms2(5,7,35);
Atoms2:=proc(a,b,c) local gu,mu,i:
option remember:
 gu:=ExtractLogs(Sabc2(a,b,c,0)):

 if gu=FAIL then
   RETURN(FAIL):
 fi:

for i from 1 to 2 do
  mu:=ExtractLogs(Sabc2(a,b,c,i)):

 
 if mu=FAIL then
   RETURN(FAIL):
 else
 gu:=gu union mu:
 fi:
od:
gu:
end:



#RaE1(a,b,c,N):abs(Sabc1f(a,b,c,N)/Sabc1f(a,b,c,N-1)):
RaE1:=proc(a,b,c,N) local gu:
gu:=evalf(abs(Sabc1f(a,b,c,N)/Sabc1f(a,b,c,N-1))):
evalf(gu,20):
end:


#RaE2(a,b,c,N):abs(Sabc2f(a,b,c,N)/Sabc2f(a,b,c,N-1)):
RaE2:=proc(a,b,c,N) local gu:
gu:=evalf(abs(Sabc2f(a,b,c,N)/Sabc2f(a,b,c,N-1))):
evalf(gu,20):
end:


#CoeffsL(L,X): Given f a linear combination of 1, X[1],.., X[nops(X)] and 1, outputs the respective coefficients
#Try(CoeffsL(Sabc1(5,7,35,100),[log(2),log(3)]);
Coeffs:=proc(L,X) local i:

if not (type(X,list) and nops(X)>0) then
 print(`X must be a non-emptylist`):
 RETURN(FAIL):
fi:

[seq(coeff(L,X[i],1),i=1..nops(X)), L-add(coeff(L,X[i],1)*X[i],i=1..nops(X))]:

end:


#Ra1(a,b,c,N): The rations of the sequence of coefficients of the atoms of Sabc1(a,b,c) with N. Try:
#Ra1(5,7,35,100); Ra1(5,17,85,100);
Ra1:=proc(a,b,c,N) local gu,mu1,mu2,lu,i:
gu:=Atoms1(a,b,c):
if gu=FAIL then
  RETURN(FAIL):
fi:

gu:=[seq(log(gu[i]),i=1..nops(gu))];

mu1:=Coeffs(Sabc1f(a,b,c,N),gu):
mu2:=Coeffs(Sabc1f(a,b,c,N-1),gu):

lu:=[seq(mu1[i]/mu2[i],i=1..nops(mu1))]:
evalf(lu,20):


end:



#Shors(a,b,c,N): The indicial polynomial and three roots of the limiting indicial equation of Ope(a,b,c,n,N). Try:
#Shors(5,7,35,N);
Shors:=proc(a,b,c,N) local gu,n,i,c1,mu:

gu:=Ope(a,b,c,n,N):

gu:=lcoeff(gu,n):

c1:=igcd(seq(coeff(gu,N,i),i=0..degree(gu,N))):
gu:=sort(normal(gu/c1)):



mu:=evalf([solve(gu,N)]):

mu:=evalf(mu,20):

[gu,sort([seq(abs(coeff(mu[i],I,0)),i=1..nops(mu))])]:



end:


#GrC(a,b,c): The growth constants of Sabc1(a,b,c,n) and Sabc2(a,b,c,n)
#It also gives the locations (or rather their squares)
GrC:=proc(a,b,c) local t,g,gu,mu,i,gu1,gu2,gu3:
g:=t*(t-a^2)*(t-b^2)/(t-c^2)^2:

gu:=[solve(numer(diff(g,t)))]:


gu:=evalf(gu):

gu:=[seq(coeff(gu[i],I,0),i=1..nops(gu))]:



for i from 1 to nops(gu) do
 if gu[i]>0 and gu[i]<a^2 then
    gu1:=gu[i]:
  elif gu[i]>0 and gu[i]>a^2 and gu[i]<b^2 then
     gu2:=gu[i]:
  elif  gu[i]>b^2 then
  gu3:=gu[i]:
 fi:
od:

mu:=[abs(subs(t=gu1,g)),abs(subs(t=gu2,g)),abs(subs(t=gu3,g))]:
mu:=evalf(mu,20):

[[evalf(gu1,20),evalf(gu2,20),evalf(gu3,20)],mu]:
end:


#Info1(a,b,c): inputs positive integers a,b,c, and outputs a list consisting of
#(i) the set of integers such that Sabc1(a,b,c,n) are linear combinations (with rational number coefficients) of their logs
#(ii) The conjectured integer K such that K^n*dn(2*n)*Sabc1(a,b,c,n) have integer coefficients
#(iii) The theoretical delta such that
#T(a,b,c,n):=Sabc1(a,b,c,n)*K^n*dn(2*n), has the property, that for some universal constant CONST
#|T(a,b,c,n)|<=CONST/(max(abs(coeff.s))^delta .
#(iv): The smallest empirical delta between 100 and 150. Try:
#Info1(5,7,35):
Info1:=proc(a,b,c) local gu1,gu2,gu3,gu4,i,X,n1:
gu1:=Atoms1(a,b,c):

if gu1=FAIL then
  RETURN(FAIL):
fi:

gu2:=GuessK1(a,b,c,200):

if gu2=FAIL then
  RETURN(FAIL):
fi:

gu2:=gu2[1]:

gu3:=deltT1(a,b,c):

X:=[seq(log(gu1[i]),i=1..nops(gu1))]:


gu4:=min(seq(deltE(Sabc1f(a,b,c,n1),X),n1=100..150)): 
gu4:=evalf(gu4,20):


[gu1,gu2,gu3,gu4]:

end:


#Info2(a,b,c): inputs positive integers a,b,c, and outputs a list consisting of
#(i) the set of integers such that Sabc2(a,b,c,n) are linear combinations (with rational number coefficients) of their logs
#(ii) The conjectured integer K such that K^n*dn(2*n)*Sabc1(a,b,c,n) have integer coefficients
#(iii) The theoretical delta such that
#T(a,b,c,n):=Sabc2(a,b,c,n)*K^n*dn(2*n), has the property, that for some universal constant CONST
#|T(a,b,c,n)|<=CONST/(max(abs(coeff.s))^delta .
#(iv): The smallest empirical delta between 100 and 150. Try:
#Info2(5,7,35):
Info2:=proc(a,b,c) local gu1,gu2,gu3,gu4,i,X,n1:
gu1:=Atoms2(a,b,c):

if gu1=FAIL then
  RETURN(FAIL):
fi:

gu2:=GuessK2(a,b,c,200):

if gu2=FAIL then
  RETURN(FAIL):
fi:

gu2:=gu2[1]:

gu3:=deltT2(a,b,c):

X:=[seq(log(gu1[i]),i=1..nops(gu1))]:


gu4:=min(seq(deltE(Sabc2f(a,b,c,n1),X),n1=100..150)): 
gu4:=evalf(gu4,20):


[gu1,gu2,gu3,gu4]:

end:

#SalPairsSlow(N): all the pairs of primes [p,q] p<q, up to the N-th prime such that Atoms1(p,q,p*q)=Atoms2(p,q,p*q)
#Try:
#SalPairsSlow(10);
SalPairsSlow:=proc(N) local i,j,p,q,gu:

gu:=[]:

for i from 1 to N do
p:=ithprime(i):
 for j from i+1 to N do
   q:=ithprime(j):

    if Atoms1(p,q,p*q)=Atoms2(p,q,p*q) then
     gu:=[op(gu),[[p,q], Atoms1(p,q,p*q)]]:
    fi:
 od:
od:

gu:
end:




#SalPairs(N): all the pairs of primes [p,q] p<q, up to the N-th prime such that Atoms1(p,q,p*q)=Atoms2(p,q,p*q)
#Try:
#SalPairs(10);
SalPairs:=proc(N) local i,j,p,q,gu:

gu:=[]:

for i from 1 to N do
p:=ithprime(i):
 for j from i+1 to N do
   q:=ithprime(j):

    if ExtractLogs(Sabc1(p,q,p*q,0))=ExtractLogs(Sabc2(p,q,p*q,0)) then
     gu:=[op(gu),[[p,q], ExtractLogs(Sabc1(p,q,p*q,0)) ]]:
    fi:
 od:
od:

gu:
end:





SalPairsPC:=proc():

[[[5, 7], {2, 3}], [[5, 17], {2, 3}], [[7, 17], {2, 3}], [[11, 19], {2, 3, 5}],
[[11, 31], {2, 3, 5}], [[13, 97], {2, 3, 7}], [[13, 127], {2, 3, 7}], [[19, 31]
, {2, 3, 5}], [[29, 41], {2, 3, 5, 7}], [[29, 71], {2, 3, 5, 7}], [[29, 251], {
2, 3, 5, 7}], [[29, 449], {2, 3, 5, 7}], [[41, 71], {2, 3, 5, 7}], [[41, 251],
{2, 3, 5, 7}], [[41, 449], {2, 3, 5, 7}], [[43, 197], {2, 3, 7, 11}], [[71, 251
], {2, 3, 5, 7}], [[71, 449], {2, 3, 5, 7}], [[89, 109], {2, 3, 5, 11}], [[89,
199], {2, 3, 5, 11}], [[89, 241], {2, 3, 5, 11}], [[97, 127], {2, 3, 7}], [[101
, 271], {2, 3, 5, 17}], [[109, 199], {2, 3, 5, 11}], [[109, 241], {2, 3, 5, 11}
], [[151, 191], {2, 3, 5, 19}], [[181, 701], {2, 3, 5, 7, 13}], [[199, 241], {2
, 3, 5, 11}], [[251, 449], {2, 3, 5, 7}], [[269, 401], {2, 3, 5, 67}], [[349, 
811], {2, 3, 5, 7, 29}], [[769, 881], {2, 3, 5, 7, 11}]]
end:



#Zug1(a,b,n1): The  pair of rational number [A,B] such that Sabc1(2*a+1,2*b+1,(2*a+1)*(2*b+1),n1)=A+B*log(b/(b+1))
#Try:
#Zug1(2,3,10);
Zug1:=proc(a,b,n1) local gu,A,B,lu,i:
gu:=Sabc1f(2*a+1,2*b+1,(2*a+1)*(2*b+1),n1):
if n1=0 then
 lu:=gu/(log(b)-log(b+1)):
 lu:=simplify(lu):
  if not type(lu,fraction) then
   RETURN(FAIL):
  else
   RETURN([0,lu]):
  fi:
fi:

if not  type(gu,`+`) then
 RETURN(FAIL):
fi:

for i from 1 to nops(gu) do
 if type(op(i,gu),fraction) then
   A:=op(i,gu):
   break:
 fi:
od:


if i=nops(gu)+1 then
  RETURN(FAIL):
fi:

B:=simplify((gu-A)/(log(b)-log(b+1))):
if not type(B,fraction) then
 RETURN(FAIL):
fi:
[A,B]:

end:


#Zug2((a,b,n1): The  pair of rational number [A,B] such that Sabc2((2*a+1),(2*b+1),(2*a+1)*(2*b+1),n1)
Zug2:=proc(a,b,n1) local gu,A,B,lu,i:
gu:=Sabc2f(2*a+1,2*b+1,(2*a+1)*(2*b+1),n1):
if n1=0 then
 lu:=gu/(log(a)-log(a+1)):
 lu:=simplify(lu):
  if not type(lu,fraction) then
   RETURN(FAIL):
  else
   RETURN([0,lu]):
  fi:
fi:


if not  type(gu,`+`) then
 RETURN(FAIL):
fi:

for i from 1 to nops(gu) do
 if type(op(i,gu),fraction) then
   A:=op(i,gu):
   break:
 fi:
od:

if i=nops(gu)+1 then
  RETURN(FAIL):
fi:


B:=simplify((gu-A)/log((a)/(a+1))):
if not type(B,fraction) then
 RETURN(FAIL):
fi:
[A,B]:

end:


#INFOabNoGood(a,b,x): inputs two  integers 1<=a<b and outputs the information about a Salikhov type theorem 
#of simultaneous approximations to log(b/(b+1)) and log(a/(a+1)). Try:
#It outputs 
#[b/(b+1), a/(a+1)], [K1,K2,delta1, delta2, nu], where K1 is such that
#Zug1(a,b,n)*lcm(1..2*n)*K1^n and Zug2(a,b,n)*lcm(1..2*n)*K2^n are integers, delta1 is
#delta1 is such that
#|A(n)+B1(n)*log((b-1)/(b+1))|<=CONST/max(A(n),B1(n))^delta1
#|A(n)+B2(n)*log((a-1)/(a+1))|<=CONST/max(A(n),B1(n))^delta2
#and nu is 1/min(delta1,delta2) the implied mu such that the corollary is
#the theorem that if Q=max(|q|,|p1|,|p2|), Q>=Q0 is sufficiently large, then
#|q+p1*log((b-1)/(b+1))+p2*log((a-1)/(a+1)))|>1/Q^nu
#It also outputs empirical delta1 and delta2. Try
#INFOabNoGood(2,3,x);

INFOabNoGood:=proc(a,b,x) local gu1,gu2:
if not (a>0 and b>0 and type(a,integer) and type(b,integer)) then
 print(a,b, `must be positive integers larger than 1`):
 RETURN(FAIL):
fi:

if not a<b then
 print(a, `must be less than `, b):
 RETURN(FAIL):
fi:

gu1:=Info1(2*a+1,2*b+1,(2*a+1)*(2*b+1)):

if gu1=FAIL then
 RETURN(FAIL):
fi:

gu2:=Info2(2*a+1,2*b+1,(2*a+1)*(2*b+1)):

if gu2=FAIL then
 RETURN(FAIL):
fi:

[[b/(b+1),a/(a+1)],[gu1[2],gu2[2],gu1[3],gu2[3], evalf(1/min(gu1[3],gu2[3]),20)],[gu1[4],gu2[4]], Shors(2*a+1,2*b+1,(2*a+1)*(2*b+1),x)]:
end:


#INFOab(a,b,x): inputs two  integers 1<=a<b and outputs the information about a Salikhov type theorem 
#of simultaneous approximations to log(b/(b+1)) and log(a/(a+1)). Try:
#It outputs 
#[b/(b+1), a/(a+1)], [K1,K2,nu], where K1 is such that
#Zug1(a,b,n)*lcm(1..2*n)*K1^n and Zug2(a,b,n)*lcm(1..2*n)*K2^n are integers, and nu is the upper bound for the measure of independence
#the theorem that if Q=max(|q|,|p1|,|p2|), Q>=Q0 is sufficiently large, then
#|q+p1*log((b-1)/(b+1))+p2*log((a-1)/(a+1)))|>1/Q^nu
#It also outputs empirical delta1 and delta2. Try
#INFOab(2,3,x);

INFOab:=proc(a,b,x) local gu1,gu2,de:
if not (a>0 and b>0 and type(a,integer) and type(b,integer)) then
 print(a,b, `must be positive integers larger than 1`):
 RETURN(FAIL):
fi:

if not a<b then
 print(a, `must be less than `, b):
 RETURN(FAIL):
fi:

de:=deltT(2*a+1,2*b+1,(2*a+1)*(2*b+1)):



gu1:=Info1(2*a+1,2*b+1,(2*a+1)*(2*b+1)):

if gu1=FAIL then
 RETURN(FAIL):
fi:

gu2:=Info2(2*a+1,2*b+1,(2*a+1)*(2*b+1)):

if gu2=FAIL then
 RETURN(FAIL):
fi:

[[b/(b+1),a/(a+1)],[gu1[2],gu2[2],evalf(1/de,20)],[gu1[4],gu2[4]], Shors(2*a+1,2*b+1,(2*a+1)*(2*b+1),x)]:
end:

#Paper(M): States  and sketches the proof of simultaneous linear diophantine approximation
#for log((b-1)/(b+1)) and log((a-1)/(a+1)) for all 3<=a<b<M a and b odd
Paper:=proc(M) local x,n,ope,i,N,gu,a,b,theta1,theta2,K1,K2,nu,pol,t0:
t0:=time():
print(`Using V. Kh. Salikhov's method to determine upper bounds for `):
print(` the linear independence measure of 1 ,log(a/(a+1)), log(b/(b+1)) for many pairs of  integers 1<=a<b<`, M ):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`In a beautiful paper, "On the Irrationality measure of log(3)",  that appeared in Doklady 2007, v. 417, no. 6, pp. 753-755`):
print(`the great number-theorist, V. Kh. Salikhov, proved the following theorem`):
print(``):
print(`Theorem (Salikhov): Suppose that q,p1,p2 are integers and Q=max(|p|,|q1|,|q2|) Q>=Q0, where Q0 is a sufficiently large number`):
print(`then |q+p1*log(2)+p2*log(3)|>1/Q^4.125 `):
print(``):
print(`Let us make the following definition`):
print(``):
print(`Definition: Given two real constants, theta1, and theta2, nu is a LINEAR INDEPENDENCE MEASURE of 1, theta1, theta2, if the following holds.`):
print(` Suppose that q,p1,p2 are integers and Q=max(|p|,|q1|,|q2|) Q>=Q0, where Q0 is a sufficiently large number`):
print(`then |q+p1*theta1+p2*theta2|>1/Q^nu `):
print(`and nu is the smallest such number. `):
print(``):
print(`Using this definition, Salikhov's theorem can be phrased as`):
print(``):
print(` nu(log(2),log(3))<=4.125 `):
print(``):
print(`In fact, Salikhov proved the equivalent statement `):
print(``):
print(` nu(log(2/3),log(3/4))<=4.125 `):
print(``):
print(`In this article we will prove many such results for the pairs of irrational numbers [log(a/(a+1)),log(b/(b+1))]`):
print(`For 2<=a<b<=`, N):
print(``):
print(`A note about proofs:`):

print(`Salikhov introduced the integrals`):
print(``):
print(E1(n)=Int(SInt(5,7,35,n,x),x=0..5)):
print(``):
print(E2(n)=Int(SInt(5,7,35,n,x),x=0..7)):
print(``):
print(``):
print(`and then introduced THREE sequences of rational numbers A1(n), A2(n), B(n), defined by `):
print(E1(n)=A1(n)+B(n)*log(3/4)):
print(``):
print(E2(n)=A2(n)+B(n)*log(2/3)):
print(``):
print(`He then used calculus and saddle point to find constants , let's call them C1, C2, C3, such that`):
print(``):
print(`|A1(n)|,|A2(n)|,|B(n)| are all OMEGA(C1^n)`):
print(``):
print(`E1(n)= OMEGA(1/C2^n)`):
print(``):
print(`E2(n)= OMEGA(1/C3^n)`):
print(``):
print(`We will use instead the fact that E1(n), E2(n) both satisfy a certain linear third order linear recurrence with polynomial coefficients`):
print(`obtained by the Almkvist-Zeilberger algorithm`):
print(``):
  ope:=Ope(5,7,35,n,N):
print(``):
print(add(coeff(ope,N,i)*E1(n+i),i=0..degree(ope,N))=0):
print(``):
print(`subject to the initial conditions`):
print(``):
print(seq(E1(i)=Sabc1(5,7,35,i),i=0..2)):
print(``):
print(`and E2(n) satisfies the same recurrence`):
print(``):
print(add(coeff(ope,N,i)*E2(n+i),i=0..degree(ope,N))=0):
print(``):
print(`but with initial conditions`):
print(``):
print(seq(E2(i)=Sabc2(5,7,35,i),i=0..2)):
print(``):
print(`Using the Poincare lemma, we get that up to polynomial corrections that disappear at the wash, C1,C2,C3 can be gotten from`):
print(` the roots of the polynomial in N (the inidicial polynomial of the constant-coefficients recurrence that is the leading term`):
print(`of the above recurrence`):
gu:=Shors(5,7,35,N)[2]:
print(gu):
print(`So far, A1(n), A2(n), B(n) are sequences of rational numbers`):
print(``):
print(`It was easy to prove that`):
print(``):
print(`B(n)*lcm(1..2*n) are all integers (or rational numbers with bounded denominators)`):
print(``):
print(`Salikhov also also found integers K1, K2, such that`):
print(``):
print(`A1(n)*K1^n*lcm(1..2*n) are all integers (or rational numbers with bounded denominators)`):
print(``):
print(`A2(n)*K2^n*lcm(1..2*n) are all integers (or rational numbers with bounded denominators)`):
print(``):
print(`in his case it so happened`):
print(`that K1=2, and K2=1`):
print(``):
print(`Since we need coefficients of the two constants log(3/4) and log(2/3) to be the SAME, and we need to make all the sequences INTEGERS`):

print(`Let K=lcm(K1,K2), in this case K=2`):
print(``):
print(`Defining A1a(n)=A1(n)*lcm(1..2n)*K^n,  Ba(n)=B(n)*lcm(1..2n)*K^n`): 
print(``):
print(`It follows that A1a(n),Ba(n) are two sequences of INTEGER sequences such that `):
print(``):
print(`|A1a(n)+Ba(n)*log(2/3)|=OMEGA(K^n*lcm(1..2n)*E1(n))`):
print(``):
print(`Since lcm(1..n)=OMEGA(e^n), and hence lcm(1..2n)=OMEGA(e^(2n)) it follows that`):
print(``):
print(`|A1a(n)+Ba(n)*log(2/3)|=OMEGA(1/max(|A1a(n)|, |Ba(n)|)^delta1 `):
print(``):
print(`where delta1=-(log(C2)+log(K)+2)/(log(C1])+log(K)+2) `):
print(``):

print(` Similarly `):
print(``):
print(`Defining A2b(n)=A2(n)*lcm(1..2n)*K^n,  Ba(n)=B(n)*lcm(1..2n)*K^n`): 
print(`note that `):
print(``):
print(`It follows that A2b(n),(and the same) Ba(n) are two sequences of INTEGER sequences such that `):
print(``):
print(`|A2b(n)+Bb(n)*log(3/4)|=K^n*lcm(1..2n)*E2(n)`):
print(``):
print(`Since lcm(1..n)=OMEGA(e^n), and hence lcm(1..2n)=OMEGA(e^(2n)) it follows that`):
print(``):
print(`|A2b(n)+Bb(n)*log(3/4)|=OMEGA(1/max(|A2b(n)|, |Ba(n)|)^delta2) `):
print(``):
print(`where delta2=-(log(C3)+log(K)+2)/(log(C1)+log(K)+2) `):
print(``):
gu:=INFOab(2,3,x):
print(``):
print(`Finally, using Lemma 3 in Salikhov's paper, that is due to M. Hata (Acta Arith. 63(4), 335-349 (1993)),`):
print(` the desired nu(log(2/3),log(3/4)) is >= then 1/min(delta1,delta2)=`, gu[2][3]):
print(``):
print(`Now comes the GENERALIZATION`):
print(``):
print(`Let a and b be positive integers larger than 1 and such that a<b`):
print(` Define `):
print(``):
print(E1(n)=Int(SInt(2*a+1,2*b+1,(2*a+1)*(2*b+1),n,x),x=0..a)):
print(``):
print(E2(n)=Int(SInt(2*a+1,2*b+1,(2*a+1)*(2*b+1),n,x),x=0..b)):
print(``):

print(`and let	introduce THREE sequences of rational numbers A1(n), A2(n), B(n), defined by `):
print(E1(n)=A1(n)+B(n)*log(b/(b+1))):
print(``):
print(E2(n)=A2(n)+B(n)*log(a/(a+1))):
print(``):
print(`Everything goes verbatim, we get find the third-order operator annihilating them, and find C1,C2,C3`):
print(``):
print(`the only thing that changes each time are the integers K1 and K2 that we confess we found empirically, but we are sure`):
print(`can be left as exercises to the readers`):
print(`then we make everything into integer sequences by mutiplying by lcm(1..2n)*lcm(K1,K2)^n`):
print(``):
print(`For the sake of brevity, we will only list a,b, theta1=log(b/(b+1)), theta2=log(a/(a+1)), K1, K2 (so the reader can confirm it)`):
print(`followed by nu`):
print(`and the polynomial in x whose roots give you C1, C2, C3`):
print(`This is enough for the reader to verify the claim. `):

for a from 2 to M do
 for b from a+1 to M do


 
gu:=INFOab(a,b,x):
 
  if gu<>FAIL and gu[2][3]>0 then
     print(`a=`, a, `b=`, b):
 print(theta1=log(b/(b+1)),theta2=log(a/(a+1)), K1=gu[2][1], K2=gu[2][2], nu=gu[2][3], pol=gu[4][1]):
 fi:
od:
od:
print(``):
print(`--------------------------------------------------------------`):
print(``):
print(`This concludes this paper, that took`, time()-t0, `seconds to generate. `):
print(``):
end:



############Start infinite families
#NuEven(a,N): a closed form expression for the exponent nu in the linear independce of {1,log(a/(a+1)), log((a+1)/(a+2))}. Try:
#when a is even.
#It also returns the cubic in N whose largest and smallest roots feature. Try:
#NuEven(a,N);
NuEven:=proc(a,N) local ope,n,gad,kat,si,hal,gu,i,K,lu:
ope:=Ope(4*a+1,4*a+3,(4*a+1)*(4*a+3),n,N):

ope:=lcoeff(ope,n):

gu:=[solve(ope,N)]:


gad:=gu[1]:

si:=coeff(evalf(subs(a=3,gad)),I,0):

for i from 2 to nops(gu) do

 hal:=coeff(evalf(subs(a=3,gu[i])),I,0):

if hal>si then
 gad:=gu[i]:
fi:
od:


kat:=gu[1]:

si:=coeff(evalf(subs(a=3,gad)),I,0):

for i from 2 to nops(gu) do

 hal:=coeff(evalf(subs(a=3,gu[i])),I,0):

if hal<si then
 kat:=gu[i]:
fi:
od:


K:=a*(a+1):

lu:=-(log(kat)+log(K)+2)/(log(gad)+log(K)+2):
lu:=1/lu:
lu:=simplify(lu):

ope:=add(factor(coeff(ope,N,i))*N^i,i=0..degree(ope,N)):


lu:=simplify(subs(a=a/2,lu)):

lu,ope:

end:




#NuOdd(a,N): a closed form expression for the exponent nu in the linear independence of {1,log((a)/(a+1)), log(( a+1)/(+2))}. 
#when a is odd and >=3
#It also returns the cubic in N whose largest and smallest roots feature. Try:
#NuOdd (a,N);
NuOdd:=proc(a,N) local ope,n,gad,kat,si,hal,gu,i,K,lu:
ope:=Ope(4*a+3,4*a+5,(4*a+3)*(4*a+5),n,N):

ope:=lcoeff(ope,n):

gu:=[solve(ope,N)]:


gad:=gu[1]:

si:=coeff(evalf(subs(a=3,gad)),I,0):

for i from 2 to nops(gu) do

 hal:=coeff(evalf(subs(a=3,gu[i])),I,0):

if hal>si then
 gad:=gu[i]:
fi:
od:


kat:=gu[1]:

si:=coeff(evalf(subs(a=3,gad)),I,0):

for i from 2 to nops(gu) do

 hal:=coeff(evalf(subs(a=3,gu[i])),I,0):

if hal<si then
 kat:=gu[i]:
fi:
od:


K:=(2*a+1)*(2*a+3):

lu:=-(log(kat)+log(K)+2)/(log(gad)+log(K)+2):
lu:=1/lu:
lu:=simplify(lu):

ope:=add(factor(coeff(ope,N,i))*N^i,i=0..degree(ope,N)):

lu:=simplify(subs(a=(a-1)/2,lu)):

lu,ope:

end:




#PaperG(a,N):inputs symbols a and N and outputs an article with a proof 
#(modulo divisibility lemmas left the reader) of the linear independence 
#of {1,log(a/(a+1)),log((a+1)/(a+2))}, together with explicit upper bounds
#for the measure of independence (different for a even and odd). Try:
#PaperG(a,N);
PaperG:=proc(a,N) local gue,guo,A1,A2,B,n,x,ope,i,lu:

print(`Proof of the Linear Independence of {1, log(a/(a+1)) and log((a+1)/(a+2))} for Every Positive Integer a`):
print(`and Corresponding Explicit Linear Independence Measures`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):

gue:=NuEven(a,N)[1]:

guo:=NuOdd(a,N)[1]:


print(`Theorem : For any positive integer a, the three numbers {1, log(a/(a+1)), log((a+1)/(a+2))} are  linearly independent`):
print(`In other words there exists a positive number, nu, such that `):
print(` if q,p1,p2 are integers and Q=max(|q|,|p1|,|p2|) Q>=Q0, where Q0 is a sufficiently large number`):
print(`then |q+p1*log(a/(a+1))+p2*log((a+1)/(a+2))|>1/Q^nu .`):
print(``):
print(`Furthermore, an upper bound for nu is as follows`):
print(``):
print(`if a is even then nu can be taken as`):
print(``):
print(gue):
print(``):
print(`and in Maple format `):
print(``):
lprint(gue):
print(``):
print(`On the other hand, if  a is odd and >=3 then nu can be taken as`):
print(``):
print(guo):
print(``):
print(`and in Maple format `):
print(``):
lprint(guo):
print(``):
print(`When a=1, then nu can be taken as`):
print(``):
lprint(INFOab(1,2,x)[2][3]):
print(``):

print(`--------------------------------------`):
print(``):
print(`Proof: Consider the two integrals`):
print(``):
print(E1(n)=Int(SInt(2*a+1,2*a+3,(2*a+1)*(2*a+3),n,x),x=0..2*a+1)):
print(``):
print(E2(n)=Int(SInt(2*a+1,2*a+3,(2*a+1)*(2*a+3),n,x),x=0..2*a+3 ) ):
print(``):
print(`It is readily seen that we can write`):
print(``):
print(E1(n)=A1(n)+B(n)*log((a+1)/(a+2))):
print(``):
print(E2(n)=A2(n)+B(n)*log((a/(a+1)))):
print(``):
print(`For THREE sequences of RATIONAL numbers A1(n), A2(n), B(n) `):
print(``):
print(`Note the crucial fact that the coefficient of log((a+1)/(a+2)) in E1(n) and the coefficient of log(a/(a+1)) in E2(n) are identical`):
print(``):
print(`but of course, A1(n), A2(n) are different. `):
print(``):
ope:=Ope(2*a+1,2*a+3,(2*a+1)*(2*a+3),n,N):
print(``):
print(`Using the amazing Almkvist-Zeilberger algorithm (using a different Maple package called EKHAD.txt), it emerges as that`):
print(`both E1(n) and E2(n), and hence also the sequences of rational numbers A1(n), A2(n), B(n) satisfy the third order linear recurrence`):
print(`equation`):
print(``):
print(add(factor(coeff(ope,N,i))*F(n+i),i=0..degree(ope,N))=0):
print(``):
print(`but of course with different initial conditions. `):
print(``):
print(`Using the Poincare lemma, we can use it to deduce the growth of E1(n), E2(n), and the common growth of A1(n),A2(n), B(n). `):
print(``):
print(`We approximate the above linear recurrence equation by a CONSTANT coefficients linear recurrence obtained by taking the leading term in n`):
print(``):
print(`The indicial polynomial, in N,  (after normalizing the constant term to be 1) `):
print(``):
ope:=lcoeff(ope,n):
ope:=normal(ope/coeff(ope,N,0)):
print(``):
ope:=add(factor(coeff(ope,N,i))*N^i,i=0..degree(ope,N)):
print(``):
print(ope):
print(``):
print(`and in Maple format`):
print(``):
lprint(ope):
print(``):
print(`Let C1(a) be the largest root, C2(a), the second largest root, and C3(a) the smallest rule of the above equation`):
print(``):
print(`that can be expressed explicitly thanks to Cardano`):
print(``):
print(`It follows from the Poincare lemma that, ignoring polynomial corrections that disappear at the end`):
print(``):
print(`A1(n),A2(n),B(n)=OMEGA(C1(a)^n)`):
print(``):
print(`E1(n)=OMEGA(C2(a)^n), E2(n)=OMEGA(C3(a)^n)`):
print(` or the other way`):
print(``):
print(`E2(n)=OMEGA(C3(a)^n), E1(n)=OMEGA(C2(a)^n)`):
print(``):
print(`Note that C1(a) is big and C2(a), C3(a) are small`):
print(``):
print(`We now need three divisibility lemmas that we leave to the reader`):
print(``):
print(`Lemma 1: B(n)*lcm(1..2*n) is always an integer `):
print(``):
print(`Lemma 2: A1(n)*K1(a)^n*lcm(1..2*n) is always an integer (or at worst has bounded denominator) `):
print(``):
print(`where K1(a)= a+2 if a is odd, and K1(a)=a/2+1 if a is even`):
print(``):
print(`Lemma 3: A2(n)*K2(a)^n*lcm(1..2*n) is always an integer (or at worst has bounded denominator) `):
print(``):
print(`where K2(a)= a if a is odd, and K2(a)=a/2 if a is even`):
print(``):
print(`Let K(a)=K1(a)*K2(a) `):
print(``):
print(`In other words, K(a)=a(a+2) if a is odd, and K(a)=(a/2)(a/2+1) if a is even `):
print(``):
print(`Multiplying the equations`):
print(``):
print( E1(n)=A1(n)+B(n)*log((a+1)/(a+2))):
print(``):
print(E2(n)=A2(n)+B(n)*log((a/(a+1)))):
print(``):
print(` by K(a)^n*lcm(1..2*n), and letting E1new(n)=K(a)^n*lcm(1..2*n)*E1(n),  E2new(n)=K(a)^n*lcm(1..2*n)*E2(n) `):
print(``):
print(` We have `):
print(`E1new(n)=A1new(n)+Bnew(n)*log((a+1)/(a+2))`):
print(``):
print(`E2new(n)=A2new(n)+Bnew(n)*log((a/(a+1)))`):
print(``):
print(`where A1new(n), A2new(n), Bnew(n) are INTEGER sequences `):
print(``):
print(`Recall the fact that lcm(1..n) is asymptotic to e^n and hence lcm(1..2*n) is asymptotic to e^(2n) .`):
print(`We now that A1new(n), A2new(n), Bnew(n) are, up to polynomial factors that disappear at the wash `):
print(``):
print(`A1new(n),A2new(n),Bnew(n)=OMEGA((C1(a)*e^2*K(a))^n )`):
print(``):
print(`E1new(n), E2new(n) =OMEGA((C3(a)*e^2*K(a))^n `):
print(``):
print(`Finally, using Lemma 3 in Salikhov's paper, that is due to M. Hata (Acta Arith. 63(4), 335-349 (1993)),`):
print(``):
print(`We get that one can take nu to be`):
print(``):
lu:=-(log(C2(a))+log(K(a))+2)/(log(C1(a))+log(K(a))+2):
print(``):
lu:=1/lu:
print(``):
print(lu):
print(``):
print(`plugging for C1(a), C3(a) their values given by the Cardano formulas and K(a)=a(a+2) if a is odd, and K(a)=(a/2)(a/2+1) if a is even `):
print(``):
print(`we get the expressions in the statement of the theorem. QED. `):
print(``):
print(`Finally let's plot the independce measure in terms of a, from a=1 to a=1000. The higher curve is for the odd case`):
print(``):
print(plot([gue,guo], a=1..1000));
end: