Coverage for /usr/lib/python3/dist-packages/scipy/special/_basic.py: 13%
583 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1#
2# Author: Travis Oliphant, 2002
3#
5import operator
6import numpy as np
7import math
8import warnings
9from numpy import (pi, asarray, floor, isscalar, iscomplex, real,
10 imag, sqrt, where, mgrid, sin, place, issubdtype,
11 extract, inexact, nan, zeros, sinc)
12from . import _ufuncs
13from ._ufuncs import (mathieu_a, mathieu_b, iv, jv, gamma,
14 psi, hankel1, hankel2, yv, kv, poch, binom)
15from . import _specfun
16from ._comb import _comb_int
19__all__ = [
20 'ai_zeros',
21 'assoc_laguerre',
22 'bei_zeros',
23 'beip_zeros',
24 'ber_zeros',
25 'bernoulli',
26 'berp_zeros',
27 'bi_zeros',
28 'clpmn',
29 'comb',
30 'digamma',
31 'diric',
32 'erf_zeros',
33 'euler',
34 'factorial',
35 'factorial2',
36 'factorialk',
37 'fresnel_zeros',
38 'fresnelc_zeros',
39 'fresnels_zeros',
40 'h1vp',
41 'h2vp',
42 'ivp',
43 'jn_zeros',
44 'jnjnp_zeros',
45 'jnp_zeros',
46 'jnyn_zeros',
47 'jvp',
48 'kei_zeros',
49 'keip_zeros',
50 'kelvin_zeros',
51 'ker_zeros',
52 'kerp_zeros',
53 'kvp',
54 'lmbda',
55 'lpmn',
56 'lpn',
57 'lqmn',
58 'lqn',
59 'mathieu_even_coef',
60 'mathieu_odd_coef',
61 'obl_cv_seq',
62 'pbdn_seq',
63 'pbdv_seq',
64 'pbvv_seq',
65 'perm',
66 'polygamma',
67 'pro_cv_seq',
68 'riccati_jn',
69 'riccati_yn',
70 'sinc',
71 'y0_zeros',
72 'y1_zeros',
73 'y1p_zeros',
74 'yn_zeros',
75 'ynp_zeros',
76 'yvp',
77 'zeta'
78]
81# mapping k to last n such that factorialk(n, k) < np.iinfo(np.int64).max
82_FACTORIALK_LIMITS_64BITS = {1: 20, 2: 33, 3: 44, 4: 54, 5: 65,
83 6: 74, 7: 84, 8: 93, 9: 101}
84# mapping k to last n such that factorialk(n, k) < np.iinfo(np.int32).max
85_FACTORIALK_LIMITS_32BITS = {1: 12, 2: 19, 3: 25, 4: 31, 5: 37,
86 6: 43, 7: 47, 8: 51, 9: 56}
89def _nonneg_int_or_fail(n, var_name, strict=True):
90 try:
91 if strict:
92 # Raises an exception if float
93 n = operator.index(n)
94 elif n == floor(n):
95 n = int(n)
96 else:
97 raise ValueError()
98 if n < 0:
99 raise ValueError()
100 except (ValueError, TypeError) as err:
101 raise err.__class__(f"{var_name} must be a non-negative integer") from err
102 return n
105def diric(x, n):
106 """Periodic sinc function, also called the Dirichlet function.
108 The Dirichlet function is defined as::
110 diric(x, n) = sin(x * n/2) / (n * sin(x / 2)),
112 where `n` is a positive integer.
114 Parameters
115 ----------
116 x : array_like
117 Input data
118 n : int
119 Integer defining the periodicity.
121 Returns
122 -------
123 diric : ndarray
125 Examples
126 --------
127 >>> import numpy as np
128 >>> from scipy import special
129 >>> import matplotlib.pyplot as plt
131 >>> x = np.linspace(-8*np.pi, 8*np.pi, num=201)
132 >>> plt.figure(figsize=(8, 8));
133 >>> for idx, n in enumerate([2, 3, 4, 9]):
134 ... plt.subplot(2, 2, idx+1)
135 ... plt.plot(x, special.diric(x, n))
136 ... plt.title('diric, n={}'.format(n))
137 >>> plt.show()
139 The following example demonstrates that `diric` gives the magnitudes
140 (modulo the sign and scaling) of the Fourier coefficients of a
141 rectangular pulse.
143 Suppress output of values that are effectively 0:
145 >>> np.set_printoptions(suppress=True)
147 Create a signal `x` of length `m` with `k` ones:
149 >>> m = 8
150 >>> k = 3
151 >>> x = np.zeros(m)
152 >>> x[:k] = 1
154 Use the FFT to compute the Fourier transform of `x`, and
155 inspect the magnitudes of the coefficients:
157 >>> np.abs(np.fft.fft(x))
158 array([ 3. , 2.41421356, 1. , 0.41421356, 1. ,
159 0.41421356, 1. , 2.41421356])
161 Now find the same values (up to sign) using `diric`. We multiply
162 by `k` to account for the different scaling conventions of
163 `numpy.fft.fft` and `diric`:
165 >>> theta = np.linspace(0, 2*np.pi, m, endpoint=False)
166 >>> k * special.diric(theta, k)
167 array([ 3. , 2.41421356, 1. , -0.41421356, -1. ,
168 -0.41421356, 1. , 2.41421356])
169 """
170 x, n = asarray(x), asarray(n)
171 n = asarray(n + (x-x))
172 x = asarray(x + (n-n))
173 if issubdtype(x.dtype, inexact):
174 ytype = x.dtype
175 else:
176 ytype = float
177 y = zeros(x.shape, ytype)
179 # empirical minval for 32, 64 or 128 bit float computations
180 # where sin(x/2) < minval, result is fixed at +1 or -1
181 if np.finfo(ytype).eps < 1e-18:
182 minval = 1e-11
183 elif np.finfo(ytype).eps < 1e-15:
184 minval = 1e-7
185 else:
186 minval = 1e-3
188 mask1 = (n <= 0) | (n != floor(n))
189 place(y, mask1, nan)
191 x = x / 2
192 denom = sin(x)
193 mask2 = (1-mask1) & (abs(denom) < minval)
194 xsub = extract(mask2, x)
195 nsub = extract(mask2, n)
196 zsub = xsub / pi
197 place(y, mask2, pow(-1, np.round(zsub)*(nsub-1)))
199 mask = (1-mask1) & (1-mask2)
200 xsub = extract(mask, x)
201 nsub = extract(mask, n)
202 dsub = extract(mask, denom)
203 place(y, mask, sin(nsub*xsub)/(nsub*dsub))
204 return y
207def jnjnp_zeros(nt):
208 """Compute zeros of integer-order Bessel functions Jn and Jn'.
210 Results are arranged in order of the magnitudes of the zeros.
212 Parameters
213 ----------
214 nt : int
215 Number (<=1200) of zeros to compute
217 Returns
218 -------
219 zo[l-1] : ndarray
220 Value of the lth zero of Jn(x) and Jn'(x). Of length `nt`.
221 n[l-1] : ndarray
222 Order of the Jn(x) or Jn'(x) associated with lth zero. Of length `nt`.
223 m[l-1] : ndarray
224 Serial number of the zeros of Jn(x) or Jn'(x) associated
225 with lth zero. Of length `nt`.
226 t[l-1] : ndarray
227 0 if lth zero in zo is zero of Jn(x), 1 if it is a zero of Jn'(x). Of
228 length `nt`.
230 See Also
231 --------
232 jn_zeros, jnp_zeros : to get separated arrays of zeros.
234 References
235 ----------
236 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
237 Functions", John Wiley and Sons, 1996, chapter 5.
238 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
240 """
241 if not isscalar(nt) or (floor(nt) != nt) or (nt > 1200):
242 raise ValueError("Number must be integer <= 1200.")
243 nt = int(nt)
244 n, m, t, zo = _specfun.jdzo(nt)
245 return zo[1:nt+1], n[:nt], m[:nt], t[:nt]
248def jnyn_zeros(n, nt):
249 """Compute nt zeros of Bessel functions Jn(x), Jn'(x), Yn(x), and Yn'(x).
251 Returns 4 arrays of length `nt`, corresponding to the first `nt`
252 zeros of Jn(x), Jn'(x), Yn(x), and Yn'(x), respectively. The zeros
253 are returned in ascending order.
255 Parameters
256 ----------
257 n : int
258 Order of the Bessel functions
259 nt : int
260 Number (<=1200) of zeros to compute
262 Returns
263 -------
264 Jn : ndarray
265 First `nt` zeros of Jn
266 Jnp : ndarray
267 First `nt` zeros of Jn'
268 Yn : ndarray
269 First `nt` zeros of Yn
270 Ynp : ndarray
271 First `nt` zeros of Yn'
273 See Also
274 --------
275 jn_zeros, jnp_zeros, yn_zeros, ynp_zeros
277 References
278 ----------
279 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
280 Functions", John Wiley and Sons, 1996, chapter 5.
281 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
283 Examples
284 --------
285 Compute the first three roots of :math:`J_1`, :math:`J_1'`,
286 :math:`Y_1` and :math:`Y_1'`.
288 >>> from scipy.special import jnyn_zeros
289 >>> jn_roots, jnp_roots, yn_roots, ynp_roots = jnyn_zeros(1, 3)
290 >>> jn_roots, yn_roots
291 (array([ 3.83170597, 7.01558667, 10.17346814]),
292 array([2.19714133, 5.42968104, 8.59600587]))
294 Plot :math:`J_1`, :math:`J_1'`, :math:`Y_1`, :math:`Y_1'` and their roots.
296 >>> import numpy as np
297 >>> import matplotlib.pyplot as plt
298 >>> from scipy.special import jnyn_zeros, jvp, jn, yvp, yn
299 >>> jn_roots, jnp_roots, yn_roots, ynp_roots = jnyn_zeros(1, 3)
300 >>> fig, ax = plt.subplots()
301 >>> xmax= 11
302 >>> x = np.linspace(0, xmax)
303 >>> x[0] += 1e-15
304 >>> ax.plot(x, jn(1, x), label=r"$J_1$", c='r')
305 >>> ax.plot(x, jvp(1, x, 1), label=r"$J_1'$", c='b')
306 >>> ax.plot(x, yn(1, x), label=r"$Y_1$", c='y')
307 >>> ax.plot(x, yvp(1, x, 1), label=r"$Y_1'$", c='c')
308 >>> zeros = np.zeros((3, ))
309 >>> ax.scatter(jn_roots, zeros, s=30, c='r', zorder=5,
310 ... label=r"$J_1$ roots")
311 >>> ax.scatter(jnp_roots, zeros, s=30, c='b', zorder=5,
312 ... label=r"$J_1'$ roots")
313 >>> ax.scatter(yn_roots, zeros, s=30, c='y', zorder=5,
314 ... label=r"$Y_1$ roots")
315 >>> ax.scatter(ynp_roots, zeros, s=30, c='c', zorder=5,
316 ... label=r"$Y_1'$ roots")
317 >>> ax.hlines(0, 0, xmax, color='k')
318 >>> ax.set_ylim(-0.6, 0.6)
319 >>> ax.set_xlim(0, xmax)
320 >>> ax.legend(ncol=2, bbox_to_anchor=(1., 0.75))
321 >>> plt.tight_layout()
322 >>> plt.show()
323 """
324 if not (isscalar(nt) and isscalar(n)):
325 raise ValueError("Arguments must be scalars.")
326 if (floor(n) != n) or (floor(nt) != nt):
327 raise ValueError("Arguments must be integers.")
328 if (nt <= 0):
329 raise ValueError("nt > 0")
330 return _specfun.jyzo(abs(n), nt)
333def jn_zeros(n, nt):
334 r"""Compute zeros of integer-order Bessel functions Jn.
336 Compute `nt` zeros of the Bessel functions :math:`J_n(x)` on the
337 interval :math:`(0, \infty)`. The zeros are returned in ascending
338 order. Note that this interval excludes the zero at :math:`x = 0`
339 that exists for :math:`n > 0`.
341 Parameters
342 ----------
343 n : int
344 Order of Bessel function
345 nt : int
346 Number of zeros to return
348 Returns
349 -------
350 ndarray
351 First `nt` zeros of the Bessel function.
353 See Also
354 --------
355 jv: Real-order Bessel functions of the first kind
356 jnp_zeros: Zeros of :math:`Jn'`
358 References
359 ----------
360 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
361 Functions", John Wiley and Sons, 1996, chapter 5.
362 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
364 Examples
365 --------
366 Compute the first four positive roots of :math:`J_3`.
368 >>> from scipy.special import jn_zeros
369 >>> jn_zeros(3, 4)
370 array([ 6.3801619 , 9.76102313, 13.01520072, 16.22346616])
372 Plot :math:`J_3` and its first four positive roots. Note
373 that the root located at 0 is not returned by `jn_zeros`.
375 >>> import numpy as np
376 >>> import matplotlib.pyplot as plt
377 >>> from scipy.special import jn, jn_zeros
378 >>> j3_roots = jn_zeros(3, 4)
379 >>> xmax = 18
380 >>> xmin = -1
381 >>> x = np.linspace(xmin, xmax, 500)
382 >>> fig, ax = plt.subplots()
383 >>> ax.plot(x, jn(3, x), label=r'$J_3$')
384 >>> ax.scatter(j3_roots, np.zeros((4, )), s=30, c='r',
385 ... label=r"$J_3$_Zeros", zorder=5)
386 >>> ax.scatter(0, 0, s=30, c='k',
387 ... label=r"Root at 0", zorder=5)
388 >>> ax.hlines(0, 0, xmax, color='k')
389 >>> ax.set_xlim(xmin, xmax)
390 >>> plt.legend()
391 >>> plt.show()
392 """
393 return jnyn_zeros(n, nt)[0]
396def jnp_zeros(n, nt):
397 r"""Compute zeros of integer-order Bessel function derivatives Jn'.
399 Compute `nt` zeros of the functions :math:`J_n'(x)` on the
400 interval :math:`(0, \infty)`. The zeros are returned in ascending
401 order. Note that this interval excludes the zero at :math:`x = 0`
402 that exists for :math:`n > 1`.
404 Parameters
405 ----------
406 n : int
407 Order of Bessel function
408 nt : int
409 Number of zeros to return
411 Returns
412 -------
413 ndarray
414 First `nt` zeros of the Bessel function.
416 See Also
417 --------
418 jvp: Derivatives of integer-order Bessel functions of the first kind
419 jv: Float-order Bessel functions of the first kind
421 References
422 ----------
423 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
424 Functions", John Wiley and Sons, 1996, chapter 5.
425 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
427 Examples
428 --------
429 Compute the first four roots of :math:`J_2'`.
431 >>> from scipy.special import jnp_zeros
432 >>> jnp_zeros(2, 4)
433 array([ 3.05423693, 6.70613319, 9.96946782, 13.17037086])
435 As `jnp_zeros` yields the roots of :math:`J_n'`, it can be used to
436 compute the locations of the peaks of :math:`J_n`. Plot
437 :math:`J_2`, :math:`J_2'` and the locations of the roots of :math:`J_2'`.
439 >>> import numpy as np
440 >>> import matplotlib.pyplot as plt
441 >>> from scipy.special import jn, jnp_zeros, jvp
442 >>> j2_roots = jnp_zeros(2, 4)
443 >>> xmax = 15
444 >>> x = np.linspace(0, xmax, 500)
445 >>> fig, ax = plt.subplots()
446 >>> ax.plot(x, jn(2, x), label=r'$J_2$')
447 >>> ax.plot(x, jvp(2, x, 1), label=r"$J_2'$")
448 >>> ax.hlines(0, 0, xmax, color='k')
449 >>> ax.scatter(j2_roots, np.zeros((4, )), s=30, c='r',
450 ... label=r"Roots of $J_2'$", zorder=5)
451 >>> ax.set_ylim(-0.4, 0.8)
452 >>> ax.set_xlim(0, xmax)
453 >>> plt.legend()
454 >>> plt.show()
455 """
456 return jnyn_zeros(n, nt)[1]
459def yn_zeros(n, nt):
460 r"""Compute zeros of integer-order Bessel function Yn(x).
462 Compute `nt` zeros of the functions :math:`Y_n(x)` on the interval
463 :math:`(0, \infty)`. The zeros are returned in ascending order.
465 Parameters
466 ----------
467 n : int
468 Order of Bessel function
469 nt : int
470 Number of zeros to return
472 Returns
473 -------
474 ndarray
475 First `nt` zeros of the Bessel function.
477 See Also
478 --------
479 yn: Bessel function of the second kind for integer order
480 yv: Bessel function of the second kind for real order
482 References
483 ----------
484 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
485 Functions", John Wiley and Sons, 1996, chapter 5.
486 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
488 Examples
489 --------
490 Compute the first four roots of :math:`Y_2`.
492 >>> from scipy.special import yn_zeros
493 >>> yn_zeros(2, 4)
494 array([ 3.38424177, 6.79380751, 10.02347798, 13.20998671])
496 Plot :math:`Y_2` and its first four roots.
498 >>> import numpy as np
499 >>> import matplotlib.pyplot as plt
500 >>> from scipy.special import yn, yn_zeros
501 >>> xmin = 2
502 >>> xmax = 15
503 >>> x = np.linspace(xmin, xmax, 500)
504 >>> fig, ax = plt.subplots()
505 >>> ax.hlines(0, xmin, xmax, color='k')
506 >>> ax.plot(x, yn(2, x), label=r'$Y_2$')
507 >>> ax.scatter(yn_zeros(2, 4), np.zeros((4, )), s=30, c='r',
508 ... label='Roots', zorder=5)
509 >>> ax.set_ylim(-0.4, 0.4)
510 >>> ax.set_xlim(xmin, xmax)
511 >>> plt.legend()
512 >>> plt.show()
513 """
514 return jnyn_zeros(n, nt)[2]
517def ynp_zeros(n, nt):
518 r"""Compute zeros of integer-order Bessel function derivatives Yn'(x).
520 Compute `nt` zeros of the functions :math:`Y_n'(x)` on the
521 interval :math:`(0, \infty)`. The zeros are returned in ascending
522 order.
524 Parameters
525 ----------
526 n : int
527 Order of Bessel function
528 nt : int
529 Number of zeros to return
531 Returns
532 -------
533 ndarray
534 First `nt` zeros of the Bessel derivative function.
537 See Also
538 --------
539 yvp
541 References
542 ----------
543 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
544 Functions", John Wiley and Sons, 1996, chapter 5.
545 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
547 Examples
548 --------
549 Compute the first four roots of the first derivative of the
550 Bessel function of second kind for order 0 :math:`Y_0'`.
552 >>> from scipy.special import ynp_zeros
553 >>> ynp_zeros(0, 4)
554 array([ 2.19714133, 5.42968104, 8.59600587, 11.74915483])
556 Plot :math:`Y_0`, :math:`Y_0'` and confirm visually that the roots of
557 :math:`Y_0'` are located at local extrema of :math:`Y_0`.
559 >>> import numpy as np
560 >>> import matplotlib.pyplot as plt
561 >>> from scipy.special import yn, ynp_zeros, yvp
562 >>> zeros = ynp_zeros(0, 4)
563 >>> xmax = 13
564 >>> x = np.linspace(0, xmax, 500)
565 >>> fig, ax = plt.subplots()
566 >>> ax.plot(x, yn(0, x), label=r'$Y_0$')
567 >>> ax.plot(x, yvp(0, x, 1), label=r"$Y_0'$")
568 >>> ax.scatter(zeros, np.zeros((4, )), s=30, c='r',
569 ... label=r"Roots of $Y_0'$", zorder=5)
570 >>> for root in zeros:
571 ... y0_extremum = yn(0, root)
572 ... lower = min(0, y0_extremum)
573 ... upper = max(0, y0_extremum)
574 ... ax.vlines(root, lower, upper, color='r')
575 >>> ax.hlines(0, 0, xmax, color='k')
576 >>> ax.set_ylim(-0.6, 0.6)
577 >>> ax.set_xlim(0, xmax)
578 >>> plt.legend()
579 >>> plt.show()
580 """
581 return jnyn_zeros(n, nt)[3]
584def y0_zeros(nt, complex=False):
585 """Compute nt zeros of Bessel function Y0(z), and derivative at each zero.
587 The derivatives are given by Y0'(z0) = -Y1(z0) at each zero z0.
589 Parameters
590 ----------
591 nt : int
592 Number of zeros to return
593 complex : bool, default False
594 Set to False to return only the real zeros; set to True to return only
595 the complex zeros with negative real part and positive imaginary part.
596 Note that the complex conjugates of the latter are also zeros of the
597 function, but are not returned by this routine.
599 Returns
600 -------
601 z0n : ndarray
602 Location of nth zero of Y0(z)
603 y0pz0n : ndarray
604 Value of derivative Y0'(z0) for nth zero
606 References
607 ----------
608 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
609 Functions", John Wiley and Sons, 1996, chapter 5.
610 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
612 Examples
613 --------
614 Compute the first 4 real roots and the derivatives at the roots of
615 :math:`Y_0`:
617 >>> import numpy as np
618 >>> from scipy.special import y0_zeros
619 >>> zeros, grads = y0_zeros(4)
620 >>> with np.printoptions(precision=5):
621 ... print(f"Roots: {zeros}")
622 ... print(f"Gradients: {grads}")
623 Roots: [ 0.89358+0.j 3.95768+0.j 7.08605+0.j 10.22235+0.j]
624 Gradients: [-0.87942+0.j 0.40254+0.j -0.3001 +0.j 0.2497 +0.j]
626 Plot the real part of :math:`Y_0` and the first four computed roots.
628 >>> import matplotlib.pyplot as plt
629 >>> from scipy.special import y0
630 >>> xmin = 0
631 >>> xmax = 11
632 >>> x = np.linspace(xmin, xmax, 500)
633 >>> fig, ax = plt.subplots()
634 >>> ax.hlines(0, xmin, xmax, color='k')
635 >>> ax.plot(x, y0(x), label=r'$Y_0$')
636 >>> zeros, grads = y0_zeros(4)
637 >>> ax.scatter(zeros.real, np.zeros((4, )), s=30, c='r',
638 ... label=r'$Y_0$_zeros', zorder=5)
639 >>> ax.set_ylim(-0.5, 0.6)
640 >>> ax.set_xlim(xmin, xmax)
641 >>> plt.legend(ncol=2)
642 >>> plt.show()
644 Compute the first 4 complex roots and the derivatives at the roots of
645 :math:`Y_0` by setting ``complex=True``:
647 >>> y0_zeros(4, True)
648 (array([ -2.40301663+0.53988231j, -5.5198767 +0.54718001j,
649 -8.6536724 +0.54841207j, -11.79151203+0.54881912j]),
650 array([ 0.10074769-0.88196771j, -0.02924642+0.5871695j ,
651 0.01490806-0.46945875j, -0.00937368+0.40230454j]))
652 """
653 if not isscalar(nt) or (floor(nt) != nt) or (nt <= 0):
654 raise ValueError("Arguments must be scalar positive integer.")
655 kf = 0
656 kc = not complex
657 return _specfun.cyzo(nt, kf, kc)
660def y1_zeros(nt, complex=False):
661 """Compute nt zeros of Bessel function Y1(z), and derivative at each zero.
663 The derivatives are given by Y1'(z1) = Y0(z1) at each zero z1.
665 Parameters
666 ----------
667 nt : int
668 Number of zeros to return
669 complex : bool, default False
670 Set to False to return only the real zeros; set to True to return only
671 the complex zeros with negative real part and positive imaginary part.
672 Note that the complex conjugates of the latter are also zeros of the
673 function, but are not returned by this routine.
675 Returns
676 -------
677 z1n : ndarray
678 Location of nth zero of Y1(z)
679 y1pz1n : ndarray
680 Value of derivative Y1'(z1) for nth zero
682 References
683 ----------
684 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
685 Functions", John Wiley and Sons, 1996, chapter 5.
686 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
688 Examples
689 --------
690 Compute the first 4 real roots and the derivatives at the roots of
691 :math:`Y_1`:
693 >>> import numpy as np
694 >>> from scipy.special import y1_zeros
695 >>> zeros, grads = y1_zeros(4)
696 >>> with np.printoptions(precision=5):
697 ... print(f"Roots: {zeros}")
698 ... print(f"Gradients: {grads}")
699 Roots: [ 2.19714+0.j 5.42968+0.j 8.59601+0.j 11.74915+0.j]
700 Gradients: [ 0.52079+0.j -0.34032+0.j 0.27146+0.j -0.23246+0.j]
702 Extract the real parts:
704 >>> realzeros = zeros.real
705 >>> realzeros
706 array([ 2.19714133, 5.42968104, 8.59600587, 11.74915483])
708 Plot :math:`Y_1` and the first four computed roots.
710 >>> import matplotlib.pyplot as plt
711 >>> from scipy.special import y1
712 >>> xmin = 0
713 >>> xmax = 13
714 >>> x = np.linspace(xmin, xmax, 500)
715 >>> zeros, grads = y1_zeros(4)
716 >>> fig, ax = plt.subplots()
717 >>> ax.hlines(0, xmin, xmax, color='k')
718 >>> ax.plot(x, y1(x), label=r'$Y_1$')
719 >>> ax.scatter(zeros.real, np.zeros((4, )), s=30, c='r',
720 ... label=r'$Y_1$_zeros', zorder=5)
721 >>> ax.set_ylim(-0.5, 0.5)
722 >>> ax.set_xlim(xmin, xmax)
723 >>> plt.legend()
724 >>> plt.show()
726 Compute the first 4 complex roots and the derivatives at the roots of
727 :math:`Y_1` by setting ``complex=True``:
729 >>> y1_zeros(4, True)
730 (array([ -0.50274327+0.78624371j, -3.83353519+0.56235654j,
731 -7.01590368+0.55339305j, -10.17357383+0.55127339j]),
732 array([-0.45952768+1.31710194j, 0.04830191-0.69251288j,
733 -0.02012695+0.51864253j, 0.011614 -0.43203296j]))
734 """
735 if not isscalar(nt) or (floor(nt) != nt) or (nt <= 0):
736 raise ValueError("Arguments must be scalar positive integer.")
737 kf = 1
738 kc = not complex
739 return _specfun.cyzo(nt, kf, kc)
742def y1p_zeros(nt, complex=False):
743 """Compute nt zeros of Bessel derivative Y1'(z), and value at each zero.
745 The values are given by Y1(z1) at each z1 where Y1'(z1)=0.
747 Parameters
748 ----------
749 nt : int
750 Number of zeros to return
751 complex : bool, default False
752 Set to False to return only the real zeros; set to True to return only
753 the complex zeros with negative real part and positive imaginary part.
754 Note that the complex conjugates of the latter are also zeros of the
755 function, but are not returned by this routine.
757 Returns
758 -------
759 z1pn : ndarray
760 Location of nth zero of Y1'(z)
761 y1z1pn : ndarray
762 Value of derivative Y1(z1) for nth zero
764 References
765 ----------
766 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
767 Functions", John Wiley and Sons, 1996, chapter 5.
768 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
770 Examples
771 --------
772 Compute the first four roots of :math:`Y_1'` and the values of
773 :math:`Y_1` at these roots.
775 >>> import numpy as np
776 >>> from scipy.special import y1p_zeros
777 >>> y1grad_roots, y1_values = y1p_zeros(4)
778 >>> with np.printoptions(precision=5):
779 ... print(f"Y1' Roots: {y1grad_roots}")
780 ... print(f"Y1 values: {y1_values}")
781 Y1' Roots: [ 3.68302+0.j 6.9415 +0.j 10.1234 +0.j 13.28576+0.j]
782 Y1 values: [ 0.41673+0.j -0.30317+0.j 0.25091+0.j -0.21897+0.j]
784 `y1p_zeros` can be used to calculate the extremal points of :math:`Y_1`
785 directly. Here we plot :math:`Y_1` and the first four extrema.
787 >>> import matplotlib.pyplot as plt
788 >>> from scipy.special import y1, yvp
789 >>> y1_roots, y1_values_at_roots = y1p_zeros(4)
790 >>> real_roots = y1_roots.real
791 >>> xmax = 15
792 >>> x = np.linspace(0, xmax, 500)
793 >>> x[0] += 1e-15
794 >>> fig, ax = plt.subplots()
795 >>> ax.plot(x, y1(x), label=r'$Y_1$')
796 >>> ax.plot(x, yvp(1, x, 1), label=r"$Y_1'$")
797 >>> ax.scatter(real_roots, np.zeros((4, )), s=30, c='r',
798 ... label=r"Roots of $Y_1'$", zorder=5)
799 >>> ax.scatter(real_roots, y1_values_at_roots.real, s=30, c='k',
800 ... label=r"Extrema of $Y_1$", zorder=5)
801 >>> ax.hlines(0, 0, xmax, color='k')
802 >>> ax.set_ylim(-0.5, 0.5)
803 >>> ax.set_xlim(0, xmax)
804 >>> ax.legend(ncol=2, bbox_to_anchor=(1., 0.75))
805 >>> plt.tight_layout()
806 >>> plt.show()
807 """
808 if not isscalar(nt) or (floor(nt) != nt) or (nt <= 0):
809 raise ValueError("Arguments must be scalar positive integer.")
810 kf = 2
811 kc = not complex
812 return _specfun.cyzo(nt, kf, kc)
815def _bessel_diff_formula(v, z, n, L, phase):
816 # from AMS55.
817 # L(v, z) = J(v, z), Y(v, z), H1(v, z), H2(v, z), phase = -1
818 # L(v, z) = I(v, z) or exp(v*pi*i)K(v, z), phase = 1
819 # For K, you can pull out the exp((v-k)*pi*i) into the caller
820 v = asarray(v)
821 p = 1.0
822 s = L(v-n, z)
823 for i in range(1, n+1):
824 p = phase * (p * (n-i+1)) / i # = choose(k, i)
825 s += p*L(v-n + i*2, z)
826 return s / (2.**n)
829def jvp(v, z, n=1):
830 """Compute derivatives of Bessel functions of the first kind.
832 Compute the nth derivative of the Bessel function `Jv` with
833 respect to `z`.
835 Parameters
836 ----------
837 v : array_like or float
838 Order of Bessel function
839 z : complex
840 Argument at which to evaluate the derivative; can be real or
841 complex.
842 n : int, default 1
843 Order of derivative. For 0 returns the Bessel function `jv` itself.
845 Returns
846 -------
847 scalar or ndarray
848 Values of the derivative of the Bessel function.
850 Notes
851 -----
852 The derivative is computed using the relation DLFM 10.6.7 [2]_.
854 References
855 ----------
856 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
857 Functions", John Wiley and Sons, 1996, chapter 5.
858 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
860 .. [2] NIST Digital Library of Mathematical Functions.
861 https://dlmf.nist.gov/10.6.E7
863 Examples
864 --------
866 Compute the Bessel function of the first kind of order 0 and
867 its first two derivatives at 1.
869 >>> from scipy.special import jvp
870 >>> jvp(0, 1, 0), jvp(0, 1, 1), jvp(0, 1, 2)
871 (0.7651976865579666, -0.44005058574493355, -0.3251471008130331)
873 Compute the first derivative of the Bessel function of the first
874 kind for several orders at 1 by providing an array for `v`.
876 >>> jvp([0, 1, 2], 1, 1)
877 array([-0.44005059, 0.3251471 , 0.21024362])
879 Compute the first derivative of the Bessel function of the first
880 kind of order 0 at several points by providing an array for `z`.
882 >>> import numpy as np
883 >>> points = np.array([0., 1.5, 3.])
884 >>> jvp(0, points, 1)
885 array([-0. , -0.55793651, -0.33905896])
887 Plot the Bessel function of the first kind of order 1 and its
888 first three derivatives.
890 >>> import matplotlib.pyplot as plt
891 >>> x = np.linspace(-10, 10, 1000)
892 >>> fig, ax = plt.subplots()
893 >>> ax.plot(x, jvp(1, x, 0), label=r"$J_1$")
894 >>> ax.plot(x, jvp(1, x, 1), label=r"$J_1'$")
895 >>> ax.plot(x, jvp(1, x, 2), label=r"$J_1''$")
896 >>> ax.plot(x, jvp(1, x, 3), label=r"$J_1'''$")
897 >>> plt.legend()
898 >>> plt.show()
899 """
900 n = _nonneg_int_or_fail(n, 'n')
901 if n == 0:
902 return jv(v, z)
903 else:
904 return _bessel_diff_formula(v, z, n, jv, -1)
907def yvp(v, z, n=1):
908 """Compute derivatives of Bessel functions of the second kind.
910 Compute the nth derivative of the Bessel function `Yv` with
911 respect to `z`.
913 Parameters
914 ----------
915 v : array_like of float
916 Order of Bessel function
917 z : complex
918 Argument at which to evaluate the derivative
919 n : int, default 1
920 Order of derivative. For 0 returns the BEssel function `yv`
922 See Also
923 --------
924 yv
926 Returns
927 -------
928 scalar or ndarray
929 nth derivative of the Bessel function.
931 Notes
932 -----
933 The derivative is computed using the relation DLFM 10.6.7 [2]_.
935 References
936 ----------
937 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
938 Functions", John Wiley and Sons, 1996, chapter 5.
939 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
941 .. [2] NIST Digital Library of Mathematical Functions.
942 https://dlmf.nist.gov/10.6.E7
944 Examples
945 --------
946 Compute the Bessel function of the second kind of order 0 and
947 its first two derivatives at 1.
949 >>> from scipy.special import yvp
950 >>> yvp(0, 1, 0), yvp(0, 1, 1), yvp(0, 1, 2)
951 (0.088256964215677, 0.7812128213002889, -0.8694697855159659)
953 Compute the first derivative of the Bessel function of the second
954 kind for several orders at 1 by providing an array for `v`.
956 >>> yvp([0, 1, 2], 1, 1)
957 array([0.78121282, 0.86946979, 2.52015239])
959 Compute the first derivative of the Bessel function of the
960 second kind of order 0 at several points by providing an array for `z`.
962 >>> import numpy as np
963 >>> points = np.array([0.5, 1.5, 3.])
964 >>> yvp(0, points, 1)
965 array([ 1.47147239, 0.41230863, -0.32467442])
967 Plot the Bessel function of the second kind of order 1 and its
968 first three derivatives.
970 >>> import matplotlib.pyplot as plt
971 >>> x = np.linspace(0, 5, 1000)
972 >>> x[0] += 1e-15
973 >>> fig, ax = plt.subplots()
974 >>> ax.plot(x, yvp(1, x, 0), label=r"$Y_1$")
975 >>> ax.plot(x, yvp(1, x, 1), label=r"$Y_1'$")
976 >>> ax.plot(x, yvp(1, x, 2), label=r"$Y_1''$")
977 >>> ax.plot(x, yvp(1, x, 3), label=r"$Y_1'''$")
978 >>> ax.set_ylim(-10, 10)
979 >>> plt.legend()
980 >>> plt.show()
981 """
982 n = _nonneg_int_or_fail(n, 'n')
983 if n == 0:
984 return yv(v, z)
985 else:
986 return _bessel_diff_formula(v, z, n, yv, -1)
989def kvp(v, z, n=1):
990 """Compute derivatives of real-order modified Bessel function Kv(z)
992 Kv(z) is the modified Bessel function of the second kind.
993 Derivative is calculated with respect to `z`.
995 Parameters
996 ----------
997 v : array_like of float
998 Order of Bessel function
999 z : array_like of complex
1000 Argument at which to evaluate the derivative
1001 n : int, default 1
1002 Order of derivative. For 0 returns the Bessel function `kv` itself.
1004 Returns
1005 -------
1006 out : ndarray
1007 The results
1009 See Also
1010 --------
1011 kv
1013 Notes
1014 -----
1015 The derivative is computed using the relation DLFM 10.29.5 [2]_.
1017 References
1018 ----------
1019 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
1020 Functions", John Wiley and Sons, 1996, chapter 6.
1021 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
1023 .. [2] NIST Digital Library of Mathematical Functions.
1024 https://dlmf.nist.gov/10.29.E5
1026 Examples
1027 --------
1028 Compute the modified bessel function of the second kind of order 0 and
1029 its first two derivatives at 1.
1031 >>> from scipy.special import kvp
1032 >>> kvp(0, 1, 0), kvp(0, 1, 1), kvp(0, 1, 2)
1033 (0.42102443824070834, -0.6019072301972346, 1.0229316684379428)
1035 Compute the first derivative of the modified Bessel function of the second
1036 kind for several orders at 1 by providing an array for `v`.
1038 >>> kvp([0, 1, 2], 1, 1)
1039 array([-0.60190723, -1.02293167, -3.85158503])
1041 Compute the first derivative of the modified Bessel function of the
1042 second kind of order 0 at several points by providing an array for `z`.
1044 >>> import numpy as np
1045 >>> points = np.array([0.5, 1.5, 3.])
1046 >>> kvp(0, points, 1)
1047 array([-1.65644112, -0.2773878 , -0.04015643])
1049 Plot the modified bessel function of the second kind and its
1050 first three derivatives.
1052 >>> import matplotlib.pyplot as plt
1053 >>> x = np.linspace(0, 5, 1000)
1054 >>> fig, ax = plt.subplots()
1055 >>> ax.plot(x, kvp(1, x, 0), label=r"$K_1$")
1056 >>> ax.plot(x, kvp(1, x, 1), label=r"$K_1'$")
1057 >>> ax.plot(x, kvp(1, x, 2), label=r"$K_1''$")
1058 >>> ax.plot(x, kvp(1, x, 3), label=r"$K_1'''$")
1059 >>> ax.set_ylim(-2.5, 2.5)
1060 >>> plt.legend()
1061 >>> plt.show()
1062 """
1063 n = _nonneg_int_or_fail(n, 'n')
1064 if n == 0:
1065 return kv(v, z)
1066 else:
1067 return (-1)**n * _bessel_diff_formula(v, z, n, kv, 1)
1070def ivp(v, z, n=1):
1071 """Compute derivatives of modified Bessel functions of the first kind.
1073 Compute the nth derivative of the modified Bessel function `Iv`
1074 with respect to `z`.
1076 Parameters
1077 ----------
1078 v : array_like or float
1079 Order of Bessel function
1080 z : array_like
1081 Argument at which to evaluate the derivative; can be real or
1082 complex.
1083 n : int, default 1
1084 Order of derivative. For 0, returns the Bessel function `iv` itself.
1086 Returns
1087 -------
1088 scalar or ndarray
1089 nth derivative of the modified Bessel function.
1091 See Also
1092 --------
1093 iv
1095 Notes
1096 -----
1097 The derivative is computed using the relation DLFM 10.29.5 [2]_.
1099 References
1100 ----------
1101 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
1102 Functions", John Wiley and Sons, 1996, chapter 6.
1103 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
1105 .. [2] NIST Digital Library of Mathematical Functions.
1106 https://dlmf.nist.gov/10.29.E5
1108 Examples
1109 --------
1110 Compute the modified Bessel function of the first kind of order 0 and
1111 its first two derivatives at 1.
1113 >>> from scipy.special import ivp
1114 >>> ivp(0, 1, 0), ivp(0, 1, 1), ivp(0, 1, 2)
1115 (1.2660658777520084, 0.565159103992485, 0.7009067737595233)
1117 Compute the first derivative of the modified Bessel function of the first
1118 kind for several orders at 1 by providing an array for `v`.
1120 >>> ivp([0, 1, 2], 1, 1)
1121 array([0.5651591 , 0.70090677, 0.29366376])
1123 Compute the first derivative of the modified Bessel function of the
1124 first kind of order 0 at several points by providing an array for `z`.
1126 >>> import numpy as np
1127 >>> points = np.array([0., 1.5, 3.])
1128 >>> ivp(0, points, 1)
1129 array([0. , 0.98166643, 3.95337022])
1131 Plot the modified Bessel function of the first kind of order 1 and its
1132 first three derivatives.
1134 >>> import matplotlib.pyplot as plt
1135 >>> x = np.linspace(-5, 5, 1000)
1136 >>> fig, ax = plt.subplots()
1137 >>> ax.plot(x, ivp(1, x, 0), label=r"$I_1$")
1138 >>> ax.plot(x, ivp(1, x, 1), label=r"$I_1'$")
1139 >>> ax.plot(x, ivp(1, x, 2), label=r"$I_1''$")
1140 >>> ax.plot(x, ivp(1, x, 3), label=r"$I_1'''$")
1141 >>> plt.legend()
1142 >>> plt.show()
1143 """
1144 n = _nonneg_int_or_fail(n, 'n')
1145 if n == 0:
1146 return iv(v, z)
1147 else:
1148 return _bessel_diff_formula(v, z, n, iv, 1)
1151def h1vp(v, z, n=1):
1152 """Compute derivatives of Hankel function H1v(z) with respect to `z`.
1154 Parameters
1155 ----------
1156 v : array_like
1157 Order of Hankel function
1158 z : array_like
1159 Argument at which to evaluate the derivative. Can be real or
1160 complex.
1161 n : int, default 1
1162 Order of derivative. For 0 returns the Hankel function `h1v` itself.
1164 Returns
1165 -------
1166 scalar or ndarray
1167 Values of the derivative of the Hankel function.
1169 See Also
1170 --------
1171 hankel1
1173 Notes
1174 -----
1175 The derivative is computed using the relation DLFM 10.6.7 [2]_.
1177 References
1178 ----------
1179 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
1180 Functions", John Wiley and Sons, 1996, chapter 5.
1181 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
1183 .. [2] NIST Digital Library of Mathematical Functions.
1184 https://dlmf.nist.gov/10.6.E7
1186 Examples
1187 --------
1188 Compute the Hankel function of the first kind of order 0 and
1189 its first two derivatives at 1.
1191 >>> from scipy.special import h1vp
1192 >>> h1vp(0, 1, 0), h1vp(0, 1, 1), h1vp(0, 1, 2)
1193 ((0.7651976865579664+0.088256964215677j),
1194 (-0.44005058574493355+0.7812128213002889j),
1195 (-0.3251471008130329-0.8694697855159659j))
1197 Compute the first derivative of the Hankel function of the first kind
1198 for several orders at 1 by providing an array for `v`.
1200 >>> h1vp([0, 1, 2], 1, 1)
1201 array([-0.44005059+0.78121282j, 0.3251471 +0.86946979j,
1202 0.21024362+2.52015239j])
1204 Compute the first derivative of the Hankel function of the first kind
1205 of order 0 at several points by providing an array for `z`.
1207 >>> import numpy as np
1208 >>> points = np.array([0.5, 1.5, 3.])
1209 >>> h1vp(0, points, 1)
1210 array([-0.24226846+1.47147239j, -0.55793651+0.41230863j,
1211 -0.33905896-0.32467442j])
1212 """
1213 n = _nonneg_int_or_fail(n, 'n')
1214 if n == 0:
1215 return hankel1(v, z)
1216 else:
1217 return _bessel_diff_formula(v, z, n, hankel1, -1)
1220def h2vp(v, z, n=1):
1221 """Compute derivatives of Hankel function H2v(z) with respect to `z`.
1223 Parameters
1224 ----------
1225 v : array_like
1226 Order of Hankel function
1227 z : array_like
1228 Argument at which to evaluate the derivative. Can be real or
1229 complex.
1230 n : int, default 1
1231 Order of derivative. For 0 returns the Hankel function `h2v` itself.
1233 Returns
1234 -------
1235 scalar or ndarray
1236 Values of the derivative of the Hankel function.
1238 See Also
1239 --------
1240 hankel2
1242 Notes
1243 -----
1244 The derivative is computed using the relation DLFM 10.6.7 [2]_.
1246 References
1247 ----------
1248 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
1249 Functions", John Wiley and Sons, 1996, chapter 5.
1250 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
1252 .. [2] NIST Digital Library of Mathematical Functions.
1253 https://dlmf.nist.gov/10.6.E7
1255 Examples
1256 --------
1257 Compute the Hankel function of the second kind of order 0 and
1258 its first two derivatives at 1.
1260 >>> from scipy.special import h2vp
1261 >>> h2vp(0, 1, 0), h2vp(0, 1, 1), h2vp(0, 1, 2)
1262 ((0.7651976865579664-0.088256964215677j),
1263 (-0.44005058574493355-0.7812128213002889j),
1264 (-0.3251471008130329+0.8694697855159659j))
1266 Compute the first derivative of the Hankel function of the second kind
1267 for several orders at 1 by providing an array for `v`.
1269 >>> h2vp([0, 1, 2], 1, 1)
1270 array([-0.44005059-0.78121282j, 0.3251471 -0.86946979j,
1271 0.21024362-2.52015239j])
1273 Compute the first derivative of the Hankel function of the second kind
1274 of order 0 at several points by providing an array for `z`.
1276 >>> import numpy as np
1277 >>> points = np.array([0.5, 1.5, 3.])
1278 >>> h2vp(0, points, 1)
1279 array([-0.24226846-1.47147239j, -0.55793651-0.41230863j,
1280 -0.33905896+0.32467442j])
1281 """
1282 n = _nonneg_int_or_fail(n, 'n')
1283 if n == 0:
1284 return hankel2(v, z)
1285 else:
1286 return _bessel_diff_formula(v, z, n, hankel2, -1)
1289def riccati_jn(n, x):
1290 r"""Compute Ricatti-Bessel function of the first kind and its derivative.
1292 The Ricatti-Bessel function of the first kind is defined as :math:`x
1293 j_n(x)`, where :math:`j_n` is the spherical Bessel function of the first
1294 kind of order :math:`n`.
1296 This function computes the value and first derivative of the
1297 Ricatti-Bessel function for all orders up to and including `n`.
1299 Parameters
1300 ----------
1301 n : int
1302 Maximum order of function to compute
1303 x : float
1304 Argument at which to evaluate
1306 Returns
1307 -------
1308 jn : ndarray
1309 Value of j0(x), ..., jn(x)
1310 jnp : ndarray
1311 First derivative j0'(x), ..., jn'(x)
1313 Notes
1314 -----
1315 The computation is carried out via backward recurrence, using the
1316 relation DLMF 10.51.1 [2]_.
1318 Wrapper for a Fortran routine created by Shanjie Zhang and Jianming
1319 Jin [1]_.
1321 References
1322 ----------
1323 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
1324 Functions", John Wiley and Sons, 1996.
1325 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
1326 .. [2] NIST Digital Library of Mathematical Functions.
1327 https://dlmf.nist.gov/10.51.E1
1329 """
1330 if not (isscalar(n) and isscalar(x)):
1331 raise ValueError("arguments must be scalars.")
1332 n = _nonneg_int_or_fail(n, 'n', strict=False)
1333 if (n == 0):
1334 n1 = 1
1335 else:
1336 n1 = n
1337 nm, jn, jnp = _specfun.rctj(n1, x)
1338 return jn[:(n+1)], jnp[:(n+1)]
1341def riccati_yn(n, x):
1342 """Compute Ricatti-Bessel function of the second kind and its derivative.
1344 The Ricatti-Bessel function of the second kind is defined as :math:`x
1345 y_n(x)`, where :math:`y_n` is the spherical Bessel function of the second
1346 kind of order :math:`n`.
1348 This function computes the value and first derivative of the function for
1349 all orders up to and including `n`.
1351 Parameters
1352 ----------
1353 n : int
1354 Maximum order of function to compute
1355 x : float
1356 Argument at which to evaluate
1358 Returns
1359 -------
1360 yn : ndarray
1361 Value of y0(x), ..., yn(x)
1362 ynp : ndarray
1363 First derivative y0'(x), ..., yn'(x)
1365 Notes
1366 -----
1367 The computation is carried out via ascending recurrence, using the
1368 relation DLMF 10.51.1 [2]_.
1370 Wrapper for a Fortran routine created by Shanjie Zhang and Jianming
1371 Jin [1]_.
1373 References
1374 ----------
1375 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
1376 Functions", John Wiley and Sons, 1996.
1377 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
1378 .. [2] NIST Digital Library of Mathematical Functions.
1379 https://dlmf.nist.gov/10.51.E1
1381 """
1382 if not (isscalar(n) and isscalar(x)):
1383 raise ValueError("arguments must be scalars.")
1384 n = _nonneg_int_or_fail(n, 'n', strict=False)
1385 if (n == 0):
1386 n1 = 1
1387 else:
1388 n1 = n
1389 nm, jn, jnp = _specfun.rcty(n1, x)
1390 return jn[:(n+1)], jnp[:(n+1)]
1393def erf_zeros(nt):
1394 """Compute the first nt zero in the first quadrant, ordered by absolute value.
1396 Zeros in the other quadrants can be obtained by using the symmetries erf(-z) = erf(z) and
1397 erf(conj(z)) = conj(erf(z)).
1400 Parameters
1401 ----------
1402 nt : int
1403 The number of zeros to compute
1405 Returns
1406 -------
1407 The locations of the zeros of erf : ndarray (complex)
1408 Complex values at which zeros of erf(z)
1410 Examples
1411 --------
1412 >>> from scipy import special
1413 >>> special.erf_zeros(1)
1414 array([1.45061616+1.880943j])
1416 Check that erf is (close to) zero for the value returned by erf_zeros
1418 >>> special.erf(special.erf_zeros(1))
1419 array([4.95159469e-14-1.16407394e-16j])
1421 References
1422 ----------
1423 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
1424 Functions", John Wiley and Sons, 1996.
1425 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
1427 """
1428 if (floor(nt) != nt) or (nt <= 0) or not isscalar(nt):
1429 raise ValueError("Argument must be positive scalar integer.")
1430 return _specfun.cerzo(nt)
1433def fresnelc_zeros(nt):
1434 """Compute nt complex zeros of cosine Fresnel integral C(z).
1436 References
1437 ----------
1438 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
1439 Functions", John Wiley and Sons, 1996.
1440 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
1442 """
1443 if (floor(nt) != nt) or (nt <= 0) or not isscalar(nt):
1444 raise ValueError("Argument must be positive scalar integer.")
1445 return _specfun.fcszo(1, nt)
1448def fresnels_zeros(nt):
1449 """Compute nt complex zeros of sine Fresnel integral S(z).
1451 References
1452 ----------
1453 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
1454 Functions", John Wiley and Sons, 1996.
1455 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
1457 """
1458 if (floor(nt) != nt) or (nt <= 0) or not isscalar(nt):
1459 raise ValueError("Argument must be positive scalar integer.")
1460 return _specfun.fcszo(2, nt)
1463def fresnel_zeros(nt):
1464 """Compute nt complex zeros of sine and cosine Fresnel integrals S(z) and C(z).
1466 References
1467 ----------
1468 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
1469 Functions", John Wiley and Sons, 1996.
1470 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
1472 """
1473 if (floor(nt) != nt) or (nt <= 0) or not isscalar(nt):
1474 raise ValueError("Argument must be positive scalar integer.")
1475 return _specfun.fcszo(2, nt), _specfun.fcszo(1, nt)
1478def assoc_laguerre(x, n, k=0.0):
1479 """Compute the generalized (associated) Laguerre polynomial of degree n and order k.
1481 The polynomial :math:`L^{(k)}_n(x)` is orthogonal over ``[0, inf)``,
1482 with weighting function ``exp(-x) * x**k`` with ``k > -1``.
1484 Notes
1485 -----
1486 `assoc_laguerre` is a simple wrapper around `eval_genlaguerre`, with
1487 reversed argument order ``(x, n, k=0.0) --> (n, k, x)``.
1489 """
1490 return _ufuncs.eval_genlaguerre(n, k, x)
1493digamma = psi
1496def polygamma(n, x):
1497 r"""Polygamma functions.
1499 Defined as :math:`\psi^{(n)}(x)` where :math:`\psi` is the
1500 `digamma` function. See [dlmf]_ for details.
1502 Parameters
1503 ----------
1504 n : array_like
1505 The order of the derivative of the digamma function; must be
1506 integral
1507 x : array_like
1508 Real valued input
1510 Returns
1511 -------
1512 ndarray
1513 Function results
1515 See Also
1516 --------
1517 digamma
1519 References
1520 ----------
1521 .. [dlmf] NIST, Digital Library of Mathematical Functions,
1522 https://dlmf.nist.gov/5.15
1524 Examples
1525 --------
1526 >>> from scipy import special
1527 >>> x = [2, 3, 25.5]
1528 >>> special.polygamma(1, x)
1529 array([ 0.64493407, 0.39493407, 0.03999467])
1530 >>> special.polygamma(0, x) == special.psi(x)
1531 array([ True, True, True], dtype=bool)
1533 """
1534 n, x = asarray(n), asarray(x)
1535 fac2 = (-1.0)**(n+1) * gamma(n+1.0) * zeta(n+1, x)
1536 return where(n == 0, psi(x), fac2)
1539def mathieu_even_coef(m, q):
1540 r"""Fourier coefficients for even Mathieu and modified Mathieu functions.
1542 The Fourier series of the even solutions of the Mathieu differential
1543 equation are of the form
1545 .. math:: \mathrm{ce}_{2n}(z, q) = \sum_{k=0}^{\infty} A_{(2n)}^{(2k)} \cos 2kz
1547 .. math:: \mathrm{ce}_{2n+1}(z, q) = \sum_{k=0}^{\infty} A_{(2n+1)}^{(2k+1)} \cos (2k+1)z
1549 This function returns the coefficients :math:`A_{(2n)}^{(2k)}` for even
1550 input m=2n, and the coefficients :math:`A_{(2n+1)}^{(2k+1)}` for odd input
1551 m=2n+1.
1553 Parameters
1554 ----------
1555 m : int
1556 Order of Mathieu functions. Must be non-negative.
1557 q : float (>=0)
1558 Parameter of Mathieu functions. Must be non-negative.
1560 Returns
1561 -------
1562 Ak : ndarray
1563 Even or odd Fourier coefficients, corresponding to even or odd m.
1565 References
1566 ----------
1567 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
1568 Functions", John Wiley and Sons, 1996.
1569 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
1570 .. [2] NIST Digital Library of Mathematical Functions
1571 https://dlmf.nist.gov/28.4#i
1573 """
1574 if not (isscalar(m) and isscalar(q)):
1575 raise ValueError("m and q must be scalars.")
1576 if (q < 0):
1577 raise ValueError("q >=0")
1578 if (m != floor(m)) or (m < 0):
1579 raise ValueError("m must be an integer >=0.")
1581 if (q <= 1):
1582 qm = 7.5 + 56.1*sqrt(q) - 134.7*q + 90.7*sqrt(q)*q
1583 else:
1584 qm = 17.0 + 3.1*sqrt(q) - .126*q + .0037*sqrt(q)*q
1585 km = int(qm + 0.5*m)
1586 if km > 251:
1587 warnings.warn("Too many predicted coefficients.", RuntimeWarning, 2)
1588 kd = 1
1589 m = int(floor(m))
1590 if m % 2:
1591 kd = 2
1593 a = mathieu_a(m, q)
1594 fc = _specfun.fcoef(kd, m, q, a)
1595 return fc[:km]
1598def mathieu_odd_coef(m, q):
1599 r"""Fourier coefficients for even Mathieu and modified Mathieu functions.
1601 The Fourier series of the odd solutions of the Mathieu differential
1602 equation are of the form
1604 .. math:: \mathrm{se}_{2n+1}(z, q) = \sum_{k=0}^{\infty} B_{(2n+1)}^{(2k+1)} \sin (2k+1)z
1606 .. math:: \mathrm{se}_{2n+2}(z, q) = \sum_{k=0}^{\infty} B_{(2n+2)}^{(2k+2)} \sin (2k+2)z
1608 This function returns the coefficients :math:`B_{(2n+2)}^{(2k+2)}` for even
1609 input m=2n+2, and the coefficients :math:`B_{(2n+1)}^{(2k+1)}` for odd
1610 input m=2n+1.
1612 Parameters
1613 ----------
1614 m : int
1615 Order of Mathieu functions. Must be non-negative.
1616 q : float (>=0)
1617 Parameter of Mathieu functions. Must be non-negative.
1619 Returns
1620 -------
1621 Bk : ndarray
1622 Even or odd Fourier coefficients, corresponding to even or odd m.
1624 References
1625 ----------
1626 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
1627 Functions", John Wiley and Sons, 1996.
1628 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
1630 """
1631 if not (isscalar(m) and isscalar(q)):
1632 raise ValueError("m and q must be scalars.")
1633 if (q < 0):
1634 raise ValueError("q >=0")
1635 if (m != floor(m)) or (m <= 0):
1636 raise ValueError("m must be an integer > 0")
1638 if (q <= 1):
1639 qm = 7.5 + 56.1*sqrt(q) - 134.7*q + 90.7*sqrt(q)*q
1640 else:
1641 qm = 17.0 + 3.1*sqrt(q) - .126*q + .0037*sqrt(q)*q
1642 km = int(qm + 0.5*m)
1643 if km > 251:
1644 warnings.warn("Too many predicted coefficients.", RuntimeWarning, 2)
1645 kd = 4
1646 m = int(floor(m))
1647 if m % 2:
1648 kd = 3
1650 b = mathieu_b(m, q)
1651 fc = _specfun.fcoef(kd, m, q, b)
1652 return fc[:km]
1655def lpmn(m, n, z):
1656 """Sequence of associated Legendre functions of the first kind.
1658 Computes the associated Legendre function of the first kind of order m and
1659 degree n, ``Pmn(z)`` = :math:`P_n^m(z)`, and its derivative, ``Pmn'(z)``.
1660 Returns two arrays of size ``(m+1, n+1)`` containing ``Pmn(z)`` and
1661 ``Pmn'(z)`` for all orders from ``0..m`` and degrees from ``0..n``.
1663 This function takes a real argument ``z``. For complex arguments ``z``
1664 use clpmn instead.
1666 Parameters
1667 ----------
1668 m : int
1669 ``|m| <= n``; the order of the Legendre function.
1670 n : int
1671 where ``n >= 0``; the degree of the Legendre function. Often
1672 called ``l`` (lower case L) in descriptions of the associated
1673 Legendre function
1674 z : float
1675 Input value.
1677 Returns
1678 -------
1679 Pmn_z : (m+1, n+1) array
1680 Values for all orders 0..m and degrees 0..n
1681 Pmn_d_z : (m+1, n+1) array
1682 Derivatives for all orders 0..m and degrees 0..n
1684 See Also
1685 --------
1686 clpmn: associated Legendre functions of the first kind for complex z
1688 Notes
1689 -----
1690 In the interval (-1, 1), Ferrer's function of the first kind is
1691 returned. The phase convention used for the intervals (1, inf)
1692 and (-inf, -1) is such that the result is always real.
1694 References
1695 ----------
1696 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
1697 Functions", John Wiley and Sons, 1996.
1698 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
1699 .. [2] NIST Digital Library of Mathematical Functions
1700 https://dlmf.nist.gov/14.3
1702 """
1703 if not isscalar(m) or (abs(m) > n):
1704 raise ValueError("m must be <= n.")
1705 if not isscalar(n) or (n < 0):
1706 raise ValueError("n must be a non-negative integer.")
1707 if not isscalar(z):
1708 raise ValueError("z must be scalar.")
1709 if iscomplex(z):
1710 raise ValueError("Argument must be real. Use clpmn instead.")
1711 if (m < 0):
1712 mp = -m
1713 mf, nf = mgrid[0:mp+1, 0:n+1]
1714 with _ufuncs.errstate(all='ignore'):
1715 if abs(z) < 1:
1716 # Ferrer function; DLMF 14.9.3
1717 fixarr = where(mf > nf, 0.0,
1718 (-1)**mf * gamma(nf-mf+1) / gamma(nf+mf+1))
1719 else:
1720 # Match to clpmn; DLMF 14.9.13
1721 fixarr = where(mf > nf, 0.0, gamma(nf-mf+1) / gamma(nf+mf+1))
1722 else:
1723 mp = m
1724 p, pd = _specfun.lpmn(mp, n, z)
1725 if (m < 0):
1726 p = p * fixarr
1727 pd = pd * fixarr
1728 return p, pd
1731def clpmn(m, n, z, type=3):
1732 """Associated Legendre function of the first kind for complex arguments.
1734 Computes the associated Legendre function of the first kind of order m and
1735 degree n, ``Pmn(z)`` = :math:`P_n^m(z)`, and its derivative, ``Pmn'(z)``.
1736 Returns two arrays of size ``(m+1, n+1)`` containing ``Pmn(z)`` and
1737 ``Pmn'(z)`` for all orders from ``0..m`` and degrees from ``0..n``.
1739 Parameters
1740 ----------
1741 m : int
1742 ``|m| <= n``; the order of the Legendre function.
1743 n : int
1744 where ``n >= 0``; the degree of the Legendre function. Often
1745 called ``l`` (lower case L) in descriptions of the associated
1746 Legendre function
1747 z : float or complex
1748 Input value.
1749 type : int, optional
1750 takes values 2 or 3
1751 2: cut on the real axis ``|x| > 1``
1752 3: cut on the real axis ``-1 < x < 1`` (default)
1754 Returns
1755 -------
1756 Pmn_z : (m+1, n+1) array
1757 Values for all orders ``0..m`` and degrees ``0..n``
1758 Pmn_d_z : (m+1, n+1) array
1759 Derivatives for all orders ``0..m`` and degrees ``0..n``
1761 See Also
1762 --------
1763 lpmn: associated Legendre functions of the first kind for real z
1765 Notes
1766 -----
1767 By default, i.e. for ``type=3``, phase conventions are chosen according
1768 to [1]_ such that the function is analytic. The cut lies on the interval
1769 (-1, 1). Approaching the cut from above or below in general yields a phase
1770 factor with respect to Ferrer's function of the first kind
1771 (cf. `lpmn`).
1773 For ``type=2`` a cut at ``|x| > 1`` is chosen. Approaching the real values
1774 on the interval (-1, 1) in the complex plane yields Ferrer's function
1775 of the first kind.
1777 References
1778 ----------
1779 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
1780 Functions", John Wiley and Sons, 1996.
1781 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
1782 .. [2] NIST Digital Library of Mathematical Functions
1783 https://dlmf.nist.gov/14.21
1785 """
1786 if not isscalar(m) or (abs(m) > n):
1787 raise ValueError("m must be <= n.")
1788 if not isscalar(n) or (n < 0):
1789 raise ValueError("n must be a non-negative integer.")
1790 if not isscalar(z):
1791 raise ValueError("z must be scalar.")
1792 if not (type == 2 or type == 3):
1793 raise ValueError("type must be either 2 or 3.")
1794 if (m < 0):
1795 mp = -m
1796 mf, nf = mgrid[0:mp+1, 0:n+1]
1797 with _ufuncs.errstate(all='ignore'):
1798 if type == 2:
1799 fixarr = where(mf > nf, 0.0,
1800 (-1)**mf * gamma(nf-mf+1) / gamma(nf+mf+1))
1801 else:
1802 fixarr = where(mf > nf, 0.0, gamma(nf-mf+1) / gamma(nf+mf+1))
1803 else:
1804 mp = m
1805 p, pd = _specfun.clpmn(mp, n, real(z), imag(z), type)
1806 if (m < 0):
1807 p = p * fixarr
1808 pd = pd * fixarr
1809 return p, pd
1812def lqmn(m, n, z):
1813 """Sequence of associated Legendre functions of the second kind.
1815 Computes the associated Legendre function of the second kind of order m and
1816 degree n, ``Qmn(z)`` = :math:`Q_n^m(z)`, and its derivative, ``Qmn'(z)``.
1817 Returns two arrays of size ``(m+1, n+1)`` containing ``Qmn(z)`` and
1818 ``Qmn'(z)`` for all orders from ``0..m`` and degrees from ``0..n``.
1820 Parameters
1821 ----------
1822 m : int
1823 ``|m| <= n``; the order of the Legendre function.
1824 n : int
1825 where ``n >= 0``; the degree of the Legendre function. Often
1826 called ``l`` (lower case L) in descriptions of the associated
1827 Legendre function
1828 z : complex
1829 Input value.
1831 Returns
1832 -------
1833 Qmn_z : (m+1, n+1) array
1834 Values for all orders 0..m and degrees 0..n
1835 Qmn_d_z : (m+1, n+1) array
1836 Derivatives for all orders 0..m and degrees 0..n
1838 References
1839 ----------
1840 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
1841 Functions", John Wiley and Sons, 1996.
1842 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
1844 """
1845 if not isscalar(m) or (m < 0):
1846 raise ValueError("m must be a non-negative integer.")
1847 if not isscalar(n) or (n < 0):
1848 raise ValueError("n must be a non-negative integer.")
1849 if not isscalar(z):
1850 raise ValueError("z must be scalar.")
1851 m = int(m)
1852 n = int(n)
1854 # Ensure neither m nor n == 0
1855 mm = max(1, m)
1856 nn = max(1, n)
1858 if iscomplex(z):
1859 q, qd = _specfun.clqmn(mm, nn, z)
1860 else:
1861 q, qd = _specfun.lqmn(mm, nn, z)
1862 return q[:(m+1), :(n+1)], qd[:(m+1), :(n+1)]
1865def bernoulli(n):
1866 """Bernoulli numbers B0..Bn (inclusive).
1868 Parameters
1869 ----------
1870 n : int
1871 Indicated the number of terms in the Bernoulli series to generate.
1873 Returns
1874 -------
1875 ndarray
1876 The Bernoulli numbers ``[B(0), B(1), ..., B(n)]``.
1878 References
1879 ----------
1880 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
1881 Functions", John Wiley and Sons, 1996.
1882 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
1883 .. [2] "Bernoulli number", Wikipedia, https://en.wikipedia.org/wiki/Bernoulli_number
1885 Examples
1886 --------
1887 >>> import numpy as np
1888 >>> from scipy.special import bernoulli, zeta
1889 >>> bernoulli(4)
1890 array([ 1. , -0.5 , 0.16666667, 0. , -0.03333333])
1892 The Wikipedia article ([2]_) points out the relationship between the
1893 Bernoulli numbers and the zeta function, ``B_n^+ = -n * zeta(1 - n)``
1894 for ``n > 0``:
1896 >>> n = np.arange(1, 5)
1897 >>> -n * zeta(1 - n)
1898 array([ 0.5 , 0.16666667, -0. , -0.03333333])
1900 Note that, in the notation used in the wikipedia article,
1901 `bernoulli` computes ``B_n^-`` (i.e. it used the convention that
1902 ``B_1`` is -1/2). The relation given above is for ``B_n^+``, so the
1903 sign of 0.5 does not match the output of ``bernoulli(4)``.
1905 """
1906 if not isscalar(n) or (n < 0):
1907 raise ValueError("n must be a non-negative integer.")
1908 n = int(n)
1909 if (n < 2):
1910 n1 = 2
1911 else:
1912 n1 = n
1913 return _specfun.bernob(int(n1))[:(n+1)]
1916def euler(n):
1917 """Euler numbers E(0), E(1), ..., E(n).
1919 The Euler numbers [1]_ are also known as the secant numbers.
1921 Because ``euler(n)`` returns floating point values, it does not give
1922 exact values for large `n`. The first inexact value is E(22).
1924 Parameters
1925 ----------
1926 n : int
1927 The highest index of the Euler number to be returned.
1929 Returns
1930 -------
1931 ndarray
1932 The Euler numbers [E(0), E(1), ..., E(n)].
1933 The odd Euler numbers, which are all zero, are included.
1935 References
1936 ----------
1937 .. [1] Sequence A122045, The On-Line Encyclopedia of Integer Sequences,
1938 https://oeis.org/A122045
1939 .. [2] Zhang, Shanjie and Jin, Jianming. "Computation of Special
1940 Functions", John Wiley and Sons, 1996.
1941 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
1943 Examples
1944 --------
1945 >>> import numpy as np
1946 >>> from scipy.special import euler
1947 >>> euler(6)
1948 array([ 1., 0., -1., 0., 5., 0., -61.])
1950 >>> euler(13).astype(np.int64)
1951 array([ 1, 0, -1, 0, 5, 0, -61,
1952 0, 1385, 0, -50521, 0, 2702765, 0])
1954 >>> euler(22)[-1] # Exact value of E(22) is -69348874393137901.
1955 -69348874393137976.0
1957 """
1958 if not isscalar(n) or (n < 0):
1959 raise ValueError("n must be a non-negative integer.")
1960 n = int(n)
1961 if (n < 2):
1962 n1 = 2
1963 else:
1964 n1 = n
1965 return _specfun.eulerb(n1)[:(n+1)]
1968def lpn(n, z):
1969 """Legendre function of the first kind.
1971 Compute sequence of Legendre functions of the first kind (polynomials),
1972 Pn(z) and derivatives for all degrees from 0 to n (inclusive).
1974 See also special.legendre for polynomial class.
1976 References
1977 ----------
1978 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
1979 Functions", John Wiley and Sons, 1996.
1980 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
1982 """
1983 if not (isscalar(n) and isscalar(z)):
1984 raise ValueError("arguments must be scalars.")
1985 n = _nonneg_int_or_fail(n, 'n', strict=False)
1986 if (n < 1):
1987 n1 = 1
1988 else:
1989 n1 = n
1990 if iscomplex(z):
1991 pn, pd = _specfun.clpn(n1, z)
1992 else:
1993 pn, pd = _specfun.lpn(n1, z)
1994 return pn[:(n+1)], pd[:(n+1)]
1997def lqn(n, z):
1998 """Legendre function of the second kind.
2000 Compute sequence of Legendre functions of the second kind, Qn(z) and
2001 derivatives for all degrees from 0 to n (inclusive).
2003 References
2004 ----------
2005 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
2006 Functions", John Wiley and Sons, 1996.
2007 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
2009 """
2010 if not (isscalar(n) and isscalar(z)):
2011 raise ValueError("arguments must be scalars.")
2012 n = _nonneg_int_or_fail(n, 'n', strict=False)
2013 if (n < 1):
2014 n1 = 1
2015 else:
2016 n1 = n
2017 if iscomplex(z):
2018 qn, qd = _specfun.clqn(n1, z)
2019 else:
2020 qn, qd = _specfun.lqnb(n1, z)
2021 return qn[:(n+1)], qd[:(n+1)]
2024def ai_zeros(nt):
2025 """
2026 Compute `nt` zeros and values of the Airy function Ai and its derivative.
2028 Computes the first `nt` zeros, `a`, of the Airy function Ai(x);
2029 first `nt` zeros, `ap`, of the derivative of the Airy function Ai'(x);
2030 the corresponding values Ai(a');
2031 and the corresponding values Ai'(a).
2033 Parameters
2034 ----------
2035 nt : int
2036 Number of zeros to compute
2038 Returns
2039 -------
2040 a : ndarray
2041 First `nt` zeros of Ai(x)
2042 ap : ndarray
2043 First `nt` zeros of Ai'(x)
2044 ai : ndarray
2045 Values of Ai(x) evaluated at first `nt` zeros of Ai'(x)
2046 aip : ndarray
2047 Values of Ai'(x) evaluated at first `nt` zeros of Ai(x)
2049 Examples
2050 --------
2051 >>> from scipy import special
2052 >>> a, ap, ai, aip = special.ai_zeros(3)
2053 >>> a
2054 array([-2.33810741, -4.08794944, -5.52055983])
2055 >>> ap
2056 array([-1.01879297, -3.24819758, -4.82009921])
2057 >>> ai
2058 array([ 0.53565666, -0.41901548, 0.38040647])
2059 >>> aip
2060 array([ 0.70121082, -0.80311137, 0.86520403])
2062 References
2063 ----------
2064 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
2065 Functions", John Wiley and Sons, 1996.
2066 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
2068 """
2069 kf = 1
2070 if not isscalar(nt) or (floor(nt) != nt) or (nt <= 0):
2071 raise ValueError("nt must be a positive integer scalar.")
2072 return _specfun.airyzo(nt, kf)
2075def bi_zeros(nt):
2076 """
2077 Compute `nt` zeros and values of the Airy function Bi and its derivative.
2079 Computes the first `nt` zeros, b, of the Airy function Bi(x);
2080 first `nt` zeros, b', of the derivative of the Airy function Bi'(x);
2081 the corresponding values Bi(b');
2082 and the corresponding values Bi'(b).
2084 Parameters
2085 ----------
2086 nt : int
2087 Number of zeros to compute
2089 Returns
2090 -------
2091 b : ndarray
2092 First `nt` zeros of Bi(x)
2093 bp : ndarray
2094 First `nt` zeros of Bi'(x)
2095 bi : ndarray
2096 Values of Bi(x) evaluated at first `nt` zeros of Bi'(x)
2097 bip : ndarray
2098 Values of Bi'(x) evaluated at first `nt` zeros of Bi(x)
2100 Examples
2101 --------
2102 >>> from scipy import special
2103 >>> b, bp, bi, bip = special.bi_zeros(3)
2104 >>> b
2105 array([-1.17371322, -3.2710933 , -4.83073784])
2106 >>> bp
2107 array([-2.29443968, -4.07315509, -5.51239573])
2108 >>> bi
2109 array([-0.45494438, 0.39652284, -0.36796916])
2110 >>> bip
2111 array([ 0.60195789, -0.76031014, 0.83699101])
2113 References
2114 ----------
2115 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
2116 Functions", John Wiley and Sons, 1996.
2117 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
2119 """
2120 kf = 2
2121 if not isscalar(nt) or (floor(nt) != nt) or (nt <= 0):
2122 raise ValueError("nt must be a positive integer scalar.")
2123 return _specfun.airyzo(nt, kf)
2126def lmbda(v, x):
2127 r"""Jahnke-Emden Lambda function, Lambdav(x).
2129 This function is defined as [2]_,
2131 .. math:: \Lambda_v(x) = \Gamma(v+1) \frac{J_v(x)}{(x/2)^v},
2133 where :math:`\Gamma` is the gamma function and :math:`J_v` is the
2134 Bessel function of the first kind.
2136 Parameters
2137 ----------
2138 v : float
2139 Order of the Lambda function
2140 x : float
2141 Value at which to evaluate the function and derivatives
2143 Returns
2144 -------
2145 vl : ndarray
2146 Values of Lambda_vi(x), for vi=v-int(v), vi=1+v-int(v), ..., vi=v.
2147 dl : ndarray
2148 Derivatives Lambda_vi'(x), for vi=v-int(v), vi=1+v-int(v), ..., vi=v.
2150 References
2151 ----------
2152 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
2153 Functions", John Wiley and Sons, 1996.
2154 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
2155 .. [2] Jahnke, E. and Emde, F. "Tables of Functions with Formulae and
2156 Curves" (4th ed.), Dover, 1945
2157 """
2158 if not (isscalar(v) and isscalar(x)):
2159 raise ValueError("arguments must be scalars.")
2160 if (v < 0):
2161 raise ValueError("argument must be > 0.")
2162 n = int(v)
2163 v0 = v - n
2164 if (n < 1):
2165 n1 = 1
2166 else:
2167 n1 = n
2168 v1 = n1 + v0
2169 if (v != floor(v)):
2170 vm, vl, dl = _specfun.lamv(v1, x)
2171 else:
2172 vm, vl, dl = _specfun.lamn(v1, x)
2173 return vl[:(n+1)], dl[:(n+1)]
2176def pbdv_seq(v, x):
2177 """Parabolic cylinder functions Dv(x) and derivatives.
2179 Parameters
2180 ----------
2181 v : float
2182 Order of the parabolic cylinder function
2183 x : float
2184 Value at which to evaluate the function and derivatives
2186 Returns
2187 -------
2188 dv : ndarray
2189 Values of D_vi(x), for vi=v-int(v), vi=1+v-int(v), ..., vi=v.
2190 dp : ndarray
2191 Derivatives D_vi'(x), for vi=v-int(v), vi=1+v-int(v), ..., vi=v.
2193 References
2194 ----------
2195 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
2196 Functions", John Wiley and Sons, 1996, chapter 13.
2197 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
2199 """
2200 if not (isscalar(v) and isscalar(x)):
2201 raise ValueError("arguments must be scalars.")
2202 n = int(v)
2203 v0 = v-n
2204 if (n < 1):
2205 n1 = 1
2206 else:
2207 n1 = n
2208 v1 = n1 + v0
2209 dv, dp, pdf, pdd = _specfun.pbdv(v1, x)
2210 return dv[:n1+1], dp[:n1+1]
2213def pbvv_seq(v, x):
2214 """Parabolic cylinder functions Vv(x) and derivatives.
2216 Parameters
2217 ----------
2218 v : float
2219 Order of the parabolic cylinder function
2220 x : float
2221 Value at which to evaluate the function and derivatives
2223 Returns
2224 -------
2225 dv : ndarray
2226 Values of V_vi(x), for vi=v-int(v), vi=1+v-int(v), ..., vi=v.
2227 dp : ndarray
2228 Derivatives V_vi'(x), for vi=v-int(v), vi=1+v-int(v), ..., vi=v.
2230 References
2231 ----------
2232 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
2233 Functions", John Wiley and Sons, 1996, chapter 13.
2234 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
2236 """
2237 if not (isscalar(v) and isscalar(x)):
2238 raise ValueError("arguments must be scalars.")
2239 n = int(v)
2240 v0 = v-n
2241 if (n <= 1):
2242 n1 = 1
2243 else:
2244 n1 = n
2245 v1 = n1 + v0
2246 dv, dp, pdf, pdd = _specfun.pbvv(v1, x)
2247 return dv[:n1+1], dp[:n1+1]
2250def pbdn_seq(n, z):
2251 """Parabolic cylinder functions Dn(z) and derivatives.
2253 Parameters
2254 ----------
2255 n : int
2256 Order of the parabolic cylinder function
2257 z : complex
2258 Value at which to evaluate the function and derivatives
2260 Returns
2261 -------
2262 dv : ndarray
2263 Values of D_i(z), for i=0, ..., i=n.
2264 dp : ndarray
2265 Derivatives D_i'(z), for i=0, ..., i=n.
2267 References
2268 ----------
2269 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
2270 Functions", John Wiley and Sons, 1996, chapter 13.
2271 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
2273 """
2274 if not (isscalar(n) and isscalar(z)):
2275 raise ValueError("arguments must be scalars.")
2276 if (floor(n) != n):
2277 raise ValueError("n must be an integer.")
2278 if (abs(n) <= 1):
2279 n1 = 1
2280 else:
2281 n1 = n
2282 cpb, cpd = _specfun.cpbdn(n1, z)
2283 return cpb[:n1+1], cpd[:n1+1]
2286def ber_zeros(nt):
2287 """Compute nt zeros of the Kelvin function ber.
2289 Parameters
2290 ----------
2291 nt : int
2292 Number of zeros to compute. Must be positive.
2294 Returns
2295 -------
2296 ndarray
2297 First `nt` zeros of the Kelvin function.
2299 See Also
2300 --------
2301 ber
2303 References
2304 ----------
2305 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
2306 Functions", John Wiley and Sons, 1996.
2307 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
2309 """
2310 if not isscalar(nt) or (floor(nt) != nt) or (nt <= 0):
2311 raise ValueError("nt must be positive integer scalar.")
2312 return _specfun.klvnzo(nt, 1)
2315def bei_zeros(nt):
2316 """Compute nt zeros of the Kelvin function bei.
2318 Parameters
2319 ----------
2320 nt : int
2321 Number of zeros to compute. Must be positive.
2323 Returns
2324 -------
2325 ndarray
2326 First `nt` zeros of the Kelvin function.
2328 See Also
2329 --------
2330 bei
2332 References
2333 ----------
2334 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
2335 Functions", John Wiley and Sons, 1996.
2336 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
2338 """
2339 if not isscalar(nt) or (floor(nt) != nt) or (nt <= 0):
2340 raise ValueError("nt must be positive integer scalar.")
2341 return _specfun.klvnzo(nt, 2)
2344def ker_zeros(nt):
2345 """Compute nt zeros of the Kelvin function ker.
2347 Parameters
2348 ----------
2349 nt : int
2350 Number of zeros to compute. Must be positive.
2352 Returns
2353 -------
2354 ndarray
2355 First `nt` zeros of the Kelvin function.
2357 See Also
2358 --------
2359 ker
2361 References
2362 ----------
2363 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
2364 Functions", John Wiley and Sons, 1996.
2365 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
2367 """
2368 if not isscalar(nt) or (floor(nt) != nt) or (nt <= 0):
2369 raise ValueError("nt must be positive integer scalar.")
2370 return _specfun.klvnzo(nt, 3)
2373def kei_zeros(nt):
2374 """Compute nt zeros of the Kelvin function kei.
2376 Parameters
2377 ----------
2378 nt : int
2379 Number of zeros to compute. Must be positive.
2381 Returns
2382 -------
2383 ndarray
2384 First `nt` zeros of the Kelvin function.
2386 See Also
2387 --------
2388 kei
2390 References
2391 ----------
2392 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
2393 Functions", John Wiley and Sons, 1996.
2394 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
2396 """
2397 if not isscalar(nt) or (floor(nt) != nt) or (nt <= 0):
2398 raise ValueError("nt must be positive integer scalar.")
2399 return _specfun.klvnzo(nt, 4)
2402def berp_zeros(nt):
2403 """Compute nt zeros of the derivative of the Kelvin function ber.
2405 Parameters
2406 ----------
2407 nt : int
2408 Number of zeros to compute. Must be positive.
2410 Returns
2411 -------
2412 ndarray
2413 First `nt` zeros of the derivative of the Kelvin function.
2415 See Also
2416 --------
2417 ber, berp
2419 References
2420 ----------
2421 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
2422 Functions", John Wiley and Sons, 1996.
2423 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
2425 """
2426 if not isscalar(nt) or (floor(nt) != nt) or (nt <= 0):
2427 raise ValueError("nt must be positive integer scalar.")
2428 return _specfun.klvnzo(nt, 5)
2431def beip_zeros(nt):
2432 """Compute nt zeros of the derivative of the Kelvin function bei.
2434 Parameters
2435 ----------
2436 nt : int
2437 Number of zeros to compute. Must be positive.
2439 Returns
2440 -------
2441 ndarray
2442 First `nt` zeros of the derivative of the Kelvin function.
2444 See Also
2445 --------
2446 bei, beip
2448 References
2449 ----------
2450 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
2451 Functions", John Wiley and Sons, 1996.
2452 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
2454 """
2455 if not isscalar(nt) or (floor(nt) != nt) or (nt <= 0):
2456 raise ValueError("nt must be positive integer scalar.")
2457 return _specfun.klvnzo(nt, 6)
2460def kerp_zeros(nt):
2461 """Compute nt zeros of the derivative of the Kelvin function ker.
2463 Parameters
2464 ----------
2465 nt : int
2466 Number of zeros to compute. Must be positive.
2468 Returns
2469 -------
2470 ndarray
2471 First `nt` zeros of the derivative of the Kelvin function.
2473 See Also
2474 --------
2475 ker, kerp
2477 References
2478 ----------
2479 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
2480 Functions", John Wiley and Sons, 1996.
2481 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
2483 """
2484 if not isscalar(nt) or (floor(nt) != nt) or (nt <= 0):
2485 raise ValueError("nt must be positive integer scalar.")
2486 return _specfun.klvnzo(nt, 7)
2489def keip_zeros(nt):
2490 """Compute nt zeros of the derivative of the Kelvin function kei.
2492 Parameters
2493 ----------
2494 nt : int
2495 Number of zeros to compute. Must be positive.
2497 Returns
2498 -------
2499 ndarray
2500 First `nt` zeros of the derivative of the Kelvin function.
2502 See Also
2503 --------
2504 kei, keip
2506 References
2507 ----------
2508 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
2509 Functions", John Wiley and Sons, 1996.
2510 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
2512 """
2513 if not isscalar(nt) or (floor(nt) != nt) or (nt <= 0):
2514 raise ValueError("nt must be positive integer scalar.")
2515 return _specfun.klvnzo(nt, 8)
2518def kelvin_zeros(nt):
2519 """Compute nt zeros of all Kelvin functions.
2521 Returned in a length-8 tuple of arrays of length nt. The tuple contains
2522 the arrays of zeros of (ber, bei, ker, kei, ber', bei', ker', kei').
2524 References
2525 ----------
2526 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
2527 Functions", John Wiley and Sons, 1996.
2528 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
2530 """
2531 if not isscalar(nt) or (floor(nt) != nt) or (nt <= 0):
2532 raise ValueError("nt must be positive integer scalar.")
2533 return (_specfun.klvnzo(nt, 1),
2534 _specfun.klvnzo(nt, 2),
2535 _specfun.klvnzo(nt, 3),
2536 _specfun.klvnzo(nt, 4),
2537 _specfun.klvnzo(nt, 5),
2538 _specfun.klvnzo(nt, 6),
2539 _specfun.klvnzo(nt, 7),
2540 _specfun.klvnzo(nt, 8))
2543def pro_cv_seq(m, n, c):
2544 """Characteristic values for prolate spheroidal wave functions.
2546 Compute a sequence of characteristic values for the prolate
2547 spheroidal wave functions for mode m and n'=m..n and spheroidal
2548 parameter c.
2550 References
2551 ----------
2552 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
2553 Functions", John Wiley and Sons, 1996.
2554 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
2556 """
2557 if not (isscalar(m) and isscalar(n) and isscalar(c)):
2558 raise ValueError("Arguments must be scalars.")
2559 if (n != floor(n)) or (m != floor(m)):
2560 raise ValueError("Modes must be integers.")
2561 if (n-m > 199):
2562 raise ValueError("Difference between n and m is too large.")
2563 maxL = n-m+1
2564 return _specfun.segv(m, n, c, 1)[1][:maxL]
2567def obl_cv_seq(m, n, c):
2568 """Characteristic values for oblate spheroidal wave functions.
2570 Compute a sequence of characteristic values for the oblate
2571 spheroidal wave functions for mode m and n'=m..n and spheroidal
2572 parameter c.
2574 References
2575 ----------
2576 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
2577 Functions", John Wiley and Sons, 1996.
2578 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
2580 """
2581 if not (isscalar(m) and isscalar(n) and isscalar(c)):
2582 raise ValueError("Arguments must be scalars.")
2583 if (n != floor(n)) or (m != floor(m)):
2584 raise ValueError("Modes must be integers.")
2585 if (n-m > 199):
2586 raise ValueError("Difference between n and m is too large.")
2587 maxL = n-m+1
2588 return _specfun.segv(m, n, c, -1)[1][:maxL]
2591def comb(N, k, exact=False, repetition=False, legacy=None):
2592 """The number of combinations of N things taken k at a time.
2594 This is often expressed as "N choose k".
2596 Parameters
2597 ----------
2598 N : int, ndarray
2599 Number of things.
2600 k : int, ndarray
2601 Number of elements taken.
2602 exact : bool, optional
2603 For integers, if `exact` is False, then floating point precision is
2604 used, otherwise the result is computed exactly. For non-integers, if
2605 `exact` is True, is disregarded.
2606 repetition : bool, optional
2607 If `repetition` is True, then the number of combinations with
2608 repetition is computed.
2609 legacy : bool, optional
2610 If `legacy` is True and `exact` is True, then non-integral arguments
2611 are cast to ints; if `legacy` is False, the result for non-integral
2612 arguments is unaffected by the value of `exact`.
2614 .. deprecated:: 1.9.0
2615 Using `legacy` is deprecated and will removed by
2616 Scipy 1.13.0. If you want to keep the legacy behaviour, cast
2617 your inputs directly, e.g.
2618 ``comb(int(your_N), int(your_k), exact=True)``.
2620 Returns
2621 -------
2622 val : int, float, ndarray
2623 The total number of combinations.
2625 See Also
2626 --------
2627 binom : Binomial coefficient considered as a function of two real
2628 variables.
2630 Notes
2631 -----
2632 - Array arguments accepted only for exact=False case.
2633 - If N < 0, or k < 0, then 0 is returned.
2634 - If k > N and repetition=False, then 0 is returned.
2636 Examples
2637 --------
2638 >>> import numpy as np
2639 >>> from scipy.special import comb
2640 >>> k = np.array([3, 4])
2641 >>> n = np.array([10, 10])
2642 >>> comb(n, k, exact=False)
2643 array([ 120., 210.])
2644 >>> comb(10, 3, exact=True)
2645 120
2646 >>> comb(10, 3, exact=True, repetition=True)
2647 220
2649 """
2650 if legacy is not None:
2651 warnings.warn(
2652 "Using 'legacy' keyword is deprecated and will be removed by "
2653 "Scipy 1.13.0. If you want to keep the legacy behaviour, cast "
2654 "your inputs directly, e.g. "
2655 "'comb(int(your_N), int(your_k), exact=True)'.",
2656 DeprecationWarning,
2657 stacklevel=2
2658 )
2659 if repetition:
2660 return comb(N + k - 1, k, exact, legacy=legacy)
2661 if exact:
2662 if int(N) == N and int(k) == k:
2663 # _comb_int casts inputs to integers, which is safe & intended here
2664 return _comb_int(N, k)
2665 elif legacy:
2666 # here at least one number is not an integer; legacy behavior uses
2667 # lossy casts to int
2668 return _comb_int(N, k)
2669 # otherwise, we disregard `exact=True`; it makes no sense for
2670 # non-integral arguments
2671 return comb(N, k)
2672 else:
2673 k, N = asarray(k), asarray(N)
2674 cond = (k <= N) & (N >= 0) & (k >= 0)
2675 vals = binom(N, k)
2676 if isinstance(vals, np.ndarray):
2677 vals[~cond] = 0
2678 elif not cond:
2679 vals = np.float64(0)
2680 return vals
2683def perm(N, k, exact=False):
2684 """Permutations of N things taken k at a time, i.e., k-permutations of N.
2686 It's also known as "partial permutations".
2688 Parameters
2689 ----------
2690 N : int, ndarray
2691 Number of things.
2692 k : int, ndarray
2693 Number of elements taken.
2694 exact : bool, optional
2695 If `exact` is False, then floating point precision is used, otherwise
2696 exact long integer is computed.
2698 Returns
2699 -------
2700 val : int, ndarray
2701 The number of k-permutations of N.
2703 Notes
2704 -----
2705 - Array arguments accepted only for exact=False case.
2706 - If k > N, N < 0, or k < 0, then a 0 is returned.
2708 Examples
2709 --------
2710 >>> import numpy as np
2711 >>> from scipy.special import perm
2712 >>> k = np.array([3, 4])
2713 >>> n = np.array([10, 10])
2714 >>> perm(n, k)
2715 array([ 720., 5040.])
2716 >>> perm(10, 3, exact=True)
2717 720
2719 """
2720 if exact:
2721 if (k > N) or (N < 0) or (k < 0):
2722 return 0
2723 val = 1
2724 for i in range(N - k + 1, N + 1):
2725 val *= i
2726 return val
2727 else:
2728 k, N = asarray(k), asarray(N)
2729 cond = (k <= N) & (N >= 0) & (k >= 0)
2730 vals = poch(N - k + 1, k)
2731 if isinstance(vals, np.ndarray):
2732 vals[~cond] = 0
2733 elif not cond:
2734 vals = np.float64(0)
2735 return vals
2738# https://stackoverflow.com/a/16327037
2739def _range_prod(lo, hi, k=1):
2740 """
2741 Product of a range of numbers spaced k apart (from hi).
2743 For k=1, this returns the product of
2744 lo * (lo+1) * (lo+2) * ... * (hi-2) * (hi-1) * hi
2745 = hi! / (lo-1)!
2747 For k>1, it correspond to taking only every k'th number when
2748 counting down from hi - e.g. 18!!!! = _range_prod(1, 18, 4).
2750 Breaks into smaller products first for speed:
2751 _range_prod(2, 9) = ((2*3)*(4*5))*((6*7)*(8*9))
2752 """
2753 if lo + k < hi:
2754 mid = (hi + lo) // 2
2755 if k > 1:
2756 # make sure mid is a multiple of k away from hi
2757 mid = mid - ((mid - hi) % k)
2758 return _range_prod(lo, mid, k) * _range_prod(mid + k, hi, k)
2759 elif lo + k == hi:
2760 return lo * hi
2761 else:
2762 return hi
2765def _exact_factorialx_array(n, k=1):
2766 """
2767 Exact computation of factorial for an array.
2769 The factorials are computed in incremental fashion, by taking
2770 the sorted unique values of n and multiplying the intervening
2771 numbers between the different unique values.
2773 In other words, the factorial for the largest input is only
2774 computed once, with each other result computed in the process.
2776 k > 1 corresponds to the multifactorial.
2777 """
2778 un = np.unique(n)
2779 # numpy changed nan-sorting behaviour with 1.21, see numpy/numpy#18070;
2780 # to unify the behaviour, we remove the nan's here; the respective
2781 # values will be set separately at the end
2782 un = un[~np.isnan(un)]
2784 # Convert to object array if np.int64 can't handle size
2785 if np.isnan(n).any():
2786 dt = float
2787 elif k in _FACTORIALK_LIMITS_64BITS.keys():
2788 if un[-1] > _FACTORIALK_LIMITS_64BITS[k]:
2789 # e.g. k=1: 21! > np.iinfo(np.int64).max
2790 dt = object
2791 elif un[-1] > _FACTORIALK_LIMITS_32BITS[k]:
2792 # e.g. k=3: 26!!! > np.iinfo(np.int32).max
2793 dt = np.int64
2794 else:
2795 dt = np.int_
2796 else:
2797 # for k >= 10, we always use object
2798 dt = object
2800 out = np.empty_like(n, dtype=dt)
2802 # Handle invalid/trivial values
2803 un = un[un > 1]
2804 out[n < 2] = 1
2805 out[n < 0] = 0
2807 # Calculate products of each range of numbers
2808 # we can only multiply incrementally if the values are k apart;
2809 # therefore we partition `un` into "lanes", i.e. its residues modulo k
2810 for lane in range(0, k):
2811 ul = un[(un % k) == lane] if k > 1 else un
2812 if ul.size:
2813 # after np.unique, un resp. ul are sorted, ul[0] is the smallest;
2814 # cast to python ints to avoid overflow with np.int-types
2815 val = _range_prod(1, int(ul[0]), k=k)
2816 out[n == ul[0]] = val
2817 for i in range(len(ul) - 1):
2818 # by the filtering above, we have ensured that prev & current
2819 # are a multiple of k apart
2820 prev = ul[i]
2821 current = ul[i + 1]
2822 # we already multiplied all factors until prev; continue
2823 # building the full factorial from the following (`prev + 1`);
2824 # use int() for the same reason as above
2825 val *= _range_prod(int(prev + 1), int(current), k=k)
2826 out[n == current] = val
2828 if np.isnan(n).any():
2829 out = out.astype(np.float64)
2830 out[np.isnan(n)] = np.nan
2831 return out
2834def factorial(n, exact=False):
2835 """
2836 The factorial of a number or array of numbers.
2838 The factorial of non-negative integer `n` is the product of all
2839 positive integers less than or equal to `n`::
2841 n! = n * (n - 1) * (n - 2) * ... * 1
2843 Parameters
2844 ----------
2845 n : int or array_like of ints
2846 Input values. If ``n < 0``, the return value is 0.
2847 exact : bool, optional
2848 If True, calculate the answer exactly using long integer arithmetic.
2849 If False, result is approximated in floating point rapidly using the
2850 `gamma` function.
2851 Default is False.
2853 Returns
2854 -------
2855 nf : float or int or ndarray
2856 Factorial of `n`, as integer or float depending on `exact`.
2858 Notes
2859 -----
2860 For arrays with ``exact=True``, the factorial is computed only once, for
2861 the largest input, with each other result computed in the process.
2862 The output dtype is increased to ``int64`` or ``object`` if necessary.
2864 With ``exact=False`` the factorial is approximated using the gamma
2865 function:
2867 .. math:: n! = \\Gamma(n+1)
2869 Examples
2870 --------
2871 >>> import numpy as np
2872 >>> from scipy.special import factorial
2873 >>> arr = np.array([3, 4, 5])
2874 >>> factorial(arr, exact=False)
2875 array([ 6., 24., 120.])
2876 >>> factorial(arr, exact=True)
2877 array([ 6, 24, 120])
2878 >>> factorial(5, exact=True)
2879 120
2881 """
2882 # don't use isscalar due to numpy/numpy#23574; 0-dim arrays treated below
2883 if np.ndim(n) == 0 and not isinstance(n, np.ndarray):
2884 # scalar cases
2885 if n is None or np.isnan(n):
2886 return np.nan
2887 elif not (np.issubdtype(type(n), np.integer)
2888 or np.issubdtype(type(n), np.floating)):
2889 raise ValueError(
2890 f"Unsupported datatype for factorial: {type(n)}\n"
2891 "Permitted data types are integers and floating point numbers"
2892 )
2893 elif n < 0:
2894 return 0
2895 elif exact and np.issubdtype(type(n), np.integer):
2896 return math.factorial(n)
2897 # we do not raise for non-integers with exact=True due to
2898 # historical reasons, though deprecation would be possible
2899 return _ufuncs._factorial(n)
2901 # arrays & array-likes
2902 n = asarray(n)
2903 if n.size == 0:
2904 # return empty arrays unchanged
2905 return n
2906 if not (np.issubdtype(n.dtype, np.integer)
2907 or np.issubdtype(n.dtype, np.floating)):
2908 raise ValueError(
2909 f"Unsupported datatype for factorial: {n.dtype}\n"
2910 "Permitted data types are integers and floating point numbers"
2911 )
2912 if exact and not np.issubdtype(n.dtype, np.integer):
2913 # legacy behaviour is to support mixed integers/NaNs;
2914 # deprecate this for exact=True
2915 n_flt = n[~np.isnan(n)]
2916 if np.allclose(n_flt, n_flt.astype(np.int64)):
2917 warnings.warn(
2918 "Non-integer arrays (e.g. due to presence of NaNs) "
2919 "together with exact=True are deprecated. Either ensure "
2920 "that the the array has integer dtype or use exact=False.",
2921 DeprecationWarning,
2922 stacklevel=2
2923 )
2924 else:
2925 msg = ("factorial with exact=True does not "
2926 "support non-integral arrays")
2927 raise ValueError(msg)
2929 if exact:
2930 return _exact_factorialx_array(n)
2931 # we do not raise for non-integers with exact=True due to
2932 # historical reasons, though deprecation would be possible
2933 res = _ufuncs._factorial(n)
2934 if isinstance(n, np.ndarray):
2935 # _ufuncs._factorial does not maintain 0-dim arrays
2936 return np.array(res)
2937 return res
2940def factorial2(n, exact=False):
2941 """Double factorial.
2943 This is the factorial with every second value skipped. E.g., ``7!! = 7 * 5
2944 * 3 * 1``. It can be approximated numerically as::
2946 n!! = 2 ** (n / 2) * gamma(n / 2 + 1) * sqrt(2 / pi) n odd
2947 = 2 ** (n / 2) * gamma(n / 2 + 1) n even
2948 = 2 ** (n / 2) * (n / 2)! n even
2950 Parameters
2951 ----------
2952 n : int or array_like
2953 Calculate ``n!!``. If ``n < 0``, the return value is 0.
2954 exact : bool, optional
2955 The result can be approximated rapidly using the gamma-formula
2956 above (default). If `exact` is set to True, calculate the
2957 answer exactly using integer arithmetic.
2959 Returns
2960 -------
2961 nff : float or int
2962 Double factorial of `n`, as an int or a float depending on
2963 `exact`.
2965 Examples
2966 --------
2967 >>> from scipy.special import factorial2
2968 >>> factorial2(7, exact=False)
2969 array(105.00000000000001)
2970 >>> factorial2(7, exact=True)
2971 105
2973 """
2974 def _approx(n):
2975 # main factor that both even/odd approximations share
2976 val = np.power(2, n / 2) * gamma(n / 2 + 1)
2977 mask = np.ones_like(n, dtype=np.float64)
2978 mask[n % 2 == 1] = sqrt(2 / pi)
2979 # analytical continuation (based on odd integers)
2980 # is scaled down by a factor of sqrt(2 / pi)
2981 # compared to the value of even integers.
2982 return val * mask
2984 # don't use isscalar due to numpy/numpy#23574; 0-dim arrays treated below
2985 if np.ndim(n) == 0 and not isinstance(n, np.ndarray):
2986 # scalar cases
2987 if n is None or np.isnan(n):
2988 return np.nan
2989 elif not np.issubdtype(type(n), np.integer):
2990 msg = "factorial2 does not support non-integral scalar arguments"
2991 raise ValueError(msg)
2992 elif n < 0:
2993 return 0
2994 elif n in {0, 1}:
2995 return 1
2996 # general integer case
2997 if exact:
2998 return _range_prod(1, n, k=2)
2999 return _approx(n)
3000 # arrays & array-likes
3001 n = asarray(n)
3002 if n.size == 0:
3003 # return empty arrays unchanged
3004 return n
3005 if not np.issubdtype(n.dtype, np.integer):
3006 raise ValueError("factorial2 does not support non-integral arrays")
3007 if exact:
3008 return _exact_factorialx_array(n, k=2)
3009 # approximation
3010 vals = zeros(n.shape)
3011 cond = (n >= 0)
3012 n_to_compute = extract(cond, n)
3013 place(vals, cond, _approx(n_to_compute))
3014 return vals
3017def factorialk(n, k, exact=True):
3018 """Multifactorial of n of order k, n(!!...!).
3020 This is the multifactorial of n skipping k values. For example,
3022 factorialk(17, 4) = 17!!!! = 17 * 13 * 9 * 5 * 1
3024 In particular, for any integer ``n``, we have
3026 factorialk(n, 1) = factorial(n)
3028 factorialk(n, 2) = factorial2(n)
3030 Parameters
3031 ----------
3032 n : int or array_like
3033 Calculate multifactorial. If `n` < 0, the return value is 0.
3034 k : int
3035 Order of multifactorial.
3036 exact : bool, optional
3037 If exact is set to True, calculate the answer exactly using
3038 integer arithmetic.
3040 Returns
3041 -------
3042 val : int
3043 Multifactorial of `n`.
3045 Raises
3046 ------
3047 NotImplementedError
3048 Raises when exact is False
3050 Examples
3051 --------
3052 >>> from scipy.special import factorialk
3053 >>> factorialk(5, 1, exact=True)
3054 120
3055 >>> factorialk(5, 3, exact=True)
3056 10
3058 """
3059 if not np.issubdtype(type(k), np.integer) or k < 1:
3060 raise ValueError(f"k must be a positive integer, received: {k}")
3061 if not exact:
3062 raise NotImplementedError
3064 helpmsg = ""
3065 if k in {1, 2}:
3066 func = "factorial" if k == 1 else "factorial2"
3067 helpmsg = f"\nYou can try to use {func} instead"
3069 # don't use isscalar due to numpy/numpy#23574; 0-dim arrays treated below
3070 if np.ndim(n) == 0 and not isinstance(n, np.ndarray):
3071 # scalar cases
3072 if n is None or np.isnan(n):
3073 return np.nan
3074 elif not np.issubdtype(type(n), np.integer):
3075 msg = "factorialk does not support non-integral scalar arguments!"
3076 raise ValueError(msg + helpmsg)
3077 elif n < 0:
3078 return 0
3079 elif n in {0, 1}:
3080 return 1
3081 return _range_prod(1, n, k=k)
3082 # arrays & array-likes
3083 n = asarray(n)
3084 if n.size == 0:
3085 # return empty arrays unchanged
3086 return n
3087 if not np.issubdtype(n.dtype, np.integer):
3088 msg = "factorialk does not support non-integral arrays!"
3089 raise ValueError(msg + helpmsg)
3090 return _exact_factorialx_array(n, k=k)
3093def zeta(x, q=None, out=None):
3094 r"""
3095 Riemann or Hurwitz zeta function.
3097 Parameters
3098 ----------
3099 x : array_like of float
3100 Input data, must be real
3101 q : array_like of float, optional
3102 Input data, must be real. Defaults to Riemann zeta.
3103 out : ndarray, optional
3104 Output array for the computed values.
3106 Returns
3107 -------
3108 out : array_like
3109 Values of zeta(x).
3111 Notes
3112 -----
3113 The two-argument version is the Hurwitz zeta function
3115 .. math::
3117 \zeta(x, q) = \sum_{k=0}^{\infty} \frac{1}{(k + q)^x};
3119 see [dlmf]_ for details. The Riemann zeta function corresponds to
3120 the case when ``q = 1``.
3122 See Also
3123 --------
3124 zetac
3126 References
3127 ----------
3128 .. [dlmf] NIST, Digital Library of Mathematical Functions,
3129 https://dlmf.nist.gov/25.11#i
3131 Examples
3132 --------
3133 >>> import numpy as np
3134 >>> from scipy.special import zeta, polygamma, factorial
3136 Some specific values:
3138 >>> zeta(2), np.pi**2/6
3139 (1.6449340668482266, 1.6449340668482264)
3141 >>> zeta(4), np.pi**4/90
3142 (1.0823232337111381, 1.082323233711138)
3144 Relation to the `polygamma` function:
3146 >>> m = 3
3147 >>> x = 1.25
3148 >>> polygamma(m, x)
3149 array(2.782144009188397)
3150 >>> (-1)**(m+1) * factorial(m) * zeta(m+1, x)
3151 2.7821440091883969
3153 """
3154 if q is None:
3155 return _ufuncs._riemann_zeta(x, out)
3156 else:
3157 return _ufuncs._zeta(x, q, out)