Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/special/_basic.py: 13%
527 statements
« prev ^ index » next coverage.py v7.3.2, created at 2023-12-12 06:31 +0000
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]
81def _nonneg_int_or_fail(n, var_name, strict=True):
82 try:
83 if strict:
84 # Raises an exception if float
85 n = operator.index(n)
86 elif n == floor(n):
87 n = int(n)
88 else:
89 raise ValueError()
90 if n < 0:
91 raise ValueError()
92 except (ValueError, TypeError) as err:
93 raise err.__class__("{} must be a non-negative integer".format(var_name)) from err
94 return n
97def diric(x, n):
98 """Periodic sinc function, also called the Dirichlet function.
100 The Dirichlet function is defined as::
102 diric(x, n) = sin(x * n/2) / (n * sin(x / 2)),
104 where `n` is a positive integer.
106 Parameters
107 ----------
108 x : array_like
109 Input data
110 n : int
111 Integer defining the periodicity.
113 Returns
114 -------
115 diric : ndarray
117 Examples
118 --------
119 >>> import numpy as np
120 >>> from scipy import special
121 >>> import matplotlib.pyplot as plt
123 >>> x = np.linspace(-8*np.pi, 8*np.pi, num=201)
124 >>> plt.figure(figsize=(8, 8));
125 >>> for idx, n in enumerate([2, 3, 4, 9]):
126 ... plt.subplot(2, 2, idx+1)
127 ... plt.plot(x, special.diric(x, n))
128 ... plt.title('diric, n={}'.format(n))
129 >>> plt.show()
131 The following example demonstrates that `diric` gives the magnitudes
132 (modulo the sign and scaling) of the Fourier coefficients of a
133 rectangular pulse.
135 Suppress output of values that are effectively 0:
137 >>> np.set_printoptions(suppress=True)
139 Create a signal `x` of length `m` with `k` ones:
141 >>> m = 8
142 >>> k = 3
143 >>> x = np.zeros(m)
144 >>> x[:k] = 1
146 Use the FFT to compute the Fourier transform of `x`, and
147 inspect the magnitudes of the coefficients:
149 >>> np.abs(np.fft.fft(x))
150 array([ 3. , 2.41421356, 1. , 0.41421356, 1. ,
151 0.41421356, 1. , 2.41421356])
153 Now find the same values (up to sign) using `diric`. We multiply
154 by `k` to account for the different scaling conventions of
155 `numpy.fft.fft` and `diric`:
157 >>> theta = np.linspace(0, 2*np.pi, m, endpoint=False)
158 >>> k * special.diric(theta, k)
159 array([ 3. , 2.41421356, 1. , -0.41421356, -1. ,
160 -0.41421356, 1. , 2.41421356])
161 """
162 x, n = asarray(x), asarray(n)
163 n = asarray(n + (x-x))
164 x = asarray(x + (n-n))
165 if issubdtype(x.dtype, inexact):
166 ytype = x.dtype
167 else:
168 ytype = float
169 y = zeros(x.shape, ytype)
171 # empirical minval for 32, 64 or 128 bit float computations
172 # where sin(x/2) < minval, result is fixed at +1 or -1
173 if np.finfo(ytype).eps < 1e-18:
174 minval = 1e-11
175 elif np.finfo(ytype).eps < 1e-15:
176 minval = 1e-7
177 else:
178 minval = 1e-3
180 mask1 = (n <= 0) | (n != floor(n))
181 place(y, mask1, nan)
183 x = x / 2
184 denom = sin(x)
185 mask2 = (1-mask1) & (abs(denom) < minval)
186 xsub = extract(mask2, x)
187 nsub = extract(mask2, n)
188 zsub = xsub / pi
189 place(y, mask2, pow(-1, np.round(zsub)*(nsub-1)))
191 mask = (1-mask1) & (1-mask2)
192 xsub = extract(mask, x)
193 nsub = extract(mask, n)
194 dsub = extract(mask, denom)
195 place(y, mask, sin(nsub*xsub)/(nsub*dsub))
196 return y
199def jnjnp_zeros(nt):
200 """Compute zeros of integer-order Bessel functions Jn and Jn'.
202 Results are arranged in order of the magnitudes of the zeros.
204 Parameters
205 ----------
206 nt : int
207 Number (<=1200) of zeros to compute
209 Returns
210 -------
211 zo[l-1] : ndarray
212 Value of the lth zero of Jn(x) and Jn'(x). Of length `nt`.
213 n[l-1] : ndarray
214 Order of the Jn(x) or Jn'(x) associated with lth zero. Of length `nt`.
215 m[l-1] : ndarray
216 Serial number of the zeros of Jn(x) or Jn'(x) associated
217 with lth zero. Of length `nt`.
218 t[l-1] : ndarray
219 0 if lth zero in zo is zero of Jn(x), 1 if it is a zero of Jn'(x). Of
220 length `nt`.
222 See Also
223 --------
224 jn_zeros, jnp_zeros : to get separated arrays of zeros.
226 References
227 ----------
228 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
229 Functions", John Wiley and Sons, 1996, chapter 5.
230 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
232 """
233 if not isscalar(nt) or (floor(nt) != nt) or (nt > 1200):
234 raise ValueError("Number must be integer <= 1200.")
235 nt = int(nt)
236 n, m, t, zo = _specfun.jdzo(nt)
237 return zo[1:nt+1], n[:nt], m[:nt], t[:nt]
240def jnyn_zeros(n, nt):
241 """Compute nt zeros of Bessel functions Jn(x), Jn'(x), Yn(x), and Yn'(x).
243 Returns 4 arrays of length `nt`, corresponding to the first `nt`
244 zeros of Jn(x), Jn'(x), Yn(x), and Yn'(x), respectively. The zeros
245 are returned in ascending order.
247 Parameters
248 ----------
249 n : int
250 Order of the Bessel functions
251 nt : int
252 Number (<=1200) of zeros to compute
254 Returns
255 -------
256 Jn : ndarray
257 First `nt` zeros of Jn
258 Jnp : ndarray
259 First `nt` zeros of Jn'
260 Yn : ndarray
261 First `nt` zeros of Yn
262 Ynp : ndarray
263 First `nt` zeros of Yn'
265 See Also
266 --------
267 jn_zeros, jnp_zeros, yn_zeros, ynp_zeros
269 References
270 ----------
271 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
272 Functions", John Wiley and Sons, 1996, chapter 5.
273 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
275 Examples
276 --------
277 Compute the first three roots of :math:`J_1`, :math:`J_1'`,
278 :math:`Y_1` and :math:`Y_1'`.
280 >>> from scipy.special import jnyn_zeros
281 >>> jn_roots, jnp_roots, yn_roots, ynp_roots = jnyn_zeros(1, 3)
282 >>> jn_roots, yn_roots
283 (array([ 3.83170597, 7.01558667, 10.17346814]),
284 array([2.19714133, 5.42968104, 8.59600587]))
286 Plot :math:`J_1`, :math:`J_1'`, :math:`Y_1`, :math:`Y_1'` and their roots.
288 >>> import numpy as np
289 >>> import matplotlib.pyplot as plt
290 >>> from scipy.special import jnyn_zeros, jvp, jn, yvp, yn
291 >>> jn_roots, jnp_roots, yn_roots, ynp_roots = jnyn_zeros(1, 3)
292 >>> fig, ax = plt.subplots()
293 >>> xmax= 11
294 >>> x = np.linspace(0, xmax)
295 >>> ax.plot(x, jn(1, x), label=r"$J_1$", c='r')
296 >>> ax.plot(x, jvp(1, x, 1), label=r"$J_1'$", c='b')
297 >>> ax.plot(x, yn(1, x), label=r"$Y_1$", c='y')
298 >>> ax.plot(x, yvp(1, x, 1), label=r"$Y_1'$", c='c')
299 >>> zeros = np.zeros((3, ))
300 >>> ax.scatter(jn_roots, zeros, s=30, c='r', zorder=5,
301 ... label=r"$J_1$ roots")
302 >>> ax.scatter(jnp_roots, zeros, s=30, c='b', zorder=5,
303 ... label=r"$J_1'$ roots")
304 >>> ax.scatter(yn_roots, zeros, s=30, c='y', zorder=5,
305 ... label=r"$Y_1$ roots")
306 >>> ax.scatter(ynp_roots, zeros, s=30, c='c', zorder=5,
307 ... label=r"$Y_1'$ roots")
308 >>> ax.hlines(0, 0, xmax, color='k')
309 >>> ax.set_ylim(-0.6, 0.6)
310 >>> ax.set_xlim(0, xmax)
311 >>> ax.legend(ncol=2, bbox_to_anchor=(1., 0.75))
312 >>> plt.tight_layout()
313 >>> plt.show()
314 """
315 if not (isscalar(nt) and isscalar(n)):
316 raise ValueError("Arguments must be scalars.")
317 if (floor(n) != n) or (floor(nt) != nt):
318 raise ValueError("Arguments must be integers.")
319 if (nt <= 0):
320 raise ValueError("nt > 0")
321 return _specfun.jyzo(abs(n), nt)
324def jn_zeros(n, nt):
325 r"""Compute zeros of integer-order Bessel functions Jn.
327 Compute `nt` zeros of the Bessel functions :math:`J_n(x)` on the
328 interval :math:`(0, \infty)`. The zeros are returned in ascending
329 order. Note that this interval excludes the zero at :math:`x = 0`
330 that exists for :math:`n > 0`.
332 Parameters
333 ----------
334 n : int
335 Order of Bessel function
336 nt : int
337 Number of zeros to return
339 Returns
340 -------
341 ndarray
342 First `nt` zeros of the Bessel function.
344 See Also
345 --------
346 jv: Real-order Bessel functions of the first kind
347 jnp_zeros: Zeros of :math:`Jn'`
349 References
350 ----------
351 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
352 Functions", John Wiley and Sons, 1996, chapter 5.
353 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
355 Examples
356 --------
357 Compute the first four positive roots of :math:`J_3`.
359 >>> from scipy.special import jn_zeros
360 >>> jn_zeros(3, 4)
361 array([ 6.3801619 , 9.76102313, 13.01520072, 16.22346616])
363 Plot :math:`J_3` and its first four positive roots. Note
364 that the root located at 0 is not returned by `jn_zeros`.
366 >>> import numpy as np
367 >>> import matplotlib.pyplot as plt
368 >>> from scipy.special import jn, jn_zeros
369 >>> j3_roots = jn_zeros(3, 4)
370 >>> xmax = 18
371 >>> xmin = -1
372 >>> x = np.linspace(xmin, xmax, 500)
373 >>> fig, ax = plt.subplots()
374 >>> ax.plot(x, jn(3, x), label=r'$J_3$')
375 >>> ax.scatter(j3_roots, np.zeros((4, )), s=30, c='r',
376 ... label=r"$J_3$_Zeros", zorder=5)
377 >>> ax.scatter(0, 0, s=30, c='k',
378 ... label=r"Root at 0", zorder=5)
379 >>> ax.hlines(0, 0, xmax, color='k')
380 >>> ax.set_xlim(xmin, xmax)
381 >>> plt.legend()
382 >>> plt.show()
383 """
384 return jnyn_zeros(n, nt)[0]
387def jnp_zeros(n, nt):
388 r"""Compute zeros of integer-order Bessel function derivatives Jn'.
390 Compute `nt` zeros of the functions :math:`J_n'(x)` on the
391 interval :math:`(0, \infty)`. The zeros are returned in ascending
392 order. Note that this interval excludes the zero at :math:`x = 0`
393 that exists for :math:`n > 1`.
395 Parameters
396 ----------
397 n : int
398 Order of Bessel function
399 nt : int
400 Number of zeros to return
402 Returns
403 -------
404 ndarray
405 First `nt` zeros of the Bessel function.
407 See Also
408 --------
409 jvp: Derivatives of integer-order Bessel functions of the first kind
410 jv: Float-order Bessel functions of the first kind
412 References
413 ----------
414 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
415 Functions", John Wiley and Sons, 1996, chapter 5.
416 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
418 Examples
419 --------
420 Compute the first four roots of :math:`J_2'`.
422 >>> from scipy.special import jnp_zeros
423 >>> jnp_zeros(2, 4)
424 array([ 3.05423693, 6.70613319, 9.96946782, 13.17037086])
426 As `jnp_zeros` yields the roots of :math:`J_n'`, it can be used to
427 compute the locations of the peaks of :math:`J_n`. Plot
428 :math:`J_2`, :math:`J_2'` and the locations of the roots of :math:`J_2'`.
430 >>> import numpy as np
431 >>> import matplotlib.pyplot as plt
432 >>> from scipy.special import jn, jnp_zeros, jvp
433 >>> j2_roots = jnp_zeros(2, 4)
434 >>> xmax = 15
435 >>> x = np.linspace(0, xmax, 500)
436 >>> fig, ax = plt.subplots()
437 >>> ax.plot(x, jn(2, x), label=r'$J_2$')
438 >>> ax.plot(x, jvp(2, x, 1), label=r"$J_2'$")
439 >>> ax.hlines(0, 0, xmax, color='k')
440 >>> ax.scatter(j2_roots, np.zeros((4, )), s=30, c='r',
441 ... label=r"Roots of $J_2'$", zorder=5)
442 >>> ax.set_ylim(-0.4, 0.8)
443 >>> ax.set_xlim(0, xmax)
444 >>> plt.legend()
445 >>> plt.show()
446 """
447 return jnyn_zeros(n, nt)[1]
450def yn_zeros(n, nt):
451 r"""Compute zeros of integer-order Bessel function Yn(x).
453 Compute `nt` zeros of the functions :math:`Y_n(x)` on the interval
454 :math:`(0, \infty)`. The zeros are returned in ascending order.
456 Parameters
457 ----------
458 n : int
459 Order of Bessel function
460 nt : int
461 Number of zeros to return
463 Returns
464 -------
465 ndarray
466 First `nt` zeros of the Bessel function.
468 See Also
469 --------
470 yn: Bessel function of the second kind for integer order
471 yv: Bessel function of the second kind for real order
473 References
474 ----------
475 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
476 Functions", John Wiley and Sons, 1996, chapter 5.
477 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
479 Examples
480 --------
481 Compute the first four roots of :math:`Y_2`.
483 >>> from scipy.special import yn_zeros
484 >>> yn_zeros(2, 4)
485 array([ 3.38424177, 6.79380751, 10.02347798, 13.20998671])
487 Plot :math:`Y_2` and its first four roots.
489 >>> import numpy as np
490 >>> import matplotlib.pyplot as plt
491 >>> from scipy.special import yn, yn_zeros
492 >>> xmin = 2
493 >>> xmax = 15
494 >>> x = np.linspace(xmin, xmax, 500)
495 >>> fig, ax = plt.subplots()
496 >>> ax.hlines(0, xmin, xmax, color='k')
497 >>> ax.plot(x, yn(2, x), label=r'$Y_2$')
498 >>> ax.scatter(yn_zeros(2, 4), np.zeros((4, )), s=30, c='r',
499 ... label='Roots', zorder=5)
500 >>> ax.set_ylim(-0.4, 0.4)
501 >>> ax.set_xlim(xmin, xmax)
502 >>> plt.legend()
503 >>> plt.show()
504 """
505 return jnyn_zeros(n, nt)[2]
508def ynp_zeros(n, nt):
509 r"""Compute zeros of integer-order Bessel function derivatives Yn'(x).
511 Compute `nt` zeros of the functions :math:`Y_n'(x)` on the
512 interval :math:`(0, \infty)`. The zeros are returned in ascending
513 order.
515 Parameters
516 ----------
517 n : int
518 Order of Bessel function
519 nt : int
520 Number of zeros to return
522 Returns
523 -------
524 ndarray
525 First `nt` zeros of the Bessel derivative function.
528 See Also
529 --------
530 yvp
532 References
533 ----------
534 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
535 Functions", John Wiley and Sons, 1996, chapter 5.
536 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
538 Examples
539 --------
540 Compute the first four roots of the first derivative of the
541 Bessel function of second kind for order 0 :math:`Y_0'`.
543 >>> from scipy.special import ynp_zeros
544 >>> ynp_zeros(0, 4)
545 array([ 2.19714133, 5.42968104, 8.59600587, 11.74915483])
547 Plot :math:`Y_0`, :math:`Y_0'` and confirm visually that the roots of
548 :math:`Y_0'` are located at local extrema of :math:`Y_0`.
550 >>> import numpy as np
551 >>> import matplotlib.pyplot as plt
552 >>> from scipy.special import yn, ynp_zeros, yvp
553 >>> zeros = ynp_zeros(0, 4)
554 >>> xmax = 13
555 >>> x = np.linspace(0, xmax, 500)
556 >>> fig, ax = plt.subplots()
557 >>> ax.plot(x, yn(0, x), label=r'$Y_0$')
558 >>> ax.plot(x, yvp(0, x, 1), label=r"$Y_0'$")
559 >>> ax.scatter(zeros, np.zeros((4, )), s=30, c='r',
560 ... label=r"Roots of $Y_0'$", zorder=5)
561 >>> for root in zeros:
562 ... y0_extremum = yn(0, root)
563 ... lower = min(0, y0_extremum)
564 ... upper = max(0, y0_extremum)
565 ... ax.vlines(root, lower, upper, color='r')
566 >>> ax.hlines(0, 0, xmax, color='k')
567 >>> ax.set_ylim(-0.6, 0.6)
568 >>> ax.set_xlim(0, xmax)
569 >>> plt.legend()
570 >>> plt.show()
571 """
572 return jnyn_zeros(n, nt)[3]
575def y0_zeros(nt, complex=False):
576 """Compute nt zeros of Bessel function Y0(z), and derivative at each zero.
578 The derivatives are given by Y0'(z0) = -Y1(z0) at each zero z0.
580 Parameters
581 ----------
582 nt : int
583 Number of zeros to return
584 complex : bool, default False
585 Set to False to return only the real zeros; set to True to return only
586 the complex zeros with negative real part and positive imaginary part.
587 Note that the complex conjugates of the latter are also zeros of the
588 function, but are not returned by this routine.
590 Returns
591 -------
592 z0n : ndarray
593 Location of nth zero of Y0(z)
594 y0pz0n : ndarray
595 Value of derivative Y0'(z0) for nth zero
597 References
598 ----------
599 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
600 Functions", John Wiley and Sons, 1996, chapter 5.
601 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
603 Examples
604 --------
605 Compute the first 4 real roots and the derivatives at the roots of
606 :math:`Y_0`:
608 >>> import numpy as np
609 >>> from scipy.special import y0_zeros
610 >>> zeros, grads = y0_zeros(4)
611 >>> with np.printoptions(precision=5):
612 ... print(f"Roots: {zeros}")
613 ... print(f"Gradients: {grads}")
614 Roots: [ 0.89358+0.j 3.95768+0.j 7.08605+0.j 10.22235+0.j]
615 Gradients: [-0.87942+0.j 0.40254+0.j -0.3001 +0.j 0.2497 +0.j]
617 Plot the real part of :math:`Y_0` and the first four computed roots.
619 >>> import matplotlib.pyplot as plt
620 >>> from scipy.special import y0
621 >>> xmin = 0
622 >>> xmax = 11
623 >>> x = np.linspace(xmin, xmax, 500)
624 >>> fig, ax = plt.subplots()
625 >>> ax.hlines(0, xmin, xmax, color='k')
626 >>> ax.plot(x, y0(x), label=r'$Y_0$')
627 >>> zeros, grads = y0_zeros(4)
628 >>> ax.scatter(zeros.real, np.zeros((4, )), s=30, c='r',
629 ... label=r'$Y_0$_zeros', zorder=5)
630 >>> ax.set_ylim(-0.5, 0.6)
631 >>> ax.set_xlim(xmin, xmax)
632 >>> plt.legend(ncol=2)
633 >>> plt.show()
635 Compute the first 4 complex roots and the derivatives at the roots of
636 :math:`Y_0` by setting ``complex=True``:
638 >>> y0_zeros(4, True)
639 (array([ -2.40301663+0.53988231j, -5.5198767 +0.54718001j,
640 -8.6536724 +0.54841207j, -11.79151203+0.54881912j]),
641 array([ 0.10074769-0.88196771j, -0.02924642+0.5871695j ,
642 0.01490806-0.46945875j, -0.00937368+0.40230454j]))
643 """
644 if not isscalar(nt) or (floor(nt) != nt) or (nt <= 0):
645 raise ValueError("Arguments must be scalar positive integer.")
646 kf = 0
647 kc = not complex
648 return _specfun.cyzo(nt, kf, kc)
651def y1_zeros(nt, complex=False):
652 """Compute nt zeros of Bessel function Y1(z), and derivative at each zero.
654 The derivatives are given by Y1'(z1) = Y0(z1) at each zero z1.
656 Parameters
657 ----------
658 nt : int
659 Number of zeros to return
660 complex : bool, default False
661 Set to False to return only the real zeros; set to True to return only
662 the complex zeros with negative real part and positive imaginary part.
663 Note that the complex conjugates of the latter are also zeros of the
664 function, but are not returned by this routine.
666 Returns
667 -------
668 z1n : ndarray
669 Location of nth zero of Y1(z)
670 y1pz1n : ndarray
671 Value of derivative Y1'(z1) for nth zero
673 References
674 ----------
675 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
676 Functions", John Wiley and Sons, 1996, chapter 5.
677 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
679 Examples
680 --------
681 Compute the first 4 real roots and the derivatives at the roots of
682 :math:`Y_1`:
684 >>> import numpy as np
685 >>> from scipy.special import y1_zeros
686 >>> zeros, grads = y1_zeros(4)
687 >>> with np.printoptions(precision=5):
688 ... print(f"Roots: {zeros}")
689 ... print(f"Gradients: {grads}")
690 Roots: [ 2.19714+0.j 5.42968+0.j 8.59601+0.j 11.74915+0.j]
691 Gradients: [ 0.52079+0.j -0.34032+0.j 0.27146+0.j -0.23246+0.j]
693 Extract the real parts:
695 >>> realzeros = zeros.real
696 >>> realzeros
697 array([ 2.19714133, 5.42968104, 8.59600587, 11.74915483])
699 Plot :math:`Y_1` and the first four computed roots.
701 >>> import matplotlib.pyplot as plt
702 >>> from scipy.special import y1
703 >>> xmin = 0
704 >>> xmax = 13
705 >>> x = np.linspace(xmin, xmax, 500)
706 >>> zeros, grads = y1_zeros(4)
707 >>> fig, ax = plt.subplots()
708 >>> ax.hlines(0, xmin, xmax, color='k')
709 >>> ax.plot(x, y1(x), label=r'$Y_1$')
710 >>> ax.scatter(zeros.real, np.zeros((4, )), s=30, c='r',
711 ... label=r'$Y_1$_zeros', zorder=5)
712 >>> ax.set_ylim(-0.5, 0.5)
713 >>> ax.set_xlim(xmin, xmax)
714 >>> plt.legend()
715 >>> plt.show()
717 Compute the first 4 complex roots and the derivatives at the roots of
718 :math:`Y_1` by setting ``complex=True``:
720 >>> y1_zeros(4, True)
721 (array([ -0.50274327+0.78624371j, -3.83353519+0.56235654j,
722 -7.01590368+0.55339305j, -10.17357383+0.55127339j]),
723 array([-0.45952768+1.31710194j, 0.04830191-0.69251288j,
724 -0.02012695+0.51864253j, 0.011614 -0.43203296j]))
725 """
726 if not isscalar(nt) or (floor(nt) != nt) or (nt <= 0):
727 raise ValueError("Arguments must be scalar positive integer.")
728 kf = 1
729 kc = not complex
730 return _specfun.cyzo(nt, kf, kc)
733def y1p_zeros(nt, complex=False):
734 """Compute nt zeros of Bessel derivative Y1'(z), and value at each zero.
736 The values are given by Y1(z1) at each z1 where Y1'(z1)=0.
738 Parameters
739 ----------
740 nt : int
741 Number of zeros to return
742 complex : bool, default False
743 Set to False to return only the real zeros; set to True to return only
744 the complex zeros with negative real part and positive imaginary part.
745 Note that the complex conjugates of the latter are also zeros of the
746 function, but are not returned by this routine.
748 Returns
749 -------
750 z1pn : ndarray
751 Location of nth zero of Y1'(z)
752 y1z1pn : ndarray
753 Value of derivative Y1(z1) for nth zero
755 References
756 ----------
757 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
758 Functions", John Wiley and Sons, 1996, chapter 5.
759 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
761 Examples
762 --------
763 Compute the first four roots of :math:`Y_1'` and the values of
764 :math:`Y_1` at these roots.
766 >>> import numpy as np
767 >>> from scipy.special import y1p_zeros
768 >>> y1grad_roots, y1_values = y1p_zeros(4)
769 >>> with np.printoptions(precision=5):
770 ... print(f"Y1' Roots: {y1grad_roots}")
771 ... print(f"Y1 values: {y1_values}")
772 Y1' Roots: [ 3.68302+0.j 6.9415 +0.j 10.1234 +0.j 13.28576+0.j]
773 Y1 values: [ 0.41673+0.j -0.30317+0.j 0.25091+0.j -0.21897+0.j]
775 `y1p_zeros` can be used to calculate the extremal points of :math:`Y_1`
776 directly. Here we plot :math:`Y_1` and the first four extrema.
778 >>> import matplotlib.pyplot as plt
779 >>> from scipy.special import y1, yvp
780 >>> y1_roots, y1_values_at_roots = y1p_zeros(4)
781 >>> real_roots = y1_roots.real
782 >>> xmax = 15
783 >>> x = np.linspace(0, xmax, 500)
784 >>> fig, ax = plt.subplots()
785 >>> ax.plot(x, y1(x), label=r'$Y_1$')
786 >>> ax.plot(x, yvp(1, x, 1), label=r"$Y_1'$")
787 >>> ax.scatter(real_roots, np.zeros((4, )), s=30, c='r',
788 ... label=r"Roots of $Y_1'$", zorder=5)
789 >>> ax.scatter(real_roots, y1_values_at_roots.real, s=30, c='k',
790 ... label=r"Extrema of $Y_1$", zorder=5)
791 >>> ax.hlines(0, 0, xmax, color='k')
792 >>> ax.set_ylim(-0.5, 0.5)
793 >>> ax.set_xlim(0, xmax)
794 >>> ax.legend(ncol=2, bbox_to_anchor=(1., 0.75))
795 >>> plt.tight_layout()
796 >>> plt.show()
797 """
798 if not isscalar(nt) or (floor(nt) != nt) or (nt <= 0):
799 raise ValueError("Arguments must be scalar positive integer.")
800 kf = 2
801 kc = not complex
802 return _specfun.cyzo(nt, kf, kc)
805def _bessel_diff_formula(v, z, n, L, phase):
806 # from AMS55.
807 # L(v, z) = J(v, z), Y(v, z), H1(v, z), H2(v, z), phase = -1
808 # L(v, z) = I(v, z) or exp(v*pi*i)K(v, z), phase = 1
809 # For K, you can pull out the exp((v-k)*pi*i) into the caller
810 v = asarray(v)
811 p = 1.0
812 s = L(v-n, z)
813 for i in range(1, n+1):
814 p = phase * (p * (n-i+1)) / i # = choose(k, i)
815 s += p*L(v-n + i*2, z)
816 return s / (2.**n)
819def jvp(v, z, n=1):
820 """Compute derivatives of Bessel functions of the first kind.
822 Compute the nth derivative of the Bessel function `Jv` with
823 respect to `z`.
825 Parameters
826 ----------
827 v : array_like or float
828 Order of Bessel function
829 z : complex
830 Argument at which to evaluate the derivative; can be real or
831 complex.
832 n : int, default 1
833 Order of derivative. For 0 returns the Bessel function `jv` itself.
835 Returns
836 -------
837 scalar or ndarray
838 Values of the derivative of the Bessel function.
840 Notes
841 -----
842 The derivative is computed using the relation DLFM 10.6.7 [2]_.
844 References
845 ----------
846 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
847 Functions", John Wiley and Sons, 1996, chapter 5.
848 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
850 .. [2] NIST Digital Library of Mathematical Functions.
851 https://dlmf.nist.gov/10.6.E7
853 Examples
854 --------
856 Compute the Bessel function of the first kind of order 0 and
857 its first two derivatives at 1.
859 >>> from scipy.special import jvp
860 >>> jvp(0, 1, 0), jvp(0, 1, 1), jvp(0, 1, 2)
861 (0.7651976865579666, -0.44005058574493355, -0.3251471008130331)
863 Compute the first derivative of the Bessel function of the first
864 kind for several orders at 1 by providing an array for `v`.
866 >>> jvp([0, 1, 2], 1, 1)
867 array([-0.44005059, 0.3251471 , 0.21024362])
869 Compute the first derivative of the Bessel function of the first
870 kind of order 0 at several points by providing an array for `z`.
872 >>> import numpy as np
873 >>> points = np.array([0., 1.5, 3.])
874 >>> jvp(0, points, 1)
875 array([-0. , -0.55793651, -0.33905896])
877 Plot the Bessel function of the first kind of order 1 and its
878 first three derivatives.
880 >>> import matplotlib.pyplot as plt
881 >>> x = np.linspace(-10, 10, 1000)
882 >>> fig, ax = plt.subplots()
883 >>> ax.plot(x, jvp(1, x, 0), label=r"$J_1$")
884 >>> ax.plot(x, jvp(1, x, 1), label=r"$J_1'$")
885 >>> ax.plot(x, jvp(1, x, 2), label=r"$J_1''$")
886 >>> ax.plot(x, jvp(1, x, 3), label=r"$J_1'''$")
887 >>> plt.legend()
888 >>> plt.show()
889 """
890 n = _nonneg_int_or_fail(n, 'n')
891 if n == 0:
892 return jv(v, z)
893 else:
894 return _bessel_diff_formula(v, z, n, jv, -1)
897def yvp(v, z, n=1):
898 """Compute derivatives of Bessel functions of the second kind.
900 Compute the nth derivative of the Bessel function `Yv` with
901 respect to `z`.
903 Parameters
904 ----------
905 v : array_like of float
906 Order of Bessel function
907 z : complex
908 Argument at which to evaluate the derivative
909 n : int, default 1
910 Order of derivative. For 0 returns the BEssel function `yv`
912 See Also
913 --------
914 yv
916 Returns
917 -------
918 scalar or ndarray
919 nth derivative of the Bessel function.
921 Notes
922 -----
923 The derivative is computed using the relation DLFM 10.6.7 [2]_.
925 References
926 ----------
927 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
928 Functions", John Wiley and Sons, 1996, chapter 5.
929 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
931 .. [2] NIST Digital Library of Mathematical Functions.
932 https://dlmf.nist.gov/10.6.E7
934 Examples
935 --------
936 Compute the Bessel function of the second kind of order 0 and
937 its first two derivatives at 1.
939 >>> from scipy.special import yvp
940 >>> yvp(0, 1, 0), yvp(0, 1, 1), yvp(0, 1, 2)
941 (0.088256964215677, 0.7812128213002889, -0.8694697855159659)
943 Compute the first derivative of the Bessel function of the second
944 kind for several orders at 1 by providing an array for `v`.
946 >>> yvp([0, 1, 2], 1, 1)
947 array([0.78121282, 0.86946979, 2.52015239])
949 Compute the first derivative of the Bessel function of the
950 second kind of order 0 at several points by providing an array for `z`.
952 >>> import numpy as np
953 >>> points = np.array([0.5, 1.5, 3.])
954 >>> yvp(0, points, 1)
955 array([ 1.47147239, 0.41230863, -0.32467442])
957 Plot the Bessel function of the second kind of order 1 and its
958 first three derivatives.
960 >>> import matplotlib.pyplot as plt
961 >>> x = np.linspace(0, 5, 1000)
962 >>> fig, ax = plt.subplots()
963 >>> ax.plot(x, yvp(1, x, 0), label=r"$Y_1$")
964 >>> ax.plot(x, yvp(1, x, 1), label=r"$Y_1'$")
965 >>> ax.plot(x, yvp(1, x, 2), label=r"$Y_1''$")
966 >>> ax.plot(x, yvp(1, x, 3), label=r"$Y_1'''$")
967 >>> ax.set_ylim(-10, 10)
968 >>> plt.legend()
969 >>> plt.show()
970 """
971 n = _nonneg_int_or_fail(n, 'n')
972 if n == 0:
973 return yv(v, z)
974 else:
975 return _bessel_diff_formula(v, z, n, yv, -1)
978def kvp(v, z, n=1):
979 """Compute derivatives of real-order modified Bessel function Kv(z)
981 Kv(z) is the modified Bessel function of the second kind.
982 Derivative is calculated with respect to `z`.
984 Parameters
985 ----------
986 v : array_like of float
987 Order of Bessel function
988 z : array_like of complex
989 Argument at which to evaluate the derivative
990 n : int, default 1
991 Order of derivative. For 0 returns the Bessel function `kv` itself.
993 Returns
994 -------
995 out : ndarray
996 The results
998 See Also
999 --------
1000 kv
1002 Notes
1003 -----
1004 The derivative is computed using the relation DLFM 10.29.5 [2]_.
1006 References
1007 ----------
1008 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
1009 Functions", John Wiley and Sons, 1996, chapter 6.
1010 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
1012 .. [2] NIST Digital Library of Mathematical Functions.
1013 https://dlmf.nist.gov/10.29.E5
1015 Examples
1016 --------
1017 Compute the modified bessel function of the second kind of order 0 and
1018 its first two derivatives at 1.
1020 >>> from scipy.special import kvp
1021 >>> kvp(0, 1, 0), kvp(0, 1, 1), kvp(0, 1, 2)
1022 (0.42102443824070834, -0.6019072301972346, 1.0229316684379428)
1024 Compute the first derivative of the modified Bessel function of the second
1025 kind for several orders at 1 by providing an array for `v`.
1027 >>> kvp([0, 1, 2], 1, 1)
1028 array([-0.60190723, -1.02293167, -3.85158503])
1030 Compute the first derivative of the modified Bessel function of the
1031 second kind of order 0 at several points by providing an array for `z`.
1033 >>> import numpy as np
1034 >>> points = np.array([0.5, 1.5, 3.])
1035 >>> kvp(0, points, 1)
1036 array([-1.65644112, -0.2773878 , -0.04015643])
1038 Plot the modified bessel function of the second kind and its
1039 first three derivatives.
1041 >>> import matplotlib.pyplot as plt
1042 >>> x = np.linspace(0, 5, 1000)
1043 >>> fig, ax = plt.subplots()
1044 >>> ax.plot(x, kvp(1, x, 0), label=r"$K_1$")
1045 >>> ax.plot(x, kvp(1, x, 1), label=r"$K_1'$")
1046 >>> ax.plot(x, kvp(1, x, 2), label=r"$K_1''$")
1047 >>> ax.plot(x, kvp(1, x, 3), label=r"$K_1'''$")
1048 >>> ax.set_ylim(-2.5, 2.5)
1049 >>> plt.legend()
1050 >>> plt.show()
1051 """
1052 n = _nonneg_int_or_fail(n, 'n')
1053 if n == 0:
1054 return kv(v, z)
1055 else:
1056 return (-1)**n * _bessel_diff_formula(v, z, n, kv, 1)
1059def ivp(v, z, n=1):
1060 """Compute derivatives of modified Bessel functions of the first kind.
1062 Compute the nth derivative of the modified Bessel function `Iv`
1063 with respect to `z`.
1065 Parameters
1066 ----------
1067 v : array_like or float
1068 Order of Bessel function
1069 z : array_like
1070 Argument at which to evaluate the derivative; can be real or
1071 complex.
1072 n : int, default 1
1073 Order of derivative. For 0, returns the Bessel function `iv` itself.
1075 Returns
1076 -------
1077 scalar or ndarray
1078 nth derivative of the modified Bessel function.
1080 See Also
1081 --------
1082 iv
1084 Notes
1085 -----
1086 The derivative is computed using the relation DLFM 10.29.5 [2]_.
1088 References
1089 ----------
1090 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
1091 Functions", John Wiley and Sons, 1996, chapter 6.
1092 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
1094 .. [2] NIST Digital Library of Mathematical Functions.
1095 https://dlmf.nist.gov/10.29.E5
1097 Examples
1098 --------
1099 Compute the modified Bessel function of the first kind of order 0 and
1100 its first two derivatives at 1.
1102 >>> from scipy.special import ivp
1103 >>> ivp(0, 1, 0), ivp(0, 1, 1), ivp(0, 1, 2)
1104 (1.2660658777520084, 0.565159103992485, 0.7009067737595233)
1106 Compute the first derivative of the modified Bessel function of the first
1107 kind for several orders at 1 by providing an array for `v`.
1109 >>> ivp([0, 1, 2], 1, 1)
1110 array([0.5651591 , 0.70090677, 0.29366376])
1112 Compute the first derivative of the modified Bessel function of the
1113 first kind of order 0 at several points by providing an array for `z`.
1115 >>> import numpy as np
1116 >>> points = np.array([0., 1.5, 3.])
1117 >>> ivp(0, points, 1)
1118 array([0. , 0.98166643, 3.95337022])
1120 Plot the modified Bessel function of the first kind of order 1 and its
1121 first three derivatives.
1123 >>> import matplotlib.pyplot as plt
1124 >>> x = np.linspace(-5, 5, 1000)
1125 >>> fig, ax = plt.subplots()
1126 >>> ax.plot(x, ivp(1, x, 0), label=r"$I_1$")
1127 >>> ax.plot(x, ivp(1, x, 1), label=r"$I_1'$")
1128 >>> ax.plot(x, ivp(1, x, 2), label=r"$I_1''$")
1129 >>> ax.plot(x, ivp(1, x, 3), label=r"$I_1'''$")
1130 >>> plt.legend()
1131 >>> plt.show()
1132 """
1133 n = _nonneg_int_or_fail(n, 'n')
1134 if n == 0:
1135 return iv(v, z)
1136 else:
1137 return _bessel_diff_formula(v, z, n, iv, 1)
1140def h1vp(v, z, n=1):
1141 """Compute derivatives of Hankel function H1v(z) with respect to `z`.
1143 Parameters
1144 ----------
1145 v : array_like
1146 Order of Hankel function
1147 z : array_like
1148 Argument at which to evaluate the derivative. Can be real or
1149 complex.
1150 n : int, default 1
1151 Order of derivative. For 0 returns the Hankel function `h1v` itself.
1153 Returns
1154 -------
1155 scalar or ndarray
1156 Values of the derivative of the Hankel function.
1158 See Also
1159 --------
1160 hankel1
1162 Notes
1163 -----
1164 The derivative is computed using the relation DLFM 10.6.7 [2]_.
1166 References
1167 ----------
1168 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
1169 Functions", John Wiley and Sons, 1996, chapter 5.
1170 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
1172 .. [2] NIST Digital Library of Mathematical Functions.
1173 https://dlmf.nist.gov/10.6.E7
1175 Examples
1176 --------
1177 Compute the Hankel function of the first kind of order 0 and
1178 its first two derivatives at 1.
1180 >>> from scipy.special import h1vp
1181 >>> h1vp(0, 1, 0), h1vp(0, 1, 1), h1vp(0, 1, 2)
1182 ((0.7651976865579664+0.088256964215677j),
1183 (-0.44005058574493355+0.7812128213002889j),
1184 (-0.3251471008130329-0.8694697855159659j))
1186 Compute the first derivative of the Hankel function of the first kind
1187 for several orders at 1 by providing an array for `v`.
1189 >>> h1vp([0, 1, 2], 1, 1)
1190 array([-0.44005059+0.78121282j, 0.3251471 +0.86946979j,
1191 0.21024362+2.52015239j])
1193 Compute the first derivative of the Hankel function of the first kind
1194 of order 0 at several points by providing an array for `z`.
1196 >>> import numpy as np
1197 >>> points = np.array([0.5, 1.5, 3.])
1198 >>> h1vp(0, points, 1)
1199 array([-0.24226846+1.47147239j, -0.55793651+0.41230863j,
1200 -0.33905896-0.32467442j])
1201 """
1202 n = _nonneg_int_or_fail(n, 'n')
1203 if n == 0:
1204 return hankel1(v, z)
1205 else:
1206 return _bessel_diff_formula(v, z, n, hankel1, -1)
1209def h2vp(v, z, n=1):
1210 """Compute derivatives of Hankel function H2v(z) with respect to `z`.
1212 Parameters
1213 ----------
1214 v : array_like
1215 Order of Hankel function
1216 z : array_like
1217 Argument at which to evaluate the derivative. Can be real or
1218 complex.
1219 n : int, default 1
1220 Order of derivative. For 0 returns the Hankel function `h2v` itself.
1222 Returns
1223 -------
1224 scalar or ndarray
1225 Values of the derivative of the Hankel function.
1227 See Also
1228 --------
1229 hankel2
1231 Notes
1232 -----
1233 The derivative is computed using the relation DLFM 10.6.7 [2]_.
1235 References
1236 ----------
1237 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
1238 Functions", John Wiley and Sons, 1996, chapter 5.
1239 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
1241 .. [2] NIST Digital Library of Mathematical Functions.
1242 https://dlmf.nist.gov/10.6.E7
1244 Examples
1245 --------
1246 Compute the Hankel function of the second kind of order 0 and
1247 its first two derivatives at 1.
1249 >>> from scipy.special import h2vp
1250 >>> h2vp(0, 1, 0), h2vp(0, 1, 1), h2vp(0, 1, 2)
1251 ((0.7651976865579664-0.088256964215677j),
1252 (-0.44005058574493355-0.7812128213002889j),
1253 (-0.3251471008130329+0.8694697855159659j))
1255 Compute the first derivative of the Hankel function of the second kind
1256 for several orders at 1 by providing an array for `v`.
1258 >>> h2vp([0, 1, 2], 1, 1)
1259 array([-0.44005059-0.78121282j, 0.3251471 -0.86946979j,
1260 0.21024362-2.52015239j])
1262 Compute the first derivative of the Hankel function of the second kind
1263 of order 0 at several points by providing an array for `z`.
1265 >>> import numpy as np
1266 >>> points = np.array([0.5, 1.5, 3.])
1267 >>> h2vp(0, points, 1)
1268 array([-0.24226846-1.47147239j, -0.55793651-0.41230863j,
1269 -0.33905896+0.32467442j])
1270 """
1271 n = _nonneg_int_or_fail(n, 'n')
1272 if n == 0:
1273 return hankel2(v, z)
1274 else:
1275 return _bessel_diff_formula(v, z, n, hankel2, -1)
1278def riccati_jn(n, x):
1279 r"""Compute Ricatti-Bessel function of the first kind and its derivative.
1281 The Ricatti-Bessel function of the first kind is defined as :math:`x
1282 j_n(x)`, where :math:`j_n` is the spherical Bessel function of the first
1283 kind of order :math:`n`.
1285 This function computes the value and first derivative of the
1286 Ricatti-Bessel function for all orders up to and including `n`.
1288 Parameters
1289 ----------
1290 n : int
1291 Maximum order of function to compute
1292 x : float
1293 Argument at which to evaluate
1295 Returns
1296 -------
1297 jn : ndarray
1298 Value of j0(x), ..., jn(x)
1299 jnp : ndarray
1300 First derivative j0'(x), ..., jn'(x)
1302 Notes
1303 -----
1304 The computation is carried out via backward recurrence, using the
1305 relation DLMF 10.51.1 [2]_.
1307 Wrapper for a Fortran routine created by Shanjie Zhang and Jianming
1308 Jin [1]_.
1310 References
1311 ----------
1312 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
1313 Functions", John Wiley and Sons, 1996.
1314 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
1315 .. [2] NIST Digital Library of Mathematical Functions.
1316 https://dlmf.nist.gov/10.51.E1
1318 """
1319 if not (isscalar(n) and isscalar(x)):
1320 raise ValueError("arguments must be scalars.")
1321 n = _nonneg_int_or_fail(n, 'n', strict=False)
1322 if (n == 0):
1323 n1 = 1
1324 else:
1325 n1 = n
1326 nm, jn, jnp = _specfun.rctj(n1, x)
1327 return jn[:(n+1)], jnp[:(n+1)]
1330def riccati_yn(n, x):
1331 """Compute Ricatti-Bessel function of the second kind and its derivative.
1333 The Ricatti-Bessel function of the second kind is defined as :math:`x
1334 y_n(x)`, where :math:`y_n` is the spherical Bessel function of the second
1335 kind of order :math:`n`.
1337 This function computes the value and first derivative of the function for
1338 all orders up to and including `n`.
1340 Parameters
1341 ----------
1342 n : int
1343 Maximum order of function to compute
1344 x : float
1345 Argument at which to evaluate
1347 Returns
1348 -------
1349 yn : ndarray
1350 Value of y0(x), ..., yn(x)
1351 ynp : ndarray
1352 First derivative y0'(x), ..., yn'(x)
1354 Notes
1355 -----
1356 The computation is carried out via ascending recurrence, using the
1357 relation DLMF 10.51.1 [2]_.
1359 Wrapper for a Fortran routine created by Shanjie Zhang and Jianming
1360 Jin [1]_.
1362 References
1363 ----------
1364 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
1365 Functions", John Wiley and Sons, 1996.
1366 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
1367 .. [2] NIST Digital Library of Mathematical Functions.
1368 https://dlmf.nist.gov/10.51.E1
1370 """
1371 if not (isscalar(n) and isscalar(x)):
1372 raise ValueError("arguments must be scalars.")
1373 n = _nonneg_int_or_fail(n, 'n', strict=False)
1374 if (n == 0):
1375 n1 = 1
1376 else:
1377 n1 = n
1378 nm, jn, jnp = _specfun.rcty(n1, x)
1379 return jn[:(n+1)], jnp[:(n+1)]
1382def erf_zeros(nt):
1383 """Compute the first nt zero in the first quadrant, ordered by absolute value.
1385 Zeros in the other quadrants can be obtained by using the symmetries erf(-z) = erf(z) and
1386 erf(conj(z)) = conj(erf(z)).
1389 Parameters
1390 ----------
1391 nt : int
1392 The number of zeros to compute
1394 Returns
1395 -------
1396 The locations of the zeros of erf : ndarray (complex)
1397 Complex values at which zeros of erf(z)
1399 Examples
1400 --------
1401 >>> from scipy import special
1402 >>> special.erf_zeros(1)
1403 array([1.45061616+1.880943j])
1405 Check that erf is (close to) zero for the value returned by erf_zeros
1407 >>> special.erf(special.erf_zeros(1))
1408 array([4.95159469e-14-1.16407394e-16j])
1410 References
1411 ----------
1412 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
1413 Functions", John Wiley and Sons, 1996.
1414 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
1416 """
1417 if (floor(nt) != nt) or (nt <= 0) or not isscalar(nt):
1418 raise ValueError("Argument must be positive scalar integer.")
1419 return _specfun.cerzo(nt)
1422def fresnelc_zeros(nt):
1423 """Compute nt complex zeros of cosine Fresnel integral C(z).
1425 References
1426 ----------
1427 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
1428 Functions", John Wiley and Sons, 1996.
1429 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
1431 """
1432 if (floor(nt) != nt) or (nt <= 0) or not isscalar(nt):
1433 raise ValueError("Argument must be positive scalar integer.")
1434 return _specfun.fcszo(1, nt)
1437def fresnels_zeros(nt):
1438 """Compute nt complex zeros of sine Fresnel integral S(z).
1440 References
1441 ----------
1442 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
1443 Functions", John Wiley and Sons, 1996.
1444 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
1446 """
1447 if (floor(nt) != nt) or (nt <= 0) or not isscalar(nt):
1448 raise ValueError("Argument must be positive scalar integer.")
1449 return _specfun.fcszo(2, nt)
1452def fresnel_zeros(nt):
1453 """Compute nt complex zeros of sine and cosine Fresnel integrals S(z) and C(z).
1455 References
1456 ----------
1457 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
1458 Functions", John Wiley and Sons, 1996.
1459 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
1461 """
1462 if (floor(nt) != nt) or (nt <= 0) or not isscalar(nt):
1463 raise ValueError("Argument must be positive scalar integer.")
1464 return _specfun.fcszo(2, nt), _specfun.fcszo(1, nt)
1467def assoc_laguerre(x, n, k=0.0):
1468 """Compute the generalized (associated) Laguerre polynomial of degree n and order k.
1470 The polynomial :math:`L^{(k)}_n(x)` is orthogonal over ``[0, inf)``,
1471 with weighting function ``exp(-x) * x**k`` with ``k > -1``.
1473 Notes
1474 -----
1475 `assoc_laguerre` is a simple wrapper around `eval_genlaguerre`, with
1476 reversed argument order ``(x, n, k=0.0) --> (n, k, x)``.
1478 """
1479 return _ufuncs.eval_genlaguerre(n, k, x)
1482digamma = psi
1485def polygamma(n, x):
1486 r"""Polygamma functions.
1488 Defined as :math:`\psi^{(n)}(x)` where :math:`\psi` is the
1489 `digamma` function. See [dlmf]_ for details.
1491 Parameters
1492 ----------
1493 n : array_like
1494 The order of the derivative of the digamma function; must be
1495 integral
1496 x : array_like
1497 Real valued input
1499 Returns
1500 -------
1501 ndarray
1502 Function results
1504 See Also
1505 --------
1506 digamma
1508 References
1509 ----------
1510 .. [dlmf] NIST, Digital Library of Mathematical Functions,
1511 https://dlmf.nist.gov/5.15
1513 Examples
1514 --------
1515 >>> from scipy import special
1516 >>> x = [2, 3, 25.5]
1517 >>> special.polygamma(1, x)
1518 array([ 0.64493407, 0.39493407, 0.03999467])
1519 >>> special.polygamma(0, x) == special.psi(x)
1520 array([ True, True, True], dtype=bool)
1522 """
1523 n, x = asarray(n), asarray(x)
1524 fac2 = (-1.0)**(n+1) * gamma(n+1.0) * zeta(n+1, x)
1525 return where(n == 0, psi(x), fac2)
1528def mathieu_even_coef(m, q):
1529 r"""Fourier coefficients for even Mathieu and modified Mathieu functions.
1531 The Fourier series of the even solutions of the Mathieu differential
1532 equation are of the form
1534 .. math:: \mathrm{ce}_{2n}(z, q) = \sum_{k=0}^{\infty} A_{(2n)}^{(2k)} \cos 2kz
1536 .. math:: \mathrm{ce}_{2n+1}(z, q) = \sum_{k=0}^{\infty} A_{(2n+1)}^{(2k+1)} \cos (2k+1)z
1538 This function returns the coefficients :math:`A_{(2n)}^{(2k)}` for even
1539 input m=2n, and the coefficients :math:`A_{(2n+1)}^{(2k+1)}` for odd input
1540 m=2n+1.
1542 Parameters
1543 ----------
1544 m : int
1545 Order of Mathieu functions. Must be non-negative.
1546 q : float (>=0)
1547 Parameter of Mathieu functions. Must be non-negative.
1549 Returns
1550 -------
1551 Ak : ndarray
1552 Even or odd Fourier coefficients, corresponding to even or odd m.
1554 References
1555 ----------
1556 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
1557 Functions", John Wiley and Sons, 1996.
1558 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
1559 .. [2] NIST Digital Library of Mathematical Functions
1560 https://dlmf.nist.gov/28.4#i
1562 """
1563 if not (isscalar(m) and isscalar(q)):
1564 raise ValueError("m and q must be scalars.")
1565 if (q < 0):
1566 raise ValueError("q >=0")
1567 if (m != floor(m)) or (m < 0):
1568 raise ValueError("m must be an integer >=0.")
1570 if (q <= 1):
1571 qm = 7.5 + 56.1*sqrt(q) - 134.7*q + 90.7*sqrt(q)*q
1572 else:
1573 qm = 17.0 + 3.1*sqrt(q) - .126*q + .0037*sqrt(q)*q
1574 km = int(qm + 0.5*m)
1575 if km > 251:
1576 warnings.warn("Too many predicted coefficients.", RuntimeWarning, 2)
1577 kd = 1
1578 m = int(floor(m))
1579 if m % 2:
1580 kd = 2
1582 a = mathieu_a(m, q)
1583 fc = _specfun.fcoef(kd, m, q, a)
1584 return fc[:km]
1587def mathieu_odd_coef(m, q):
1588 r"""Fourier coefficients for even Mathieu and modified Mathieu functions.
1590 The Fourier series of the odd solutions of the Mathieu differential
1591 equation are of the form
1593 .. math:: \mathrm{se}_{2n+1}(z, q) = \sum_{k=0}^{\infty} B_{(2n+1)}^{(2k+1)} \sin (2k+1)z
1595 .. math:: \mathrm{se}_{2n+2}(z, q) = \sum_{k=0}^{\infty} B_{(2n+2)}^{(2k+2)} \sin (2k+2)z
1597 This function returns the coefficients :math:`B_{(2n+2)}^{(2k+2)}` for even
1598 input m=2n+2, and the coefficients :math:`B_{(2n+1)}^{(2k+1)}` for odd
1599 input m=2n+1.
1601 Parameters
1602 ----------
1603 m : int
1604 Order of Mathieu functions. Must be non-negative.
1605 q : float (>=0)
1606 Parameter of Mathieu functions. Must be non-negative.
1608 Returns
1609 -------
1610 Bk : ndarray
1611 Even or odd Fourier coefficients, corresponding to even or odd m.
1613 References
1614 ----------
1615 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
1616 Functions", John Wiley and Sons, 1996.
1617 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
1619 """
1620 if not (isscalar(m) and isscalar(q)):
1621 raise ValueError("m and q must be scalars.")
1622 if (q < 0):
1623 raise ValueError("q >=0")
1624 if (m != floor(m)) or (m <= 0):
1625 raise ValueError("m must be an integer > 0")
1627 if (q <= 1):
1628 qm = 7.5 + 56.1*sqrt(q) - 134.7*q + 90.7*sqrt(q)*q
1629 else:
1630 qm = 17.0 + 3.1*sqrt(q) - .126*q + .0037*sqrt(q)*q
1631 km = int(qm + 0.5*m)
1632 if km > 251:
1633 warnings.warn("Too many predicted coefficients.", RuntimeWarning, 2)
1634 kd = 4
1635 m = int(floor(m))
1636 if m % 2:
1637 kd = 3
1639 b = mathieu_b(m, q)
1640 fc = _specfun.fcoef(kd, m, q, b)
1641 return fc[:km]
1644def lpmn(m, n, z):
1645 """Sequence of associated Legendre functions of the first kind.
1647 Computes the associated Legendre function of the first kind of order m and
1648 degree n, ``Pmn(z)`` = :math:`P_n^m(z)`, and its derivative, ``Pmn'(z)``.
1649 Returns two arrays of size ``(m+1, n+1)`` containing ``Pmn(z)`` and
1650 ``Pmn'(z)`` for all orders from ``0..m`` and degrees from ``0..n``.
1652 This function takes a real argument ``z``. For complex arguments ``z``
1653 use clpmn instead.
1655 Parameters
1656 ----------
1657 m : int
1658 ``|m| <= n``; the order of the Legendre function.
1659 n : int
1660 where ``n >= 0``; the degree of the Legendre function. Often
1661 called ``l`` (lower case L) in descriptions of the associated
1662 Legendre function
1663 z : float
1664 Input value.
1666 Returns
1667 -------
1668 Pmn_z : (m+1, n+1) array
1669 Values for all orders 0..m and degrees 0..n
1670 Pmn_d_z : (m+1, n+1) array
1671 Derivatives for all orders 0..m and degrees 0..n
1673 See Also
1674 --------
1675 clpmn: associated Legendre functions of the first kind for complex z
1677 Notes
1678 -----
1679 In the interval (-1, 1), Ferrer's function of the first kind is
1680 returned. The phase convention used for the intervals (1, inf)
1681 and (-inf, -1) is such that the result is always real.
1683 References
1684 ----------
1685 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
1686 Functions", John Wiley and Sons, 1996.
1687 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
1688 .. [2] NIST Digital Library of Mathematical Functions
1689 https://dlmf.nist.gov/14.3
1691 """
1692 if not isscalar(m) or (abs(m) > n):
1693 raise ValueError("m must be <= n.")
1694 if not isscalar(n) or (n < 0):
1695 raise ValueError("n must be a non-negative integer.")
1696 if not isscalar(z):
1697 raise ValueError("z must be scalar.")
1698 if iscomplex(z):
1699 raise ValueError("Argument must be real. Use clpmn instead.")
1700 if (m < 0):
1701 mp = -m
1702 mf, nf = mgrid[0:mp+1, 0:n+1]
1703 with _ufuncs.errstate(all='ignore'):
1704 if abs(z) < 1:
1705 # Ferrer function; DLMF 14.9.3
1706 fixarr = where(mf > nf, 0.0,
1707 (-1)**mf * gamma(nf-mf+1) / gamma(nf+mf+1))
1708 else:
1709 # Match to clpmn; DLMF 14.9.13
1710 fixarr = where(mf > nf, 0.0, gamma(nf-mf+1) / gamma(nf+mf+1))
1711 else:
1712 mp = m
1713 p, pd = _specfun.lpmn(mp, n, z)
1714 if (m < 0):
1715 p = p * fixarr
1716 pd = pd * fixarr
1717 return p, pd
1720def clpmn(m, n, z, type=3):
1721 """Associated Legendre function of the first kind for complex arguments.
1723 Computes the associated Legendre function of the first kind of order m and
1724 degree n, ``Pmn(z)`` = :math:`P_n^m(z)`, and its derivative, ``Pmn'(z)``.
1725 Returns two arrays of size ``(m+1, n+1)`` containing ``Pmn(z)`` and
1726 ``Pmn'(z)`` for all orders from ``0..m`` and degrees from ``0..n``.
1728 Parameters
1729 ----------
1730 m : int
1731 ``|m| <= n``; the order of the Legendre function.
1732 n : int
1733 where ``n >= 0``; the degree of the Legendre function. Often
1734 called ``l`` (lower case L) in descriptions of the associated
1735 Legendre function
1736 z : float or complex
1737 Input value.
1738 type : int, optional
1739 takes values 2 or 3
1740 2: cut on the real axis ``|x| > 1``
1741 3: cut on the real axis ``-1 < x < 1`` (default)
1743 Returns
1744 -------
1745 Pmn_z : (m+1, n+1) array
1746 Values for all orders ``0..m`` and degrees ``0..n``
1747 Pmn_d_z : (m+1, n+1) array
1748 Derivatives for all orders ``0..m`` and degrees ``0..n``
1750 See Also
1751 --------
1752 lpmn: associated Legendre functions of the first kind for real z
1754 Notes
1755 -----
1756 By default, i.e. for ``type=3``, phase conventions are chosen according
1757 to [1]_ such that the function is analytic. The cut lies on the interval
1758 (-1, 1). Approaching the cut from above or below in general yields a phase
1759 factor with respect to Ferrer's function of the first kind
1760 (cf. `lpmn`).
1762 For ``type=2`` a cut at ``|x| > 1`` is chosen. Approaching the real values
1763 on the interval (-1, 1) in the complex plane yields Ferrer's function
1764 of the first kind.
1766 References
1767 ----------
1768 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
1769 Functions", John Wiley and Sons, 1996.
1770 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
1771 .. [2] NIST Digital Library of Mathematical Functions
1772 https://dlmf.nist.gov/14.21
1774 """
1775 if not isscalar(m) or (abs(m) > n):
1776 raise ValueError("m must be <= n.")
1777 if not isscalar(n) or (n < 0):
1778 raise ValueError("n must be a non-negative integer.")
1779 if not isscalar(z):
1780 raise ValueError("z must be scalar.")
1781 if not (type == 2 or type == 3):
1782 raise ValueError("type must be either 2 or 3.")
1783 if (m < 0):
1784 mp = -m
1785 mf, nf = mgrid[0:mp+1, 0:n+1]
1786 with _ufuncs.errstate(all='ignore'):
1787 if type == 2:
1788 fixarr = where(mf > nf, 0.0,
1789 (-1)**mf * gamma(nf-mf+1) / gamma(nf+mf+1))
1790 else:
1791 fixarr = where(mf > nf, 0.0, gamma(nf-mf+1) / gamma(nf+mf+1))
1792 else:
1793 mp = m
1794 p, pd = _specfun.clpmn(mp, n, real(z), imag(z), type)
1795 if (m < 0):
1796 p = p * fixarr
1797 pd = pd * fixarr
1798 return p, pd
1801def lqmn(m, n, z):
1802 """Sequence of associated Legendre functions of the second kind.
1804 Computes the associated Legendre function of the second kind of order m and
1805 degree n, ``Qmn(z)`` = :math:`Q_n^m(z)`, and its derivative, ``Qmn'(z)``.
1806 Returns two arrays of size ``(m+1, n+1)`` containing ``Qmn(z)`` and
1807 ``Qmn'(z)`` for all orders from ``0..m`` and degrees from ``0..n``.
1809 Parameters
1810 ----------
1811 m : int
1812 ``|m| <= n``; the order of the Legendre function.
1813 n : int
1814 where ``n >= 0``; the degree of the Legendre function. Often
1815 called ``l`` (lower case L) in descriptions of the associated
1816 Legendre function
1817 z : complex
1818 Input value.
1820 Returns
1821 -------
1822 Qmn_z : (m+1, n+1) array
1823 Values for all orders 0..m and degrees 0..n
1824 Qmn_d_z : (m+1, n+1) array
1825 Derivatives for all orders 0..m and degrees 0..n
1827 References
1828 ----------
1829 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
1830 Functions", John Wiley and Sons, 1996.
1831 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
1833 """
1834 if not isscalar(m) or (m < 0):
1835 raise ValueError("m must be a non-negative integer.")
1836 if not isscalar(n) or (n < 0):
1837 raise ValueError("n must be a non-negative integer.")
1838 if not isscalar(z):
1839 raise ValueError("z must be scalar.")
1840 m = int(m)
1841 n = int(n)
1843 # Ensure neither m nor n == 0
1844 mm = max(1, m)
1845 nn = max(1, n)
1847 if iscomplex(z):
1848 q, qd = _specfun.clqmn(mm, nn, z)
1849 else:
1850 q, qd = _specfun.lqmn(mm, nn, z)
1851 return q[:(m+1), :(n+1)], qd[:(m+1), :(n+1)]
1854def bernoulli(n):
1855 """Bernoulli numbers B0..Bn (inclusive).
1857 Parameters
1858 ----------
1859 n : int
1860 Indicated the number of terms in the Bernoulli series to generate.
1862 Returns
1863 -------
1864 ndarray
1865 The Bernoulli numbers ``[B(0), B(1), ..., B(n)]``.
1867 References
1868 ----------
1869 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
1870 Functions", John Wiley and Sons, 1996.
1871 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
1872 .. [2] "Bernoulli number", Wikipedia, https://en.wikipedia.org/wiki/Bernoulli_number
1874 Examples
1875 --------
1876 >>> import numpy as np
1877 >>> from scipy.special import bernoulli, zeta
1878 >>> bernoulli(4)
1879 array([ 1. , -0.5 , 0.16666667, 0. , -0.03333333])
1881 The Wikipedia article ([2]_) points out the relationship between the
1882 Bernoulli numbers and the zeta function, ``B_n^+ = -n * zeta(1 - n)``
1883 for ``n > 0``:
1885 >>> n = np.arange(1, 5)
1886 >>> -n * zeta(1 - n)
1887 array([ 0.5 , 0.16666667, -0. , -0.03333333])
1889 Note that, in the notation used in the wikipedia article,
1890 `bernoulli` computes ``B_n^-`` (i.e. it used the convention that
1891 ``B_1`` is -1/2). The relation given above is for ``B_n^+``, so the
1892 sign of 0.5 does not match the output of ``bernoulli(4)``.
1894 """
1895 if not isscalar(n) or (n < 0):
1896 raise ValueError("n must be a non-negative integer.")
1897 n = int(n)
1898 if (n < 2):
1899 n1 = 2
1900 else:
1901 n1 = n
1902 return _specfun.bernob(int(n1))[:(n+1)]
1905def euler(n):
1906 """Euler numbers E(0), E(1), ..., E(n).
1908 The Euler numbers [1]_ are also known as the secant numbers.
1910 Because ``euler(n)`` returns floating point values, it does not give
1911 exact values for large `n`. The first inexact value is E(22).
1913 Parameters
1914 ----------
1915 n : int
1916 The highest index of the Euler number to be returned.
1918 Returns
1919 -------
1920 ndarray
1921 The Euler numbers [E(0), E(1), ..., E(n)].
1922 The odd Euler numbers, which are all zero, are included.
1924 References
1925 ----------
1926 .. [1] Sequence A122045, The On-Line Encyclopedia of Integer Sequences,
1927 https://oeis.org/A122045
1928 .. [2] Zhang, Shanjie and Jin, Jianming. "Computation of Special
1929 Functions", John Wiley and Sons, 1996.
1930 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
1932 Examples
1933 --------
1934 >>> import numpy as np
1935 >>> from scipy.special import euler
1936 >>> euler(6)
1937 array([ 1., 0., -1., 0., 5., 0., -61.])
1939 >>> euler(13).astype(np.int64)
1940 array([ 1, 0, -1, 0, 5, 0, -61,
1941 0, 1385, 0, -50521, 0, 2702765, 0])
1943 >>> euler(22)[-1] # Exact value of E(22) is -69348874393137901.
1944 -69348874393137976.0
1946 """
1947 if not isscalar(n) or (n < 0):
1948 raise ValueError("n must be a non-negative integer.")
1949 n = int(n)
1950 if (n < 2):
1951 n1 = 2
1952 else:
1953 n1 = n
1954 return _specfun.eulerb(n1)[:(n+1)]
1957def lpn(n, z):
1958 """Legendre function of the first kind.
1960 Compute sequence of Legendre functions of the first kind (polynomials),
1961 Pn(z) and derivatives for all degrees from 0 to n (inclusive).
1963 See also special.legendre for polynomial class.
1965 References
1966 ----------
1967 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
1968 Functions", John Wiley and Sons, 1996.
1969 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
1971 """
1972 if not (isscalar(n) and isscalar(z)):
1973 raise ValueError("arguments must be scalars.")
1974 n = _nonneg_int_or_fail(n, 'n', strict=False)
1975 if (n < 1):
1976 n1 = 1
1977 else:
1978 n1 = n
1979 if iscomplex(z):
1980 pn, pd = _specfun.clpn(n1, z)
1981 else:
1982 pn, pd = _specfun.lpn(n1, z)
1983 return pn[:(n+1)], pd[:(n+1)]
1986def lqn(n, z):
1987 """Legendre function of the second kind.
1989 Compute sequence of Legendre functions of the second kind, Qn(z) and
1990 derivatives for all degrees from 0 to n (inclusive).
1992 References
1993 ----------
1994 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
1995 Functions", John Wiley and Sons, 1996.
1996 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
1998 """
1999 if not (isscalar(n) and isscalar(z)):
2000 raise ValueError("arguments must be scalars.")
2001 n = _nonneg_int_or_fail(n, 'n', strict=False)
2002 if (n < 1):
2003 n1 = 1
2004 else:
2005 n1 = n
2006 if iscomplex(z):
2007 qn, qd = _specfun.clqn(n1, z)
2008 else:
2009 qn, qd = _specfun.lqnb(n1, z)
2010 return qn[:(n+1)], qd[:(n+1)]
2013def ai_zeros(nt):
2014 """
2015 Compute `nt` zeros and values of the Airy function Ai and its derivative.
2017 Computes the first `nt` zeros, `a`, of the Airy function Ai(x);
2018 first `nt` zeros, `ap`, of the derivative of the Airy function Ai'(x);
2019 the corresponding values Ai(a');
2020 and the corresponding values Ai'(a).
2022 Parameters
2023 ----------
2024 nt : int
2025 Number of zeros to compute
2027 Returns
2028 -------
2029 a : ndarray
2030 First `nt` zeros of Ai(x)
2031 ap : ndarray
2032 First `nt` zeros of Ai'(x)
2033 ai : ndarray
2034 Values of Ai(x) evaluated at first `nt` zeros of Ai'(x)
2035 aip : ndarray
2036 Values of Ai'(x) evaluated at first `nt` zeros of Ai(x)
2038 Examples
2039 --------
2040 >>> from scipy import special
2041 >>> a, ap, ai, aip = special.ai_zeros(3)
2042 >>> a
2043 array([-2.33810741, -4.08794944, -5.52055983])
2044 >>> ap
2045 array([-1.01879297, -3.24819758, -4.82009921])
2046 >>> ai
2047 array([ 0.53565666, -0.41901548, 0.38040647])
2048 >>> aip
2049 array([ 0.70121082, -0.80311137, 0.86520403])
2051 References
2052 ----------
2053 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
2054 Functions", John Wiley and Sons, 1996.
2055 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
2057 """
2058 kf = 1
2059 if not isscalar(nt) or (floor(nt) != nt) or (nt <= 0):
2060 raise ValueError("nt must be a positive integer scalar.")
2061 return _specfun.airyzo(nt, kf)
2064def bi_zeros(nt):
2065 """
2066 Compute `nt` zeros and values of the Airy function Bi and its derivative.
2068 Computes the first `nt` zeros, b, of the Airy function Bi(x);
2069 first `nt` zeros, b', of the derivative of the Airy function Bi'(x);
2070 the corresponding values Bi(b');
2071 and the corresponding values Bi'(b).
2073 Parameters
2074 ----------
2075 nt : int
2076 Number of zeros to compute
2078 Returns
2079 -------
2080 b : ndarray
2081 First `nt` zeros of Bi(x)
2082 bp : ndarray
2083 First `nt` zeros of Bi'(x)
2084 bi : ndarray
2085 Values of Bi(x) evaluated at first `nt` zeros of Bi'(x)
2086 bip : ndarray
2087 Values of Bi'(x) evaluated at first `nt` zeros of Bi(x)
2089 Examples
2090 --------
2091 >>> from scipy import special
2092 >>> b, bp, bi, bip = special.bi_zeros(3)
2093 >>> b
2094 array([-1.17371322, -3.2710933 , -4.83073784])
2095 >>> bp
2096 array([-2.29443968, -4.07315509, -5.51239573])
2097 >>> bi
2098 array([-0.45494438, 0.39652284, -0.36796916])
2099 >>> bip
2100 array([ 0.60195789, -0.76031014, 0.83699101])
2102 References
2103 ----------
2104 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
2105 Functions", John Wiley and Sons, 1996.
2106 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
2108 """
2109 kf = 2
2110 if not isscalar(nt) or (floor(nt) != nt) or (nt <= 0):
2111 raise ValueError("nt must be a positive integer scalar.")
2112 return _specfun.airyzo(nt, kf)
2115def lmbda(v, x):
2116 r"""Jahnke-Emden Lambda function, Lambdav(x).
2118 This function is defined as [2]_,
2120 .. math:: \Lambda_v(x) = \Gamma(v+1) \frac{J_v(x)}{(x/2)^v},
2122 where :math:`\Gamma` is the gamma function and :math:`J_v` is the
2123 Bessel function of the first kind.
2125 Parameters
2126 ----------
2127 v : float
2128 Order of the Lambda function
2129 x : float
2130 Value at which to evaluate the function and derivatives
2132 Returns
2133 -------
2134 vl : ndarray
2135 Values of Lambda_vi(x), for vi=v-int(v), vi=1+v-int(v), ..., vi=v.
2136 dl : ndarray
2137 Derivatives Lambda_vi'(x), for vi=v-int(v), vi=1+v-int(v), ..., vi=v.
2139 References
2140 ----------
2141 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
2142 Functions", John Wiley and Sons, 1996.
2143 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
2144 .. [2] Jahnke, E. and Emde, F. "Tables of Functions with Formulae and
2145 Curves" (4th ed.), Dover, 1945
2146 """
2147 if not (isscalar(v) and isscalar(x)):
2148 raise ValueError("arguments must be scalars.")
2149 if (v < 0):
2150 raise ValueError("argument must be > 0.")
2151 n = int(v)
2152 v0 = v - n
2153 if (n < 1):
2154 n1 = 1
2155 else:
2156 n1 = n
2157 v1 = n1 + v0
2158 if (v != floor(v)):
2159 vm, vl, dl = _specfun.lamv(v1, x)
2160 else:
2161 vm, vl, dl = _specfun.lamn(v1, x)
2162 return vl[:(n+1)], dl[:(n+1)]
2165def pbdv_seq(v, x):
2166 """Parabolic cylinder functions Dv(x) and derivatives.
2168 Parameters
2169 ----------
2170 v : float
2171 Order of the parabolic cylinder function
2172 x : float
2173 Value at which to evaluate the function and derivatives
2175 Returns
2176 -------
2177 dv : ndarray
2178 Values of D_vi(x), for vi=v-int(v), vi=1+v-int(v), ..., vi=v.
2179 dp : ndarray
2180 Derivatives D_vi'(x), for vi=v-int(v), vi=1+v-int(v), ..., vi=v.
2182 References
2183 ----------
2184 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
2185 Functions", John Wiley and Sons, 1996, chapter 13.
2186 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
2188 """
2189 if not (isscalar(v) and isscalar(x)):
2190 raise ValueError("arguments must be scalars.")
2191 n = int(v)
2192 v0 = v-n
2193 if (n < 1):
2194 n1 = 1
2195 else:
2196 n1 = n
2197 v1 = n1 + v0
2198 dv, dp, pdf, pdd = _specfun.pbdv(v1, x)
2199 return dv[:n1+1], dp[:n1+1]
2202def pbvv_seq(v, x):
2203 """Parabolic cylinder functions Vv(x) and derivatives.
2205 Parameters
2206 ----------
2207 v : float
2208 Order of the parabolic cylinder function
2209 x : float
2210 Value at which to evaluate the function and derivatives
2212 Returns
2213 -------
2214 dv : ndarray
2215 Values of V_vi(x), for vi=v-int(v), vi=1+v-int(v), ..., vi=v.
2216 dp : ndarray
2217 Derivatives V_vi'(x), for vi=v-int(v), vi=1+v-int(v), ..., vi=v.
2219 References
2220 ----------
2221 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
2222 Functions", John Wiley and Sons, 1996, chapter 13.
2223 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
2225 """
2226 if not (isscalar(v) and isscalar(x)):
2227 raise ValueError("arguments must be scalars.")
2228 n = int(v)
2229 v0 = v-n
2230 if (n <= 1):
2231 n1 = 1
2232 else:
2233 n1 = n
2234 v1 = n1 + v0
2235 dv, dp, pdf, pdd = _specfun.pbvv(v1, x)
2236 return dv[:n1+1], dp[:n1+1]
2239def pbdn_seq(n, z):
2240 """Parabolic cylinder functions Dn(z) and derivatives.
2242 Parameters
2243 ----------
2244 n : int
2245 Order of the parabolic cylinder function
2246 z : complex
2247 Value at which to evaluate the function and derivatives
2249 Returns
2250 -------
2251 dv : ndarray
2252 Values of D_i(z), for i=0, ..., i=n.
2253 dp : ndarray
2254 Derivatives D_i'(z), for i=0, ..., i=n.
2256 References
2257 ----------
2258 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
2259 Functions", John Wiley and Sons, 1996, chapter 13.
2260 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
2262 """
2263 if not (isscalar(n) and isscalar(z)):
2264 raise ValueError("arguments must be scalars.")
2265 if (floor(n) != n):
2266 raise ValueError("n must be an integer.")
2267 if (abs(n) <= 1):
2268 n1 = 1
2269 else:
2270 n1 = n
2271 cpb, cpd = _specfun.cpbdn(n1, z)
2272 return cpb[:n1+1], cpd[:n1+1]
2275def ber_zeros(nt):
2276 """Compute nt zeros of the Kelvin function ber.
2278 Parameters
2279 ----------
2280 nt : int
2281 Number of zeros to compute. Must be positive.
2283 Returns
2284 -------
2285 ndarray
2286 First `nt` zeros of the Kelvin function.
2288 See Also
2289 --------
2290 ber
2292 References
2293 ----------
2294 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
2295 Functions", John Wiley and Sons, 1996.
2296 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
2298 """
2299 if not isscalar(nt) or (floor(nt) != nt) or (nt <= 0):
2300 raise ValueError("nt must be positive integer scalar.")
2301 return _specfun.klvnzo(nt, 1)
2304def bei_zeros(nt):
2305 """Compute nt zeros of the Kelvin function bei.
2307 Parameters
2308 ----------
2309 nt : int
2310 Number of zeros to compute. Must be positive.
2312 Returns
2313 -------
2314 ndarray
2315 First `nt` zeros of the Kelvin function.
2317 See Also
2318 --------
2319 bei
2321 References
2322 ----------
2323 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
2324 Functions", John Wiley and Sons, 1996.
2325 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
2327 """
2328 if not isscalar(nt) or (floor(nt) != nt) or (nt <= 0):
2329 raise ValueError("nt must be positive integer scalar.")
2330 return _specfun.klvnzo(nt, 2)
2333def ker_zeros(nt):
2334 """Compute nt zeros of the Kelvin function ker.
2336 Parameters
2337 ----------
2338 nt : int
2339 Number of zeros to compute. Must be positive.
2341 Returns
2342 -------
2343 ndarray
2344 First `nt` zeros of the Kelvin function.
2346 See Also
2347 --------
2348 ker
2350 References
2351 ----------
2352 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
2353 Functions", John Wiley and Sons, 1996.
2354 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
2356 """
2357 if not isscalar(nt) or (floor(nt) != nt) or (nt <= 0):
2358 raise ValueError("nt must be positive integer scalar.")
2359 return _specfun.klvnzo(nt, 3)
2362def kei_zeros(nt):
2363 """Compute nt zeros of the Kelvin function kei.
2365 Parameters
2366 ----------
2367 nt : int
2368 Number of zeros to compute. Must be positive.
2370 Returns
2371 -------
2372 ndarray
2373 First `nt` zeros of the Kelvin function.
2375 See Also
2376 --------
2377 kei
2379 References
2380 ----------
2381 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
2382 Functions", John Wiley and Sons, 1996.
2383 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
2385 """
2386 if not isscalar(nt) or (floor(nt) != nt) or (nt <= 0):
2387 raise ValueError("nt must be positive integer scalar.")
2388 return _specfun.klvnzo(nt, 4)
2391def berp_zeros(nt):
2392 """Compute nt zeros of the derivative of the Kelvin function ber.
2394 Parameters
2395 ----------
2396 nt : int
2397 Number of zeros to compute. Must be positive.
2399 Returns
2400 -------
2401 ndarray
2402 First `nt` zeros of the derivative of the Kelvin function.
2404 See Also
2405 --------
2406 ber, berp
2408 References
2409 ----------
2410 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
2411 Functions", John Wiley and Sons, 1996.
2412 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
2414 """
2415 if not isscalar(nt) or (floor(nt) != nt) or (nt <= 0):
2416 raise ValueError("nt must be positive integer scalar.")
2417 return _specfun.klvnzo(nt, 5)
2420def beip_zeros(nt):
2421 """Compute nt zeros of the derivative of the Kelvin function bei.
2423 Parameters
2424 ----------
2425 nt : int
2426 Number of zeros to compute. Must be positive.
2428 Returns
2429 -------
2430 ndarray
2431 First `nt` zeros of the derivative of the Kelvin function.
2433 See Also
2434 --------
2435 bei, beip
2437 References
2438 ----------
2439 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
2440 Functions", John Wiley and Sons, 1996.
2441 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
2443 """
2444 if not isscalar(nt) or (floor(nt) != nt) or (nt <= 0):
2445 raise ValueError("nt must be positive integer scalar.")
2446 return _specfun.klvnzo(nt, 6)
2449def kerp_zeros(nt):
2450 """Compute nt zeros of the derivative of the Kelvin function ker.
2452 Parameters
2453 ----------
2454 nt : int
2455 Number of zeros to compute. Must be positive.
2457 Returns
2458 -------
2459 ndarray
2460 First `nt` zeros of the derivative of the Kelvin function.
2462 See Also
2463 --------
2464 ker, kerp
2466 References
2467 ----------
2468 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
2469 Functions", John Wiley and Sons, 1996.
2470 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
2472 """
2473 if not isscalar(nt) or (floor(nt) != nt) or (nt <= 0):
2474 raise ValueError("nt must be positive integer scalar.")
2475 return _specfun.klvnzo(nt, 7)
2478def keip_zeros(nt):
2479 """Compute nt zeros of the derivative of the Kelvin function kei.
2481 Parameters
2482 ----------
2483 nt : int
2484 Number of zeros to compute. Must be positive.
2486 Returns
2487 -------
2488 ndarray
2489 First `nt` zeros of the derivative of the Kelvin function.
2491 See Also
2492 --------
2493 kei, keip
2495 References
2496 ----------
2497 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
2498 Functions", John Wiley and Sons, 1996.
2499 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
2501 """
2502 if not isscalar(nt) or (floor(nt) != nt) or (nt <= 0):
2503 raise ValueError("nt must be positive integer scalar.")
2504 return _specfun.klvnzo(nt, 8)
2507def kelvin_zeros(nt):
2508 """Compute nt zeros of all Kelvin functions.
2510 Returned in a length-8 tuple of arrays of length nt. The tuple contains
2511 the arrays of zeros of (ber, bei, ker, kei, ber', bei', ker', kei').
2513 References
2514 ----------
2515 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
2516 Functions", John Wiley and Sons, 1996.
2517 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
2519 """
2520 if not isscalar(nt) or (floor(nt) != nt) or (nt <= 0):
2521 raise ValueError("nt must be positive integer scalar.")
2522 return (_specfun.klvnzo(nt, 1),
2523 _specfun.klvnzo(nt, 2),
2524 _specfun.klvnzo(nt, 3),
2525 _specfun.klvnzo(nt, 4),
2526 _specfun.klvnzo(nt, 5),
2527 _specfun.klvnzo(nt, 6),
2528 _specfun.klvnzo(nt, 7),
2529 _specfun.klvnzo(nt, 8))
2532def pro_cv_seq(m, n, c):
2533 """Characteristic values for prolate spheroidal wave functions.
2535 Compute a sequence of characteristic values for the prolate
2536 spheroidal wave functions for mode m and n'=m..n and spheroidal
2537 parameter c.
2539 References
2540 ----------
2541 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
2542 Functions", John Wiley and Sons, 1996.
2543 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
2545 """
2546 if not (isscalar(m) and isscalar(n) and isscalar(c)):
2547 raise ValueError("Arguments must be scalars.")
2548 if (n != floor(n)) or (m != floor(m)):
2549 raise ValueError("Modes must be integers.")
2550 if (n-m > 199):
2551 raise ValueError("Difference between n and m is too large.")
2552 maxL = n-m+1
2553 return _specfun.segv(m, n, c, 1)[1][:maxL]
2556def obl_cv_seq(m, n, c):
2557 """Characteristic values for oblate spheroidal wave functions.
2559 Compute a sequence of characteristic values for the oblate
2560 spheroidal wave functions for mode m and n'=m..n and spheroidal
2561 parameter c.
2563 References
2564 ----------
2565 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
2566 Functions", John Wiley and Sons, 1996.
2567 https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html
2569 """
2570 if not (isscalar(m) and isscalar(n) and isscalar(c)):
2571 raise ValueError("Arguments must be scalars.")
2572 if (n != floor(n)) or (m != floor(m)):
2573 raise ValueError("Modes must be integers.")
2574 if (n-m > 199):
2575 raise ValueError("Difference between n and m is too large.")
2576 maxL = n-m+1
2577 return _specfun.segv(m, n, c, -1)[1][:maxL]
2580def comb(N, k, exact=False, repetition=False, legacy=True):
2581 """The number of combinations of N things taken k at a time.
2583 This is often expressed as "N choose k".
2585 Parameters
2586 ----------
2587 N : int, ndarray
2588 Number of things.
2589 k : int, ndarray
2590 Number of elements taken.
2591 exact : bool, optional
2592 For integers, if `exact` is False, then floating point precision is
2593 used, otherwise the result is computed exactly. For non-integers, if
2594 `exact` is True, the inputs are currently cast to integers, though
2595 this behavior is deprecated (see below).
2596 repetition : bool, optional
2597 If `repetition` is True, then the number of combinations with
2598 repetition is computed.
2599 legacy : bool, optional
2600 If `legacy` is True and `exact` is True, then non-integral arguments
2601 are cast to ints; if `legacy` is False, the result for non-integral
2602 arguments is unaffected by the value of `exact`.
2604 .. deprecated:: 1.9.0
2605 Non-integer arguments are currently being cast to integers when
2606 `exact=True`. This behaviour is deprecated and the default will
2607 change to avoid the cast in SciPy 1.11.0. To opt into the future
2608 behavior set `legacy=False`. If you want to keep the
2609 argument-casting but silence this warning, cast your inputs
2610 directly, e.g. ``comb(int(your_N), int(your_k), exact=True)``.
2612 Returns
2613 -------
2614 val : int, float, ndarray
2615 The total number of combinations.
2617 See Also
2618 --------
2619 binom : Binomial coefficient considered as a function of two real
2620 variables.
2622 Notes
2623 -----
2624 - Array arguments accepted only for exact=False case.
2625 - If N < 0, or k < 0, then 0 is returned.
2626 - If k > N and repetition=False, then 0 is returned.
2628 Examples
2629 --------
2630 >>> import numpy as np
2631 >>> from scipy.special import comb
2632 >>> k = np.array([3, 4])
2633 >>> n = np.array([10, 10])
2634 >>> comb(n, k, exact=False)
2635 array([ 120., 210.])
2636 >>> comb(10, 3, exact=True)
2637 120
2638 >>> comb(10, 3, exact=True, repetition=True)
2639 220
2641 """
2642 if repetition:
2643 return comb(N + k - 1, k, exact, legacy=legacy)
2644 if exact:
2645 if int(N) != N or int(k) != k:
2646 if legacy:
2647 warnings.warn(
2648 "Non-integer arguments are currently being cast to "
2649 "integers when exact=True. This behaviour is "
2650 "deprecated and the default will change to avoid the cast "
2651 "in SciPy 1.11.0. To opt into the future behavior set "
2652 "legacy=False. If you want to keep the argument-casting "
2653 "but silence this warning, cast your inputs directly, "
2654 "e.g. comb(int(your_N), int(your_k), exact=True).",
2655 DeprecationWarning, stacklevel=2
2656 )
2657 else:
2658 return comb(N, k)
2659 # _comb_int casts inputs to integers
2660 return _comb_int(N, k)
2661 else:
2662 k, N = asarray(k), asarray(N)
2663 cond = (k <= N) & (N >= 0) & (k >= 0)
2664 vals = binom(N, k)
2665 if isinstance(vals, np.ndarray):
2666 vals[~cond] = 0
2667 elif not cond:
2668 vals = np.float64(0)
2669 return vals
2672def perm(N, k, exact=False):
2673 """Permutations of N things taken k at a time, i.e., k-permutations of N.
2675 It's also known as "partial permutations".
2677 Parameters
2678 ----------
2679 N : int, ndarray
2680 Number of things.
2681 k : int, ndarray
2682 Number of elements taken.
2683 exact : bool, optional
2684 If `exact` is False, then floating point precision is used, otherwise
2685 exact long integer is computed.
2687 Returns
2688 -------
2689 val : int, ndarray
2690 The number of k-permutations of N.
2692 Notes
2693 -----
2694 - Array arguments accepted only for exact=False case.
2695 - If k > N, N < 0, or k < 0, then a 0 is returned.
2697 Examples
2698 --------
2699 >>> import numpy as np
2700 >>> from scipy.special import perm
2701 >>> k = np.array([3, 4])
2702 >>> n = np.array([10, 10])
2703 >>> perm(n, k)
2704 array([ 720., 5040.])
2705 >>> perm(10, 3, exact=True)
2706 720
2708 """
2709 if exact:
2710 if (k > N) or (N < 0) or (k < 0):
2711 return 0
2712 val = 1
2713 for i in range(N - k + 1, N + 1):
2714 val *= i
2715 return val
2716 else:
2717 k, N = asarray(k), asarray(N)
2718 cond = (k <= N) & (N >= 0) & (k >= 0)
2719 vals = poch(N - k + 1, k)
2720 if isinstance(vals, np.ndarray):
2721 vals[~cond] = 0
2722 elif not cond:
2723 vals = np.float64(0)
2724 return vals
2727# https://stackoverflow.com/a/16327037
2728def _range_prod(lo, hi):
2729 """
2730 Product of a range of numbers.
2732 Returns the product of
2733 lo * (lo+1) * (lo+2) * ... * (hi-2) * (hi-1) * hi
2734 = hi! / (lo-1)!
2736 Breaks into smaller products first for speed:
2737 _range_prod(2, 9) = ((2*3)*(4*5))*((6*7)*(8*9))
2738 """
2739 if lo + 1 < hi:
2740 mid = (hi + lo) // 2
2741 return _range_prod(lo, mid) * _range_prod(mid + 1, hi)
2742 if lo == hi:
2743 return lo
2744 return lo * hi
2747def factorial(n, exact=False):
2748 """
2749 The factorial of a number or array of numbers.
2751 The factorial of non-negative integer `n` is the product of all
2752 positive integers less than or equal to `n`::
2754 n! = n * (n - 1) * (n - 2) * ... * 1
2756 Parameters
2757 ----------
2758 n : int or array_like of ints
2759 Input values. If ``n < 0``, the return value is 0.
2760 exact : bool, optional
2761 If True, calculate the answer exactly using long integer arithmetic.
2762 If False, result is approximated in floating point rapidly using the
2763 `gamma` function.
2764 Default is False.
2766 Returns
2767 -------
2768 nf : float or int or ndarray
2769 Factorial of `n`, as integer or float depending on `exact`.
2771 Notes
2772 -----
2773 For arrays with ``exact=True``, the factorial is computed only once, for
2774 the largest input, with each other result computed in the process.
2775 The output dtype is increased to ``int64`` or ``object`` if necessary.
2777 With ``exact=False`` the factorial is approximated using the gamma
2778 function:
2780 .. math:: n! = \\Gamma(n+1)
2782 Examples
2783 --------
2784 >>> import numpy as np
2785 >>> from scipy.special import factorial
2786 >>> arr = np.array([3, 4, 5])
2787 >>> factorial(arr, exact=False)
2788 array([ 6., 24., 120.])
2789 >>> factorial(arr, exact=True)
2790 array([ 6, 24, 120])
2791 >>> factorial(5, exact=True)
2792 120
2794 """
2795 if exact:
2796 if np.ndim(n) == 0:
2797 if np.isnan(n):
2798 return n
2799 return 0 if n < 0 else math.factorial(n)
2800 else:
2801 n = asarray(n)
2802 un = np.unique(n).astype(object)
2804 # Convert to object array of long ints if np.int_ can't handle size
2805 if np.isnan(n).any():
2806 dt = float
2807 elif un[-1] > 20:
2808 dt = object
2809 elif un[-1] > 12:
2810 dt = np.int64
2811 else:
2812 dt = np.int_
2814 out = np.empty_like(n, dtype=dt)
2816 # Handle invalid/trivial values
2817 # Ignore runtime warning when less operator used w/np.nan
2818 with np.errstate(all='ignore'):
2819 un = un[un > 1]
2820 out[n < 2] = 1
2821 out[n < 0] = 0
2823 # Calculate products of each range of numbers
2824 if un.size:
2825 val = math.factorial(un[0])
2826 out[n == un[0]] = val
2827 for i in range(len(un) - 1):
2828 prev = un[i] + 1
2829 current = un[i + 1]
2830 val *= _range_prod(prev, current)
2831 out[n == current] = val
2833 if np.isnan(n).any():
2834 out = out.astype(np.float64)
2835 out[np.isnan(n)] = n[np.isnan(n)]
2836 return out
2837 else:
2838 out = _ufuncs._factorial(n)
2839 return out
2842def factorial2(n, exact=False):
2843 """Double factorial.
2845 This is the factorial with every second value skipped. E.g., ``7!! = 7 * 5
2846 * 3 * 1``. It can be approximated numerically as::
2848 n!! = special.gamma(n/2+1)*2**((m+1)/2)/sqrt(pi) n odd
2849 = 2**(n/2) * (n/2)! n even
2851 Parameters
2852 ----------
2853 n : int or array_like
2854 Calculate ``n!!``. Arrays are only supported with `exact` set
2855 to False. If ``n < -1``, the return value is 0.
2856 Otherwise if ``n <= 0``, the return value is 1.
2857 exact : bool, optional
2858 The result can be approximated rapidly using the gamma-formula
2859 above (default). If `exact` is set to True, calculate the
2860 answer exactly using integer arithmetic.
2862 Returns
2863 -------
2864 nff : float or int
2865 Double factorial of `n`, as an int or a float depending on
2866 `exact`.
2868 Examples
2869 --------
2870 >>> from scipy.special import factorial2
2871 >>> factorial2(7, exact=False)
2872 array(105.00000000000001)
2873 >>> factorial2(7, exact=True)
2874 105
2876 """
2877 if exact:
2878 if n < -1:
2879 return 0
2880 if n <= 0:
2881 return 1
2882 val = 1
2883 for k in range(n, 0, -2):
2884 val *= k
2885 return val
2886 else:
2887 n = asarray(n)
2888 vals = zeros(n.shape, 'd')
2889 cond1 = (n % 2) & (n >= -1)
2890 cond2 = (1-(n % 2)) & (n >= -1)
2891 oddn = extract(cond1, n)
2892 evenn = extract(cond2, n)
2893 nd2o = oddn / 2.0
2894 nd2e = evenn / 2.0
2895 place(vals, cond1, gamma(nd2o + 1) / sqrt(pi) * pow(2.0, nd2o + 0.5))
2896 place(vals, cond2, gamma(nd2e + 1) * pow(2.0, nd2e))
2897 return vals
2900def factorialk(n, k, exact=True):
2901 """Multifactorial of n of order k, n(!!...!).
2903 This is the multifactorial of n skipping k values. For example,
2905 factorialk(17, 4) = 17!!!! = 17 * 13 * 9 * 5 * 1
2907 In particular, for any integer ``n``, we have
2909 factorialk(n, 1) = factorial(n)
2911 factorialk(n, 2) = factorial2(n)
2913 Parameters
2914 ----------
2915 n : int
2916 Calculate multifactorial. If ``n < 1 - k``, the return value is 0.
2917 Otherwise if ``n <= 0``, the return value is 1.
2918 k : int
2919 Order of multifactorial.
2920 exact : bool, optional
2921 If exact is set to True, calculate the answer exactly using
2922 integer arithmetic.
2924 Returns
2925 -------
2926 val : int
2927 Multifactorial of `n`.
2929 Raises
2930 ------
2931 NotImplementedError
2932 Raises when exact is False
2934 Examples
2935 --------
2936 >>> from scipy.special import factorialk
2937 >>> factorialk(5, 1, exact=True)
2938 120
2939 >>> factorialk(5, 3, exact=True)
2940 10
2942 """
2943 if exact:
2944 if n < 1-k:
2945 return 0
2946 if n <= 0:
2947 return 1
2948 val = 1
2949 for j in range(n, 0, -k):
2950 val = val*j
2951 return val
2952 else:
2953 raise NotImplementedError
2956def zeta(x, q=None, out=None):
2957 r"""
2958 Riemann or Hurwitz zeta function.
2960 Parameters
2961 ----------
2962 x : array_like of float
2963 Input data, must be real
2964 q : array_like of float, optional
2965 Input data, must be real. Defaults to Riemann zeta.
2966 out : ndarray, optional
2967 Output array for the computed values.
2969 Returns
2970 -------
2971 out : array_like
2972 Values of zeta(x).
2974 Notes
2975 -----
2976 The two-argument version is the Hurwitz zeta function
2978 .. math::
2980 \zeta(x, q) = \sum_{k=0}^{\infty} \frac{1}{(k + q)^x};
2982 see [dlmf]_ for details. The Riemann zeta function corresponds to
2983 the case when ``q = 1``.
2985 See Also
2986 --------
2987 zetac
2989 References
2990 ----------
2991 .. [dlmf] NIST, Digital Library of Mathematical Functions,
2992 https://dlmf.nist.gov/25.11#i
2994 Examples
2995 --------
2996 >>> import numpy as np
2997 >>> from scipy.special import zeta, polygamma, factorial
2999 Some specific values:
3001 >>> zeta(2), np.pi**2/6
3002 (1.6449340668482266, 1.6449340668482264)
3004 >>> zeta(4), np.pi**4/90
3005 (1.0823232337111381, 1.082323233711138)
3007 Relation to the `polygamma` function:
3009 >>> m = 3
3010 >>> x = 1.25
3011 >>> polygamma(m, x)
3012 array(2.782144009188397)
3013 >>> (-1)**(m+1) * factorial(m) * zeta(m+1, x)
3014 2.7821440091883969
3016 """
3017 if q is None:
3018 return _ufuncs._riemann_zeta(x, out)
3019 else:
3020 return _ufuncs._zeta(x, q, out)