# Symbolic Computation of 
Formal Power Series in Maxima

February 2021

Bertrand Teguia Tabuguia

This worksheet demonstrates some capacities of a new package FPS (Formal Power Series) that can be used to recover power series formulas and sometimes discover unexpected ones. The package implements symbolic algorithms from

Teguia Tabuguia, Bertrand and Koepf, Wolfram. Symbolic computation of hypergeometric type and non-holonomic power series. arXiv preprint arXiv:2102.04157. 2021 -- This paper was divided and rewritten into two parts for better descriptions of the methods:

Teguia Tabuguia, Bertrand and Koepf, Wolfram. Symbolic conversion of holonomic functions to hypergeometric type power series. To appear in the Computer Algebra issue of the journal of Programming and Computer Software.

Teguia Tabuguia, Bertrand and Koepf, Wolfram. On the representation of non-holonomic power series. To be presented at the Maple Conference 2021.

Teguia Tabuguia, Bertrand. Power Series Representations of Hypergeometric Type and Non-Holonomic Functions in Computer Algebra. Ph.D. thesis, University of Kassel, Germany. https://kobra.uni-kassel.de/handle/123456789/11598. May 2020. 

The results go far beyond what Maxima's powerseries command, which uses a pattern matching approach, currently does.

batchload(FPS);

## Overview of the method

Consider for example f(z):=arctan(z)+cos(z). 
First, the code computes a differential equation of a special type, called holonomic, satisfied by f(z). This is done in the following way.

f:atan(z)+cos(z);

DE:HolonomicDE(f,F(z));

Second, the differential equation is converted into a recurrence equation.

RE:DEtoRE(DE,F(z),a[n]);

This can also be done directly.

FindRE(f,z,a[n]);

Then the coefficient of the power series sought is computed as (m-fold hypergeometric term) solutions of the obtained RE. The algorithm used is new.

mfoldHyper(RE,a[n]);

These are the coefficients of the even part of the expansion (There are only 2-fold hypergeometric term solutions over the rationals). For the odd part, mfoldHyper can be used as follows.

mfoldHyper(RE,a[n],2,1);

Finally a linear combination of the involved power series is sought with verification of the initial values, and we get

FPS(f,z,n);

At the current time, this complete approach for hypergeometric type power series is new for all computer algebra systems. The same computations can be done for series with negative and fractional powers.

FPS(sin(z)/z^8,z,n);

FPS(cos(sqrt(z))+sin(z^(1/3)),z,n);

Apart from series related to hypergeometric type series, there is another approach used wich mainly deals with expressions that are not holonomic like log(1+sin(z)), (1+tan(z))/(1-tan(z)), etc. We find it a little inconvenient to introduce Bernoulli or Euler numbers in the expansion since it could not give full information to the user. Instead we return recursive formulas that the user can use to recover truncated series.

f:log(1+sin(z));

DE:QDE(f,F(z));

FindQRE(f,z,a[n]);

FPS(f,z,n);

FPS((1+tan(z))/(1-tan(z)),z,n);

QNF(tan(z),z,n);

The formulas for expressions like tan(z) and 1/(exp(z)-1) are easily related to hypergeometric type series. Therefore FPS yields

FPS(tan(z),z,n);

FPS(1/(exp(z)-1),z,n);

Note that the coefficient of the latter generates what is called Bernoulli numbers.

## Provided commands

The reader should be able to guess the meaning of the symbols used in the syntax below. The main commands provided by the package are:
1. FPS(f,z,n): to compute power series. By default the point of expansion is zero. A different point of expansion can be given as fourth argument (FPS(f,z,n,z0)).
2. mfoldHyper(RE,a[n]): to find m-fold hypergeometric term solutions of holonomic recurrence equations, for all positive integer m. A third argument C (mfoldHyper(RE,a[n],C)) can be spefified for computations involving algebraic numbers (Complex numbers, irrational numbers, parameters).
3. HolonomicDE(f,F(z)): to compute a holonomic differential equation satisfied by a given expression.
4. QDE(f,F(z)): to compute a quadratic differential equation satisfied by a given expression.
5. NLDE(f,F(z)): to compute an algebraic differential equation satisfied by a given expression. Sometimes this is used by FPS for non-holonomic expressions.
6. FindRE(f,z,a[n]): to compute a holonomic recurrence equation satisfied by a given expression
7. FindQRE(f,z,a[n]): to compute a non-holonomic recurrence equation (from QDE) satisfied by a given expression.
8. Taylor(f,z,z0,N): to compute truncated taylor expansions of holonomic functions. This is usefull sometimes since it is faster for computing larger-order Taylor expansions.
9. QTaylor(f,z,z0,N): to compute truncated series for non holonomic expressions. However this is not efficient in general since the recurrence equations used are not linear.
10. QNF(f,z,n): direct approach to compute power series of non-holonomic expressions. This is also used by FPS(f,z,n), but can be used to see a different result for tan(z) for example.
11. CompatibleDE(DE1,DE2,F(z)): to check the compatibility of two differential equations. This can be used for holonomic expressions in general.

Some variables of the package: 
1. Nmax with default value 10, is the maximum order of differential equations computed using HolonomicDE and internally in FPS.
2. QNmax with default value 21, is the analogue for non-holonomic differential equations. Note that the differential operator used in the code does not have same order with diff(f,z,d). The corresponding function is delta2diff(f,z,d).
3. NLmax with default value 30, is the analogue for the general algebraic differential equation case. The derivative operator is deltadiff(f,z,d).
4. NLDEflag with default value false, can be useful for some demonstration with some non-holonomic functions. However, this will not be presented in this worksheet.

## Classical examples

FPS(exp(z),z,n);

FPS(log(1+z),z,n);

FPS(1/(1-z),z,n);

FPS(cos(z),z,n);

FPS(atan(z),z,n);

FPS(atanh(z),z,n);

FPS(cosh(z),z,n);

FPS(sinh(z),z,n);

FPS(cos(z),z,n);

FPS(sin(z),z,n);

FPS(sin(z)+cos(z),z,n);

FPS(asin(z),z,n);

FPS(asinh(z),z,n);

declare(p,constant)$

FPS(1/(1-z)^p,z,n);

## More thoughtfull series

FPS(erf(z),z,n);

FPS(asin(z)^2/z^2,z,n);

FPS(exp(asinh(z)),z,n);

FPS(exp(asin(z)),z,n);

FPS(sqrt((1-sqrt(1-z))/z),z,n);

FPS((z+sqrt(1+z^2))^p,z,n);

FPS(exp(z)-2*exp(-z/2)*cos(sqrt(3)*z/2+%pi/3),z,n);

FPS(1/2*log((1+z)/(1-z))-atan(z),z,n);

FPS(log(1+z+z^2+z^3),z,n);

FPS(cos(4*acos(z)),z,n);

FPS(asin(z)+acos(z),z,n);

FPS(atan((z+p)/(1-z*p)),z,n);

FPS(atan((z+p)/(1-z*p))-atan(z),z,n);

FPS(sin(3*acos(z))/sqrt(1-z^2),z,n);

FPS(1/(sqrt(1-4*z))*((1-sqrt(1-4*z))/2*z)^2,z,n);

FPS(sqrt(sqrt(8*z^3+1)-1),z,n);

FPS(1/(z+z^2),z,n);

FPS(cos(z)/z^1000,z,n);

FPS(log(1+z)+cos(z),z,n);

FPS(cos(z)+z*sin(z),z,n);

FPS(sin(z)+z*cos(z),z,n);

FPS(log(1+sqrt(z)+z+z^(3/2)),z,n);

FPS(atan(z)+exp(z^2),z,n);

FPS(exp(z)*cos(z),z,n);

FPS(exp(z)*sin(z),z,n);

FPS((cos(z)+sin(z))^2,z,n);

FPS((cos(z)+sin(z))^3,z,n);

FPS(asech(z),z,n);

FPS(exp(z)+log(1+z),z,n,%e);

FPS(sin(z)+exp(z),z,n,%pi);

FPS(atan(z),z,n,inf);

FPS(exp(1/z),z,n,inf);

FPS(z*exp(-z)*expintegral_ei(z),z,n,inf);

FPS(sqrt(%pi)*exp(z)*(1-erf(sqrt(z))),z,n,inf);

FPS(exp(-1/z)*expintegral_ei(1/z)/z,z,n);

FPS(acos(z^(1/2))+exp(z^2),z,n);

FPS(exp(asinh(z^2))+1/(1-z^(2/3)),z,n);

FPS(z*log(z)^2+asin(z),z,n);

FPS(log(1+sqrt(z))+atan(z^(1/3)),z,n);

FPS(sin(2*z)+cos(z),z,n,%pi/2);

FPS(asin(z)*cos(z),z,n);

FPS(atan(z)*log(1+z^2),z,n);

FPS(log(1+z)^2,z,n);

FPS(exp(z^(3/2))*asin(z^(3/4)),z,n);

FPS(log(1+z)^5,z,n);

FPS(atan(z)^2,z,n);

FPS(exp(z+z^2),z,n);

FPS(exp(z^(3/2)+sqrt(z)),z,n);

FPS(sec(sqrt(z)),z,n);

FPS(tanh(z),z,n);

FPS(exp(z)*log(1+z),z,n);

FPS(sinh(log(1+z)),z,n);

FPS(cosh(log(1+z)),z,n);

FPS(cos(log(1+z)),z,n);

FPS(sin(log(1+z)),z,n);

FPS(acos(z)^2,z,n);

FPS(asin(z)^2,z,n);

FPS(acos(z)^2-asin(z)^2,z,n);

FPS(sin(z)+sinh(z),z,n);

FPS(cos(2*z)+sin(3*z),z,n);

FPS(sin(z)*cos(z),z,n);

FPS(sinh(z)*sin(z),z,n);

FPS(cosh(z)*cos(z),z,n);

FPS(exp(z)*cosh(z),z,n);

FPS(atanh(z)+sqrt(1+z)+1/sqrt(1+z),z,n);

FPS(sinh(z)^2+cosh(z)^2,z,n);

FPS(-cosh(z)^2+sinh(z)^2,z,n);

declare(q,constant)$

FPS(1/(p-z^2)/(q-z^3),z,n);

FPS(exp(2*atanh(sin(2*z)/(1+cos(2*z)))),z,n);

FPS((1+tan(z))/(1-tan(z))-exp(2*atanh(sin(2*z)/(1+cos(2*z)))),z,n);