1Function: sumnummonieninit
2Section: sums
3C-Name: sumnummonieninit
4Prototype: DGDGDGp
5Help: sumnummonieninit({asymp},{w},{n0 = 1}): initialize tables for Monien summation of a series with positive terms.
6Doc: initialize tables for Monien summation of a series $\sum_{n\geq n_0}
7 f(n)$ where $f(1/z)$ has a complex analytic continuation in a (complex)
8 neighbourhood of the segment $[0,1]$.
9
10 By default, assume that $f(n) = O(n^{-2})$ and has a nonzero asymptotic
11 expansion
12 $$f(n) = \sum_{i\geq 2} a_i / n^i$$
13 at infinity. Note that the sum starts at $i = 2$! The argument \kbd{asymp}
14 allows to specify different expansions:
15
16 \item a real number $\beta > 0$ means
17  $$f(n) = \sum_{i\geq 1} a_i / n^{i + \beta}$$
18 (Now the summation starts at $1$.)
19
20 \item a vector $[\alpha,\beta]$ of reals, where we must have $\alpha > 0$
21 and $\alpha + \beta > 1$ to ensure convergence, means that
22  $$f(n) = \sum_{i\geq 1} a_i / n^{\alpha i + \beta}$$
23 Note that $\kbd{asymp} = [1, \beta]$ is equivalent to
24 $\kbd{asymp}=\beta$.
25
26 \bprog
27 ? \p57
28 ? s = sumnum(n = 1, sin(1/sqrt(n)) / n); \\ reference point
29
30 ? \p38
31 ? sumnummonien(n = 1, sin(1/sqrt(n)) / n) - s
32 %2 = -0.001[...] \\ completely wrong
33
34 ? t = sumnummonieninit(1/2);  \\ f(n) = sum_i 1 / n^(i+1/2)
35 ? sumnummonien(n = 1, sin(1/sqrt(n)) / n, t) - s
36 %3 = 0.E-37 \\ now correct
37 @eprog\noindent (As a matter of fact, in the above summation, the
38 result given by \kbd{sumnum} at \kbd{\bs p38} is slighly incorrect,
39 so we had to increase the accuracy to \kbd{\bs p57}.)
40
41 The argument $w$ is used to sum expressions of the form
42 $$ \sum_{n\geq n_0} f(n) w(n),$$
43 for varying $f$ \emph{as above}, and fixed weight function $w$, where we
44 further assume that the auxiliary sums
45 $$g_w(m) = \sum_{n\geq n_0} w(n) / n^{\alpha m + \beta} $$
46 converge for all $m\geq 1$. Note that for nonnegative integers $k$,
47 and weight $w(n) = (\log n)^k$, the function $g_w(m) = \zeta^{(k)}(\alpha m +
48 \beta)$ has a simple expression; for general weights, $g_w$ is
49 computed using \kbd{sumnum}. The following variants are available
50
51 \item an integer $k \geq 0$, to code $w(n) = (\log n)^k$;
52
53 \item a \typ{CLOSURE} computing the values $w(n)$, where we
54 assume that $w(n) = O(n^\epsilon)$ for all $\epsilon > 0$;
55
56 \item a vector $[w, \kbd{fast}]$, where $w$ is a closure as above
57 and \kbd{fast} is a scalar;
58 we assume that $w(n) = O(n^{\kbd{fast}+\epsilon})$; note that
59 $\kbd{w} = [w, 0]$ is equivalent to $\kbd{w} = w$. Note that if
60 $w$ decreases exponentially, \kbd{suminf} should be used instead.
61
62 The subsequent calls to \kbd{sumnummonien} \emph{must} use the same value
63 of $n_0$ as was used here.
64 \bprog
65 ? \p300
66 ? sumnummonien(n = 1, n^-2*log(n)) + zeta'(2)
67 time = 328 ms.
68 %1 = -1.323[...]E-6 \\ completely wrong, f does not satisfy hypotheses !
69 ? tab = sumnummonieninit(, 1); \\ codes w(n) = log(n)
70 time = 3,993 ms.
71 ? sumnummonien(n = 1, n^-2, tab) + zeta'(2)
72 time = 41 ms.
73 %3 = -5.562684646268003458 E-309  \\ now perfect
74
75 ? tab = sumnummonieninit(, n->log(n)); \\ generic, slower
76 time = 9,808 ms.
77 ? sumnummonien(n = 1, n^-2, tab) + zeta'(2)
78 time = 40 ms.
79 %5 = -5.562684646268003458 E-309  \\ identical result
80 @eprog
81