About implementing special functions

Into

On this page I'd like to describe the process of implementing new special functions. There should be a list of all point one usually has to take care of.

Warning: This information is still work in progress. Some arguments may be incomplete or wrong.

Getting started

The class Function is the base class for all functions. It is located in sympy.core. First we start a new file with the name of the new function. Then we import the base class with

from sympy.core.function import Function

Now we can subclass from it. The example is based on the Gaussian error function erf(z). We start a new python class with:

class erf(Function):
    """
    The Gauss error function.
    """

I added a first small docstring. (See below for information how a good docstring should look like.)

TODO:

• How to export the new function (add in __init__)
• Better to append in files where similar functions are
• Common baseclass for families of functions (trigonometric, bessel, ...)

Basic properties

The new function erf(z) takes only one (complex) argument. We tell sympy this with the following varibale assignment.

    nargs = 1

This has to be done on the toplevel in the new class. (On the same level where you usually put new def ...: statements.)

At this point we can go into a python shell and call the new function.

>>> z = Symbol("z")
>>> erf(z)
erf(z)

However there is not much more we can do with it yet.

The arguments (in this case z) are stored in the member variable self.args which is always a tuple, even if there is only one element in it.

We can always retrieve the argument z by self.args[0] If the function would take two arguments f(x,y) then we would use self.args[0] and self.args[1] for getteing back x and y.

TODO:

• Are the args sorted? (Remebers me of a atan2 fix to finish)

Function evaluation

To allow for evaluation of our function for given arguments (numbers or symbolic expressions) we have to add the following function to the class:

    @classmethod
    def eval(cls, arg):
        ...

Excerpt from the docstring:

The eval() method is called when the class cls is about to be instantiated and it should return either some simplified instance (possible of some other class) [...]

Now we can put all the knowledge about erf(z) inside this function. For example we know the value erf(0) to be 0, hence we write:

        # Value at zero
        if z is S.Zero:
            return S.Zero

The code in the eval method will usually take the form of a big nested if tree.

Special values

We can now add further well known function values for:

• Values at single points in C
• Values at (real) infinities oo and -oo

In our example of erf we know:

        elif arg is S.Infinity:
            return S.One
        elif arg is S.NegativeInfinity:
            return S.NegativeOne

this works well for finitely many special points z \in C. If we have a function with a more complicated structure the the code will become more involved. For example the evaluation of cos(k*2*pi) with k an integer should always return 0.

Symmetries

• Even/odd functions
• Complex mirror symmetry

Complex mirror symmetry

The definition of what is called "mirror symmetry" is given by the following transformation rule: Assume z to be a complex number, then it holds that conjugate(f(z)) = f(conjugate(z)).

We can implement this rule (read from left to right) with the code snippet:

    def _eval_conjugate(self):
        return self.func(self.args[0].conjugate())

It's better to implement the symmetry rule this way than the other way around (i.e., in .eval()) because it can then work its way down recursively. For example, f(2 + I).conjugate() would fully evaluate to f(2 - I).

Differentiation

Multiple arguments

Series expansions

• Which methods to implement?

Taylor series expansions

If we know a closed form expression for the n-th term in the Taylor series expansion we can implement the method taylor_term(n, x, ...). The example shows the implementation of from a Fresnel function:

    @staticmethod
    def taylor_term(n, x, *previous_terms):
        if n < 0:
            return S.Zero
        else:
            x = sympify(x)                                                                               
            return x*(-x**4)**n*(S(2)**(-2*n)*pi**(2*n))/((4*n+1)*C.factorial(2*n))

This is then helpful for expressing Taylor expansions like:

z = Symbol("z")
n = Symbol("n", integer=True)
Sum(fresnelc(z).taylor_term(n, z), (n,0,oo))

Integration

• Enabling function in the integration routines.

• Extending the tables for the meijerG routines

• See sphinx doc: [Extending the Hypergeometric Tables] (http://docs.sympy.org/0.7.1/modules/simplify/hyperexpand.html#extending-the-hypergeometric-tables)

Limit calculation

• Simple cases

• Extensions to allow a function to be tractable by the Gruntz algorithm