Numerical derivation of *expr* with respect to X at X = a.

? derivnum(x=0,sin(exp(x))) - cos(1) %1 = -1.262177448 E-29

A clumsier approach, which would not work in library mode, is

? f(x) = sin(exp(x)) ? f'(0) - cos(1) %1 = -1.262177448 E-29

When a is a power series, compute `derivnum(t = a,f)`

as f'(a) =
(f(a))'/a'.

The library syntax is

. Also
available is **derivnum**(void *E, GEN (*eval)(void*,GEN), GEN a, long prec)`GEN `

, which also allows power series for a.**derivfun**(void *E, GEN (*eval)(void *, GEN), GEN a, long prec)

Numerical
integration of (2iPi)^{-1}*expr* with respect to X on the circle
|X-a |= R.
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 %1 = -2.398082982 E-104 - 7.94487211 E-107*I

The library syntax is

.**intcirc**(void *E, GEN (*eval)(void*,GEN), 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*. Endpoints 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)(void*,GEN), GEN a, GEN b, GEN z, GEN tab, long prec)

Numerical
integration of *expr*(X)exp(-2iPi zX) from a to b, in other words
Fourier transform (from a to b) of the function represented by
*expr*. Note the minus sign. Endpoints 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)(void*,GEN), 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*. Endpoints 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)(void*,GEN), GEN a, GEN b, GEN z, GEN tab, long prec)

Initialize 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

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

syntax.

Numerical integration of (2iPi)^{-1}*expr*(X)e^{Xz} with respect
to X on the line Re(X) = sig. 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
abscissa of integration, and then the integrand 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 abscissa 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)(void*,GEN), GEN sig,GEN z, GEN tab, long prec)

Numerical
integration of (2iPi)^{-1}*expr*(X)z^{-X} with respect to X on the
line Re(X) = sig, 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
abscissa of integration, and then the 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)(void*,GEN), GEN sig, GEN z, GEN tab, long prec)

Numerical integration
of (2iPi)^{-1}s(X)z^{-X} with respect to X on the line Re(X) = sig.
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 `GEN `

.**intmellininvshort**(GEN sig, GEN z, GEN tab, long prec)

Numerical integration
of *expr* on ]a,b[ with respect to X. The integrand may have values
belonging to a vector space over the real numbers; in particular, it can be
complex-valued or vector-valued. But it is assumed that the function is regular
on ]a,b[. If the endpoints a and b are finite and the function is regular
there, the situation is simple:

? intnum(x = 0,1, x^2) %1 = 0.3333333333333333333333333333 ? intnum(x = 0,Pi/2, [cos(x), sin(x)]) %2 = [1.000000000000000000000000000, 1.000000000000000000000000000]

An endpoint equal to ± oo is coded as the single-component vector
[±1]. You are welcome to set, e.g `oo = [1]`

or `INFINITY = [1]`

,
then using `+oo`

, `-oo`

, `-INFINITY`

, etc. will have the expected
behavior.

? oo = [1]; \\ for clarity ? intnum(x = 1,+oo, 1/x^2) %2 = 1.000000000000000000000000000

In basic usage, it is assumed that the function does not decrease exponentially fast at infinity:

? intnum(x=0,+oo, exp(-x)) *** at top-level: intnum(x=0,+oo,exp(- *** ^-------------------- *** exp: exponent (expo) overflow

We shall see in a moment how to avoid the last problem, after describing
the last argument *tab*, which is both optional and technical. The
routine uses weights, which are mostly independent of the function being
integrated, evaluated at many sampling points. If *tab* is

***** a positive integer m, we use 2^m sampling points, hopefully
increasing accuracy. But note that the running 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. If *tab* is omitted, the algorithm guesses
a reasonable value for m depending on the current precision only, which
should be sufficient for regular functions. That value may be obtained from
`intnumstep`

, and increased in case of difficulties.

***** a set of integration tables as output by `intnuminit`

,
they are used directly. This is useful if several integrations of the same
type are performed (on the same kind of interval and functions, for a given
accuracy), in particular for multivariate integrals, since we then skip
expensive precomputations.

**Specifying the behavior at endpoints.**
This is done as follows. An endpoint a is either given as such (a scalar,
real or complex, or [±1] for ± oo ), or as a two component vector
[a,alpha], to indicate the behavior of the integrand in a neighborhood
of a.

If a is finite, the code [a,alpha] means the function has a singularity of the form (x-a)^{alpha}, up to logarithms. (If alpha \ge 0, we only assume the function is regular, which is the default assumption.) If a wrong singularity exponent is used, the result will lose a catastrophic number of decimals:

? intnum(x=0, 1, x^(-1/2)) \\ assume x^{-1/2} is regular at 0 %1 = 1.999999999999999999990291881 ? intnum(x=[0,-1/2], 1, x^(-1/2)) \\ no, it's not %2 = 2.000000000000000000000000000 ? intnum(x=[0,-1/10], 1, x^(-1/2)) %3 = 1.999999999999999999999946438 \\ using a wrong exponent is bad

If a is ± oo , which is coded as [± 1], the situation is more complicated, and [[±1],alpha] means:

***** alpha = 0 (or no alpha at all, i.e. simply [±1]) assumes that the
integrand 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). This includes oscillating but quickly
decreasing functions such as exp(-x)sin(x).

? oo = [1]; ? intnum(x=0, +oo, exp(-2*x)) *** at top-level: intnum(x=0,+oo,exp(- *** ^-------------------- *** exp: exponent (expo) overflow ? intnum(x=0, [+oo, 2], exp(-2*x)) %1 = 0.5000000000000000000000000000 \\ OK! ? intnum(x=0, [+oo, 4], exp(-2*x)) %2 = 0.4999999999999999999961990984 \\ wrong exponent ==> imprecise result ? intnum(x=0, [+oo, 20], exp(-2*x)) %2 = 0.4999524997739071283804510227 \\ disaster

***** 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 non-oscillating function tending slowly to 0 (e.g. like a negative power of x), then

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

***** alpha = -k* I 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).

**Note.** If 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],k* I], but a good rule of thumb is that
if the oscillations are much weaker than the exponential decrease, choose
[[±1],alpha], otherwise choose [[±1],k* I], 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 integrand is almost always a
product of an exponentially decreasing and an oscillating factor. 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.

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 67 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 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

,
where an omitted **intnum**(void *E, GEN (*eval)(void*,GEN), GEN a,GEN b,GEN tab, long prec)*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 abscissa and weight
corresponding to t = 0 (phi(0) and phi'(0)). *tab*[4] and
*tab*[5] contain the abscissas 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 abscissas 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 `GEN `

.**intnuminit**(GEN a, GEN b, long m, long prec)

Initialize tables for integrations from a to b using abscissas
ph(t) and weights ph'(t). Note that there is no equal sign after the
variable name t since t always goes from - oo to + oo , but it
is ph(t) which goes from a to b, and this is not checked. If *flag* = 1
or 2, multiply the reserved table length by 4^{*flag*}, to avoid corresponding
error.

The library syntax is **intnuminitgen**(void *E, GEN (*eval)(void*,GEN), GEN a, GEN b, long m, long flag, long prec)

Numerical integration of *expr* (smooth in ]a,b[), with respect to
X. Suitable for low accuracy; if *expr* is very regular (e.g. analytic
in a large region) and high accuracy is desired, try `intnum`

first.

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 an appropriate huge number. However beware of real underflow
when dealing with rapidly decreasing functions. For example, in order to
compute the int_0^ oo e^{-x^2}dx to 28 decimal digits, then one can
set infinity equal to 10 for example, and certainly not to `1e1000`

.

The library syntax is

,
where **intnumromb**(void *E, GEN (*eval)(void*,GEN), GEN a, GEN b, long flag, long prec)`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 `long `

.**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 = 128 ms. ? prod(i=1, 100, 1 - X^i, 1 + O(X^101)) time = 8 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)

Of course, in this specific case, it is faster to use `eta`

,
which is computed using Euler's formula.

? prod(i=1, 1000, 1 - X^i, 1 + O(X^1001)); time = 589 ms. ? \ps1000 seriesprecision = 1000 significant terms ? eta(X) - % time = 8ms. %4 = O(X^1001)

The library syntax is

.**produit**(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)(void*,GEN), 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. In particular, non-convergent products result in infinite
loops. 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)(void*,GEN), 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. (You will get an error message
`roots must be bracketed in solve`

if this does not hold.)
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)(void*,GEN),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, 10^4, 1/i); \\ rational number: denominator has 4345 digits. time = 236 ms. ? sum(i=1, 5000, 1/i, 0.) time = 8 ms. %2 = 9.787606036044382264178477904

The library syntax is

.**somme**(GEN a, GEN b, char *expr, GEN x)

Numerical summation of the series *expr*, which should be an
alternating series, the formal variable X starting at a. Use an
algorithm of Cohen, Villegas and Zagier (*Experiment. Math.* **9**
(2000), no. 1, 3--12).

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

The routine is heuristic and a rigorous proof assumes that the values of
*expr* are the moments of a positive measure on [0,1]. 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**]). It should be used to try and guess the value of
an infinite sum. (However, see the example at the end of
Section [**Label: se:userfundef**].)

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) \\ Had to hit < C-C > *** at top-level: 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

. Also
available is **sumalt**(void *E, GEN (*eval)(void*,GEN),GEN a,long prec)`sumalt2`

with the same arguments (*flag* = 1).

Sum of expression *expr* over the positive divisors of n.
This function is a trivial wrapper essentially equivalent to

D = divisors(n); for (i = 1, #D, X = D[i]; eval(expr))(except that

`X`

is lexically scoped to the `sumdiv`

loop). If `sumdivmult`

.
Sum of *multiplicative* expression *expr* over the positive
divisors d of n. Assume that *expr* evaluates to f(d)
where f is multiplicative: f(1) = 1 and f(ab) = f(a)f(b) for coprime
a and b.

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) \\ Had to hit < C-C > *** at top-level: 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)(void*,GEN), 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 %3 = -1.321115060 E320 + 0.E311*I \\ nonsense ? sumnum(n=1, [2,log(2)], 2^(-n), /*omitted*/, 1) - 1 \\ of real type time = 5,860 ms. %4 = -1.5... E-236 \\ slow and lost 70 decimals ? m = intnumstep() %5 = 9 ? sumnum(n=1,[2,log(2)], 2^(-n), m+1, 1) - 1 time = 11,770 ms. %6 = -1.9... E-305 \\ now perfect, but slow.

The library syntax is

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

Numerical
summation of (-1)^X*expr*(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.

**Warning 2.** 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); \\ abscissa 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)(void*,GEN), 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
abscissa of integration coded as in `sumnum`

, and m is as in
`intnuminit`

.

The library syntax is `GEN `

.**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, then we use `sumalt`

. For regular functions, the
function `sumnum`

is in general much faster once the initializations
have been made using `sumnuminit`

.

The routine is heuristic and assumes that *expr* is more or less a
decreasing function of X. In particular, the result will be completely
wrong if *expr* is 0 too often. We do not check either that all terms
have the same sign. As `sumalt`

, this function should be used to
try and guess the value of an infinite sum.

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

The library syntax is

. Also
available is **sumpos**(void *E, GEN (*eval)(void*,GEN),GEN a,long prec)`sumpos2`

with the same arguments (*flag* = 1).