Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add TaylorSeriesRationalFunction and improve Derivative for univariate rational functions #4104

Merged
merged 2 commits into from
Sep 4, 2020
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions doc/ref/ratfun.xml
Original file line number Diff line number Diff line change
Expand Up @@ -430,6 +430,8 @@ extended versions for Laurent polynomials exist).

<#Include Label="UnivariateRationalFunctionByCoefficients">

<#Include Label="TaylorSeriesRationalFunction">

</Section>


Expand Down
21 changes: 21 additions & 0 deletions lib/ratfun.gd
Original file line number Diff line number Diff line change
Expand Up @@ -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");


#############################################################################
##
Expand Down
50 changes: 47 additions & 3 deletions lib/ratfunul.gi
Original file line number Diff line number Diff line change
Expand Up @@ -1069,18 +1069,62 @@ 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,DenominatorOfRationalFunction(cf),
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;
hulpke marked this conversation as resolved.
Show resolved Hide resolved
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;
Comment on lines +1118 to +1124
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Assuming we require deg >= 0 (which is not clear to me from the description of the function), could also do this

Suggested change
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;
t:=Value(f,at)*x^0;
for i in [1..deg] do
f:=Derivative(f);
t:=t+Value(f,at)*(x-at)^i/Factorial(i);
od;

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this is a wash in terms of cost or code readability, so I left it out. (Added test for deg\ge0)

return t;
end);

#############################################################################
##
#F Discriminant( <f> ) . . . . . . . . . . . . discriminant of polynomial f
Expand Down