Q:=RationalField();
P<x>:=PolynomialRing(Q);
(x^2+1)^3*(x+2)^2;

// Integration de la partie polynomiale.
function intPoly(g,f)
  A:=Parent(g);
  X:=A.1;
//  q:=g div f;
//  r:=g mod f;
  q,r:=Quotrem(g,f);
//  s:=A!0;
//  if q ne 0 then
//    s:=&+([Coefficient(q,i)/(i+1)*X^(i+1):i in [0..Degree(q)]]);
//  end if;
//  return s,r,f;
  return Integral(q),r,f;
end function;
// s est la partie polynomiale integree ; reste a integrer r/f.

intPoly(x^2+1,x+2);

// Integration par la methode de Hermite.
function intHermite(r,f)
  A:=Parent(f);
  S:=SquarefreePartialFractionDecomposition(r/f);
  // Pour accumuler les fractions rationnelles toutes integrees.
  R:=[];
  // Pour accumuler la partie logarithmique.
  L:=[];
  // On pourra utiliser la decomposition de bas en haut.
  i:=#S;
  // Pour eviter un bug.
  S cat:=[<1,1,0>];
  // Variable auxilliaire pour reperer quand on change de g.
  g2:=1;
  while i gt 0 do
    if S[i][2] eq 1 then
      if S[i][3] ne 0 then
        L cat:=[S[i]];
      end if;
    else
      g:=S[i][1];
      j:=S[i][2];
      p:=S[i][3];
      if g ne g2 then
        g2:=g;
        dg:=Derivative(g);
        un,s,t:=XGCD(dg,g); // un = s*dg + t*g
      end if;
      sp:=s*p;
      tp:=t*p;
//      s2:=sp mod g;
//      q:=sp div g;
      q,s2:=Quotrem(sp,g);
      t2:=tp+q*dg;
      incr:=t2+Derivative(s2)/(j-1);
      if i eq 1 or S[i-1][2] ne j-1 then S[i]:=<g,j-1,A!0>; i+:=1; end if;
      S[i-1][3]+:=incr;
      R cat:=[-s2/(j-1)/g^(j-1)];
    end if;
    i-:=1;
  end while;
  return R,[<l[1],l[3]>:l in L];
end function;

intHermite(x,(x^2+1)^2);

// Integration d'une partie logarithmique a/b.
function intRothsteinTrager(a,b)
  db:=Derivative(b);
  A:=Parent(b); // Q[X]
  B:=PolynomialRing(Parent(Coefficient(b,0)),2); // Q[X,Y]
  h:=hom<A->B|B.1>; // Q[X] -> Q[X,Y]
  proj:=hom<B->A|A!0,A.1>; // Q[X,Y] -> Q[Y] (avec Y note X).
  R:=Resultant(h(a)-B.2*h(db),h(b),B.1);
  R:=proj(R/LeadingCoefficient(R));
  F:=Factorization(R);
  if &or([Degree(f[1]) ne 1:f in F]) then
    return <"partie logarithmique non-integrable sans extension",F>;
  end if;
  return [<c,GCD(b,a-c*db)> where c is Roots(f[1])[1][1]:f in F];
end function;

intRothsteinTrager(P!1,x^2+1);
intRothsteinTrager(P!1,x^3+x);

// Integration d'une fraction rationnelle generale.
function integration(g,f)
  s,r,f:=intPoly(g,f);
  R,L:=intHermite(r,f);
  return s,R,[intRothsteinTrager(l[2],l[1]):l in L];
end function;

integration(P!1,(2*x^2+1)^2);

Q:=RationalField();
pol<A>:=PolynomialRing(Q);
F<a>:=NumberField(A^2+2);
P<x>:=PolynomialRing(F);
integration(P!1,(2*x^2+1)^2);