Symbolic definite integration: methods and open issues

Daniel Lichtblau

Abstract

The computation of definite integrals presents one with a variety of choices. There are various methods such as Newton-Leibniz or Slater's convolution method. There are issues such as whether to split or merge sums, how to search for singularities on the path of integration, when to issue conditional results, how to assess (possibly conditional) convergence, and more. These various considerations moreover interact with one another in a multitude of ways. Herein we discuss these various issues and illustrate with examples.

Introduction

In principle, algorithmic integration is fairly well understood. There are algorithms to handle various sorts of integrands and integration bounds. But, in the words of librettist Ira Gershwin, "It ain't necessarily so" [],  and actual practice is an altogether different situation. Troublesome areas include:

Determining (possibly conditional) convergence.

Recognizing and handling cancellation of additive singularities.

Intricacies of the algorithmic implementation, e.g. need for recursion or reliance on powerful transformations.

Handling of parameters.

Working with legacy code.

In this report I will describe and illustrate many of the issues I have encountered in the process of working on the Mathematica [] Integrate code base. As best I can tell, they also apply in part to other existing programs (though approaches to handling said issues may vary).

This report discusses the situation primarily as it evolved through Mathematica versions 5 and 6 development. I wish to thank Oleg Marichev (in honor of whose 60th birthday I spoke on this topic) [] for many discussions where he patiently tried to walk me through issues involved in convolution of MeijerG functions. I also thank Bhuvanesh Bhatt, a senior Software Quality Assurance engineer at Wolfram Research, for keeping me apprised of the situation regarding integration bugs in Mathematica, and for the countless hours he spent tracking and diagnosing them, and making test cases.

While we focus on particulars of Mathematica handling, I would be remiss in not pointing out that there is other work on this topic, some of which takes very different approaches to similar issues. In particular, see Davenport's synopsis of this same topic [], and various articles authored or coauthored by Jeffrey on removal of path singularities [, , ].

Basic Structure of Integrate Code in Mathematica

Indefinite Integration

The indefinite integration code consists primarily of a partial implementation of the Risch algorithm [, ], in addition to extensive table lookup methods. These latter attempt to handle integrands with exponentials, trigs (and/or hyperbolics), elliptic integrals, and integrands involving special functions, particularly after the Risch methods have given up. The actual situation is slightly more complicated in that Risch code may handle part of an integrand, sending the rest to table lookup code. That in turn may do transformations and call the integration process recursively.

In addition to what is stated above, some integrand transformations are attempted prior to all else. The idea here is to recognize cases where one integrand may be converted to something that is easier to work with, but equivalent modulo a differential constant.

Also useful are mathematically equivalent transformations such as factorization and partial fractions decomposition. An added wrinkle is that such transformations can invert one another, and it is often difficult to recognize which one will be beneficial to a given integrand. Thus one must take care to avoid attempting them blindly, and, in so doing, incurring the wrath of the gods of infinite recursion.

Definite Integration

Definite integration is done via a number of methods as indicated below.

Special case contour integration.

Newton-Leibniz code specialized for integrands of the form  rational × trig or rational × exponential.

Newton-Leibniz code specialized for integrands containing logs or polylogs.

A general case implementation of Newton-Leibniz integration.

An implementation of definite integration by convolution of MeijerG functions []. This requires that we integrate from 0 to infinity. It also requires that the integrand be represented as a power of the integration variable times one or two MeijerG functions. Since a UnitStep function may be represented as a MeijerG, we may lift the infinite range requirement whenever the integrand requires only one MeijerG. This method in particular uses several transformation tactics to handle algebraics, trigs, logs, exponentials, and so forth.

The overall implementation is best described as a polyalgorithm that calls on any or all of the above in various ways, depending on heuristics (or perhaps the mood of the little man inside the code).

Brief Descriptions and Examples of the Two Primary Methods of Definite Integration

Newton-Leibniz

First we find an antiderivative to the integrand. We then look at the integration path, finding candidate singular points. We will treat these in one of two ways (based on crude heuristics): either issue a proviso on parameters that guarantee the point is not a "bad" point, or else split the integration path at that point and use limits as we approach from either side. Once the integration path is split into suitable segments we take appropriate limits of the antiderivative as we approach the segment endpoints from appropriate directions, then sum these signed values to get the resulting definite integral.

Here is a classical example that incorrectly gave zero in a prior version of Mathematica. Note that one path, from -1+i to -1-i, must be split where it crosses the (negative) real axis.

Integrate[1/z, {z, 1 + I, -1 + I, -1 - I, 1 - I, 1 + I}]

Slater Convolution via MeijerG Products

The idea in this case is to make suitable transformations so that we have an integrand in, say, the variable z, as a product (z^k)MeijerG[(Subscript[expression, 1])]MeijerG[(Subscript[expression, 2])] where we integrate from 0 to infinity. While this path may seem rather specialized, bear in mind that we can do translations and rotations to move a path to the real axis, and we can represent UnitStep as a MeijerG function. This means we can treat a finite interval on the positive axis as though it was semi-infinite, at the cost of a MeijerG factor. The basic convolution theorem is discussed in [].

Several common functions such as trigs, various exponentials, Bessel functions, and certain common algebraics e.g. (((a+b (x^n)))^m) may be represented as MeijerG functions, while others may be brought to MeijerG form after some suitable transformation of variables. Hence we look for viable transformations and then do a table driven conversion in an effort to obtain the desired form.

Once in the factored form above (or perhaps a simpler form wherein any factor or factors are not present) we go into convolution code, at the heart of which lies a function that formulates the integrated result as a new MeijerG expression via the Slater theorem. It then performs some manipulations in an effort to recast in terms of better known functions.

Here is a simple example.

val=Integrate[BesselJ[2,x] (x^2) (E^(-x+2)),{x,0,∞},GenerateConditions->False]
nval=NIntegrate[BesselJ[2,x] (x^2) (E^(-x+2)),{x,0,∞}]
Chop[val-nval,(1/(10^8))]

This operated by first converting the integrand to an "inert" form representing the integrand product as

(x^2) MeijerG[{{},{}},{{0},{}},x] MeijerG[{{},{}},{{1},{-1}},((x^2)/4)]

Software Engineering Problems Associated with this Body of Code

Below is a list of some of the issues that any powerful symbolic integration implementation will face.

Different modules, written by different people, were often not on the best of speaking terms. Hence various pieces of functionality of potentially general use might require reinventing, in some cases with bad consequences.

Much is legacy code. It was written over a period of about 15 years, by several different people. It is spread over some several dozen source code files. No one person understands all of it, and some parts are no doubt not understood by anyone (this is one of the pitfalls of legacy code). I find it unlikely that in an area of this scope any powerful body of code will be entirely understood by the team that develops it.

Some of the steps used are quite fragile. Specifically they may be sensitive to small changes in seemingly unrelated functions such as Together (providing a "canonical" form to rational functions), Apart (a partial fractions decomposition, which is in fact stitched at the hip to Integrate), Factor, Solve, Simplify, and others. After encountering this problem in many forms I have come to conclude that it really is endemic to symbolic integration and not simply an artifact of implementation details. The entire idea of working with transformations of integrands, coupled with the mathematical impossibility of creating canonical forms for all possible expressions, and the difficulty of forcing transformations to be improvements for a given purpose, make this a thorny issue. Indeed, certain types of standard tactics such as integration by parts or use of l'Hospital's rule for limit extraction can actually lead to infinite recursion if not done with considerable caution.

There are tradeoffs to be made between speed and power. One wants to try certain transformations, for example, in order to handle certain classes of problem. But there is no plausible way in general to perfectly delimit cases that may hang in the transformation process. Hence some form of delimiting is required e.g. by time or operation count.

Many classes of integration problem can make use of assumptions regarding, say, behavior of integrand at infinity. This creates a potential for bugs when the code is not sufficiently careful to delimit what inputs enter handlers that may utilize such assumptions. Moreover even where a handler explicitly tries to determine such behavior there is the issue alluded to above, where, say, limit extraction may hang if not suitably constrained. But when constraints (on time, memory, or some measure of operation count) are used, how is one to handle an aborted intermediate result? With a warning message? By giving up? By continuing as though it had been a "good" case e.g. of convergence at infinity?

Symbolic integration is one of the most complicated pieces of machinery to be found in algorithmic mathematics. The reliance on other functions as well as intricacies of integrand transformation attempts certainly implies that it will have a nontrivial implementation. Hence it is vital that the pieces be adequately documented. This is a rule throughout the field of software development, but the importance cannot be overstated in the context of integration.

When I began work on this body of code, more than one out of every four open bugs in the Mathematica kernel (with around 2000 functions) was in the category of definite integration (a part of one function, Integrate). This mass posed several problems in and of itself. First, the function clearly receives widespread usage, and hence it is difficult to start from scratch. Second, the scope made it a bit difficult to perform adequate triage (that is to say, it is hard to see the forest for the trees). Third, the massive overhaul needed strongly implied that there would be considerable breakage, at least in the short term. Fourth is the likelihood that the scope of trouble exceeds the capabilities of any one developer (I can attest that it certainly exceeds the capabilities of this particular developer).

Issues in the Implementation of Indefinite Integration

We will discuss in brief some issues that tend to be specific to indefinite integration. As our primary focus is on definite integration we defer to the next section those that are common to both.

The Curse of Recursing

Often an integral may be broken into two parts. The technical term for the first part is "done" while that for the second is "not done". The second term might be further rewritten in ways that return one to the original problem or a variant thereof, thus leading to recursive splitting. As an example of an integrand that might elicit such behavior, try Integrate[Sin[x]/(Sqrt[a-x]*(1+x)),x] in version 4 of Mathematica. One approach to fixing this involves use of hashing to recognize integrands that have previously come our way.

Another common cause of descent into the infinite is the utilization of pairs of inverse transformations. We might, for example, convert a trig to exponentials in one handler, and convert back to trigs in another, applied later. We try to avoid this pitfall by using Block to (what else?) block the dual handler when we do one such transformation.

There is a bright side. We now tend to handle more problems than in the past. Here is an example that we now handle due to more active use of transformations of the integrand. It comes from an integration test suite in [].

∫(((Tan[x])^2)/(Sqrt[(1-(((Sin[x])^2)/2)) (1-((2 ((Sin[x])^2))/5))]))\[DifferentialD]x

Out on a Limb with a Cut Branch

Frequently we require transformations that bring into play multivalued functions. The consequence is that we may arrive at an antiderivative that is only correct up to a piecewise multiplicative constant. Hence we now attempt to restore the proper factor. This is not always trivial and often leads to a significantly enlarged form of result. Moreover it is not trivial to recognize in all cases how to correctly reverse effects of such transformations.

Transformations

In many cases it is well understood how one might obtain an antiderivative for a particular class of integrand e.g. rational functions of trigs. But the needed transformations must be applied carefully in order to avoid potential explosion in intermediate complexity or that of the final result.

To Expand or Not Expand, Indefinitely

This is a (very important) special case of a transformation. As the issues are a subset of those that arise in the context of definite integration we cover it there instead.

Issues in the Implementation of Definite Integration

In this section we arrive at the main focus of this report. Within symbolic definite integration one encounters a wide array of issues. In the subsections below I will endeavor to present and illustrate many of them.

To Expand or Not Expand, Definitely

When one has certain types of input it makes sense to expand over summands and integrate term-by-term. In other cases this can be a very bad thing to do. So the question is when should one expand the input and loop over summands. We will illustrate the issue with a few examples.

First, a "simple" integrand involving a rational function.

p=7 (x^13)+10 (x^8)+4 (x^7)-7 (x^6)-4 (x^3)-4 (x^2)+3 x+3;
q=(x^14)-2 (x^8)-2 (x^7)-2 (x^4)-4 (x^3)-(x^2)+2 x+1;
result=(∫_2^5)(p/q)\[DifferentialD]x

We do a numeric check that this is correct.

N[result]-NIntegrate[p/q,{x,2,5}]

Now we expand and integrate over each term, then try to simplify the result. It is a mess.

result2=Simplify[Map[Integrate[#,{x,2,5}]&,Expand[(p/q)]]]

N[result2-result]

We see that the mess at least is correct. But it is obvious that we are better served by avoiding such expansion in this example.

Here is an example where one must not expand. The reason is that the summands have cancelling singularities.

ii=(∫_0^∞)((E^(I/x))-1) Sin[x]\[DifferentialD]x

ii=(∫_0^∞)((E^(I/x))-1) Sin[x]\[DifferentialD]x

Again we validate this with a numeric check.

N[ii]-NIntegrate[((E^(I/x))-1) Sin[x],{x,.001,1000},WorkingPrecision->25]

Clearly we cannot expand and handle summands separately.

Map[Integrate[#,{x,0,Infinity}]&,Expand[((E^(I/x))-1) Sin[x]]]

By contrast, here is a case where preexpansion helps considerably.

integrand=((a c Cos[t])/(g s))+((b q Cos[2 t])/(c f))+((c Cos[3 t])/(a d))+((d f Cos[4 t])/(a h n))+((e q Cos[5 t])/(a g))+((f l Cos[6 t])/(m r w))+((b g Cos[7 t])/(n o x))+((h Sin[t])/(b c))+((i Sin[2 t])/(e h r))+((j y Sin[3 t])/(l p))+((d k Sin[4 t])/(c j))+((a l m Sin[5 t])/(b f h s))+((m p Sin[6 t])/(j k))+((n q Sin[7 t])/(c x));
Timing[(∫_0^(2 π))integrand\[DifferentialD]t]

Versions of Mathematica that did not expand but instead tried to work with the entire integrand at one time took about a minute to handle the same integration. The reason is that transformations involving several trigonometric summands can be costly both in time and memory.

Yet another drawback to expansion is that the individual pieces might have different ideas regarding provisos for parameters, in particular if they go through different routes in the code. If conflicting provisos emerge we might either get a useless result (if we fail to recognize that the conditions cannot all hold) or else be forced to redo the integral, thus having wasted time processing summands individually.

Heuristics Involving Method Choices and Transformations of Integrand

Quite often one encounters an integrand that might be handled in different ways. As the result can vary considerably in terms of speed and/or complexity it is a (wide) open problem to optimally dispatch based on structure of the integrand. There is one nice feature to the decision process of whether or not to try an approach based on evaluating an antiderivative at endpoints: we must first find that antiderivative. If we fail at that stage, we certainly know it is not a viable approach!

Here is an example where use of convolution approach gives a result faster.

Timing[(∫_0^∞)(1/(x^8))Sin[x] Sin[(x/3)] Sin[(x/5)] Sin[(x/7)] Sin[(x/9)] Sin[(x/11)] Sin[(x/13)] Sin[(x/15)]\[DifferentialD]x]

Here is an example where evaluation of an antiderivative gives a preferred form of result.

Integrate[Sin[x-y],{y,x,Pi},Assumptions->{0<x<Pi},GenerateConditions->False]

While the nondefault option value may make it look pathological, it is spawned from the example below which does not use any nondefault option settings.

Integrate[Abs[Sin[x-y]],{x,0,Pi},{y,0,Pi}]

One would certainly prefer that either method give the "nice" form of result shown above. It is an unfortunate fact that this state of affairs is quite difficult to achieve. Thus one must use heuristics to distinguish cases that might receive more favorable treatment for one or the other method. Noting that the integrands in these examples are both primarily trigonometric functions of linear functions of the variable, it becomes clear that formulation of good heuristics is by no means trivial.

A related issue of heuristics is that frequently an integrand may appear amenable to any of several different transformations. For example, for a quadratic radical function of the integration variable we might want to transform to a MeijerG function, or do  a linear change of coordinates in an effort to simplify or remove the radical.

Another related issue is that some transformations may not be valid for all possible values of a parameter. For example, one might wish to write (((a+b (x^d)))^n) as (a^n) (((1+((b (x^d))/a)))^n) (for integration variable x) in order to convert to a MeijerG form. This transformation is of course not valid for all values of the parameter a. We then must choose between issuing a conditional answer, or abandoning this approach and trying another that might give an unconditional result.

Still another related issue is in how to transform products into MeijerG functions. There may be several possibilities, with the quality of outcome dependent on the choices made.

Generations of Generation of Provisos (Results That Depend on Conditions)

Below I indicate ways in which generation of conditions can be problematic. Several are illustrated with older versions of Mathematica because we have made improvements that render the specific example obsolete. But the ideas behind them are general and no doubt related problems still lurk in the current implementation.

Propagation of Conditions for Multiple Integrals

For multiple integration one may want to propagate conditions from inner to outer integrations. This is problematic if the conditions generated cannot adequately be simplified. On the other hand if we fail to generate and propagate them we run a greater risk of obtaining a bad result.

Here is an example where in version 5.0 we required conditions in order to propagate singularity information.

result=Integrate[((Abs[x-y])^n),{x,0,1},{y,0,1},GenerateConditions->True]

A numeric check validates this.

(result/.n->3.2) - NIntegrate[((Abs[x-y])^3.2),{x,0,1},{y,0,1}]

Note that if we do not insist on generation of conditions we get a result that is quite obviously incorrect insofar as the correct one clearly must lie between zero and one for positive values of n.

Integrate[Abs[x - y]^n, {x, 0, 1}, {y, 0, 1}, GenerateConditions->False]

(∫_0^1)(∫_0^1)((Abs[x-y])^n)\[DifferentialD]y\[DifferentialD]x/.n->3.2

Excessive Conditions

Various algorithm implementations may force generation of unneeded conditions. For example, MeijerG convolution will require that certain values lie in "wedges" emanating from the origin, and Newton-Leibniz methods may issue conditions based on the specific form of the antiderivative.

Here is a simple example using version 5.0. Special case code for handling exponentials via convolution wants to insist that a parameter take a negative real part.

Integrate[(E^(-((a (((x-b))^2))/(sigma^2)))),{x,-∞,∞},Assumptions->{a>0,sigma∈Reals}]

The current behavior is nicer.

Integrate[(E^(-((a (((x-b))^2))/(sigma^2)))),{x,-∞,∞},Assumptions->{a>0,sigma∈Reals}]

Necessity of Condition Generation

Often one may wish to ignore conditions. This happens for example in cases where we know in advance that they will be satisfied by parameter values we may later choose. A problem is that results can depend on what path is taken in the code. If the integration is split into parts, say because the path is split, and if the parts take paths in the code that make different assumptions about parameters, then results might be entirely incorrect. One manifestation might be incorrect cancellation. We illustrate with the preceding example. In the code it handles half the integration path using one assumption about the parameter b, handles the other half range using a contrary assumption, and the end result is quite wrong.

Integrate[(E^(-((a (((x-b))^2))/(sigma^2)))),{x,-∞,∞},GenerateConditions->False,Assumptions->{a>0,sigma∈Reals}]

The upshot is that partial results cancelled, with no indication that they require conflicting conditions to hold.

Genericity of Generated Conditions

Conditions generated may be only generically correct. Here is an example that exhibits this phenomenon. Specifically, when Im[a]==Im[b] (that is, a vertical integration path) we have a problem because that their difference appears in some denominators.

(∫_a^b)Log[x]\[DifferentialD]x

Assessment of Convergence

Testing for convergence is closely related to generation of conditions insofar as the former can depend on parameter values. But convergence testing, even in the absence of symbolic parameters, is no easy matter (it is complicated by oscillatory factors, possible cancellation of singular terms, and so forth). Hence code that tests integrands for convergence tends to be of an ad hoc nature. While it usually serves well it is by no means a science at this time, and I do not know how to make it one. Were one to issue "uncertain of convergence" messages in all possible cases where we cannot verify convergence or identify conditions to generate that would ensure convergence, we would have a flood of messages. Even worse, we would lose many integrals. The reason is that summands can separately spawn conflicting conditions, and this can happen in cases when they should instead cancel singularities in pairs.

Here are some simple examples that may help to give an idea of the difficulties lurking within convergence assessment. Mathematica will evaluate them correctly but this was not always the case.

(∫_0^∞)((Cos[(x^2)])/(Log[x]))\[DifferentialD]x

(∫_0^∞)(E^x) Sech[a x] x\[DifferentialD]x

The first diverges, but not due to problems at either endpoint (there is a pole at 1). The second converges conditionally, and version 4 of Mathematica even knows this, but then gets the condition wrong.

Below is an example that gives an idea of what is involved in sorting out conditions for convergence. We need to look at possibilities of oscillatory, exponentially growing, and exponentially damped factors in order to figure out the correct conditions on parameters. In this case the original assumption suffices to guarantee convergence.

Integrate[(Sin[betty*r]*Sin[gramma*r])/E^(alf*r^2),{r,-Infinity,Infinity},Assumptions->Re[alf]>0]

Here is a slightly more difficult variant, one that requires nontrivial conditions on parameters beyond what is given in the assumptions option.

Integrate[(Sin[betty*r^2]*Sin[gramma*r^2])/E^(alf*r^2),{r,-Infinity,Infinity},Assumptions->Re[alf]>0]

As another example we show an integral that converges. The integrand is a sum of two terms, each of which separately will diverge. If the code splits this then it must recognize that there will be a cancellation of singular parts.

Integrate[(Log[x]-Log[a])/(x-a),{x,0,a},Assumptions->a>1]

Difficulties Involving Parameters and Detection of Singularities

One issue is in finding singular points on the integration path. Even something as simple as a trig function can be problematic.  There are other tar pits lurking beneath the surface of the swamp.

Parametric Singularities

The presence of elliptic functions in the antiderivative usually removes any hope of correctly detecting parameter dependent singularities.  Below is an example.

(∫_0^(π/2))(-(m MeijerG[{{1,1},{}},{{(1/2),(3/2)},{}},(1/(1-m))])+π EllipticK[m] (1+(-1+m) ((Sin[y])^2)))/(π (Sqrt[1+(-1+m) ((Sin[y])^2)]))\[DifferentialD]y

Whether the result is worth the screen real estate it occupies is open to debate.

Transcendentals and Singularity Detection

Often finding singularities depends on finding roots of expressions. These need not be algebraic, and a general purpose transcendental root finder is a nontrivial undertaking (also it might be slow). In the example below, we find the bad point at the origin, and correctly decide the integral is divergent.

(∫_(-(1/2))^(1/2))((-2+(x^2)+20 Cos[x])/(((-6 x+(x^3)+60 Sin[x]))^2))\[DifferentialD]x

But minor modifications give rise to a singularity that is not at an algebraic number nor at a "convenient" multiple thereof (e.g. π times a rational).

ii=(∫_1^5)((-2+(x^2)+20 Cos[x])/(-6 x+(x^3)+60 Sin[x]))\[DifferentialD]x
N[ii]
NIntegrate[((-2+(x^2)+20 Cos[x])/(-6 x+(x^3)+60 Sin[x])),{x,1,5}]

Clearly this diverges. The denominator has a pole in the integration path, and it is a (nonalgebraic) root of a transcendental equation.

FindRoot[-6 x+(x^3)+60 Sin[x]==0,{x,2}]

Parametrized Singularities in Multiple Integrals

Parametrized singularities can cause trouble, particularly in multidimensional integration. Here is an example from []. We begin with a small variation on the actual problem.

(∫_0^π)(∫_0^π)Abs[Sin[x-y]]\[DifferentialD]y\[DifferentialD]x

Now we double the region of integration in both directions.

(∫_0^(2 π))(∫_0^(2 π))Abs[Sin[x-y]]\[DifferentialD]y\[DifferentialD]x

The second one is off by a factor of 2. The problem is that one singularity line, x==y, is recognized. But the broken line x==y+π modulo 2π goes unrecognized. Hence we do not handle ranges correctly, splitting only into two rather than three pieces in the first level of integration. Note that this is in version 5.0 of Mathematica; later versions do get the correct value of 8 π.

Branch Cuts Intersecting the Integration Path

It is important to assess whether an integration path crosses a branch cut of an antiderivative (so that we might split the path into segments). The example below does this in order to get the correct result.

i1=(∫_(-1-I)^(-2+I))Log[x]\[DifferentialD]x
N[i1]
NIntegrate[Log[x],{x,-1-I,-2+I}]

Version 4 of Mathematica failed to catch the crossing and gave a result with imaginary part off by 3 π. Needless to say, this problem becomes vastly more difficult if the crossing might or might not exist based on parameter values. One might try the example below to get an indication of what might be a reasonable result.

Integrate[Log[x],{x,-1-I,y}]

Problems with Algebraic Manipulation

Here is an example that was quite wrong prior to version 5.1 (we show the result from version 5.0).

ii=(∫_0^∞)2 y (E^y) BesselK[0,(Sqrt[3]) y Sign[y]]\[DifferentialD]y
N[ii]
NIntegrate[2 y (E^y) BesselK[0,(Sqrt[3]) y Sign[y]],{y,0,∞}]

The problem was in the innards of the convolution code. At a key step we made use of an exponential of the form exponential((-2 i) π-log((4/3))+2 log(2)) and we subsequently required a square root. Once we replace exponential with Exp we lose track of a factor of (-1) and this gave rise to a bad result. While this example is now repaired the general problem remains of how to find and forestall this phenomenon.

Mathematica's integration code once had an interesting method for eluding this branch cut sort of problem. The tactic was to convert various exponential things to trigs in the hope that anything "bad" e.g. explicitly complex values, would pack up for vacation. Often this worked, and indeed I introduced some bugs simply by disabling some of this and making conversions in the reverse direction. But my feeling is that the problem should be addressed at a more basic level, by finding the culprits and spanking the branch cuts out of them. As an indication that this earlier method had severe limitations, I note that it too failed on the example above.

The next examples have bugs in current versions due to branch cut issues not correctly handled in convolution code. They show that not all changes can be viewed as progress, insofar as they gave correct results in version 4. Needless to say, I hope to address these in the near future.

Integrate[(E^(-I x t)) ((E^((1/3) (-I) (x^3)))-1) (x^(-(3/2))),{x,0,∞},Assumptions->t>0]

A numerical check will reveal that this is simply not correct.

The next example does not even involve a symbolic parameter. Again, the culprit deep down is in convolution of MeijerG functions.

ii=Integrate[ArcSin[v]*Log[v]/v,{v,0,1}]
N[ii]
NIntegrate[ArcSin[v]*Log[v]/v,{v,0,1}]

General Considerations

Here are some of the questions that require thought in the design and implementation of definite integration. They have emerged from study and overhaul of our existing code base, but I believe they apply more generally e.g. to definite summation and related polyalgorithm computational calculus.

Given a choice of methods, which should one attempt first? This can have serious repercussions in terms of speed. For example, what might be a fast MeijerG convolution can be very slow to evaluate as an indefinite integral followed by extraction of limiting values. And of course the opposite can happen. Or the methods might be comparable in speed but give results of vastly different complexity of form.

When or how should one presimplify the input?

When or how should one simplify the result?

Special case methods may be very helpful for speed or simplicity of result. But they also proliferate opportunities for bugs. Hence it requires care to think through what classes are important to handle in these ways.

Some methods require intrinsically "slow" technology. For example, refinement of conditions (which is sometimes essential in order that they not blow up) may require some level of CAD support behind the scenes. Even limit extraction for Newton-Leibniz methods can be slow. We are thus faced with questions of when to apply such technology and how to prevent it from causing many inputs to hang.

In regard to prevention of hanging in computationally intensive technology noted above, we have found it quite necessary to place time constraints on certain pieces of code. (Motivation: often they succeed. If they fail, so be it, and we then try other things.) This gives rise to a new set of problems. One is that asynchronous interrupt handling, required by TimeConstrained, is imperfect and in rare cases will cause a kernel crash. Another is that results now take on a platform dependent nature, and this is seriously unsettling. A possible future direction that will alleviate this: have potentially slow code stopped by some measure of operation count rather than asynchronous interrupts.

Indefinite integration is (generally) more powerful in version 5 than in the past. This has come at the price of speed: it simply tries more transformations and the like. This in turn means all Newton-Leibniz code is at risk of getting slower. We now cache some results from indefinite integration but  this is at best a partial fix to the problem.

Provisional moral, or maybe conundrum: The more things improve, the more opportunity for closely related things to deteriorate.

Summary and Directions for Future Work

Development and overhaul of Integrate has had successes and at least some setbacks. One should hope these latter may prove to be temporary.

Among the successes, virtually all infinite recursion problems from earlier versions has been fixed. Cases where a method failed and we did not move on to try another method, that is, premature bailout, have likewise been addressed. Numerous problems in convergence assessment and singularity detection have likewise been fixed. Much of the idiosyncratic ad hoc code has been removed or rewritten. Most open problems fall into the various categories described above, which, while large, is better understood than the integration swamp that we had four years ago.

This has come at a cost of speed. Moreover we now have some level of platform dependent behavior. Work remains to iron out bugs in convolution code and elsewhere. Heuristics to determine use of expansion, simplification, refinement of conditions, and the like are crude and certainly imperfect. Modularity of the code, which affects usability, maintaince and further development, has improved but is likewise still far from perfect. The same may be said for its documentation.

This presents a synopsis of the issues faced in development of a robust definite integration. To some extent will also give a snapshot of the recent status of Integrate in Mathematica. Lest the reader be left with the idea that many of the issues presented above are specific to the Mathematica implementation of Integrate, we refer to [] for a discussion that overlaps ours but is not specific to any particular body of code.

My opinion is that future work in definite integration within Mathematica should go in a few directions. I list various possibilites below.

Strengthening the MeijerG convolution code, with attention paid to form of results, range of integrands covered, and correct handling of convergence and branch cut problems.

Better assessment of convergence conditions. This will include splitting the GenerateConditions option into two or three distinct ones so that meanings are more clear. At present the code may use it to test for convergence, or to look for parameter conditions that guarantee convergence.

Figuring out a platform independent way to measure work done, and base an interruption scheme thereon. It is a fact of life that certain possibly expensive operations are needed to make an integration code work well, and they can also hang examples if not terminated within reasonable time. Basing this on something as crude as TimeConstrained is problematic, among other reasons because results are not independent of the speed of the machine on which they are run.

Continued debugging based on existing bugs in that code. The overall count of bugs in Integrate is now about 30 percent of what it was when this work was first undertaken. This brings the quantity to a point where a careful categorization of specific problem areas becomes realistic. My anticipation is that some of these specific areas will simply need point fixes. What remains as problematic should be amenable, in parts, to rewriting without deleterious implications to the overall structure and function of Integrate code.

Continued addition to the test suite so as to avoid large scale breakage during overhaul. In past a problem that led to the massive bugs list was that we did not have an adequate test suite It had many integration examples, but not enough for a function of the scope of Integrate). Hence bug fixes might well cause breakage elsewhere and no "trip wire" existed in terms of a sufficiently large test suite. At present, virtually all fixed bugs get made into new tests.

References

V. Adamchik and O. Marichev. The algorithm for calculating integrals of hypergeometric type functions and its realization in Reduce system. Proceedings of the 1990 International Symposium on Symbolic and Algebraic Computation (ISSAC 90) 212-224, S. Watanabe and M. Nagata, editors. 1990. Addison-Wesley.

M. Bronstein. Symbolic Integration I: Transcendental Functions. Algorithms in Computation and Mathematics Volume 1,  1997. Springer.

K. Charlwood. Some integration trials on computer algebra systems. Manuscript, 2002.

J. Davenport. The difficulties of definite integration, 2003. Available electronically at:
http://staff.bath.ac.uk/masjhd/Calculemus2003-paper.pdf
http://www-calfor.lip6.fr/~rr/Calculemus03/davenport.pdf

K. Geddes, S. Czapor, and G. Labahn. Algorithms for Computer Algebra. 1992, Kluwer.

I. Gershwin. It Ain't Necessarily So (song from "Porgy and Bess"). Music by George Gershwin. Lyrics by Ira Gershwin.  Gershwin Publishing Co., 1935.

D. Jeffrey. Integration to obtain expressions valid on domains of maximum extent. Proceedings of the 1993 International Symposium on Symbolic and Algebraic Computation (ISSAC 93) 34-41, M. Bronstein, editor, 1993. ACM Press.

D. Jeffrey. Rectifying transformations for the integration of rational trigonometric functions. Journal of Symbolic Computation 24:567-573, 1997.

D. Jeffrey and A. Rich. The evaluation of trigonometric integrals avoiding spurious discontinuities. ACM Transactions on Mathematical Software 20:124-135, 1994.

D. Lichtblau. "Symbolic Definite Integration". Talk presented at the 2005 Mathematica Technology Conference.
http://library.wolfram.com/infocenter/Conferences/5832/

[Math Horizons 2002]. Problem 157: A Double Integral, proposed by K. S. Murray. Problem Section, M. Klamkin and A. Liu, editors. Math Horizons 10(2):32-35, 2002.

S. Wolfram. The Mathematica Book (5th edition). Wolfram Media, 2003.