From 508c4d9c94c0610ea3a9c2f5a5d2d9aed9924df2 Mon Sep 17 00:00:00 2001
From: Alexander Hulpke <hulpke@math.colostate.edu>
Date: Thu, 20 Aug 2020 13:33:35 -0600
Subject: [PATCH 1/2] ENHANCE: Added Taylor series for rational functions

---
 doc/ref/ratfun.xml |  2 ++
 lib/ratfun.gd      | 21 ++++++++++++++++++++
 lib/ratfunul.gi    | 49 +++++++++++++++++++++++++++++++++++++++++++---
 3 files changed, 69 insertions(+), 3 deletions(-)

diff --git a/doc/ref/ratfun.xml b/doc/ref/ratfun.xml
index 900ab74d6a..5d3ff47831 100644
--- a/doc/ref/ratfun.xml
+++ b/doc/ref/ratfun.xml
@@ -430,6 +430,8 @@ extended versions for Laurent polynomials exist).
 
 <#Include Label="UnivariateRationalFunctionByCoefficients">
 
+<#Include Label="TaylorSeriesRationalFunction">
+
 </Section>
 
 
diff --git a/lib/ratfun.gd b/lib/ratfun.gd
index 5103a0bca3..316829f938 100644
--- a/lib/ratfun.gd
+++ b/lib/ratfun.gd
@@ -1340,6 +1340,27 @@ DeclareOperation("DegreeIndeterminate",[IsPolynomial,IsPosInt]);
 DeclareAttribute("Derivative",IsUnivariateRationalFunction);
 DeclareOperation("Derivative",[IsPolynomialFunction,IsPosInt]);
 
+#############################################################################
+##
+#A  TaylorSeriesRationalFunction( <ratfun>, <at>, <deg> )
+##
+##  <#GAPDoc Label="TaylorSeriesRationalFunction">
+##  <ManSection>
+##  <Attr Name="TaylorSeriesRationalFunction" Arg='ratfun, at, deg]'/>
+##
+##  <Description>
+##  Computes the taylor series up to degree <A>deg</A> of <A>ratfun</A> at
+##  <A>at</A>.
+##  <Example><![CDATA[
+##  gap> TaylorSeriesRationalFunction((x^5+3*x+7)/(x^5+x+1),0,11);
+##  -50*x^11+36*x^10-26*x^9+22*x^8-18*x^7+14*x^6-10*x^5+4*x^4-4*x^3+4*x^2-4*x+7
+##  ]]></Example>
+##  </Description>
+##  </ManSection>
+##  <#/GAPDoc>
+##
+DeclareGlobalFunction("TaylorSeriesRationalFunction");
+
 
 #############################################################################
 ##
diff --git a/lib/ratfunul.gi b/lib/ratfunul.gi
index cd4c033922..7efc505f33 100644
--- a/lib/ratfunul.gi
+++ b/lib/ratfunul.gi
@@ -1069,18 +1069,61 @@ end);
 RedispatchOnCondition(Derivative,true,
   [IsPolynomial],[IsLaurentPolynomial],0);
 
-InstallOtherMethod(Derivative,"uratfun,ind",true,
+InstallOtherMethod(Derivative,"uratfun",true,
   [IsUnivariateRationalFunction],0,
 function(ratfun)
-local num,den;
+local num,den,R,cf,pow,dr,a;
+  # try to be good for the case of iterated derivatives by using if the
+  # denominator is a power
   num:=NumeratorOfRationalFunction(ratfun);
   den:=DenominatorOfRationalFunction(ratfun);
-  return (Derivative(num)*den-num*Derivative(den))/(den^2);
+  R:=CoefficientsRing(DefaultRing([den]));
+  cf:=Collected(Factors(den));
+  pow:=Gcd(List(cf,x->x[2]));
+  if pow>1 then
+    dr:=Product(List(cf,x->x[1]^(x[2]/pow))); # root
+  else
+    dr:=den;
+  fi;
+  #cf:=(Derivative(num)*dr-num*pow*Derivative(dr))/dr^(pow+1);
+  num:=Derivative(num)*dr-num*pow*Derivative(dr);
+  den:=dr^(pow+1);
+  if IsOne(Gcd(num,dr)) then
+    # avoid cancellation
+    num:=CoefficientsOfUnivariateLaurentPolynomial(num);
+    den:=CoefficientsOfUnivariateLaurentPolynomial(den);
+    cf:=UnivariateRationalFunctionByExtRepNC(FamilyObj(ratfun),
+      num[1],den[1],num[2]-den[2],
+        IndeterminateNumberOfUnivariateRationalFunction(ratfun));
+    # maintain factorization of denominator
+    StoreFactorsPol(R,den,ListWithIdenticalEntries(pow+1,dr));
+  else
+    cf:=num/den; # cancellation
+  fi;
+  return cf;
+  #return (Derivative(num)*den-num*Derivative(den))/(den^2);
 end);
 
 RedispatchOnCondition(Derivative,true,
   [IsPolynomial],[IsUnivariateRationalFunction],0);
 
+
+InstallGlobalFunction(TaylorSeriesRationalFunction,function(f,at,deg)
+local t,i,x;
+  if not (IsUnivariateRationalFunction(f) and deg>=0) then
+    Error("function is not univariate pol, or degree negative");
+  fi;
+  x:=IndeterminateOfUnivariateRationalFunction(f);
+  t:=Zero(f);
+  for i in [0..deg] do
+    if i>0 then 
+      f:=Derivative(f);
+    fi;
+    t:=t+Value(f,at)*(x-at)^i/Factorial(i);
+  od;
+  return t;
+end);
+
 #############################################################################
 ##
 #F  Discriminant( <f> ) . . . . . . . . . . . . discriminant of polynomial f

From 8519e207fe5454a224f7d37ca0fd203c0efa3085 Mon Sep 17 00:00:00 2001
From: Alexander Hulpke <hulpke@math.colostate.edu>
Date: Fri, 4 Sep 2020 08:37:44 -0600
Subject: [PATCH 2/2] Fixed the suggested simplification (list/pol confusion)

---
 lib/ratfunul.gi | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/lib/ratfunul.gi b/lib/ratfunul.gi
index 7efc505f33..80d31a46c4 100644
--- a/lib/ratfunul.gi
+++ b/lib/ratfunul.gi
@@ -1096,7 +1096,8 @@ local num,den,R,cf,pow,dr,a;
       num[1],den[1],num[2]-den[2],
         IndeterminateNumberOfUnivariateRationalFunction(ratfun));
     # maintain factorization of denominator
-    StoreFactorsPol(R,den,ListWithIdenticalEntries(pow+1,dr));
+    StoreFactorsPol(R,DenominatorOfRationalFunction(cf),
+      ListWithIdenticalEntries(pow+1,dr));
   else
     cf:=num/den; # cancellation
   fi;