Last week I told you why rs_series doesn’t work with negative or fractional powers because of the constraints of a polynomial back-end and why we need to modify polys. The situation isn’t that hopeless actually. Let’s talk about negative and fractional powers one by one.

Negative Powers

The reason negative exponents work in ring_series is because I modified PolyElement to allow so. In hind sight, it wasn’t the right decision and needs to be replaced with something that doesn’t alter polys. It is rather surprising that I came across a possible solution so late (Now we know why good documentation is so important).

I already knew that polys allows us to create a FractionField. A fraction field over a domain R consists of elements of the form a/b where a and b belong to R. In our case we are interested in the fraction field of polynomial ring, i.e, fractions with polynomials as numerator and denominator. So a/bis not a * b**(-1) but is a / b, where a and b are polynomials.

What was new to me was that just like ring, polys also has sparse field. In effect, it allows us to create sparse rational functions without altering anything.

I modified some functions in ring_series to work with a rational function field here, and it works quite well indeed.

In [42]: from sympy.polys.fields import *

In [43]: F, a, b = field('a, b', QQ)

In [44]: p = rs_sin(a + b, a, 3)/a**5

In [45]: p*a
Out[45]: (-3*a**2*b - 3*a*b**2 + 6*a - b**3 + 6*b)/(6*a**4)

Note that all these are field operations and I haven’t modified field.py in any way. Elegant!

But then again, having a field increases the complexity as we need to evaluate the numerator and denominator separately.

Fractional Powers

Fractional powers are a much trickier case as there is no simple solution to it as above. What we can do is optimise the option I had presented in my last post, i.e, have each fractional power as a generator. But doing that opens up a Pandora’s box. Simple things such as sqrt(a)**2 == a do not hold true any more. The current rs_series treats sqrt(a) as a constant if we are expanding with respect to a:

In [22]: rs_series(sin(a**QQ(1,2)), a**QQ(1,2),5)
Out[22]: -1/6*(sqrt(a))**3 + (sqrt(a))

In [23]: rs_series(sin(a**QQ(1,2)), a,5)
Out[23]: sin(sqrt(a))

So, if we indeed decide to tread this path, we would need to replace a here with sqrt(a)**2. This really complicates the situation as we need to figure out what to replace with. In any calculation the powers change multiple times and each time we’ll need to figure out how to rewrite the series.

Next Week

It is now mostly a design decision whether we want the ring_series to be confined within the polys module. The polys environment allows efficient manipulations of laurent series (with FracField), but I doubt we can achieve the speed we want with puiseux series without modifying polys. One possible solution is to separate the modified parts of polys along with ring_series from polys. We are using polys only because it has the data structure that we want. Separating them would allow us to simultaneously make use of its back-end and not introduce unnecessary complexity in our representation.

Other than that, documentation is another priority now. I had planned to do it earlier too, but couldn’t. This week’s discovery has reminded me of its importance.

Cheers!