(* Now we start with a quadratic denominator *)

f[a_,m_]:=f[a,m]=Integrate[1/(x^2+2*a*x+1)^(m+1),{x,0,Infinity}];

(* The parameters are a and m *)

(* The first example is *)

f[a,0]

(* Mathematica assumes that a is a complex parameter *)

(* The command "Assumptions" will let Mathematica know that a is a real parameter *)

(* We take it between 0 and 1 to make the roots of the denominator not real *)

g[a_,m_]:=g[a,m]=Integrate[1/(x^2+2*a*x+1)^(m+1),{x,0,Infinity},Assumptions->(0<a < 1)];

g[a,0]

(* Observe that both numerator and denominator are NOT real numbers *)

FullSimplify[g[a,0]]

(* This does not help *)

g[a,1]

g[a,2]

FullSimplify[g[a,2]]

(* Now we introduce the command that will transform Sqrt[1/(-1+a^2)] into 1/Sqrt[-1+a^2] *)

g1[a_,m_]:=g1[a,m]=g[a,m]/.{Sqrt[1/(-1+a^2)]->1/Sqrt[-1+a^2]};

g[a,0]

g1[a,0]

(* These two are the same *)

g[a,1]

g1[a,1]

(* Here you see that g1[a,1] is simpler than g[a,1] *)

g[a,2]

g1[a,2]

(* Now compute the first four *)

(k:=1)&&While[k≤4,(Print[k,"   ",g1[a,k]];k++)]

(* And now the first twenty (why not, it is free *)

(k:=1)&&While[k≤20,(Print[k,"   ",g1[a,k]];k++)]

(* Note that the form appears clear. First of all, for m at least 14, we get the function Sqrt[1-a^2]*ArcTan[Sqrt[1-a^2]/a]. I will call that b *)

g2[a_,m_]:=g1[a,m]/.{Sqrt[1-a^2]*ArcTan[Sqrt[1-a^2]/a]->b};

(* Let's see what happens *)

g1[a,20]

g2[a,20]

(* The ArcTan function has disappeared and now IT LOOKS SIMPLER *)

(* Now we see that the denominators have a specific form *)

g3[a_,m_]:=g3[a,m]=(1-a^2)^(m+1)*g2[a,m];

g2[a,20]

g3[a,20]

(* It does not automatically cancel. We try this *)

g4[a_,m_]:=g4[a,m]=(-1+a^2)^(m+1)*g2[a,m];

g4[a,20]

(* Now expand *)

Expand[g4[a,20]]

(* We see that g4[a,m] is of the form P[a,m] + C[m]*b *)

g4[a,20]

(* This command gives the coefficient of b *)

Coefficient[g4[a,20],b]

Coefficient[g4[a,30],b]

(* Now we will try to guess the forms of these coefficients *)

g4[a,20]-Coefficient[g4[a,20],b]*b//Together

(* This is the polynomial in a *)

g5[a_,m_]:=g4[a,m]-Coefficient[g4[a,m],b]*b//Together;

g5[a,14]

g5[a,15]

g5[a,16]

den[a_,m_]:=Denominator[g5[a,m]];

Table[den[a,m],{m,14,20}]

(* Nowwe compute quotients of consecutive denominators *)

den[a,15]/den[a,14]

den[a,16]/den[a,15]

den[a,17]/den[a,16]

(* Sometimes they are integers, sometimes not *)

quo[m_]:=den[a,m+1]/den[a,m];

Table[quo[m],{m,14,20}]

Table[quo[m],{m,14,40}]

(* We will try to explain this *)