Coverage for /usr/lib/python3/dist-packages/sympy/polys/orthopolys.py: 27%
116 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1"""Efficient functions for generating orthogonal polynomials."""
2from sympy.core.symbol import Dummy
3from sympy.polys.densearith import (dup_mul, dup_mul_ground,
4 dup_lshift, dup_sub, dup_add)
5from sympy.polys.domains import ZZ, QQ
6from sympy.polys.polytools import named_poly
7from sympy.utilities import public
9def dup_jacobi(n, a, b, K):
10 """Low-level implementation of Jacobi polynomials."""
11 if n < 1:
12 return [K.one]
13 m2, m1 = [K.one], [(a+b)/K(2) + K.one, (a-b)/K(2)]
14 for i in range(2, n+1):
15 den = K(i)*(a + b + i)*(a + b + K(2)*i - K(2))
16 f0 = (a + b + K(2)*i - K.one) * (a*a - b*b) / (K(2)*den)
17 f1 = (a + b + K(2)*i - K.one) * (a + b + K(2)*i - K(2)) * (a + b + K(2)*i) / (K(2)*den)
18 f2 = (a + i - K.one)*(b + i - K.one)*(a + b + K(2)*i) / den
19 p0 = dup_mul_ground(m1, f0, K)
20 p1 = dup_mul_ground(dup_lshift(m1, 1, K), f1, K)
21 p2 = dup_mul_ground(m2, f2, K)
22 m2, m1 = m1, dup_sub(dup_add(p0, p1, K), p2, K)
23 return m1
25@public
26def jacobi_poly(n, a, b, x=None, polys=False):
27 r"""Generates the Jacobi polynomial `P_n^{(a,b)}(x)`.
29 Parameters
30 ==========
32 n : int
33 Degree of the polynomial.
34 a
35 Lower limit of minimal domain for the list of coefficients.
36 b
37 Upper limit of minimal domain for the list of coefficients.
38 x : optional
39 polys : bool, optional
40 If True, return a Poly, otherwise (default) return an expression.
41 """
42 return named_poly(n, dup_jacobi, None, "Jacobi polynomial", (x, a, b), polys)
44def dup_gegenbauer(n, a, K):
45 """Low-level implementation of Gegenbauer polynomials."""
46 if n < 1:
47 return [K.one]
48 m2, m1 = [K.one], [K(2)*a, K.zero]
49 for i in range(2, n+1):
50 p1 = dup_mul_ground(dup_lshift(m1, 1, K), K(2)*(a-K.one)/K(i) + K(2), K)
51 p2 = dup_mul_ground(m2, K(2)*(a-K.one)/K(i) + K.one, K)
52 m2, m1 = m1, dup_sub(p1, p2, K)
53 return m1
55def gegenbauer_poly(n, a, x=None, polys=False):
56 r"""Generates the Gegenbauer polynomial `C_n^{(a)}(x)`.
58 Parameters
59 ==========
61 n : int
62 Degree of the polynomial.
63 x : optional
64 a
65 Decides minimal domain for the list of coefficients.
66 polys : bool, optional
67 If True, return a Poly, otherwise (default) return an expression.
68 """
69 return named_poly(n, dup_gegenbauer, None, "Gegenbauer polynomial", (x, a), polys)
71def dup_chebyshevt(n, K):
72 """Low-level implementation of Chebyshev polynomials of the first kind."""
73 if n < 1:
74 return [K.one]
75 m2, m1 = [K.one], [K.one, K.zero]
76 for i in range(2, n+1):
77 m2, m1 = m1, dup_sub(dup_mul_ground(dup_lshift(m1, 1, K), K(2), K), m2, K)
78 return m1
80def dup_chebyshevu(n, K):
81 """Low-level implementation of Chebyshev polynomials of the second kind."""
82 if n < 1:
83 return [K.one]
84 m2, m1 = [K.one], [K(2), K.zero]
85 for i in range(2, n+1):
86 m2, m1 = m1, dup_sub(dup_mul_ground(dup_lshift(m1, 1, K), K(2), K), m2, K)
87 return m1
89@public
90def chebyshevt_poly(n, x=None, polys=False):
91 r"""Generates the Chebyshev polynomial of the first kind `T_n(x)`.
93 Parameters
94 ==========
96 n : int
97 Degree of the polynomial.
98 x : optional
99 polys : bool, optional
100 If True, return a Poly, otherwise (default) return an expression.
101 """
102 return named_poly(n, dup_chebyshevt, ZZ,
103 "Chebyshev polynomial of the first kind", (x,), polys)
105@public
106def chebyshevu_poly(n, x=None, polys=False):
107 r"""Generates the Chebyshev polynomial of the second kind `U_n(x)`.
109 Parameters
110 ==========
112 n : int
113 Degree of the polynomial.
114 x : optional
115 polys : bool, optional
116 If True, return a Poly, otherwise (default) return an expression.
117 """
118 return named_poly(n, dup_chebyshevu, ZZ,
119 "Chebyshev polynomial of the second kind", (x,), polys)
121def dup_hermite(n, K):
122 """Low-level implementation of Hermite polynomials."""
123 if n < 1:
124 return [K.one]
125 m2, m1 = [K.one], [K(2), K.zero]
126 for i in range(2, n+1):
127 a = dup_lshift(m1, 1, K)
128 b = dup_mul_ground(m2, K(i-1), K)
129 m2, m1 = m1, dup_mul_ground(dup_sub(a, b, K), K(2), K)
130 return m1
132def dup_hermite_prob(n, K):
133 """Low-level implementation of probabilist's Hermite polynomials."""
134 if n < 1:
135 return [K.one]
136 m2, m1 = [K.one], [K.one, K.zero]
137 for i in range(2, n+1):
138 a = dup_lshift(m1, 1, K)
139 b = dup_mul_ground(m2, K(i-1), K)
140 m2, m1 = m1, dup_sub(a, b, K)
141 return m1
143@public
144def hermite_poly(n, x=None, polys=False):
145 r"""Generates the Hermite polynomial `H_n(x)`.
147 Parameters
148 ==========
150 n : int
151 Degree of the polynomial.
152 x : optional
153 polys : bool, optional
154 If True, return a Poly, otherwise (default) return an expression.
155 """
156 return named_poly(n, dup_hermite, ZZ, "Hermite polynomial", (x,), polys)
158@public
159def hermite_prob_poly(n, x=None, polys=False):
160 r"""Generates the probabilist's Hermite polynomial `He_n(x)`.
162 Parameters
163 ==========
165 n : int
166 Degree of the polynomial.
167 x : optional
168 polys : bool, optional
169 If True, return a Poly, otherwise (default) return an expression.
170 """
171 return named_poly(n, dup_hermite_prob, ZZ,
172 "probabilist's Hermite polynomial", (x,), polys)
174def dup_legendre(n, K):
175 """Low-level implementation of Legendre polynomials."""
176 if n < 1:
177 return [K.one]
178 m2, m1 = [K.one], [K.one, K.zero]
179 for i in range(2, n+1):
180 a = dup_mul_ground(dup_lshift(m1, 1, K), K(2*i-1, i), K)
181 b = dup_mul_ground(m2, K(i-1, i), K)
182 m2, m1 = m1, dup_sub(a, b, K)
183 return m1
185@public
186def legendre_poly(n, x=None, polys=False):
187 r"""Generates the Legendre polynomial `P_n(x)`.
189 Parameters
190 ==========
192 n : int
193 Degree of the polynomial.
194 x : optional
195 polys : bool, optional
196 If True, return a Poly, otherwise (default) return an expression.
197 """
198 return named_poly(n, dup_legendre, QQ, "Legendre polynomial", (x,), polys)
200def dup_laguerre(n, alpha, K):
201 """Low-level implementation of Laguerre polynomials."""
202 m2, m1 = [K.zero], [K.one]
203 for i in range(1, n+1):
204 a = dup_mul(m1, [-K.one/K(i), (alpha-K.one)/K(i) + K(2)], K)
205 b = dup_mul_ground(m2, (alpha-K.one)/K(i) + K.one, K)
206 m2, m1 = m1, dup_sub(a, b, K)
207 return m1
209@public
210def laguerre_poly(n, x=None, alpha=0, polys=False):
211 r"""Generates the Laguerre polynomial `L_n^{(\alpha)}(x)`.
213 Parameters
214 ==========
216 n : int
217 Degree of the polynomial.
218 x : optional
219 alpha : optional
220 Decides minimal domain for the list of coefficients.
221 polys : bool, optional
222 If True, return a Poly, otherwise (default) return an expression.
223 """
224 return named_poly(n, dup_laguerre, None, "Laguerre polynomial", (x, alpha), polys)
226def dup_spherical_bessel_fn(n, K):
227 """Low-level implementation of fn(n, x)."""
228 if n < 1:
229 return [K.one, K.zero]
230 m2, m1 = [K.one], [K.one, K.zero]
231 for i in range(2, n+1):
232 m2, m1 = m1, dup_sub(dup_mul_ground(dup_lshift(m1, 1, K), K(2*i-1), K), m2, K)
233 return dup_lshift(m1, 1, K)
235def dup_spherical_bessel_fn_minus(n, K):
236 """Low-level implementation of fn(-n, x)."""
237 m2, m1 = [K.one, K.zero], [K.zero]
238 for i in range(2, n+1):
239 m2, m1 = m1, dup_sub(dup_mul_ground(dup_lshift(m1, 1, K), K(3-2*i), K), m2, K)
240 return m1
242def spherical_bessel_fn(n, x=None, polys=False):
243 """
244 Coefficients for the spherical Bessel functions.
246 These are only needed in the jn() function.
248 The coefficients are calculated from:
250 fn(0, z) = 1/z
251 fn(1, z) = 1/z**2
252 fn(n-1, z) + fn(n+1, z) == (2*n+1)/z * fn(n, z)
254 Parameters
255 ==========
257 n : int
258 Degree of the polynomial.
259 x : optional
260 polys : bool, optional
261 If True, return a Poly, otherwise (default) return an expression.
263 Examples
264 ========
266 >>> from sympy.polys.orthopolys import spherical_bessel_fn as fn
267 >>> from sympy import Symbol
268 >>> z = Symbol("z")
269 >>> fn(1, z)
270 z**(-2)
271 >>> fn(2, z)
272 -1/z + 3/z**3
273 >>> fn(3, z)
274 -6/z**2 + 15/z**4
275 >>> fn(4, z)
276 1/z - 45/z**3 + 105/z**5
278 """
279 if x is None:
280 x = Dummy("x")
281 f = dup_spherical_bessel_fn_minus if n < 0 else dup_spherical_bessel_fn
282 return named_poly(abs(n), f, ZZ, "", (QQ(1)/x,), polys)