There are two tasks related to preparation of plots of functions: first, to produce the numbers required for a plot, and second, to draw a plot with axes, symbols, a legend, perhaps additional illustrations and so on. Here we only concern ourselves with the first task, that of preparation of the numerical data for a plot. There are many plotting programs that can read a file with numbers and plot it in any desired manner.
Generating data for plots of functions generally does not require high-precision calculations. However, we need an algorithm that can be adjusted to produce data to different levels of precision. In some particularly ill-behaved cases, a precise plot will not be possible and we would not want to waste time producing data that is too accurate for what it is worth.
A simple approach to plotting would be to divide the interval into many equal subintervals and to evaluate the function on the resulting grid. Precision of the plot can be adjusted by choosing a larger or a smaller number of points.
However, this approach is not optimal. Sometimes a function changes rapidly near one point but slowly everywhere else. For example, f(x)=1/x changes very quickly at small x. Suppose we need to plot this function between 0 and 100. It would be wasteful to use the same subdivision interval everywhere: a finer grid is only required over a small portion of the plotting range near x=0.
The adaptive plotting routine Plot2D'adaptive uses a simple algorithm to select the optimal grid to approximate a function of one argument f(x). The algorithm repeatedly subdivides the grid intervals near points where the existing grid does not represent the function well enough. A similar algorithm for adaptive grid refinement could be used for numerical integration. The idea is that plotting and numerical integration require the same kind of detailed knowledge about the behavior of the function.
The algorithm first splits the interval into a specified initial number of equal subintervals, and then repeatedly splits each subinterval in half until the function is well enough approximated by the resulting grid. The integer parameter depth gives the maximum number of binary splittings for a given initial interval; thus, at most 2^depth additional grid points will be generated. The function Plot2D'adaptive should return a list of pairs of points {{x1,y1}, {x2,y2}, ...} to be used directly for plotting.
The adaptive plotting algorithm works like this:
This algorithm works well if the initial number of points and the depth parameter are large enough. These parameters can be adjusted to balance the available computing time and the desired level of detail in the resulting plot.
Singularities in the function are handled by the step 3. Namely, the change in the sequence a, a[1], b, b[1], c is always considered to be "too rapid" if one of these values is a non-number (e.g. Infinity or Undefined). Thus, the interval immediately adjacent to a singularity will be plotted at the highest allowed refinement level. When preparing the plotting data, the singular points are simply not printed to the data file, so that a plotting programs does not encounter any problems.
The coefficients c[k] for grids with a constant step h can be found, for example, by solving the following system of equations,
In the same way it is possible to find quadratures for the integral over a subinterval rather than over the whole interval of x. In the current implementation of the adaptive plotting algorithm, two quadratures are used: the 3-point quadrature ( n=2) and the 4-point quadrature ( n=3) for the integral over the first subinterval, Integrate(x,a[0],a[1])f(x). Their coefficients are (5/12, 2/3, -1/12) and ( 3/8, 19/24, -5/24, 1/24). An example of using the first of these subinterval quadratures would be the approximation
The task of surface plotting is to obtain a picture of a two-dimensional surface as if it were a solid object in three dimensions. A graphical representation of a surface is a complicated task. Sometimes it is required to use particular coordinates or projections, to colorize the surface, to remove hidden lines and so on. We shall only be concerned with the task of obtaining the data for a plot from a given function of two variables f(x,y). Specialized programs can take a text file with the data and let the user interactively produce a variety of surface plots.
The currently implemented algorithm in the function Plot3DS is very similar to the adaptive plotting algorithm for two-dimensional plots. A given rectangular plotting region a[1]<=x<=a[2], b[1]<=y<=b[2] is subdivided to produce an equally spaced rectangular grid of points. This is the initial grid which will be adaptively refined where necessary. The refinement algorithm will divide a given rectangle in four quarters if the available function values indicate that the function does not change smoothly enough on that rectangle.
The criterion of a "smooth enough" change is very similar to the procedure outlined in the previous section. The change is "smooth enough" if all points are finite, nonsingular values, and if the integral of the function over the rectangle is sufficiently well approximated by a certain low-order "cubature" formula.
The two-dimensional integral of the function is estimated using the following 5-point Newton-Cotes cubature:
1/12 0 1/12 0 2/3 0 1/12 0 1/12 |
An example of using this cubature would be the approximation
Similarly, an 8-point cubature with zero sum is used to estimate the error:
-1/3 2/3 1/6 -1/6 -2/3 -1/2 1/2 0 1/3 |
One minor problem with adaptive surface plotting is that the resulting set of points may not correspond to a rectangular grid in the parameter space (x, y). This is because some rectangles from the initial grid will need to be bisected more times than others. So, unless adaptive refinement is disabled, the function Plot3DS produces a somewhat disordered set of points. However, most surface plotting programs require that the set of data points be a rectangular grid in the parameter space. So a smoothing and interpolation procedure is necessary to convert a non-gridded set of data points ("scattered" data) to a gridded set.
The program gnuplot has this facility (the "set dgrid3d" command), although its implementation may not be optimal for all purposes. The current solution is to use the dgrid3d command and to generate a grid at the highest level of bisection ever used during the adaptive refinement process. This quite possibly generates many more grid points than necessary, but a more optimal solution would be more much more time-consuming or would require a specialized external program. One such set of programs is the GMT ("generic mapping tools") utility suite.
A two-dimensional parametric plot is a line in a two-dimensional space, defined by two equations such as x=f(t), y=g(t). Two functions f, g and a range of the independent variable t, for example, t[1]<=t<=t[2], need to be specified.
Parametric plots can be used to represent plots of functions in non-Euclidean coordinates. For example, to plot the function rho=Cos(4*phi)^2 in polar coordinates ( rho, phi), one can rewrite the Euclidean coordinates through the polar coordinates, x=rho*Cos(phi), y=rho*Sin(phi), and use the equivalent parametric plot with phi as the parameter: x=Cos(4*phi)^2*Cos(phi), y=Cos(4*phi)^2*Sin(phi).
Sometimes higher-dimensional parametric plots are required. A line plot in three dimensions is defined by three functions of one variable, for example, x=f(t), y=g(t), z=h(t), and a range of the parameter t. A surface plot in three dimensions is defined by three functions of two variables each, for example, x=f(u,v), y=g(u,v), z=h(u,v), and a rectangular domain in the (u, v) space.
The data for parametric plots can be generated separately using the same adaptive plotting algorithms as for ordinary function plots, as if all functions such as f(t) or g(u,v) were unrelated functions. The result would be several separate data sets for the x, y, ... coordinates. These data sets could then be combined using an interactive plotting program.
A different question is whether a CAS really needs to be able to evaluate, say, 10,000 digits of the value of a Bessel function of some 10,000-digit complex argument. It seems likely that no applied problem of natural sciences would need floating-point computations of special functions with such a high precision. However, arbitrary-precision computations are certainly useful in some mathematical applications; e.g. some mathematical identities can be first guessed by a floating-point computation with many digits and then proved.
Very high precision computations of special functions might be useful in the future. But it is already quite clear that computations with moderately high precision (say, 50 or 100 decimal digits) are useful for applied problems. For example, to obtain the leading asymptotic of an analytic function, we could expand it in series and take the first term. But we need to check that the coefficient at what we think is the leading term of the series does not vanish. This coefficient could be a certain "exact" number such as (Cos(355)+1)^2. This number is "exact" in the sense that it is made of integers and elementary functions. But we cannot say a priori that this number is nonzero. The problem of "zero determination" (finding out whether a certain "exact" number is zero) is known to be algorithmically unsolvable if we allow transcendental functions. The only practical general approach seems to be to compute the number in question with many digits. Usually a few digits are enough, but occasionally several hundred digits are needed.
Implementing an efficient algorithm that computes 100 digits of Sin(3/7) already involves many of the issues that would also be relevant for a 10,000 digit computation. Modern algorithms allow evaluations of all elementary functions in time that is asymptotically logarithmic in the number of digits P and linear in the cost of long multiplication (usually denoted M(P)). Almost all special functions can be evaluated in time that is asymptotically linear in P and in M(P). (However, this asymptotic cost sometimes applies only to very high precision, e.g., P>1000, and different algorithms need to be implemented for calculations in lower precision.)
In Yacas we strive to implement all numerical functions to arbitrary precision. All integer or rational functions return exact results, and all floating-point functions return their value with P correct decimal digits (assuming sufficient precision of the arguments). The current value of P is accessed as GetPrecision() and may be changed by Precision(...).
Implementing an arbitrary-precision floating-point computation of a function f(x), such as f(x)=Exp(x), typically needs the following:
In calculations with machine precision where the number of digits is fixed, the problem of round-off errors is quite prominent. Every arithmetic operation causes a small loss of precision; as a result, a few last digits of the final value are usually incorrect. But if we have an arbitrary precision capability, we can always increase precision by a few more digits during intermediate computations and thus eliminate all round-off error in the final result. We should, of course, take care not to increase the working precision unnecessarily, because any increase of precision means slower calculations. Taking twice as many digits as needed and hoping that the result is precise is not a good solution.
Selecting algorithms for computations is the most non-trivial part of the implementation. We want to achieve arbitrarily high precision, so we need to find either a series, or a continued fraction, or a sequence given by explicit formula, that converges to the function in a controlled way. It is not enough to use a table of precomputed values or a fixed approximation formula that has a limited precision.
In the last 30 years, the interest in arbitrary-precision computations grew and many efficient algorithms for elementary and special functions were published. Most algorithms are iterative. Almost always it is very important to know in advance how many iterations are needed for given x, P. This knowledge allows to estimate the computational cost, in terms of the required precision P and of the cost of long multiplication M(P), and choose the best algorithm.
Typically all operations will fall into one of the following categories (sorted by the increasing cost):
The cost of long multiplication M(P) is between O(P^2) for low precision and O(P*Ln(P)) for very high precision. In some cases, a different algorithm should be chosen if the precision is high enough to allow M(P) faster than O(P^2).
Some algorithms also need storage space (e.g. an efficient algorithm for summation of the Taylor series uses O(Ln(P)) temporary P-digit numbers).
Below we shall normally denote by P the required number of decimal digits. The formulae frequently contain conspicuous factors of Ln(10), so it will be clear how to obtain analogous expressions for another base. (Most implementations use a binary base rather than a decimal base since it is more convenient for many calculations.)
Suppose we truncate the series after n-th term and the series converges "well enough" after that term. Then the error will be approximately equal to the first term we dropped. (This is what we really mean by "converges well enough" and this will generally be the case in all applications, because we would not want to use a series that does not converge well enough.)
The term we dropped is x^(n+1)/(n+1)!. To estimate n! for large n, one can use the inequality
If we use the upper bound on n! from this estimate, we find that the term we dropped is bounded by
We can try to guess the result. The largest term on the LHS grows as n0*Ln(n0) and it should be approximately equal to P*Ln(10); but Ln(n0) grows very slowly, so this gives us a hint that n0 is proportional to P*Ln(10). As a first try, we set n0=P*Ln(10)-2 and compare the RHS with the LHS; we find that we have overshot by a factor Ln(P)-1+Ln(Ln(10)), which is not a large factor. We can now compensate and divide n0 by this factor, so our second try is
Our final result is that it is enough to take
Here is a simple estimate of the normal round-off error in a computation of n terms of a power series. Suppose that the sum of the series is of order 1, that the terms monotonically decrease in magnitude, and that adding one term requires two multiplications and one addition. If all calculations are performed with absolute precision epsilon=10^(-P), then the total accumulated round-off error is 3*n*epsilon. If the relative error is 3*n*epsilon, it means that our answer is something like a*(1+3*n*epsilon) where a is the correct answer. We can see that out of the total P digits of this answer, only the first k decimal digits are correct, where k= -Ln(3*n*epsilon)/Ln(10). In other words, we have lost
This estimate assumes several things about the series (basically, that the series is "well-behaved"). These assumptions must be verified in each particular case. For example, if the series begins with some large terms but converges to a very small value, this estimate is wrong (see the next subsection).
In the previous exercise we found the number of terms n for Exp(x). So now we know how many extra digits of working precision we need for this particular case.
Below we shall have to perform similar estimates of the required number of terms and of the accumulated round-off error in our analysis of the algorithms.
Consider the computation of Sin(x) by the truncated Taylor series
First, we determine the necessary number of terms N. The magnitude of the sum is never larger than 1. Therefore we need the N-th term of the series to be smaller than 10^(-P). The inequality is (2*N+1)! >10^(P+M*(2*N+1)). We obtain that 2*N+2>e*10^M is a necessary condition, and if P is large, we find approximately
However, taking enough terms does not yet guarantee a good result. The terms of the series grow at first and then start to decrease. The sum of these terms is, however, small. Therefore there is some cancellation and we need to increase the working precision to avoid the round-off. Let us estimate the required working precision.
We need to find the magnitude of the largest term of the series. The ratio of the next term to the previous term is x/(2*k*(2*k+1)) and therefore the maximum will be when this ratio becomes equal to 1, i.e. for 2*k<=>Sqrt(x). Therefore the largest term is of order x^Sqrt(x)/Sqrt(x)! and so we need about M/2*Sqrt(x) decimal digits before the decimal point to represent this term. But we also need to keep at least P digits after the decimal point, or else the round-off error will erase the significant digits of the result. In addition, we will have unavoidable round-off error due to O(P) arithmetic operations. So we should increase precision again by P+Ln(P)/Ln(10) digits plus a few guard digits.
As an example, to compute Sin(10) to P=50 decimal digits with this method, we need a working precision of about 60 digits, while to compute Sin(10000) we need to work with about 260 digits. This shows how inefficient the Taylor series for Sin(x) becomes for large arguments x. A simple transformation x=2*Pi*n+x' would have reduced x to at most 7, and the unnecessary computations with 260 digits would be avoided. The main cause of this inefficiency is that we have to add and subtract extremely large numbers to get a relatively small result of order 1.
We find that the method of Taylor series for Sin(x) at large x is highly inefficient because of round-off error and should be complemented by other methods. This situation seems to be typical for Taylor series.
The algorithms for basic arithmetic in the internal math version are currently rather slow compared with gmp. If P is the number of digits of precision, then multiplication and division take M(P)=O(P^2) operations in the internal math. (Of course, multiplication and division by a short integer takes time linear in P.) Much faster algorithms (Karatsuba, Toom-Cook, FFT multiplication, Newton-Raphson division etc.) are implemented in gmp, CLN and some other libraries. The asymptotic cost of multiplication for very large precision is M(P)<=>O(P^1.6) for the Karatsuba method and M(P)=O(P*Ln(P)*Ln(Ln(P))) for the FFT method. In the estimates of computation cost in this book we shall assume that M(P) is at least linear in P and maybe a bit slower.
The costs of multiplication may be different in various arbitrary-precision arithmetic libraries and on different computer platforms. As a rough guide, one can assume that the straightforward O(P^2) multiplication is good until 100-200 decimal digits, the asymptotically fastest method of FFT multiplication is good at the precision of about 5,000 or more decimal digits, and the Karatsuba multiplication is best in the middle range.
Warning: calculations with internal Yacas math using precision exceeding 10,000 digits are currently impractically slow.
In some algorithms it is necessary to compute the integer parts of expressions such as a*Ln(b)/Ln(10) or a*Ln(10)/Ln(2), where a, b are short integers of order O(P). Such expressions are frequently needed to estimate the number of terms in the Taylor series or similar parameters of the algorithms. In these cases, it is important that the result is not underestimated. However, it would be wasteful to compute 1000*Ln(10)/Ln(2) in great precision only to discard most of that information by taking the integer part of that number. It is more efficient to approximate such constants from above by short rational numbers, for example, Ln(10)/Ln(2)<28738/8651 and Ln(2)<7050/10171. The error of such an approximation will be small enough for practical purposes. The function BracketRational can be used to find optimal rational approximations.
The function IntLog (see below) efficiently computes the integer part of a logarithm (for an integer base, not a natural logarithm). If more precision is desired in calculating Ln(a)/Ln(b) for integer a, b, one can compute IntLog(a^k,b) for some integer k and then divide by k.
The exponent E is easy to obtain:
Once we found E, we can write x=10^(E+m) where m=Exp(1000)/Ln(10)-E is a floating-point number, 0<m<1. Then M=10^m. To find M with P (decimal) digits, we need m with also at least P digits. Therefore, we actually need to evaluate Exp(1000)/Ln(10) with 434+P decimal digits before we can find P digits of the mantissa of x. We ran into a perhaps surprising situation: one needs a high-precision calculation even to find the first digit of x, because it is necessary to find the exponent E exactly as an integer, and E is a rather large integer. A normal double-precision numerical calculation would give an overflow error at this point.
Suppose we have found the number x=Exp(Exp(1000)) with some precision. What about finding Sin(x)? Now, this is extremely difficult, because to find even the first digit of Sin(x) we have to evaluate x with absolute error of at most 0.5. We know, however, that the number x has approximately 10^434 digits before the decimal point. Therefore, we would need to calculate x with at least that many digits. Computations with 10^434 digits is clearly far beyond the capability of modern computers. It seems unlikely that even the sign of Sin(Exp(Exp(1000))) will be obtained in the near future.
Suppose that N is the largest integer that our arbitrary-precision facility can reasonably handle. (For Yacas internal math library, N is about 10^10000.) Then it follows that numbers x of order 10^N can be calculated with at most one (1) digit of floating-point precision, while larger numbers cannot be calculated with any precision at all.
It seems that very large numbers can be obtained in practice only through exponentiation or powers. It is unlikely that such numbers will arise from sums or products of reasonably-sized numbers in some formula.
If numbers larger than 10^N are created only by exponentiation operations, then special exponential notation could be used to represent them. For example, a very large number z could be stored and manipulated as an unevaluated exponential z=Exp(M*10^E) where M>=1 is a floating-point number with P digits of mantissa and E is an integer, Ln(N)<E<N. Let us call such objects "exponentially large numbers" or "exp-numbers" for short.
In practice, we should decide on a threshold value N and promote a number to an exp-number when its logarithm exceeds N.
Note that an exp-number z might be positive or negative, e.g. z= -Exp(M*10^E).
Arithmetic operations can be applied to the exp-numbers. However, exp-large arithmetic is of limited use because of an almost certainly critical loss of precision. The power and logarithm operations can be meaningfully performed on exp-numbers z. For example, if z=Exp(M*10^E) and p is a normal floating-point number, then z^p=Exp(p*M*10^E) and Ln(z)=M*10^E. We can also multiply or divide two exp-numbers. But it makes no sense to multiply an exp-number z by a normal number because we cannot represent the difference between z and say 2.52*z. Similarly, adding z to anything else would result in a total underflow, since we do not actually know a single digit of the decimal representation of z. So if z1 and z2 are exp-numbers, then z1+z2 is simply equal to either z1 or z2 depending on which of them is larger.
We find that an exp-number z acts as an effective "infinity" compared with normal numbers. But exp-numbers cannot be used as a device for computing limits: the unavoidable underflow will almost certainly produce wrong results. For example, trying to verify
Taking a logarithm of an exp-number brings it back to the realm of normal, representable numbers. However, taking an exponential of an exp-number results in a number which is not representable even as an exp-number. This is because an exp-number z needs to have its exponent E represented exactly as an integer, but Exp(z) has an exponent of order O(z) which is not a representable number. The monstrous number Exp(z) could be only written as Exp(Exp(M*10^E)), i.e. as a "doubly exponentially large" number, or "2-exp-number" for short. Thus we obtain a hierarchy of iterated exp-numbers. Each layer is "unrepresentably larger" than the previous one.
The same considerations apply to very small numbers of the order 10^(-N) or smaller. Such numbers can be manipulated as "exponentially small numbers", i.e. expressions of the form Exp(-M*10^E) with floating-point mantissa M>=1 and integer E satisfying Ln(N)<E<N. Exponentially small numbers act as an effective zero compared with normal numbers.
Taking a logarithm of an exp-small number makes it again a normal representable number. However, taking an exponential of an exp-small number produces 1 because of underflow. To obtain a "doubly exponentially small" number, we need to take a reciprocal of a doubly exponentially large number, or take the exponent of an exponentially large negative power. In other words, Exp(-M*10^E) is exp-small, while Exp(-Exp(M*10^E)) is 2-exp-small.
The practical significance of exp-numbers is rather limited. We cannot obtain even a single significant digit of an exp-number. A "computation" with exp-numbers is essentially a floating-point computation with logarithms of these exp-numbers. A practical problem that needs numbers of this magnitude can probably be restated in terms of more manageable logarithms of such numbers. In practice, exp-numbers could be useful not as a means to get a numerical answer, but as a warning sign of critical overflow or underflow.
Usually one considers infinite continued fractions, i.e. the sequences a[i], b[i] are infinite. The value of an infinite continued fraction is defined as the limit of the fraction truncated after a very large number of terms. (A continued traction can be truncated after n-th term if one replaces b[n] by 0.)
An infinite continued fraction does not always converge. Convergence depends on the values of the terms.
The representation of a number via a continued fraction is not unique because we could, for example, multiply the numerator and the denominator of any simple fraction inside it by any number. Therefore one may consider some normalized representations. A continued fraction is called "regular" if b[k]=1 for all k, all a[k] are integer and a[k]>0 for k>=1. Regular continued fractions always converge.
The algorithm for converting a rational number r=n/m into a regular continued fraction is simple. First, we determine the integer part of r, which is Div(n,m). If it is negative, we need to subtract one, so that r=n[0]+x and the remainder x is nonnegative and less than 1. The remainder x=Mod(n,m)/m is then inverted, r[1]:=1/x=m/Mod(n,m) and so we have completed the first step in the decomposition, r=n[0]+1/r[1]; now n[0] is integer but r[1] is perhaps not integer. We repeat the same procedure on r[1], obtain the next integer term n[1] and the remainder r[2] and so on, until such n that r[n] is an integer and there is no more work to do. This process will always terminate.
If r is a real number which is known by its floating-point representation at some precision, then we can use the same algorithm because all floating-point values are actually rational numbers.
Real numbers known by their exact representations can sometimes be expressed as infinite continued fractions, for example
The functions GuessRational and NearRational take a real number x and use continued fractions to find rational approximations r=p/q<=>x with "optimal" (small) numerators and denominators p, q.
Suppose we know that a certain number x is rational but we have only a floating-point representation of x with a limited precision, for example, x<=>1.5662650602409638. We would like to guess a rational form for x (in this example x=130/83). The function GuessRational solves this problem.
Consider the following example. The number 17/3 has a continued fraction expansion {5,1,2}. Evaluated as a floating point number with limited precision, it may become something like 17/3+0.00001, where the small number represents a round-off error. The continued fraction expansion of this number is {5, 1, 2, 11110, 1, 5, 1, 3, 2777, 2}. The presence of an unnaturally large term 11110 clearly signifies the place where the floating-point error was introduced; all terms following it should be discarded to recover the continued fraction {5,1,2} and from it the initial number 17/3.
If a continued fraction for a number x is cut right before an unusually large term and evaluated, the resulting rational number has a small denominator and is very close to x. This works because partial continued fractions provide "optimal" rational approximations for the final (irrational) number, and because the magnitude of the terms of the partial fraction is related to the magnitude of the denominator of the resulting rational approximation.
GuessRational(x, prec) needs to choose the place where it should cut the continued fraction. The algorithm for this is somewhat heuristic but it works well enough. The idea is to cut the continued fraction when adding one more term would change the result by less than the specified precision. To realize this in practice, we need an estimate of how much a continued fraction changes when we add one term.
The routine GuessRational uses a (somewhat weak) upper bound for the difference of continued fractions that differ only by an additional last term:
The above estimate for delta hinges on the inequality
This algorithm works well if x is computed with enough precision; namely, it must be computed to at least as many digits as there are in the numerator and the denominator of the fraction combined. Also, the parameter prec should not be too large (or else the algorithm will find another rational number with a larger denominator that approximates x "better" than the precision to which you know x).
The related function NearRational(x, prec) works somewhat differently. The goal is to find an "optimal" rational number, i.e. with smallest numerator and denominator, that is within the distance 10^(-prec) of a given value x. The function NearRational does not always give the same answer as GuessRational.
The algorithm for NearRational comes from the HAKMEM [Beeler et al. 1972], Item 101C. Their description is terse but clear:
Problem: Given an interval, find in it the rational number with the smallest numerator and denominator. Solution: Express the endpoints as continued fractions. Find the first term where they differ and add 1 to the lesser term, unless it's last. Discard the terms to the right. What's left is the continued fraction for the "smallest" rational in the interval. (If one fraction terminates but matches the other as far as it goes, append an infinity and proceed as above.) |
The HAKMEM text [Beeler et al. 1972] contains several interesting insights relevant to continued fractions and other numerical algorithms.
In this section we describe some methods for computing general continued fractions and for estimating the number of terms needed to achieve a given precision.
Let us introduce some notation. A continued fraction
This method requires one long division at each step. There may be significant round-off error if a[m] and b[m] have opposite signs, but otherwise the round-off error is very small because a convergent continued fraction is not sensitive to small changes in its terms.
The idea is that the starting value of the backward recurrence should be chosen not as a[n] but as another number that more closely approximates the infinite remainder of the fraction. The infinite remainder, which we can symbolically write as F[n][Infinity], can be sometimes estimated analytically (obviously, we are unable to compute the remainder exactly). In simple cases, F[n][Infinity] changes very slowly at large n (warning: this is not always true and needs to be verified in each particular case!). Suppose that F[n][Infinity] is approximately constant; then it must be approximately equal to F[n+1][Infinity]. Therefore, if we solve the (quadratic) equation
We may use more terms of the original continued fraction starting from a[n] and obtain a more precise estimate of the remainder. In each case we will only have to solve one quadratic equation.
The "top-down" method is slower but provides an automatic error estimate and can be used to evaluate a continued fraction with more and more terms until the desired precision is achieved. The idea
The formula for f[k] is the following. First the auxiliary sequence P[k], Q[k] for k>=1 needs to be defined by P[1]=0, Q[1]=1, and P[k+1]:=b[k]*Q[k], Q[k+1]:=P[k]+a[k]*Q[k]. Then define f[0]:=a[0] and
Evaluating the next element f[k] requires four long multiplications and one division. This is significantly slower, compared with just one long division or two long multiplications in the bottom-up method. Therefore it is desirable to have an a priori estimate of the convergence rate and to be able to choose the number of terms before the computation. Below we shall consider some examples when the formula for f[k] allows to estimate the required number of terms analytically.
The bottom-up methods are simpler and faster than the top-down methods but require to know the number of terms in advance. In many cases the required number of terms can be estimated analytically, and then the bottom-up methods are always preferable. But in some cases the convergence analysis is very complicated.
The plain bottom-up method requires one long division at each step, while the bottom-up method with two recurrences requires two long multiplications at each step. Since the time needed for a long division is usually about four times that for a long multiplication (e.g. when the division is computed by Newton's method), the second variation of the bottom-up method is normally faster.
The estimate of the remainder improves the convergence of the bottom-up method and should always be used if available.
If an estimate of the number of terms is not possible, the top-down methods should be used, looping until the running error estimate shows enough precision. This incurs a performance penalty of at least 100% and at most 300%. The top-down method with two steps at once should be used only when the formula for f[k] results in alternating signs.
Usually, a continued fraction representation of a function will converge geometrically or slower, i.e. at least O(P) terms are needed for a precision of P digits. If a geometrically convergent Taylor series representation is also available, the continued fraction method will be slower because it requires at least as many or more long multiplications per term. Also, in most cases the Taylor series can be computed much more efficiently using the rectangular scheme. (See, e.g., the section on ArcTan(x) for a more detailed consideration.)
However, there are some functions for which a Taylor series is not easily computable or does not converge but a continued fraction is available. For example, the incomplete Gamma function and related functions can be computed using continued fractions in some domains of their arguments.
So far we have reduced the difference between F[m][n+1] and F[m][n] to a similar difference on the next level m+1 instead of m; i.e. we can increment m but keep n fixed. We can apply this formula to F[0][n+1]-F[0][n], i.e. for m=0, and continue applying the same recurrence relation until m reaches n. The result is
Now the problem is to simplify the two long products in the denominator. We notice that F[1][n] has F[2][n] in the denominator, and therefore F[1][n]*F[2][n]=F[2][n]*a[1]+b[1]. The next product is F[1][n]*F[2][n]*F[3][n] and it simplifies to a linear function of F[3][n], namely F[1][n]*F[2][n]*F[3][n] = (b[1]+a[1]*a[2])*F[3][n]+b[1]*a[2]. So we can see that there is a general formula
Having found the coefficients P[k], Q[k], we can now rewrite the long products in the denominator, e.g.
For example, the continued fraction
There are some cases when a continued fraction representation is efficient. The complementary error function Erfc(x) can be computed using the continued fraction due to Laplace (e.g. [Thacher 1963]),
The error function is a particular case of the incomplete Gamma function
Suppose we are given the terms a[k], b[k] that define an infinite continued fraction, and we need to estimate its convergence rate. We have to find the number of terms n for which the error of approximation is less than a given epsilon. In our notation, we need to solve Abs(f[n+1])<epsilon for n.
The formula that we derived for f[n+1] gives an error estimate for the continued fraction truncated at the n-th term. But this formula contains the numbers Q[n] in the denominator. The main problem is to find how quickly the sequence Q[n] grows. The recurrence relation for this sequence can be rewritten as
We have used this bound to estimate the relative error of the continued fraction expansion for ArcTan(x) at small x (elsewhere in this book). However, we found that at large x this bound becomes greater than 1. This does not mean that the continued fraction does not converge and cannot be used to compute ArcTan(x) when x>1, but merely indicates that the "simple bound" is too weak. The sequence Q[n] actually grows faster than the product of all a[k] and we need a tighter bound on this growth. In many cases such a bound can be obtained by the method of generating functions.
The asymptotic growth of the sequence Q[n] can be estimated by the method of steepest descent, also known as Laplace's method. (See, e.g., [Olver 1974], ch. 3, sec. 7.5.) This method is somewhat complicated but quite powerful. The method requires that we find an integral representation for Q[n] (usually a contour integral in the complex plane). Then we can convert the integral into an asymptotic series in k^(-1).
Along with the general presentation of the method, we shall consider an example when the convergence rate can be obtained analytically. The example is the representation of the complementary error function Erfc(x),
The "simple bound" would give Abs(f[n+1])<=v^n*n! and this expression grows with n. But we know that the above continued fraction actually converges for any v, so f[n+1] must tend to zero for large n. It seems that the "simple bound" is not strong enough for any v and we need a better bound.
An integral representation for Q[n] can be obtained using the method of generating functions. Consider a function G(s) defined by the infinite series
Note that the above series for the function G(s) may or may not converge for any given s; we shall manipulate G(s) as a formal power series until we obtain an explicit representation. What we really need is an analytic continuation of G(s) to the complex s.
It is generally the case that if we know a simple linear recurrence relation for a sequence, then we can also easily find its generating function. The generating function will satisfy a linear differential equation. To guess this equation, we write down the series for G(s) and its derivative G'(s) and try to find their linear combination which is identically zero because of the recurrence relation. (There is, of course, a computer algebra algorithm for doing this automatically.)
Taking the derivative G'(s) produces the forward-shifted series
In the case of our sequence Q[n] above, the recurrence relation is
The second step is to obtain an integral representation for Q[n], so that we could use the method of steepest descents and find its asymptotic at large n.
In our notation Q[n+1] is equal to the n-th derivative of the generating function at s=0:
There are two ways to proceed. One is to obtain an integral representation for G(s), for instance
The second possibility is to express Q[n] as a contour integral in the complex plane around s=0 in the counter-clockwise direction:
In the particular case of the continued fraction for Erfc(x), the calculations are somewhat easier if Re(v)>0 (where v:=1/(2*x^2)). Full details are given in a separate section. The result for Re(v)>0 is
Note that this is not merely a bound but an actual asymptotic estimate of f[n+1]. (Stirling's formula can also be derived using the method of steepest descent from an integral representation of the Gamma function, in a similar way.)
Defined as above, the value of f[n+1] is in general a complex number. The absolute value of f[n+1] can be found using the formula
When Re(v)<=0, the same formula can be used (this can be shown by a more careful consideration of the branches of the square roots). The continued fraction does not converge when Re(v)<0 and Im(v)=0 (i.e. for pure imaginary x). This can be seen from the above formula: in that case Re(v)= -Abs(v) and Abs(f[n+1]) does not decrease when n grows.
These estimates show that the error of the continued fraction approximation to Erfc(x) (when it converges) decreases with n slower than in a geometric progression. This means that we need to take O(P^2) terms to get P digits of precision.
To use the method of steepest descent, we represent the integrand as an exponential of some function g(t,n) and find "stationary points" where this function has local maxima:
We only need to consider very large values of n, so we can neglect terms of order O(1/Sqrt(n)) or smaller. We find that, in our case, two peaks of Re(g) occur at approximately t1<=> -1/2+Sqrt(n*v) and t2<=> -1/2-Sqrt(n*v). We assume that n is large enough so that n*v>1/2. Then the first peak is at a positive t and the second peak is at a negative t. The contribution of the peaks can be computed from the Taylor approximation of g(t,n) near the peaks. We can expand, for example,
Then we obtain the estimate
Usually one of the stationary points has the largest value of Re(g(s)); this is the dominant stationary point. If s0 is the dominant stationary point and g2=(Deriv(s,2)g(s0)) is the second derivative of g at that point, then the asymptotic of the integral is
We have to choose a new contour and check the convergence of the resulting integral separately. In each case we may need to isolate the singularities of G(s) or to find the regions of infinity where G(s) quickly decays (so that the infinite parts of the contour might be moved there). There is no prescription that works for all functions G(s).
Let us return to our example. For G(s)=Exp(s+(v*s^2)/2), the function g(s) has no singularities except the pole at s=0. There are two stationary points located at the (complex) roots s1, s2 of the quadratic equation v*s^2+s-n=0. Note that v is an arbitrary (nonzero) complex number. We now need to find which of the two stationary points gives the dominant contribution. By comparing Re(g(s1)) and Re(g(s2)) we find that the point s with the largest real part gives the dominant contribution. However, if Re(s1)=Re(s2) (this happens only if v is real and v<0, i.e. if x is pure imaginary), then both stationary points contribute equally. Barring that possibility, we find (with the usual definition of the complex square root) that the dominant contribution for large n is from the stationary point at
This formula agrees with the asymptotic for Q[n+1] found above for real v>0, when we use Stirling's formula for (n-1)!:
The treatment for Re(v)<0 is similar.
Newton's method sometimes suffers from a sensitivity to the initial guess. If the initial value x[0] is not chosen sufficiently close to the root, the iterations may converge very slowly or not converge at all. To remedy this, one can combine Newton's iteration with simple bisection. Once the root is bracketed inside an interval (a, b), one checks whether (a+b)/2 is a better approximation for the root than that obtained from Newton's iteration. This guarantees at least linear convergence in the worst case.
For some equations f(x)=0, Newton's method converges faster than quadratically. For example, solving Sin(x)=0 in the neighborhood of x=3.14159 gives "cubic" convergence, i.e. the number of correct digits is tripled at each step. This happens because Sin(x) near its root x=Pi has a vanishing second derivative and thus the function is particularly well approximated by a straight line.
Halley's method can be generalized to any function f(x). A cubically convergent iteration is always obtained if we replace the equation f(x)=0 by an equivalent equation
The Halley iteration for the equation f(x)=0 can be written as
Halley's iteration, despite its faster convergence rate, may be more cumbersome to evaluate than Newton's iteration and so it may not provide a more efficient numerical method for a given function. Only in some special cases is Halley's iteration just as simple to compute as Newton's iteration.
Halley's method is sometimes less sensitive to the choice of the initial point x[0]. An extreme example of sensitivity to the initial point is the equation x^(-2)=12 for which Newton's iteration x'=3/2*x-6*x^3 converges to the root only from initial points 0<x[0]<0.5 and wildly diverges otherwise, while Halley's iteration converges to the root from any x[0]>0.
It is at any rate not true that Halley's method always converges better than Newton's method. For instance, it diverges on the equation 2*Cos(x)=x unless started at x[0] within the interval (-1/6*Pi, 7/6*Pi). Another example is the equation Ln(x)=a. This equation allows to compute x=Exp(a) if a fast method for computing Ln(x) is available (e.g. the AGM-based method). For this equation, Newton's iteration
When it converges, Halley's iteration can still converge very slowly for certain functions f(x), for example, for f(x)=x^n-a if n^n>a. For such functions that have very large and rapidly changing derivatives, no general method can converge faster than linearly. In other words, a simple bisection will generally do just as well as any sophisticated iteration, until the root is approximated very precisely. Halley's iteration combined with bisection seems to be a good choice for such problems.
In the above examples, y is a small quantity and the series represents corrections to the initial value x, therefore the order of convergence is equal to the first discarded order of y in the series.
These simple constructions are possible because the functions satisfy simple identities, such as Exp(a+b)=Exp(a)*Exp(b) or Sqrt(a*b)=Sqrt(a)*Sqrt(b). For other functions the formulae quickly become very complicated and unsuitable for practical computations.
For practical evaluation, iterations must be supplemented with "quality control". For example, if x0 and x1 are two consecutive approximations that are already very close, we can quickly compute the achieved (relative) precision by finding the number of leading zeros in the number
Suppose x is an approximation that is correct to P digits; then we expect the quantity x' to be correct to 2*P digits. Therefore we should perform calculations in the first formula with 2*P digits; this means three long multiplications, 3*M(2*P). Now consider the calculation in the second formula. First, the quantity y:=1-a*x^2 is computed using two 2*P-digit multiplications.
The advantage is even greater with higher-order methods. For example, a fourth-order iteration for the inverse square root can be written as
The asymptotic cost of finding the root x of the equation f(x)=0 with P digits of precision is usually the same as the cost of computing f(x) [Brent 1975]. The main argument can be summarized by the following simple example. To get the result to P digits, we need O(Ln(P)) Newton's iterations. At each iteration we shall have to compute the function f(x) to a certain number of digits. Suppose that we start with one correct digit and that each iteration costs us c*M(2*P) operations where c is a given constant, while the number of correct digits grows from P to 2*P. Then the total cost of k iterations is
Thus, if we have a fast method of computing, say, ArcTan(x), then we immediately obtain a method of computing Tan(x) which is asymptotically as fast (up to a constant).
Increasing the order by 1 costs us comparatively little, and we may change the order k at any time. Is there a particular order k that gives the smallest computational cost and should be used for all iterations, or the order needs to be adjusted during the computation? A natural question is to find the optimal computational strategy.
It is difficult to fully analyze this question, but it seems that choosing a particular order k for all iterations is close to the optimal strategy.
A general "strategy" is a set of orders S(P,P[0])=(k[1], k[2], ..., k[n]) to be chosen at the first, second, ..., n-th iteration, given the initial precision P[0] and the required final precision P. At each iteration, the precision will be multiplied by the factor k[i]. The optimal strategy S(P,P[0]) is a certain function of P[0] and P such that the required precision is reached, i.e.
If we assume that the cost of multiplication M(P) is proportional to some power of P, for instance M(P)=P^mu, then the cost of each iteration and the total cost are homogeneous functions of P and P[0]. Therefore the optimal strategy is a function only of the ratio P/P[0]. We can multiply both P[0] and P by a constant factor and the optimal strategy will remain the same. We can denote the optimal strategy S(P/P[0]).
We can check whether it is better to use several iterations at smaller orders instead of one iteration at a large order. Suppose that M(P)=P^mu, the initial precision is 1 digit, and the final precision P=k^n. We can use either n iterations of the order k or 1 iteration of the order P. The cost of one iteration of order P at target precision P is C(P,P), whereas the total cost of n iterations of order k is
So far we have only considered strategies that use the same order k for all iterations, and we have not yet shown that such strategies are the best ones. We now give a plausible argument (not quite a rigorous proof) to justify this claim.
Consider the optimal strategy S(P^2) for the initial precision 1 and the final precision P^2, when P is very large. Since it is better to use several iterations at lower orders, we may assume that the strategy S(P^2) contains many iterations and that one of these iterations reaches precision P. Then the strategy S(P^2) is equivalent to a sequence of the two optimal strategies to go from 1 to P and from P to P^2. However, both strategies must be the same because the optimal strategy only depends on the ratio of precisions. Therefore, the optimal strategy S(P^2) is a sequence of two identical strategies (S(P), S(P)).
Suppose that k[1] is the first element of S(P). The optimal strategy to go from precision k[1] to precision P*k[1] is also S(P). Therefore the second element of S(P) is also equal to k[1], and by extension all elements of S(P) are the same.
A similar consideration gives the optimal strategy for other iterations that compute inverses of analytic functions, such as Newton's iteration for the inverse square root or for higher roots. The difference is that the value of c should be chosen as the equivalent number of multiplications needed to compute the function. For instance, c=1 for division and c=2 for the inverse square root iteration.
The conclusion is that in each case we should compute the optimal order k in advance and use this order for all iterations.
Divisions by large integers k! and separate evaluations of powers x^k are avoided if we store the previous term. The next term can be obtained by a short division of the previous term by k and a long multiplication by x. Then we only need O(N) long multiplications to evaluate the series. Usually the required number of terms N=O(P), so the total cost is O(P*M(P)).
There is no accumulation of round-off error in this method if x is small enough (in the case of Exp(x), a sufficient condition is Abs(x)<1/2). To see this, suppose that x is known to P digits (with relative error 10^(-P)). Since Abs(x)<1/2, the n-th term x^n/n! <4^(-n) (this is a rough estimate but it is enough). Since each multiplication by x results in adding 1 significant bit of relative round-off error, the relative error of x^n/n! is about 2^n times the relative error of x, i.e. 2^n*10^(-P). So the absolute round-off error of x^n/n! is not larger than
In practice, one could truncate the precision of x^n/n! as n grows, leaving a few guard bits each time to keep the round-off error negligibly small and yet to gain some computation speed. This however does not change the asymptotic complexity of the method---it remains O(P*M(P)).
However, if x is a small rational number, then the multiplication by x is short and takes O(P) operations. In that case, the total complexity of the method is O(P^2) which is always faster than O(P*M(P)).
If the coefficients a[k] are related by a simple ratio, then Horner's scheme may be modified to simplify the calculations. For example, the Horner scheme for the Taylor series for Exp(x) may be written as
Similarly to the simple summation method, the working precision for Horner's scheme may be adjusted to reduce the computation time: for example, x*a[N-1] needs to be computed with just a few significant digits if x is small. This does not change the asymptotic complexity of the method: it requires O(N)=O(P) long multiplications by x, so for general real x the complexity is again O(P*M(P)). However, if x is a small rational number, then the multiplication by x is short and takes O(P) operations. In that case, the total complexity of the method is O(P^2).
The "rectangular" algorithm uses 2*Sqrt(N) long multiplications (assuming that the coefficients of the series are short rational numbers) and Sqrt(N) units of storage. For high-precision floating-point x, this method provides a significant advantage over Horner's scheme.
Suppose we need to evaluate Sum(k,0,N,a[k]*x^k) and we know the number of terms N in advance. Suppose also that the coefficients a[k] are rational numbers with small numerators and denominators, so a multiplication a[k]*x is not a long multiplication (usually, either a[k] or the ratio a[k]/a[k-1] is a short rational number). Then we can organize the calculation in a rectangular array with c columns and r rows like this,
The total required number of long multiplications is r+c+Ln(r)-2. The minimum number of multiplications, given that r*c>=N, is around 2*Sqrt(N) at r<=>Sqrt(N)-1/2. Therefore, by arranging the Taylor series in a rectangle with sides r and c, we obtain an algorithm which costs O(Sqrt(N)) instead of O(N) long multiplications and requires Sqrt(N) units of storage.
One might wonder if we should not try to arrange the Taylor series in a cube or another multidimensional matrix instead of a rectangle. However, calculations show that this does not save time: the optimal arrangement is the two-dimensional rectangle.
The rectangular method saves the number of long multiplications by x but increases the number of short multiplications and additions. If x is a small integer or a small rational number, multiplications by x are fast and it does not make sense to use the rectangular method. Direct evaluation schemes are more efficient in that case.
Reducing the working precision saves some computation time. (We also need to estimate M but this can usually be done quickly by bit counting.) Instead of O(Sqrt(P)) long multiplications at precision P, we now need one long multiplication at precision P, another long multiplication at precision P-M, and so on. This technique will not change the asymptotic complexity which remains O(Sqrt(P)*M(P)), but it will reduce the constant factor in front of the O.
Like the previous two methods, there is no accumulated round-off error if x is small.
In the first case, it is better to use either Horner's scheme (for small P, slow multiplication) or the binary splitting technique (for large P, fast multiplication). The rectangular method is actually slower than Horner's scheme if x and the coefficients a[k] are small rational numbers. In the second case (when x is a floating-point number), it is better to use the "rectangular" algorithm.
In both cases we need to know the number of terms in advance, as we will have to repeat the whole calculation if a few more terms are needed. The simple summation method rarely gives an advantage over Horner's scheme, because it is almost always the case that one can easily compute the number of terms required for any target precision.
Note that if the argument x is not small, round-off error will become significant and needs to be considered separately for a given series.
For example, consider the Taylor series for Sin(x),
The above series expansions are asymptotic in the following sense: if we truncate the series and then take the limit of very large x, then the difference between the two sides of the equation goes to zero.
It is important that the series be first truncated and then the limit of large x be taken. Usually, an asymptotic series, if taken as an infinite series, does not actually converge for any finite x. This can be seen in the examples above. For instance, in the asymptotic series for Erfc(x) the n-th term has (2*n-1)!! in the numerator which grows faster than the n-th power of any number. The terms of the series decrease at first but then eventually start to grow, even if we select a large value of x.
The way to use an asymptotic series for a numerical calculation is to truncate the series well before the terms start to grow.
Error estimates of the asymptotic series are sometimes difficult, but the rule of thumb seems to be that the error of the approximation is usually not greater than the first discarded term of the series. This can be understood intuitively as follows. Suppose we truncate the asymptotic series at a point where the terms still decrease, safely before they start to grow. For example, let the terms around the 100-th term be A[100], A[101], A[102], ..., each of these numbers being significantly smaller than the previous one, and suppose we retain A[100] but drop the terms after it. Then our approximation would have been a lot better if we retained A[101] as well. (This step of the argument is really an assumption about the behavior of the series; it seems that this assumption is correct in many practically important cases.) Therefore the error of the approximation is approximately equal to A[101].
The inherent limitation of the method of asymptotic series is that for any given x, there will be a certain place in the series where the term has the minimum absolute value (after that, the series is unusable), and the error of the approximation cannot be smaller than that term.
For example, take the above asymptotic series for Erfc(x). The logarithm of the absolute value of the n-th term can be estimated using Stirling's formula for the factorial as
We find that for a given finite x, no matter how large, there is a maximum precision that can be achieved with the asymptotic series; if we need more precision, we have to use a different method.
However, sometimes the function we are evaluating allows identity transformations that relate f(x) to f(y) with y>x. For example, the Gamma function satisfies x*Gamma(x)=Gamma(x+1). In this case we can transform the function so that we would need to evaluate it at large enough x for the asymptotic series to give us enough precision.
More formally, one can define the function of two arguments AGM(x,y) as the limit of the sequence a[k] where a[k+1]=1/2*(a[k]+b[k]), b[k+1]=Sqrt(a[k]*b[k]), and the initial values are a[0]=x, b[0]=y. (The limit of the sequence b[k] is the same.) This function is obviously linear, AGM(c*x,c*y)=c*AGM(x,y), so in principle it is enough to compute AGM(1,x) or arbitrarily select c for convenience.
Gauss and Legendre knew that the limit of the AGM sequence is related to the complete elliptic integral,
The AGM sequence is also defined for complex values a, b. One needs to take a square root Sqrt(a*b), which requires a branch cut to be well-defined. Selecting the natural cut along the negative real semiaxis (Re(x)<0, Im(x)=0), we obtain an AGM sequence that converges for any initial values x, y with positive real part.
Let us estimate the convergence rate of the AGM sequence starting from x, y, following the paper [Brent 1975]. Clearly the worst case is when the numbers x and y are very different (one is much larger than another). In this case the numbers a[k], b[k] become approximately equal after about k=1/Ln(2)*Ln(Abs(Ln(x/y))) iterations (note: Brent's paper online mistypes this as 1/Ln(2)*Abs(Ln(x/y))). This is easy to see: if x is much larger than y, then at each step the ratio r:=x/y is transformed into r'=1/2*Sqrt(r). When the two numbers become roughly equal to each other, one needs about Ln(n)/Ln(2) more iterations to make the first n (decimal) digits of a[k] and b[k] coincide, because the relative error epsilon=1-b/a decays approximately as epsilon[k]<=>1/8*Exp(-2^k).
Unlike Newton's iteration, the AGM sequence does not correct errors, so all numbers need to be computed with full precision. Actually, slightly more precision is needed to compensate for accumulated round-off error. Brent (in [Brent 1975]) says that O(Ln(Ln(n))) bits of accuracy are lost to round-off error if there are total of n iterations.
The AGM sequence can be used for fast computations of Pi, Ln(x) and ArcTan(x). However, currently the limitations of Yacas internal math make these methods less efficient than simpler methods based on Taylor series and Newton iterations.
If we need to take O(P) terms of the series to obtain P digits of precision, then ordinary methods would require O(P^2) arithmetic operations. (Each term needs O(P) operations because all coefficients are rational numbers with O(P) digits and we need to perform a few short multiplications or divisions.) The binary splitting method requires O(M(P*Ln(P))*Ln(P)) operations instead of the O(P^2) operations. In other words, we need to perform long multiplications of integers of size O(P*Ln(P)) digits, but we need only O(Ln(P)) such multiplications. The binary splitting method performs better than the straightforward summation method if the cost of multiplication is lower than O(P^2)/Ln(P). This is usually true only for large enough precision (at least a thousand digits).
Thus there are two main limitations of the binary splitting method:
The main advantages of the method are:
For example, the Taylor series for ArcSin(x) (when x is a short rational number) is of this form:
The goal is to compute the sum S(0,N) with a chosen number of terms N. Instead of computing the rational number S directly, the binary splitting method propose to compute the following four integers P, Q, B, and T:
Thus the range [0, N) is split in half on each step. At the base of recursion the four integers P, Q, B, and T are computed directly. At the end of the calculation (top level of recursion), one floating-point division is performed to recover S=T/(B*Q). It is clear that the four integers carry the full information needed to continue the calculation with more terms. So this algorithm is easy to checkpoint and parallelize.
The integers P, Q, B, and T grow during the calculation to O(N*Ln(N)) bits, and we need to multiply these large integers. However, there are only O(Ln(N)) steps of recursion and therefore O(Ln(N)) long multiplications are needed. If the series converges linearly, we need N=O(P) terms to obtain P digits of precision. Therefore, the total asymptotic cost of the method is O(M(P*Ln(P))*Ln(P)) operations.
A more general form of the binary splitting technique is also given in [Haible et al. 1998]. The generalization applies to series for the form
The binary splitting technique can also be used for series with complex integer coefficients, or more generally for coefficients in any finite algebraic extension of integers, e.q. Z[ Sqrt(2)] (the ring of numbers of the form p+q*Sqrt(2) where p, q are integers). Thus we may compute the Bessel function J0(Sqrt(3)) using the binary splitting method and obtain exact intermediate results of the form p+q*Sqrt(3). But this will still not help compute J0(Pi). This is a genuine limitation of the binary splitting method.