5.10 Series expansion

Expressions know how to expand themselves as a Taylor series or (more generally) a Laurent series. As in most conventional Computer Algebra Systems, no distinction is made between those two. There is a class of its own for storing such series (class pseries) and a built-in function (called Order) for storing the order term of the series. As a consequence, if you want to work with series, i.e. multiply two series, you need to call the method ex::series again to convert it to a series object with the usual structure (expansion plus order term). A sample application from special relativity could read:

#include <ginac/ginac.h>
using namespace std;
using namespace GiNaC;

int main()
{
    symbol v("v"), c("c");
    
    ex gamma = 1/sqrt(1 - pow(v/c,2));
    ex mass_nonrel = gamma.series(v==0, 10);
    
    cout << "the relativistic mass increase with v is " << endl
         << mass_nonrel << endl;
    
    cout << "the inverse square of this series is " << endl
         << pow(mass_nonrel,-2).series(v==0, 10) << endl;
}

Only calling the series method makes the last output simplify to 1-v^2/c^2+O(v^10), without that call we would just have a long series raised to the power -2.

As another instructive application, let us calculate the numerical value of Archimedes' constant Pi (for which there already exists the built-in constant Pi) using John Machin's amazing formula Pi==16*atan(1/5)-4*atan(1/239). This equation (and similar ones) were used for over 200 years for computing digits of pi (see Pi Unleashed). We may expand the arcus tangent around 0 and insert the fractions 1/5 and 1/239. However, as we have seen, a series in GiNaC carries an order term with it and the question arises what the system is supposed to do when the fractions are plugged into that order term. The solution is to use the function series_to_poly() to simply strip the order term off:

#include <ginac/ginac.h>
using namespace GiNaC;

ex machin_pi(int degr)
{
    symbol x;
    ex pi_expansion = series_to_poly(atan(x).series(x,degr));
    ex pi_approx = 16*pi_expansion.subs(x==numeric(1,5))
                   -4*pi_expansion.subs(x==numeric(1,239));
    return pi_approx;
}

int main()
{
    using std::cout;  // just for fun, another way of...
    using std::endl;  // ...dealing with this namespace std.
    ex pi_frac;
    for (int i=2; i<12; i+=2) {
        pi_frac = machin_pi(i);
        cout << i << ":\t" << pi_frac << endl
             << "\t" << pi_frac.evalf() << endl;
    }
    return 0;
}

Note how we just called .series(x,degr) instead of .series(x==0,degr). This is a simple shortcut for ex's method series(): if the first argument is a symbol the expression is expanded in that symbol around point 0. When you run this program, it will type out:

2:      3804/1195
        3.1832635983263598326
4:      5359397032/1706489875
        3.1405970293260603143
6:      38279241713339684/12184551018734375
        3.141621029325034425
8:      76528487109180192540976/24359780855939418203125
        3.141591772182177295
10:     327853873402258685803048818236/104359128170408663038552734375
        3.1415926824043995174