Coverage for /usr/lib/python3/dist-packages/sympy/polys/appellseqs.py: 37%
54 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1r"""
2Efficient functions for generating Appell sequences.
4An Appell sequence is a zero-indexed sequence of polynomials `p_i(x)`
5satisfying `p_{i+1}'(x)=(i+1)p_i(x)` for all `i`. This definition leads
6to the following iterative algorithm:
8.. math :: p_0(x) = c_0,\ p_i(x) = i \int_0^x p_{i-1}(t)\,dt + c_i
10The constant coefficients `c_i` are usually determined from the
11just-evaluated integral and `i`.
13Appell sequences satisfy the following identity from umbral calculus:
15.. math :: p_n(x+y) = \sum_{k=0}^n \binom{n}{k} p_k(x) y^{n-k}
17References
18==========
20.. [1] https://en.wikipedia.org/wiki/Appell_sequence
21.. [2] Peter Luschny, "An introduction to the Bernoulli function",
22 https://arxiv.org/abs/2009.06743
23"""
24from sympy.polys.densearith import dup_mul_ground, dup_sub_ground, dup_quo_ground
25from sympy.polys.densetools import dup_eval, dup_integrate
26from sympy.polys.domains import ZZ, QQ
27from sympy.polys.polytools import named_poly
28from sympy.utilities import public
30def dup_bernoulli(n, K):
31 """Low-level implementation of Bernoulli polynomials."""
32 if n < 1:
33 return [K.one]
34 p = [K.one, K(-1,2)]
35 for i in range(2, n+1):
36 p = dup_integrate(dup_mul_ground(p, K(i), K), 1, K)
37 if i % 2 == 0:
38 p = dup_sub_ground(p, dup_eval(p, K(1,2), K) * K(1<<(i-1), (1<<i)-1), K)
39 return p
41@public
42def bernoulli_poly(n, x=None, polys=False):
43 r"""Generates the Bernoulli polynomial `\operatorname{B}_n(x)`.
45 `\operatorname{B}_n(x)` is the unique polynomial satisfying
47 .. math :: \int_{x}^{x+1} \operatorname{B}_n(t) \,dt = x^n.
49 Based on this, we have for nonnegative integer `s` and integer
50 `a` and `b`
52 .. math :: \sum_{k=a}^{b} k^s = \frac{\operatorname{B}_{s+1}(b+1) -
53 \operatorname{B}_{s+1}(a)}{s+1}
55 which is related to Jakob Bernoulli's original motivation for introducing
56 the Bernoulli numbers, the values of these polynomials at `x = 1`.
58 Examples
59 ========
61 >>> from sympy import summation
62 >>> from sympy.abc import x
63 >>> from sympy.polys import bernoulli_poly
64 >>> bernoulli_poly(5, x)
65 x**5 - 5*x**4/2 + 5*x**3/3 - x/6
67 >>> def psum(p, a, b):
68 ... return (bernoulli_poly(p+1,b+1) - bernoulli_poly(p+1,a)) / (p+1)
69 >>> psum(4, -6, 27)
70 3144337
71 >>> summation(x**4, (x, -6, 27))
72 3144337
74 >>> psum(1, 1, x).factor()
75 x*(x + 1)/2
76 >>> psum(2, 1, x).factor()
77 x*(x + 1)*(2*x + 1)/6
78 >>> psum(3, 1, x).factor()
79 x**2*(x + 1)**2/4
81 Parameters
82 ==========
84 n : int
85 Degree of the polynomial.
86 x : optional
87 polys : bool, optional
88 If True, return a Poly, otherwise (default) return an expression.
90 See Also
91 ========
93 sympy.functions.combinatorial.numbers.bernoulli
95 References
96 ==========
98 .. [1] https://en.wikipedia.org/wiki/Bernoulli_polynomials
99 """
100 return named_poly(n, dup_bernoulli, QQ, "Bernoulli polynomial", (x,), polys)
102def dup_bernoulli_c(n, K):
103 """Low-level implementation of central Bernoulli polynomials."""
104 p = [K.one]
105 for i in range(1, n+1):
106 p = dup_integrate(dup_mul_ground(p, K(i), K), 1, K)
107 if i % 2 == 0:
108 p = dup_sub_ground(p, dup_eval(p, K.one, K) * K((1<<(i-1))-1, (1<<i)-1), K)
109 return p
111@public
112def bernoulli_c_poly(n, x=None, polys=False):
113 r"""Generates the central Bernoulli polynomial `\operatorname{B}_n^c(x)`.
115 These are scaled and shifted versions of the plain Bernoulli polynomials,
116 done in such a way that `\operatorname{B}_n^c(x)` is an even or odd function
117 for even or odd `n` respectively:
119 .. math :: \operatorname{B}_n^c(x) = 2^n \operatorname{B}_n
120 \left(\frac{x+1}{2}\right)
122 Parameters
123 ==========
125 n : int
126 Degree of the polynomial.
127 x : optional
128 polys : bool, optional
129 If True, return a Poly, otherwise (default) return an expression.
130 """
131 return named_poly(n, dup_bernoulli_c, QQ, "central Bernoulli polynomial", (x,), polys)
133def dup_genocchi(n, K):
134 """Low-level implementation of Genocchi polynomials."""
135 if n < 1:
136 return [K.zero]
137 p = [-K.one]
138 for i in range(2, n+1):
139 p = dup_integrate(dup_mul_ground(p, K(i), K), 1, K)
140 if i % 2 == 0:
141 p = dup_sub_ground(p, dup_eval(p, K.one, K) // K(2), K)
142 return p
144@public
145def genocchi_poly(n, x=None, polys=False):
146 r"""Generates the Genocchi polynomial `\operatorname{G}_n(x)`.
148 `\operatorname{G}_n(x)` is twice the difference between the plain and
149 central Bernoulli polynomials, so has degree `n-1`:
151 .. math :: \operatorname{G}_n(x) = 2 (\operatorname{B}_n(x) -
152 \operatorname{B}_n^c(x))
154 The factor of 2 in the definition endows `\operatorname{G}_n(x)` with
155 integer coefficients.
157 Parameters
158 ==========
160 n : int
161 Degree of the polynomial plus one.
162 x : optional
163 polys : bool, optional
164 If True, return a Poly, otherwise (default) return an expression.
166 See Also
167 ========
169 sympy.functions.combinatorial.numbers.genocchi
170 """
171 return named_poly(n, dup_genocchi, ZZ, "Genocchi polynomial", (x,), polys)
173def dup_euler(n, K):
174 """Low-level implementation of Euler polynomials."""
175 return dup_quo_ground(dup_genocchi(n+1, ZZ), K(-n-1), K)
177@public
178def euler_poly(n, x=None, polys=False):
179 r"""Generates the Euler polynomial `\operatorname{E}_n(x)`.
181 These are scaled and reindexed versions of the Genocchi polynomials:
183 .. math :: \operatorname{E}_n(x) = -\frac{\operatorname{G}_{n+1}(x)}{n+1}
185 Parameters
186 ==========
188 n : int
189 Degree of the polynomial.
190 x : optional
191 polys : bool, optional
192 If True, return a Poly, otherwise (default) return an expression.
194 See Also
195 ========
197 sympy.functions.combinatorial.numbers.euler
198 """
199 return named_poly(n, dup_euler, QQ, "Euler polynomial", (x,), polys)
201def dup_andre(n, K):
202 """Low-level implementation of Andre polynomials."""
203 p = [K.one]
204 for i in range(1, n+1):
205 p = dup_integrate(dup_mul_ground(p, K(i), K), 1, K)
206 if i % 2 == 0:
207 p = dup_sub_ground(p, dup_eval(p, K.one, K), K)
208 return p
210@public
211def andre_poly(n, x=None, polys=False):
212 r"""Generates the Andre polynomial `\mathcal{A}_n(x)`.
214 This is the Appell sequence where the constant coefficients form the sequence
215 of Euler numbers ``euler(n)``. As such they have integer coefficients
216 and parities matching the parity of `n`.
218 Luschny calls these the *Swiss-knife polynomials* because their values
219 at 0 and 1 can be simply transformed into both the Bernoulli and Euler
220 numbers. Here they are called the Andre polynomials because
221 `|\mathcal{A}_n(n\bmod 2)|` for `n \ge 0` generates what Luschny calls
222 the *Andre numbers*, A000111 in the OEIS.
224 Examples
225 ========
227 >>> from sympy import bernoulli, euler, genocchi
228 >>> from sympy.abc import x
229 >>> from sympy.polys import andre_poly
230 >>> andre_poly(9, x)
231 x**9 - 36*x**7 + 630*x**5 - 5124*x**3 + 12465*x
233 >>> [andre_poly(n, 0) for n in range(11)]
234 [1, 0, -1, 0, 5, 0, -61, 0, 1385, 0, -50521]
235 >>> [euler(n) for n in range(11)]
236 [1, 0, -1, 0, 5, 0, -61, 0, 1385, 0, -50521]
237 >>> [andre_poly(n-1, 1) * n / (4**n - 2**n) for n in range(1, 11)]
238 [1/2, 1/6, 0, -1/30, 0, 1/42, 0, -1/30, 0, 5/66]
239 >>> [bernoulli(n) for n in range(1, 11)]
240 [1/2, 1/6, 0, -1/30, 0, 1/42, 0, -1/30, 0, 5/66]
241 >>> [-andre_poly(n-1, -1) * n / (-2)**(n-1) for n in range(1, 11)]
242 [-1, -1, 0, 1, 0, -3, 0, 17, 0, -155]
243 >>> [genocchi(n) for n in range(1, 11)]
244 [-1, -1, 0, 1, 0, -3, 0, 17, 0, -155]
246 >>> [abs(andre_poly(n, n%2)) for n in range(11)]
247 [1, 1, 1, 2, 5, 16, 61, 272, 1385, 7936, 50521]
249 Parameters
250 ==========
252 n : int
253 Degree of the polynomial.
254 x : optional
255 polys : bool, optional
256 If True, return a Poly, otherwise (default) return an expression.
258 See Also
259 ========
261 sympy.functions.combinatorial.numbers.andre
263 References
264 ==========
266 .. [1] Peter Luschny, "An introduction to the Bernoulli function",
267 https://arxiv.org/abs/2009.06743
268 """
269 return named_poly(n, dup_andre, ZZ, "Andre polynomial", (x,), polys)