numerical integration of expr with respect to X on the circle |X-a |= R, divided by 2iPi. In other words, when expr is a meromorphic function, sum of the residues in the corresponding disk. tab is as in intnum, except that if computed with intnuminit it should be with the endpoints [-1, 1].

? \p105
 ? intcirc(s=1, 0.5, zeta(s)) - 1
 time = 3,460 ms.
 %1 = -2.40... E-104 - 2.7... E-106*I

The library syntax is intcirc(void *E, GEN (*eval)(GEN,void*), GEN a,GEN R,GEN tab, long prec).

numerical integration of expr(X)cos(2Pi zX) from a to b, in other words Fourier cosine transform (from a to b) of the function represented by expr. a and b are coded as in intnum, and are not necessarily at infinity, but if they are, oscillations (i.e. [[±1],alpha I]) are forbidden.

The library syntax is intfouriercos(void *E, GEN (*eval)(GEN,void*), GEN a, GEN b, GEN z, GEN tab, long prec).

numerical integration of expr(X)exp(-2Pi zX) from a to b, in other words Fourier transform (from a to b) of the function represented by expr. Note the minus sign. a and b are coded as in intnum, and are not necessarily at infinity but if they are, oscillations (i.e. [[±1],alpha I]) are forbidden.

The library syntax is intfourierexp(void *E, GEN (*eval)(GEN,void*), GEN a, GEN b, GEN z, GEN tab, long prec).

numerical integration of expr(X)sin(2Pi zX) from a to b, in other words Fourier sine transform (from a to b) of the function represented by expr. a and b are coded as in intnum, and are not necessarily at infinity but if they are, oscillations (i.e. [[±1],alpha I]) are forbidden.

The library syntax is intfouriersin(void *E, GEN (*eval)(GEN,void*), GEN a, GEN b, GEN z, GEN tab, long prec).

initalize tables for use with integral transforms such as intmellininv, etc., where a and b are coded as in intnum, expr is the function s(X) to which the integral transform is to be applied (which will multiply the weights of integration) and m is as in intnuminit. If flag is nonzero, assumes that s(-X) = \overline{s(X)}, which makes the computation twice as fast. See intmellininvshort for examples of the use of this function, which is particularly useful when the function s(X) is lengthy to compute, such as a gamma product.

The library syntax is intfuncinit(void *E, GEN (*eval)(GEN,void*), GEN a,GEN b,long m, long flag, long prec). Note that the order of m and flag are reversed compared to the GP syntax.

numerical integration of expr(X)e^{Xz} with respect to X on the line Re(X) = sig, divided by 2iPi, in other words, inverse Laplace transform of the function corresponding to expr at the value z.

sig is coded as follows. Either it is a real number sigma, equal to the abcissa of integration, and then the function to be integrated is assumed to be slowly decreasing when the imaginary part of the variable tends to ± oo . Or it is a two component vector [sigma,alpha], where sigma is as before, and either alpha = 0 for slowly decreasing functions, or alpha > 0 for functions decreasing like exp(-alpha t). Note that it is not necessary to choose the exact value of alpha. tab is as in intnum.

It is often a good idea to use this function with a value of m one or two higher than the one chosen by default (which can be viewed thanks to the function intnumstep), or to increase the abcissa of integration sigma. For example:

? \p 105
 ? intlaplaceinv(x=2, 1, 1/x) - 1
 time = 350 ms.
 %1 = 7.37... E-55 + 1.72... E-54*I \\ not so good

 ? m = intnumstep()
 %2 = 7
 ? intlaplaceinv(x=2, 1, 1/x, m+1) - 1
 time = 700 ms.
 %3 = 3.95... E-97 + 4.76... E-98*I \\ better

 ? intlaplaceinv(x=2, 1, 1/x, m+2) - 1
 time = 1400 ms.
 %4 = 0.E-105 + 0.E-106*I \\ perfect but slow.

 ? intlaplaceinv(x=5, 1, 1/x) - 1
 time = 340 ms.
 %5 = -5.98... E-85 + 8.08... E-85*I \\ better than %1

 ? intlaplaceinv(x=5, 1, 1/x, m+1) - 1
 time = 680 ms.
 %6 = -1.09... E-106 + 0.E-104*I \\ perfect, fast.

 ? intlaplaceinv(x=10, 1, 1/x) - 1
 time = 340 ms.
 %7 = -4.36... E-106 + 0.E-102*I \\ perfect, fastest, but why sig = 10?

 ? intlaplaceinv(x=100, 1, 1/x) - 1
 time = 330 ms.
 %7 = 1.07... E-72 + 3.2... E-72*I \\ too far now...

The library syntax is intlaplaceinv(void *E, GEN (*eval)(GEN,void*), GEN sig,GEN z, GEN tab, long prec).

numerical integration of expr(X)z^{-X} with respect to X on the line Re(X) = sig, divided by 2iPi, in other words, inverse Mellin transform of the function corresponding to expr at the value z.

sig is coded as follows. Either it is a real number sigma, equal to the abcissa of integration, and then the function to be integrated is assumed to decrease exponentially fast, of the order of exp(-t) when the imaginary part of the variable tends to ± oo . Or it is a two component vector [sigma,alpha], where sigma is as before, and either alpha = 0 for slowly decreasing functions, or alpha > 0 for functions decreasing like exp(-alpha t), such as gamma products. Note that it is not necessary to choose the exact value of alpha, and that alpha = 1 (equivalent to sig alone) is usually sufficient. tab is as in intnum.

As all similar functions, this function is provided for the convenience of the user, who could use intnum directly. However it is in general better to use intmellininvshort.

? \p 105
 ? intmellininv(s=2,4, gamma(s)^3);
 time = 1,190 ms. \\ reasonable.

 ? \p 308
 ? intmellininv(s=2,4, gamma(s)^3);
 time = 51,300 ms. \\ slow because of Gamma(s)^3.

The library syntax is intmellininv(void *E, GEN (*eval)(GEN,void*), GEN sig, GEN z, GEN tab, long prec).

numerical integration of s(X)z^{-X} with respect to X on the line Re(X) = sig, divided by 2iPi, in other words, inverse Mellin transform of s(X) at the value z. Here s(X) is implicitly contained in tab in intfuncinit format, typically

  tab = intfuncinit(T = [-1], [1], s(sig + I*T))

or similar commands. Take the example of the inverse Mellin transform of Gamma(s)^3 given in intmellininv:

? \p 105
 ? oo = [1]; \\ for clarity

 ? A = intmellininv(s=2,4, gamma(s)^3);
 time = 2,500 ms. \\ not too fast because of Gamma(s)^3.

 \\  function of real type, decreasing as exp(-3Pi/2.|t|)

 ? tab = intfuncinit(t=[-oo, 3*Pi/2],[oo, 3*Pi/2], gamma(2+I*t)^3, 1);
 time = 1,370 ms.
 ? intmellininvshort(2,4, tab) - A
 time = 50 ms.
 %4 = -1.26... - 3.25...E-109*I \\ 50 times faster than 
A and perfect.

 ? tab2 = intfuncinit(t=-oo, oo, gamma(2+I*t)^3, 1);
 ? intmellininvshort(2,4, tab2)
 %6 = -1.2...E-42 - 3.2...E-109*I  \\ 63 digits lost

In the computation of tab, it was not essential to include the exact exponential decrease of Gamma(2+it)^3. But as the last example shows, a rough indication must be given, otherwise slow decrease is assumed, resulting in catastrophic loss of accuracy.

The library syntax is intmellininvshort(GEN sig, GEN z, GEN tab, long prec).

numerical integration of expr on [a,b] (possibly infinite interval) with respect to X, where a and b are coded as explained below. The integrand may have values belonging to a vector space over the real numbers; in particular, it can be complex-valued or vector-valued.

If tab is omitted, necessary integration tables are computed using intnuminit according to the current precision. It may be a positive integer m, and tables are computed assuming the integration step is 1/2^m. Finally tab can be a table output by intnuminit, in which case it is used directly. This is important if several integrations of the same type are performed (on the same kind of interval and functions, and the same accuracy), since it saves expensive precomputations.

If tab is omitted the algorithm guesses a reasonable value for m depending on the current precision. That value may be obtained as

  intnumstep()

However this value may be off from the optimal one, and this is important since the integration time is roughly proportional to 2^m. One may try consecutive values of m until they give the same value up to an accepted error.

The endpoints a and b are coded as follows. If a is not at ± oo , it is either coded as a scalar (real or complex), or as a two component vector [a,alpha], where the function is assumed to have a singularity of the form (x-a)^{alpha+epsilon} at a, where epsilon indicates that powers of logarithms are neglected. In particular, [a,alpha] with alpha >= 0 is equivalent to a. If a wrong singularity exponent is used, the result will lose a catastrophic number of decimals, for instance approximately half the number of digits will be correct if alpha = -1/2 is omitted.

The endpoints of integration can be ± oo , which is coded as [± 1] or as [[±1],alpha]. Here alpha codes the behaviour of the function at ± oo as follows.

* alpha = 0 (or no alpha at all, i.e. simply [±1]) assumes that the function to be integrated tends to zero, but not exponentially fast, and not oscillating such as sin(x)/x.

* alpha > 0 assumes that the function tends to zero exponentially fast approximately as exp(-alpha x), including reasonably oscillating functions such as exp(-x)sin(x). The precise choice of alpha, while useful in extreme cases, is not critical, and may be off by a factor of 10 or more from the correct value.

* alpha < -1 assumes that the function tends to 0 slowly, like x^{alpha}. Here it is essential to give the correct alpha, if possible, but on the other hand alpha <= -2 is equivalent to alpha = 0, in other words to no alpha at all.

The last two codes are reserved for oscillating functions. Let k > 0 real, and g(x) a nonoscillating function tending to 0, then

* alpha = k I assumes that the function behaves like cos(kx)g(x).

* alpha = -kI assumes that the function behaves like sin(kx)g(x).

Here it is critical to give the exact value of k. If the oscillating part is not a pure sine or cosine, one must expand it into a Fourier series, use the above codings, and sum the resulting contributions. Otherwise you will get nonsense. Note that cos(kx) (and similarly sin(kx)) means that very function, and not a translated version such as cos(kx+a).

If for instance f(x) = cos(kx)g(x) where g(x) tends to zero exponentially fast as exp(-alpha x), it is up to the user to choose between [[±1],alpha] and [[±1],kI], but a good rule of thumb is that if the oscillations are much weaker than the exponential decrease, choose [[±1],alpha], otherwise choose [[±1],kI], although the latter can reasonably be used in all cases, while the former cannot. To take a specific example, in the inverse Mellin transform, the function to be integrated is almost always exponentially decreasing times oscillating. If we choose the oscillating type of integral we perhaps obtain the best results, at the expense of having to recompute our functions for a different value of the variable z giving the transform, preventing us to use a function such as intmellininvshort. On the other hand using the exponential type of integral, we obtain less accurate results, but we skip expensive recomputations. See intmellininvshort and intfuncinit for more explanations.

Note. If you do not like the code [±1] for ± oo , you are welcome to set, e.g oo = [1] or INFINITY = [1], then using +oo, -oo, -INFINITY, etc. will have the expected behaviour.

We shall now see many examples to get a feeling for what the various parameters achieve. All examples below assume precision is set to 105 decimal digits. We first type

? \p 105
 ? oo = [1]  \\ for clarity

Apparent singularities. Even if the function f(x) represented by expr has no singularities, it may be important to define the function differently near special points. For instance, if f(x) = 1 /(exp(x)-1) - exp(-x)/x, then int_0^ oo f(x)dx = gamma, Euler's constant Euler. But

? f(x) = 1/(exp(x)-1) - exp(-x)/x
 ? intnum(x = 0, [oo,1],  f(x)) - Euler
 %1 = 6.00... E-67

thus only correct to 76 decimal digits. This is because close to 0 the function f is computed with an enormous loss of accuracy. A better solution is

? f(x) = 1/(exp(x)-1)-exp(-x)/x
 ? F = truncate( f(t + O(t^7)) ); \\ expansion around t = 0

 ? g(x) = if (x > 1e-18, f(x), subst(F,t,x))  \\ note that 6.18 > 105

 ? intnum(x = 0, [oo,1],  g(x)) - Euler
 %2 = 0.E-106 \\ perfect

It is up to the user to determine constants such as the 10^{-18} and 7 used above.

True singularities. With true singularities the result is much worse. For instance

? intnum(x = 0, 1,  1/sqrt(x)) - 2
 %1 = -1.92... E-59 \\ only 59 correct decimals

 
 ? intnum(x = [0,-1/2], 1,  1/sqrt(x)) - 2
 %2 = 0.E-105 \\ better

Oscillating functions.

? intnum(x = 0, oo, sin(x) / x) - Pi/2
 %1 = 20.78.. \\ nonsense

 ? intnum(x = 0, [oo,1], sin(x)/x) - Pi/2
 %2 = 0.004.. \\ bad

 ? intnum(x = 0, [oo,-I], sin(x)/x) - Pi/2
 %3 = 0.E-105 \\ perfect

 ? intnum(x = 0, [oo,-I], sin(2*x)/x) - Pi/2  \\ oops, wrong k

 %4 = 0.07...
 ? intnum(x = 0, [oo,-2*I], sin(2*x)/x) - Pi/2
 %5 = 0.E-105 \\ perfect

 
 ? intnum(x = 0, [oo,-I], sin(x)^3/x) - Pi/4
 %6 = 0.0092... \\ bad

 ? sin(x)^3 - (3*sin(x)-sin(3*x))/4
 %7 = O(x^17)

We may use the above linearization and compute two oscillating integrals with "infinite endpoints" [oo, -I] and [oo, -3*I] respectively, or notice the obvious change of variable, and reduce to the single integral (1/2)int_0^ oo sin(x)/xdx. We finish with some more complicated examples:

? intnum(x = 0, [oo,-I], (1-cos(x))/x^2) - Pi/2
 %1 = -0.0004... \\ bad

 ? intnum(x = 0, 1, (1-cos(x))/x^2) \
 + intnum(x = 1, oo, 1/x^2) - intnum(x = 1, [oo,I], cos(x)/x^2) - Pi/2
 %2 = -2.18... E-106 \\ OK

 
 ? intnum(x = 0, [oo, 1], sin(x)^3*exp(-x)) - 0.3
 %3 = 5.45... E-107 \\ OK

 ? intnum(x = 0, [oo,-I], sin(x)^3*exp(-x)) - 0.3
 %4 = -1.33... E-89 \\ lost 16 decimals. Try higher m:

 ? m = intnumstep()
 %5 = 7 \\ the value of m actually used above.

 ? tab = intnuminit(0,[oo,-I], m+1); \\ try m one higher.

 ? intnum(x = 0, oo, sin(x)^3*exp(-x), tab) - 0.3
 %6 = 5.45... E-107 \\ OK this time.

Warning. Like sumalt, intnum often assigns a reasonable value to diverging integrals. Use these values at your own risk! For example:

? intnum(x = 0, [oo, -I], x^2*sin(x))
 %1 = -2.0000000000...

Note the formula int_0^ oo sin(x)/x^sdx = cos(Pi s/2) Gamma(1-s) , a priori valid only for 0 < Re(s) < 2, but the right hand side provides an analytic continuation which may be evaluated at s = -2...

Multivariate integration. Using successive univariate integration with respect to different formal parameters, it is immediate to do naive multivariate integration. But it is important to use a suitable intnuminit to precompute data for the internal integrations at least!

For example, to compute the double integral on the unit disc x^2+y^2 <= 1 of the function x^2+y^2, we can write

? tab = intnuminit(-1,1);
 ? intnum(x=-1,1, intnum(y=-sqrt(1-x^2),sqrt(1-x^2), x^2+y^2, tab), tab)

The first tab is essential, the second optional. Compare:

? tab = intnuminit(-1,1);
 time = 30 ms.
 ? intnum(x=-1,1, intnum(y=-sqrt(1-x^2),sqrt(1-x^2), x^2+y^2));
 time = 54,410 ms. \\ slow

 ? intnum(x=-1,1, intnum(y=-sqrt(1-x^2),sqrt(1-x^2), x^2+y^2, tab), tab);
 time = 7,210 ms.  \\ faster

However, the intnuminit program is usually pessimistic when it comes to choosing the integration step 2^{-m}. It is often possible to improve the speed by trial and error. Continuing the above example:

? test(M) =
 {
   tab = intnuminit(-1,1, M);
   intnum(x=-1,1, intnum(y=-sqrt(1-x^2),sqrt(1-x^2), x^2+y^2,tab), tab) - Pi/2
 }
 ? m = intnumstep() \\ what value of m did it take ?

 %1 = 7
 ? test(m - 1)
 time = 1,790 ms.
 %2 = -2.05... E-104 \\ 4 = 2^2 times faster and still OK.

 ? test(m - 2)
 time = 430 ms.
 %3 = -1.11... E-104 \\ 16 = 2^4 times faster and still OK.

 ? test(m - 3)
 time = 120 ms.
 %3 = -7.23... E-60 \\ 64 = 2^6 times faster, lost 45 decimals.

The library syntax is intnum(void *E, GEN (*eval)(GEN,void*), GEN a,GEN b,GEN tab, long prec), where an omitted tab is coded as NULL.

initialize tables for integration from a to b, where a and b are coded as in intnum. Only the compactness, the possible existence of singularities, the speed of decrease or the oscillations at infinity are taken into account, and not the values. For instance intnuminit(-1,1) is equivalent to intnuminit(0,Pi), and intnuminit([0,-1/2],[1]) is equivalent to {\tt intnuminit([-1],[-1,-1/2])}. If m is not given, it is computed according to the current precision. Otherwise the integration step is 1/2^m. Reasonable values of m are m = 6 or m = 7 for 100 decimal digits, and m = 9 for 1000 decimal digits.

The result is technical, but in some cases it is useful to know the output. Let x = phi(t) be the change of variable which is used. tab[1] contains the integer m as above, either given by the user or computed from the default precision, and can be recomputed directly using the function intnumstep. tab[2] and tab[3] contain respectively the abcissa and weight corresponding to t = 0 (phi(0) and phi'(0)). tab[4] and tab[5] contain the abcissas and weights corresponding to positive t = nh for 1 <= n <= N and h = 1/2^m (phi(nh) and phi'(nh)). Finally tab[6] and tab[7] contain either the abcissas and weights corresponding to negative t = nh for -N <= n <= -1, or may be empty (but not always) if phi(t) is an odd function (implicitly we would have tab[6] = -tab[4] and tab[7] = tab[5]).

The library syntax is intnuminit(GEN a, GEN b, long m, long prec).

numerical integration of expr (smooth in ]a,b[), with respect to X. This function is deprecated, use intnum instead.

Set flag = 0 (or omit it altogether) when a and b are not too large, the function is smooth, and can be evaluated exactly everywhere on the interval [a,b].

If flag = 1, uses a general driver routine for doing numerical integration, making no particular assumption (slow).

flag = 2 is tailored for being used when a or b are infinite. One must have ab > 0, and in fact if for example b = + oo , then it is preferable to have a as large as possible, at least a >= 1.

If flag = 3, the function is allowed to be undefined (but continuous) at a or b, for example the function sin(x)/x at x = 0.

The user should not require too much accuracy: 18 or 28 decimal digits is OK, but not much more. In addition, analytical cleanup of the integral must have been done: there must be no singularities in the interval or at the boundaries. In practice this can be accomplished with a simple change of variable. Furthermore, for improper integrals, where one or both of the limits of integration are plus or minus infinity, the function must decrease sufficiently rapidly at infinity. This can often be accomplished through integration by parts. Finally, the function to be integrated should not be very small (compared to the current precision) on the entire interval. This can of course be accomplished by just multiplying by an appropriate constant.

Note that infinity can be represented with essentially no loss of accuracy by 1e1000. However beware of real underflow when dealing with rapidly decreasing functions. For example, if one wants to compute the int_0^ oo e^{-x^2}dx to 28 decimal digits, then one should set infinity equal to 10 for example, and certainly not to 1e1000.

The library syntax is intnumromb(void *E, GEN (*eval)(GEN,void*), GEN a, GEN b, long flag, long prec), where eval(x, E) returns the value of the function at x. You may store any additional information required by eval in E, or set it to NULL.

give the value of m used in all the intnum and sumnum programs, hence such that the integration step is equal to 1/2^m.

The library syntax is intnumstep(long prec).

product of expression expr, initialized at x, the formal parameter X going from a to b. As for sum, the main purpose of the initialization parameter x is to force the type of the operations being performed. For example if it is set equal to the integer 1, operations will start being done exactly. If it is set equal to the real 1., they will be done using real numbers having the default precision. If it is set equal to the power series 1+O(X^k) for a certain k, they will be done using power series of precision at most k. These are the three most common initializations.

As an extreme example, compare

? prod(i=1, 100, 1 - X^i);  \\ this has degree 5050 !!

 time = 3,335 ms.
 ? prod(i=1, 100, 1 - X^i, 1 + O(X^101))
 time = 43 ms.
 %2 = 1 - X - X^2 + X^5 + X^7 - X^12 - X^15 + X^22 + X^26 - X^35 - X^40 + \
   X^51 + X^57 - X^70 - X^77 + X^92 + X^100 + O(X^101)

The library syntax is produit(entree *ep, GEN a, GEN b, char *expr, GEN x).

product of expression expr, initialized at 1. (i.e.to a real number equal to 1 to the current realprecision), the formal parameter X ranging over the prime numbers between a and b.

The library syntax is prodeuler(void *E, GEN (*eval)(GEN,void*), GEN a,GEN b, long prec).

infinite product of expression expr, the formal parameter X starting at a. The evaluation stops when the relative error of the expression minus 1 is less than the default precision. The expressions must always evaluate to an element of C.

If flag = 1, do the product of the (1+expr) instead.

The library syntax is prodinf(void *E, GEN (*eval)(GEN, void*), GEN a, long prec) (flag = 0), or prodinf1 with the same arguments (flag = 1).

find a real root of expression expr between a and b, under the condition expr(X = a) * expr(X = b) <= 0. This routine uses Brent's method and can fail miserably if expr is not defined in the whole of [a,b] (try solve(x = 1, 2, tan(x)).

The library syntax is zbrent(void *E,GEN (*eval)(GEN,void*),GEN a,GEN b,long prec).

sum of expression expr, initialized at x, the formal parameter going from a to b. As for prod, the initialization parameter x may be given to force the type of the operations being performed.

As an extreme example, compare

? sum(i=1, 5000, 1/i); \\ rational number: denominator has 2166 digits.

 time = 1,241 ms.
 ? sum(i=1, 5000, 1/i, 0.)
 time = 158 ms.
 %2 = 9.094508852984436967261245533

The library syntax is somme(entree *ep, GEN a, GEN b, char *expr, GEN x). This is to be used as follows: ep represents the dummy variable used in the expression expr

/* compute a^2 + ... + b^2 */
 {
   /* define the dummy variable "i" */
   entree *ep = is_entry("i");
   /* sum for a <= i <= b */
   return somme(ep, a, b, "i^2", gen_0);
 }

numerical summation of the series expr, which should be an alternating series, the formal variable X starting at a. Use an algorithm of F.Villegas as modified by D.Zagier (improves on Euler-Van Wijngaarden method).

If flag = 1, use a variant with slightly different polynomials. Sometimes faster.

Divergent alternating series can sometimes be summed by this method, as well as series which are not exactly alternating (see for example Section [Label: se:user_defined]). If the series already converges geometrically, suminf is often a better choice:

? \p28
 ? sumalt(i = 1, -(-1)^i / i)  - log(2)
 time = 0 ms.
 %1 = -2.524354897 E-29
 ? suminf(i = 1, -(-1)^i / i)
   *** suminf: user interrupt after 10min, 20,100 ms.
 ? \p1000
 ? sumalt(i = 1, -(-1)^i / i)  - log(2)
 time = 90 ms.
 %2 = 4.459597722 E-1002
 
 ? sumalt(i = 0, (-1)^i / i!) - exp(-1)
 time = 670 ms.
 %3 = -4.03698781490633483156497361352190615794353338591897830587 E-944
 ? suminf(i = 0, (-1)^i / i!) - exp(-1)
 time = 110 ms.
 %4 = -8.39147638 E-1000   \\  faster and more accurate

The library syntax is sumalt(void *E, GEN (*eval)(GEN,void*),GEN a,long prec). Also available is sumalt2 with the same arguments (flag = 1).

sum of expression expr over the positive divisors of n.

Arithmetic functions like sigma use the multiplicativity of the underlying expression to speed up the computation. In the present version 2.3.1, there is no way to indicate that expr is multiplicative in n, hence specialized functions should be preferred whenever possible.

The library syntax is divsum(entree *ep, GEN num, char *expr).

infinite sum of expression expr, the formal parameter X starting at a. The evaluation stops when the relative error of the expression is less than the default precision for 3 consecutive evaluations. The expressions must always evaluate to a complex number.

If the series converges slowly, make sure realprecision is low (even 28 digits may be too much). In this case, if the series is alternating or the terms have a constant sign, sumalt and sumpos should be used instead.

? \p28
 ? suminf(i = 1, -(-1)^i / i)
   *** suminf: user interrupt after 10min, 20,100 ms.
 ? sumalt(i = 1, -(-1)^i / i) - log(2)
 time = 0 ms.
 %1 = -2.524354897 E-29

The library syntax is suminf(void *E, GEN (*eval)(GEN,void*), GEN a, long prec).

numerical summation of expr, the variable X taking integer values from ceiling of a to + oo , where expr is assumed to be a holomorphic function f(X) for Re(X) >= sigma.

The parameter sigma belongs to R is coded in the argument sig as follows: it is either

* a real number sigma. Then the function f is assumed to decrease at least as 1/X^2 at infinity, but not exponentially;

* a two-component vector [sigma,alpha], where sigma is as before, alpha < -1. The function f is assumed to decrease like X^{alpha}. In particular, alpha <= -2 is equivalent to no alpha at all.

* a two-component vector [sigma,alpha], where sigma is as before, alpha > 0. The function f is assumed to decrease like exp(-alpha X). In this case it is essential that alpha be exactly the rate of exponential decrease, and it is usually a good idea to increase the default value of m used for the integration step. In practice, if the function is exponentially decreasing sumnum is slower and less accurate than sumpos or suminf, so should not be used.

The function uses the intnum routines and integration on the line Re(s) = sigma. The optional argument tab is as in intnum, except it must be initialized with sumnuminit instead of intnuminit.

When tab is not precomputed, sumnum can be slower than sumpos, when the latter is applicable. It is in general faster for slowly decreasing functions.

Finally, if flag is nonzero, we assume that the function f to be summed is of real type, i.e. satisfies \overline{f(z)} = f(\overline{z}), which speeds up the computation.

? \p 308
 ? a = sumpos(n=1, 1/(n^3+n+1));
 time = 1,410 ms.
 ? tab = sumnuminit(2);
 time = 1,620 ms. \\ slower but done once and for all.

 ? b = sumnum(n=1, 2, 1/(n^3+n+1), tab);
 time = 460 ms. \\ 3 times as fast as 
sumpos

 ? a - b
 %4 = -1.0... E-306 + 0.E-320*I \\ perfect.

 ? sumnum(n=1, 2, 1/(n^3+n+1), tab, 1) - a; \\ function of real type

 time = 240 ms.
 %2 = -1.0... E-306 \\ twice as fast, no imaginary part.

 ? c = sumnum(n=1, 2, 1/(n^2+1), tab, 1);
 time = 170 ms. \\ fast

 ? d = sumpos(n=1, 1 / (n^2+1));
 time = 2,700 ms. \\ slow.

 ? d - c
 time = 0 ms.
 %5 = 1.97... E-306 \\ perfect.

For slowly decreasing function, we must indicate singularities:

? \p 308
 ? a = sumnum(n=1, 2, n^(-4/3));
 time = 9,930 ms. \\ slow because of the computation of n^{-4/3}.

 ? a - zeta(4/3)
 time = 110 ms.
 %1 = -2.42... E-107 \\ lost 200 decimals because of singularity at  oo 

 ? b = sumnum(n=1, [2,-4/3], n^(-4/3), /*omitted*/, 1); \\ of real type

 time = 12,210 ms.
 ? b - zeta(4/3)
 %3 = 1.05... E-300 \\ better

Since the complex values of the function are used, beware of determination problems. For instance:

? \p 308
 ? tab = sumnuminit([2,-3/2]);
 time = 1,870 ms.
 ? sumnum(n=1,[2,-3/2], 1/(n*sqrt(n)), tab,1) - zeta(3/2)
 time = 690 ms.
 %1 = -1.19... E-305 \\ fast and correct

 ? sumnum(n=1,[2,-3/2], 1/sqrt(n^3), tab,1) - zeta(3/2)
 time = 730 ms.
 %2 = -1.55... \\ nonsense. However

 ? sumnum(n=1,[2,-3/2], 1/n^(3/2), tab,1) - zeta(3/2)
 time = 8,990 ms.
 %3 = -1.19... E-305 \\ perfect, as 1/(n*sqrt{n}) above but much slower

For exponentially decreasing functions, sumnum is given for completeness, but one of suminf or sumpos should always be preferred. If you experiment with such functions and sumnum anyway, indicate the exact rate of decrease and increase m by 1 or 2:

? suminf(n=1, 2^(-n)) - 1
 time = 10 ms.
 %1 = -1.11... E-308 \\ fast and perfect

 ? sumpos(n=1, 2^(-n)) - 1
 time = 10 ms.
 %2 = -2.78... E-308 \\ also fast and perfect

 ? sumnum(n=1,2, 2^(-n)) - 1
    *** sumnum: precision too low in mpsc1 \\ nonsense

 ? sumnum(n=1, [2,log(2)], 2^(-n), /*omitted*/, 1) - 1 \\ of real type

 time = 5,860 ms.
 %3 = -1.5... E-236 \\ slow and lost 70 decimals

 ? m = intnumstep()
 %4 = 9
 ? sumnum(n=1,[2,log(2)], 2^(-n), m+1, 1) - 1
 time = 11,770 ms.
 %5 = -1.9... E-305 \\ now perfect, but slow.

The library syntax is sumnum(void *E, GEN (*eval)(GEN,void*), GEN a,GEN sig,GEN tab,long flag, long prec).

numerical summation of (-1)^Xexpr(X), the variable X taking integer values from ceiling of a to + oo , where expr is assumed to be a holomorphic function for Re(X) >= sig (or sig[1]).

Warning. This function uses the intnum routines and is orders of magnitude slower than sumalt. It is only given for completeness and should not be used in practice.

Warning2. The expression expr must not include the (-1)^X coefficient. Thus sumalt(n = a,(-1)^nf(n)) is (approximately) equal to sumnumalt(n = a,sig,f(n)).

sig is coded as in sumnum. However for slowly decreasing functions (where sig is coded as [sigma,alpha] with alpha < -1), it is not really important to indicate alpha. In fact, as for sumalt, the program will often give meaningful results (usually analytic continuations) even for divergent series. On the other hand the exponential decrease must be indicated.

tab is as in intnum, but if used must be initialized with sumnuminit. If flag is nonzero, assumes that the function f to be summed is of real type, i.e. satisfies \overline{f(z)} = f(\overline{z}), and then twice faster when tab is precomputed.

? \p 308
 ? tab = sumnuminit(2, /*omitted*/, -1); \\ abcissa sigma = 2, alternating sums.

 time = 1,620 ms. \\ slow, but done once and for all.

 ? a = sumnumalt(n=1, 2, 1/(n^3+n+1), tab, 1);
 time = 230 ms. \\ similar speed to 
sumnum

 ? b = sumalt(n=1, (-1)^n/(n^3+n+1));
 time = 0 ms. \\ infinitely faster!

 ? a - b
 time = 0 ms.
 %1 = -1.66... E-308 \\ perfect

The library syntax is sumnumalt(void *E, GEN (*eval)(GEN,void*), GEN a, GEN sig, GEN tab, long flag, long prec).

initialize tables for numerical summation using sumnum (with sgn = 1) or sumnumalt (with sgn = -1), sig is the abcissa of integration coded as in sumnum, and m is as in intnuminit.

The library syntax is sumnuminit(GEN sig, long m, long sgn, long prec).

numerical summation of the series expr, which must be a series of terms having the same sign, the formal variable X starting at a. The algorithm used is Van Wijngaarden's trick for converting such a series into an alternating one, and is quite slow. For regular functions, the function sumnum is in general much faster once the initializations have been made using sumnuminit.

If flag = 1, use slightly different polynomials. Sometimes faster.

The library syntax is sumpos(void *E, GEN (*eval)(GEN,void*),GEN a,long prec). Also available is sumpos2 with the same arguments (flag = 1).