## Status Update—Week 8

In Uncategorized on July 16, 2011 by nessgrh

## Introduction

So this week I finally extended the integration lookup tables to now include trigonomic, hyperbolic and exponential functions, certain algebraic functions, error functions, the logarithm, the error function, and various bessel functions. I also did some miscellaneous fixes and extensions which I might mention in passing. Furthermore I decided to create a “TODO later” list: basically a list of things that I really really do want to fix, but that will not go into my first integrals pull request. For I have been working on this for so (comparatively) long, I feel I need to stabilise the branch. Indeed the only thing left to do before I subsmit it is a decent treatment of branched functions. After that I think progress can be incremential, and I will start working on the “TODO later” stuff in a gsoc-3 branch.

This post consists of three parts: the first shows results of this weeks work. The second describes some items on the “TODO later” list (to give an idea of what will be missing in the upcoming pull request). The third branch describes my ideas for polar numbers in detail, which I will start implementing tomorrow (I did not actually do any work today, instead I want to the Calypso Water park (near Ottawa). Hurray for flexible work hours.)

## Visible Results of the Extended Tables

Basically now the code can solve the same sort of problems as before, just on many more classes of functions. One simple example is some indefinite integration. A few comments are in order. First of all, there is no use hiding that indefinite integration is hardly the strength of the algorithm. The examples show that it can be made to work in some situations, but basically any composition of functions or similar will stop it from working. (The same is true for definite integration, of course, but there nobody expects much more *g*.) We see in equations 4, 5, 7, 9 (sorry for the non-consecutive numbering) that integration of simple combinations of trigonometric and polynomial functions can be done, and can indeed be faster than existing code. Equations 12, 14, 15 show recursive mellin transforms in action, of course this gets much slower. Equations 18, 19, 20 and 38 show integration of bessel functions. Incidentially I could not verify the last by differentiation, so I will have to look into this again—before submitting the pull request, of course. Equations 40 and 41 illustrate integration of bessel functions. Equation 42 shows a slight shortcoming: the integral can be expressed using struve functions (I think), but they are just not in the table for expanding g-functions (and neither in sympy). Equation 43 shows a more serious limitation: the integral can be done using integration by parts, but the integrand is *not* a meijer g-function (mellin transform involves digamma functions), so this is not an integral the algorithm will be able to determine.

Let’s now look at some definite integration. There is not very much to say here, except that I show these integrals because they can be done. There are many more like them that can be done, and extremely many more like them that cannot be done. Some of these might look odd (23, 36); I will explain what is going on in the “TODO later” section.

## The “TODO Later” list

One problem is that hyperexpand() sometimes is not very clever. For example in equation 28 of the definite integration examples, we see that apparently the integral cannot be done in named special functions. hyperexpand() retains the meijer g-function because the hypergeometric functions cannot be expanded; this is usually a good idea since a g-function expansion in hypergeometric functions can be very messy. However in this case the hypergeometric functions actually cancel! hyperexpand() note this and return the simplified result. However things are not quite so easy. Sometimes there is choice involved in what hypergeometric functions to expand into, and only some choices cancel. Moreover, sometimes the hypergeometric functions do not cancel, but instead can be evaluated at special arguments. Here are a few examples. Equation 3 shows some non-obvious cancellation: the gammas have to be simplified, and then the symmetry of the arguments of the hypergeometric functions has to be used. Equation five shows an example that can be done because of the special argument, indeed 2F1 can always be summed at unity, using a theorem of Gauss’. Equation 7 shows something more complicated still, I know it can be done but I am not sure how (maybe using a so-called quadratic transformation). One side-effect of this shortcoming is that currently I have had to put formulae for e.g. $J_a(x) J_b(x)$ into the table, even though they should be derivable (at a time cost, of course). This is a bad thing since it leads to an explosion of entries.

Another thing that should be fixed is that the convergence conditions for the integral of two g-functions are sometimes just too weak. For example none of the following can currently be done: $\int_0^\infty \sin{x} \cos{x} \frac{\mathrm{d}x}{\sqrt{x}}$, $\int_0^\infty e^{-x} \sinh{\frac{x}{2}} \mathrm{d}x$, $\int_0^\infty e^{-x} \sinh{x} \frac{\mathrm{d}x}{x\sqrt{x}}$, $\int_0^\infty e^x H(1-x) \mathrm{d}x$. This also stops the recursive algorithm from working.

Then the tables can and should of course be extended further. Notably missing currently are inverse trigonometric functions and powers of trigonometric/hyperbolic functions (though these should be derivable).

As another point, the code to transform an integrand and limits into a known form is not currently very clever. For example $\int_1^\infty \frac{\sin(zu)}{(u^2 - 1)^a} \mathrm{d}u$ can be done by rewriting it as $\int_0^\infty H(u^2-1) \frac{\sin(zu)}{(u^2 - 1)^a} \mathrm{d}u$, but the code cannot do this, yet.

On a related note, I should probably write a function to simplify bessel functions, e.g. equation 36 of the definite integration examples is just a single bessel I function (times a step function), but this is far from obvious. Importantly this needs proper treatment of polar numbers to get the branch factors right.

## Introducing Polar Numbers

Finally, my new plan to treat branch cuts. I think there is no way around implementing proper polar numbers, in some way. Let’s first talk about the mathematics, and then about a possible implementation.

Let $\mathcal{S}$ again denote the riemann surface of the logarithm, and let $\mathbb{C}$ denote the ordinary complex line. [We shall sometimes think of $\mathcal{S}$ as an open subset of a topological space $S$ which contains one additional point called "origin"; notably there is no complex structure near it.] We have the holomorphic map $Exp: \mathbb{C} \rightarrow \mathcal{S}$ which is indeed a biholomorphism. Its inverse is $Log: \mathcal{S} \rightarrow \mathbb{C}$. We also have a holomorphic projection map $p: \mathcal{S} \rightarrow \mathbb{C}, Exp(z) \mapsto exp(z)$, where $exp(z)$ denotes the ordinary exponential function from $\mathbb{C}$ to itself. There is a section (one-sided inverse) of $p$ related to the standard branch of the argument which we shall denote $q: \mathbb{C} \rightarrow S, exp(k + i \theta) \mapsto Exp(k + i \theta), 0 \in \mathbb{C} \mapsto 0 \in S$, where $k$ is real and $-\pi < \theta \le \pi$. This function is not continuous at the non-positive reals, but holomorphic outside thereof.

Multiplication in $\mathcal{S}$ is defined in the obvious way. If $f: D \rightarrow \mathcal{S}$ is a meromorphic (say) function, then $p \circ f: D \rightarrow \mathcal{C}$ is also meromorphic. Similarly if $g: \mathbb{C} \rightarrow D$ is any meromorphic (say) function, then so is $g \circ p: \mathcal{S} \rightarrow D$. This allows us to extend many more functions from $\mathbb{C}$ to $\mathcal{S}$, and this is often done implicitely. Let’s mention two peculiarites: the sum of two elements in $\mathcal{S}$ cannot be defined very satisfactorily (although of course their sum in $\mathbb{C}$ is well-defined and well-behaved). Also if $a \in \mathcal{S}, b \in \mathbb{C}$ we can define $a^b = Exp(b Log(a))$, and if instead $b \in \mathcal{S}$ we can define $a^b = Exp(p(b) Log(a))$. With these definitions all the usual algebraic relations (like $(xy)^z = x^z y^z$) hold for all $x, y, z \in \mathcal{S}$.

Now about the implementation. First we need a function to create/represent elements of $\mathcal{S}$. I had thought about exp(x, polar=True), polar being a new assumption. While I still think introducing an assumption is a good thing, I’m not so sure about this notation any more. The advantage is that the object looks very much like exp, so much code (printing, simplification) knows what to do with it. The disadvantage is that the assumption is easily lost, e.g. obj.func(*obj.args) will lose it, and also code might just create new exp objects, not knowing that there can be assumptions around. So probably creating a new function exp_polar is better.

Next I want to add an assumption polar which can be given to symbols. It does not interact much with other assumptions, except that positive shall imply polar. Much simplification code (notably in expand_log, powdenest, powsimp, logcombine) can now be updated to look for the polar assumption instead of the cruder positive.

Finally, for sake of usability, I propose the following functions: polarify(expr, subs=True/False) traverses the expression tree and turns any non-positive numeric constants (e.g. 1+I) into their image under $q$ (namely sqrt(2)*exp_polar(I*pi/4)). It also converts instances of exp into exp_polar. If subs=True, also all symbols that are not polar will be substituted for polar dummies; thus polarify() behaves much like posify(). A function unpolarify(expr) that essentially computes the image of expr under $p$, in as simple a form as possible. Thus it converts the outermost instances of exp_polar into exp, and recursively does the same thing for arguments of unbranched functions (e.g. powers with integral exponents, or things like sin(x)—here probably the function class should give a hint, similar to nargs). Then functions like expand, powdenest can accept a new argument polar=True/False which works like force, except that it calls polarify instead of posify. Finally for convenience Expr should get a method extract_branch_factor (a wrapper around extract_multiplicatively I imagine) that extracts factors of the form exp_polar(I*pi*n), n nt integer, which allows branched functions to expose their branching behaviour.

This is about it; I feel I have spent forever writing this post (it’s actually just 2.5 hours :D).

### 7 Responses to “Status Update—Week 8”

1. It’s better to just run with the cache off (isympy -C) then to include cache clearing in your timings.

2. This indefinite integration is interesting (to me esp.). Here, we see that it is faster than the current slow implementation (not surprizing), but reduced in power in the sense that it gives rediculuously complicated answers when very simple ones exist (like http://www.wolframalpha.com/input/?i=integrate%28sin%28x%29*cos%28x%29*exp%28-x%29%2C+x%29 vs. your [48]). I will show this to Mateusz and see what he thinks.

• Oh, yes, I forgot to comment on this. I meant to say that the answer is ridiculous.

3. I showed some of your stuff to Mateusz. We looked at trying to get your result from integrate(erf(x)*besselj(a, x), (x, 0, oo) to look like this one (http://www.wolframalpha.com/input/?_=1310834567360&i=integrate(exp(-x)*besselj(a%2c+x)%2c+(x%2c+0%2c+oo))&fp=1&incTime=true). We had to send it through several simplifications and manual substitutions in SymPy to get them equivalent. So clearly, the biggest deficiency here is the lack of powerful simplification functions in SymPy.

4. I agree with your plan for polar. It’s better to create exp_polar, as assumptions aren’t generally passed to functions (only when creating Symbols). I also agree about creating a polar assumption. I hope you aren’t bogged down by our mess of an assumptions system :)

polarify() sounds like something that will work 99% of the time, until you get some corner cases where it is not as easy as just simple substitution (more advanced mathematical relations would be involved). But you will have the best ideas on what to do here as you have been right there in the code with the stuff that doesn’t work.

By the way, Mateusz said that you should always try to verify your results against Mathematica.

• > I agree with your plan for polar.

That’s great.

> By the way, Mateusz said that you should always try to verify your results against Mathematica.

I surely do that with every new and reasonably simple result, certainly if it goes into the tests. Admittedly when I get these horrible expressions I usually make no attempt at ascertening their correctness. Trying to get these into somewhat better shape is on my todo later list, but I don’t yet know where to start working on this (I haven’t investigated it at all so far).

5. It would also be fun to start comparing your results against Sage. Mateusz told me he thinks that once our combined integration power is in SymPy (my stuff and your stuff), that Sage will start using us for integration. I think for your branch, it might already be more powerful for definite integration than any of the open source alternatives.