Coverage for /home/martinb/.local/share/virtualenvs/camcops/lib/python3.6/site-packages/scipy/special/_basic.py : 13%

Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1#
2# Author: Travis Oliphant, 2002
3#
5import operator
6import numpy as np
7import math
8from numpy import (pi, asarray, floor, isscalar, iscomplex, real,
9 imag, sqrt, where, mgrid, sin, place, issubdtype,
10 extract, inexact, nan, zeros, sinc)
11from . import _ufuncs as ufuncs
12from ._ufuncs import (mathieu_a, mathieu_b, iv, jv, gamma,
13 psi, hankel1, hankel2, yv, kv, ndtri,
14 poch, binom, hyp0f1)
15from . import specfun
16from . import orthogonal
17from ._comb import _comb_int
20__all__ = [
21 'ai_zeros',
22 'assoc_laguerre',
23 'bei_zeros',
24 'beip_zeros',
25 'ber_zeros',
26 'bernoulli',
27 'berp_zeros',
28 'bi_zeros',
29 'clpmn',
30 'comb',
31 'digamma',
32 'diric',
33 'erf_zeros',
34 'euler',
35 'factorial',
36 'factorial2',
37 'factorialk',
38 'fresnel_zeros',
39 'fresnelc_zeros',
40 'fresnels_zeros',
41 'gamma',
42 'h1vp',
43 'h2vp',
44 'hankel1',
45 'hankel2',
46 'hyp0f1',
47 'iv',
48 'ivp',
49 'jn_zeros',
50 'jnjnp_zeros',
51 'jnp_zeros',
52 'jnyn_zeros',
53 'jv',
54 'jvp',
55 'kei_zeros',
56 'keip_zeros',
57 'kelvin_zeros',
58 'ker_zeros',
59 'kerp_zeros',
60 'kv',
61 'kvp',
62 'lmbda',
63 'lpmn',
64 'lpn',
65 'lqmn',
66 'lqn',
67 'mathieu_a',
68 'mathieu_b',
69 'mathieu_even_coef',
70 'mathieu_odd_coef',
71 'ndtri',
72 'obl_cv_seq',
73 'pbdn_seq',
74 'pbdv_seq',
75 'pbvv_seq',
76 'perm',
77 'polygamma',
78 'pro_cv_seq',
79 'psi',
80 'riccati_jn',
81 'riccati_yn',
82 'sinc',
83 'y0_zeros',
84 'y1_zeros',
85 'y1p_zeros',
86 'yn_zeros',
87 'ynp_zeros',
88 'yv',
89 'yvp',
90 'zeta'
91]
94def _nonneg_int_or_fail(n, var_name, strict=True):
95 try:
96 if strict:
97 # Raises an exception if float
98 n = operator.index(n)
99 elif n == floor(n):
100 n = int(n)
101 else:
102 raise ValueError()
103 if n < 0:
104 raise ValueError()
105 except (ValueError, TypeError) as err:
106 raise err.__class__("{} must be a non-negative integer".format(var_name))
107 return n
110def diric(x, n):
111 """Periodic sinc function, also called the Dirichlet function.
113 The Dirichlet function is defined as::
115 diric(x, n) = sin(x * n/2) / (n * sin(x / 2)),
117 where `n` is a positive integer.
119 Parameters
120 ----------
121 x : array_like
122 Input data
123 n : int
124 Integer defining the periodicity.
126 Returns
127 -------
128 diric : ndarray
130 Examples
131 --------
132 >>> from scipy import special
133 >>> import matplotlib.pyplot as plt
135 >>> x = np.linspace(-8*np.pi, 8*np.pi, num=201)
136 >>> plt.figure(figsize=(8, 8));
137 >>> for idx, n in enumerate([2, 3, 4, 9]):
138 ... plt.subplot(2, 2, idx+1)
139 ... plt.plot(x, special.diric(x, n))
140 ... plt.title('diric, n={}'.format(n))
141 >>> plt.show()
143 The following example demonstrates that `diric` gives the magnitudes
144 (modulo the sign and scaling) of the Fourier coefficients of a
145 rectangular pulse.
147 Suppress output of values that are effectively 0:
149 >>> np.set_printoptions(suppress=True)
151 Create a signal `x` of length `m` with `k` ones:
153 >>> m = 8
154 >>> k = 3
155 >>> x = np.zeros(m)
156 >>> x[:k] = 1
158 Use the FFT to compute the Fourier transform of `x`, and
159 inspect the magnitudes of the coefficients:
161 >>> np.abs(np.fft.fft(x))
162 array([ 3. , 2.41421356, 1. , 0.41421356, 1. ,
163 0.41421356, 1. , 2.41421356])
165 Now find the same values (up to sign) using `diric`. We multiply
166 by `k` to account for the different scaling conventions of
167 `numpy.fft.fft` and `diric`:
169 >>> theta = np.linspace(0, 2*np.pi, m, endpoint=False)
170 >>> k * special.diric(theta, k)
171 array([ 3. , 2.41421356, 1. , -0.41421356, -1. ,
172 -0.41421356, 1. , 2.41421356])
173 """
174 x, n = asarray(x), asarray(n)
175 n = asarray(n + (x-x))
176 x = asarray(x + (n-n))
177 if issubdtype(x.dtype, inexact):
178 ytype = x.dtype
179 else:
180 ytype = float
181 y = zeros(x.shape, ytype)
183 # empirical minval for 32, 64 or 128 bit float computations
184 # where sin(x/2) < minval, result is fixed at +1 or -1
185 if np.finfo(ytype).eps < 1e-18:
186 minval = 1e-11
187 elif np.finfo(ytype).eps < 1e-15:
188 minval = 1e-7
189 else:
190 minval = 1e-3
192 mask1 = (n <= 0) | (n != floor(n))
193 place(y, mask1, nan)
195 x = x / 2
196 denom = sin(x)
197 mask2 = (1-mask1) & (abs(denom) < minval)
198 xsub = extract(mask2, x)
199 nsub = extract(mask2, n)
200 zsub = xsub / pi
201 place(y, mask2, pow(-1, np.round(zsub)*(nsub-1)))
203 mask = (1-mask1) & (1-mask2)
204 xsub = extract(mask, x)
205 nsub = extract(mask, n)
206 dsub = extract(mask, denom)
207 place(y, mask, sin(nsub*xsub)/(nsub*dsub))
208 return y
211def jnjnp_zeros(nt):
212 """Compute zeros of integer-order Bessel functions Jn and Jn'.
214 Results are arranged in order of the magnitudes of the zeros.
216 Parameters
217 ----------
218 nt : int
219 Number (<=1200) of zeros to compute
221 Returns
222 -------
223 zo[l-1] : ndarray
224 Value of the lth zero of Jn(x) and Jn'(x). Of length `nt`.
225 n[l-1] : ndarray
226 Order of the Jn(x) or Jn'(x) associated with lth zero. Of length `nt`.
227 m[l-1] : ndarray
228 Serial number of the zeros of Jn(x) or Jn'(x) associated
229 with lth zero. Of length `nt`.
230 t[l-1] : ndarray
231 0 if lth zero in zo is zero of Jn(x), 1 if it is a zero of Jn'(x). Of
232 length `nt`.
234 See Also
235 --------
236 jn_zeros, jnp_zeros : to get separated arrays of zeros.
238 References
239 ----------
240 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
241 Functions", John Wiley and Sons, 1996, chapter 5.
242 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html
244 """
245 if not isscalar(nt) or (floor(nt) != nt) or (nt > 1200):
246 raise ValueError("Number must be integer <= 1200.")
247 nt = int(nt)
248 n, m, t, zo = specfun.jdzo(nt)
249 return zo[1:nt+1], n[:nt], m[:nt], t[:nt]
252def jnyn_zeros(n, nt):
253 """Compute nt zeros of Bessel functions Jn(x), Jn'(x), Yn(x), and Yn'(x).
255 Returns 4 arrays of length `nt`, corresponding to the first `nt`
256 zeros of Jn(x), Jn'(x), Yn(x), and Yn'(x), respectively. The zeros
257 are returned in ascending order.
259 Parameters
260 ----------
261 n : int
262 Order of the Bessel functions
263 nt : int
264 Number (<=1200) of zeros to compute
266 Returns
267 -------
268 Jn : ndarray
269 First `nt` zeros of Jn
270 Jnp : ndarray
271 First `nt` zeros of Jn'
272 Yn : ndarray
273 First `nt` zeros of Yn
274 Ynp : ndarray
275 First `nt` zeros of Yn'
277 See Also
278 --------
279 jn_zeros, jnp_zeros, yn_zeros, ynp_zeros
281 References
282 ----------
283 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
284 Functions", John Wiley and Sons, 1996, chapter 5.
285 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html
287 """
288 if not (isscalar(nt) and isscalar(n)):
289 raise ValueError("Arguments must be scalars.")
290 if (floor(n) != n) or (floor(nt) != nt):
291 raise ValueError("Arguments must be integers.")
292 if (nt <= 0):
293 raise ValueError("nt > 0")
294 return specfun.jyzo(abs(n), nt)
297def jn_zeros(n, nt):
298 r"""Compute zeros of integer-order Bessel functions Jn.
300 Compute `nt` zeros of the Bessel functions :math:`J_n(x)` on the
301 interval :math:`(0, \infty)`. The zeros are returned in ascending
302 order. Note that this interval excludes the zero at :math:`x = 0`
303 that exists for :math:`n > 0`.
305 Parameters
306 ----------
307 n : int
308 Order of Bessel function
309 nt : int
310 Number of zeros to return
312 Returns
313 -------
314 ndarray
315 First `n` zeros of the Bessel function.
317 See Also
318 --------
319 jv
321 References
322 ----------
323 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
324 Functions", John Wiley and Sons, 1996, chapter 5.
325 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html
327 Examples
328 --------
329 >>> import scipy.special as sc
331 We can check that we are getting approximations of the zeros by
332 evaluating them with `jv`.
334 >>> n = 1
335 >>> x = sc.jn_zeros(n, 3)
336 >>> x
337 array([ 3.83170597, 7.01558667, 10.17346814])
338 >>> sc.jv(n, x)
339 array([-0.00000000e+00, 1.72975330e-16, 2.89157291e-16])
341 Note that the zero at ``x = 0`` for ``n > 0`` is not included.
343 >>> sc.jv(1, 0)
344 0.0
346 """
347 return jnyn_zeros(n, nt)[0]
350def jnp_zeros(n, nt):
351 r"""Compute zeros of integer-order Bessel function derivatives Jn'.
353 Compute `nt` zeros of the functions :math:`J_n'(x)` on the
354 interval :math:`(0, \infty)`. The zeros are returned in ascending
355 order. Note that this interval excludes the zero at :math:`x = 0`
356 that exists for :math:`n > 1`.
358 Parameters
359 ----------
360 n : int
361 Order of Bessel function
362 nt : int
363 Number of zeros to return
365 Returns
366 -------
367 ndarray
368 First `n` zeros of the Bessel function.
370 See Also
371 --------
372 jvp, jv
374 References
375 ----------
376 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
377 Functions", John Wiley and Sons, 1996, chapter 5.
378 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html
380 Examples
381 --------
382 >>> import scipy.special as sc
384 We can check that we are getting approximations of the zeros by
385 evaluating them with `jvp`.
387 >>> n = 2
388 >>> x = sc.jnp_zeros(n, 3)
389 >>> x
390 array([3.05423693, 6.70613319, 9.96946782])
391 >>> sc.jvp(n, x)
392 array([ 2.77555756e-17, 2.08166817e-16, -3.01841885e-16])
394 Note that the zero at ``x = 0`` for ``n > 1`` is not included.
396 >>> sc.jvp(n, 0)
397 0.0
399 """
400 return jnyn_zeros(n, nt)[1]
403def yn_zeros(n, nt):
404 r"""Compute zeros of integer-order Bessel function Yn(x).
406 Compute `nt` zeros of the functions :math:`Y_n(x)` on the interval
407 :math:`(0, \infty)`. The zeros are returned in ascending order.
409 Parameters
410 ----------
411 n : int
412 Order of Bessel function
413 nt : int
414 Number of zeros to return
416 Returns
417 -------
418 ndarray
419 First `n` zeros of the Bessel function.
421 See Also
422 --------
423 yn, yv
425 References
426 ----------
427 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
428 Functions", John Wiley and Sons, 1996, chapter 5.
429 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html
431 Examples
432 --------
433 >>> import scipy.special as sc
435 We can check that we are getting approximations of the zeros by
436 evaluating them with `yn`.
438 >>> n = 2
439 >>> x = sc.yn_zeros(n, 3)
440 >>> x
441 array([ 3.38424177, 6.79380751, 10.02347798])
442 >>> sc.yn(n, x)
443 array([-1.94289029e-16, 8.32667268e-17, -1.52655666e-16])
445 """
446 return jnyn_zeros(n, nt)[2]
449def ynp_zeros(n, nt):
450 r"""Compute zeros of integer-order Bessel function derivatives Yn'(x).
452 Compute `nt` zeros of the functions :math:`Y_n'(x)` on the
453 interval :math:`(0, \infty)`. The zeros are returned in ascending
454 order.
456 Parameters
457 ----------
458 n : int
459 Order of Bessel function
460 nt : int
461 Number of zeros to return
463 See Also
464 --------
465 yvp
467 References
468 ----------
469 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
470 Functions", John Wiley and Sons, 1996, chapter 5.
471 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html
473 Examples
474 --------
475 >>> import scipy.special as sc
477 We can check that we are getting approximations of the zeros by
478 evaluating them with `yvp`.
480 >>> n = 2
481 >>> x = sc.ynp_zeros(n, 3)
482 >>> x
483 array([ 5.00258293, 8.3507247 , 11.57419547])
484 >>> sc.yvp(n, x)
485 array([ 2.22044605e-16, -3.33066907e-16, 2.94902991e-16])
487 """
488 return jnyn_zeros(n, nt)[3]
491def y0_zeros(nt, complex=False):
492 """Compute nt zeros of Bessel function Y0(z), and derivative at each zero.
494 The derivatives are given by Y0'(z0) = -Y1(z0) at each zero z0.
496 Parameters
497 ----------
498 nt : int
499 Number of zeros to return
500 complex : bool, default False
501 Set to False to return only the real zeros; set to True to return only
502 the complex zeros with negative real part and positive imaginary part.
503 Note that the complex conjugates of the latter are also zeros of the
504 function, but are not returned by this routine.
506 Returns
507 -------
508 z0n : ndarray
509 Location of nth zero of Y0(z)
510 y0pz0n : ndarray
511 Value of derivative Y0'(z0) for nth zero
513 References
514 ----------
515 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
516 Functions", John Wiley and Sons, 1996, chapter 5.
517 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html
519 """
520 if not isscalar(nt) or (floor(nt) != nt) or (nt <= 0):
521 raise ValueError("Arguments must be scalar positive integer.")
522 kf = 0
523 kc = not complex
524 return specfun.cyzo(nt, kf, kc)
527def y1_zeros(nt, complex=False):
528 """Compute nt zeros of Bessel function Y1(z), and derivative at each zero.
530 The derivatives are given by Y1'(z1) = Y0(z1) at each zero z1.
532 Parameters
533 ----------
534 nt : int
535 Number of zeros to return
536 complex : bool, default False
537 Set to False to return only the real zeros; set to True to return only
538 the complex zeros with negative real part and positive imaginary part.
539 Note that the complex conjugates of the latter are also zeros of the
540 function, but are not returned by this routine.
542 Returns
543 -------
544 z1n : ndarray
545 Location of nth zero of Y1(z)
546 y1pz1n : ndarray
547 Value of derivative Y1'(z1) for nth zero
549 References
550 ----------
551 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
552 Functions", John Wiley and Sons, 1996, chapter 5.
553 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html
555 """
556 if not isscalar(nt) or (floor(nt) != nt) or (nt <= 0):
557 raise ValueError("Arguments must be scalar positive integer.")
558 kf = 1
559 kc = not complex
560 return specfun.cyzo(nt, kf, kc)
563def y1p_zeros(nt, complex=False):
564 """Compute nt zeros of Bessel derivative Y1'(z), and value at each zero.
566 The values are given by Y1(z1) at each z1 where Y1'(z1)=0.
568 Parameters
569 ----------
570 nt : int
571 Number of zeros to return
572 complex : bool, default False
573 Set to False to return only the real zeros; set to True to return only
574 the complex zeros with negative real part and positive imaginary part.
575 Note that the complex conjugates of the latter are also zeros of the
576 function, but are not returned by this routine.
578 Returns
579 -------
580 z1pn : ndarray
581 Location of nth zero of Y1'(z)
582 y1z1pn : ndarray
583 Value of derivative Y1(z1) for nth zero
585 References
586 ----------
587 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
588 Functions", John Wiley and Sons, 1996, chapter 5.
589 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html
591 """
592 if not isscalar(nt) or (floor(nt) != nt) or (nt <= 0):
593 raise ValueError("Arguments must be scalar positive integer.")
594 kf = 2
595 kc = not complex
596 return specfun.cyzo(nt, kf, kc)
599def _bessel_diff_formula(v, z, n, L, phase):
600 # from AMS55.
601 # L(v, z) = J(v, z), Y(v, z), H1(v, z), H2(v, z), phase = -1
602 # L(v, z) = I(v, z) or exp(v*pi*i)K(v, z), phase = 1
603 # For K, you can pull out the exp((v-k)*pi*i) into the caller
604 v = asarray(v)
605 p = 1.0
606 s = L(v-n, z)
607 for i in range(1, n+1):
608 p = phase * (p * (n-i+1)) / i # = choose(k, i)
609 s += p*L(v-n + i*2, z)
610 return s / (2.**n)
613def jvp(v, z, n=1):
614 """Compute derivatives of Bessel functions of the first kind.
616 Compute the nth derivative of the Bessel function `Jv` with
617 respect to `z`.
619 Parameters
620 ----------
621 v : float
622 Order of Bessel function
623 z : complex
624 Argument at which to evaluate the derivative; can be real or
625 complex.
626 n : int, default 1
627 Order of derivative
629 Returns
630 -------
631 scalar or ndarray
632 Values of the derivative of the Bessel function.
634 Notes
635 -----
636 The derivative is computed using the relation DLFM 10.6.7 [2]_.
638 References
639 ----------
640 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
641 Functions", John Wiley and Sons, 1996, chapter 5.
642 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html
643 .. [2] NIST Digital Library of Mathematical Functions.
644 https://dlmf.nist.gov/10.6.E7
646 """
647 n = _nonneg_int_or_fail(n, 'n')
648 if n == 0:
649 return jv(v, z)
650 else:
651 return _bessel_diff_formula(v, z, n, jv, -1)
654def yvp(v, z, n=1):
655 """Compute derivatives of Bessel functions of the second kind.
657 Compute the nth derivative of the Bessel function `Yv` with
658 respect to `z`.
660 Parameters
661 ----------
662 v : float
663 Order of Bessel function
664 z : complex
665 Argument at which to evaluate the derivative
666 n : int, default 1
667 Order of derivative
669 Returns
670 -------
671 scalar or ndarray
672 nth derivative of the Bessel function.
674 Notes
675 -----
676 The derivative is computed using the relation DLFM 10.6.7 [2]_.
678 References
679 ----------
680 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
681 Functions", John Wiley and Sons, 1996, chapter 5.
682 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html
683 .. [2] NIST Digital Library of Mathematical Functions.
684 https://dlmf.nist.gov/10.6.E7
686 """
687 n = _nonneg_int_or_fail(n, 'n')
688 if n == 0:
689 return yv(v, z)
690 else:
691 return _bessel_diff_formula(v, z, n, yv, -1)
694def kvp(v, z, n=1):
695 """Compute nth derivative of real-order modified Bessel function Kv(z)
697 Kv(z) is the modified Bessel function of the second kind.
698 Derivative is calculated with respect to `z`.
700 Parameters
701 ----------
702 v : array_like of float
703 Order of Bessel function
704 z : array_like of complex
705 Argument at which to evaluate the derivative
706 n : int
707 Order of derivative. Default is first derivative.
709 Returns
710 -------
711 out : ndarray
712 The results
714 Examples
715 --------
716 Calculate multiple values at order 5:
718 >>> from scipy.special import kvp
719 >>> kvp(5, (1, 2, 3+5j))
720 array([-1.84903536e+03+0.j , -2.57735387e+01+0.j ,
721 -3.06627741e-02+0.08750845j])
724 Calculate for a single value at multiple orders:
726 >>> kvp((4, 4.5, 5), 1)
727 array([ -184.0309, -568.9585, -1849.0354])
729 Notes
730 -----
731 The derivative is computed using the relation DLFM 10.29.5 [2]_.
733 References
734 ----------
735 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
736 Functions", John Wiley and Sons, 1996, chapter 6.
737 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html
738 .. [2] NIST Digital Library of Mathematical Functions.
739 https://dlmf.nist.gov/10.29.E5
741 """
742 n = _nonneg_int_or_fail(n, 'n')
743 if n == 0:
744 return kv(v, z)
745 else:
746 return (-1)**n * _bessel_diff_formula(v, z, n, kv, 1)
749def ivp(v, z, n=1):
750 """Compute derivatives of modified Bessel functions of the first kind.
752 Compute the nth derivative of the modified Bessel function `Iv`
753 with respect to `z`.
755 Parameters
756 ----------
757 v : array_like
758 Order of Bessel function
759 z : array_like
760 Argument at which to evaluate the derivative; can be real or
761 complex.
762 n : int, default 1
763 Order of derivative
765 Returns
766 -------
767 scalar or ndarray
768 nth derivative of the modified Bessel function.
770 See Also
771 --------
772 iv
774 Notes
775 -----
776 The derivative is computed using the relation DLFM 10.29.5 [2]_.
778 References
779 ----------
780 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
781 Functions", John Wiley and Sons, 1996, chapter 6.
782 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html
783 .. [2] NIST Digital Library of Mathematical Functions.
784 https://dlmf.nist.gov/10.29.E5
786 """
787 n = _nonneg_int_or_fail(n, 'n')
788 if n == 0:
789 return iv(v, z)
790 else:
791 return _bessel_diff_formula(v, z, n, iv, 1)
794def h1vp(v, z, n=1):
795 """Compute nth derivative of Hankel function H1v(z) with respect to `z`.
797 Parameters
798 ----------
799 v : array_like
800 Order of Hankel function
801 z : array_like
802 Argument at which to evaluate the derivative. Can be real or
803 complex.
804 n : int, default 1
805 Order of derivative
807 Returns
808 -------
809 scalar or ndarray
810 Values of the derivative of the Hankel function.
812 Notes
813 -----
814 The derivative is computed using the relation DLFM 10.6.7 [2]_.
816 References
817 ----------
818 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
819 Functions", John Wiley and Sons, 1996, chapter 5.
820 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html
821 .. [2] NIST Digital Library of Mathematical Functions.
822 https://dlmf.nist.gov/10.6.E7
824 """
825 n = _nonneg_int_or_fail(n, 'n')
826 if n == 0:
827 return hankel1(v, z)
828 else:
829 return _bessel_diff_formula(v, z, n, hankel1, -1)
832def h2vp(v, z, n=1):
833 """Compute nth derivative of Hankel function H2v(z) with respect to `z`.
835 Parameters
836 ----------
837 v : array_like
838 Order of Hankel function
839 z : array_like
840 Argument at which to evaluate the derivative. Can be real or
841 complex.
842 n : int, default 1
843 Order of derivative
845 Returns
846 -------
847 scalar or ndarray
848 Values of the derivative of the Hankel 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/f_src/special_functions/special_functions.html
859 .. [2] NIST Digital Library of Mathematical Functions.
860 https://dlmf.nist.gov/10.6.E7
862 """
863 n = _nonneg_int_or_fail(n, 'n')
864 if n == 0:
865 return hankel2(v, z)
866 else:
867 return _bessel_diff_formula(v, z, n, hankel2, -1)
870def riccati_jn(n, x):
871 r"""Compute Ricatti-Bessel function of the first kind and its derivative.
873 The Ricatti-Bessel function of the first kind is defined as :math:`x
874 j_n(x)`, where :math:`j_n` is the spherical Bessel function of the first
875 kind of order :math:`n`.
877 This function computes the value and first derivative of the
878 Ricatti-Bessel function for all orders up to and including `n`.
880 Parameters
881 ----------
882 n : int
883 Maximum order of function to compute
884 x : float
885 Argument at which to evaluate
887 Returns
888 -------
889 jn : ndarray
890 Value of j0(x), ..., jn(x)
891 jnp : ndarray
892 First derivative j0'(x), ..., jn'(x)
894 Notes
895 -----
896 The computation is carried out via backward recurrence, using the
897 relation DLMF 10.51.1 [2]_.
899 Wrapper for a Fortran routine created by Shanjie Zhang and Jianming
900 Jin [1]_.
902 References
903 ----------
904 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
905 Functions", John Wiley and Sons, 1996.
906 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html
907 .. [2] NIST Digital Library of Mathematical Functions.
908 https://dlmf.nist.gov/10.51.E1
910 """
911 if not (isscalar(n) and isscalar(x)):
912 raise ValueError("arguments must be scalars.")
913 n = _nonneg_int_or_fail(n, 'n', strict=False)
914 if (n == 0):
915 n1 = 1
916 else:
917 n1 = n
918 nm, jn, jnp = specfun.rctj(n1, x)
919 return jn[:(n+1)], jnp[:(n+1)]
922def riccati_yn(n, x):
923 """Compute Ricatti-Bessel function of the second kind and its derivative.
925 The Ricatti-Bessel function of the second kind is defined as :math:`x
926 y_n(x)`, where :math:`y_n` is the spherical Bessel function of the second
927 kind of order :math:`n`.
929 This function computes the value and first derivative of the function for
930 all orders up to and including `n`.
932 Parameters
933 ----------
934 n : int
935 Maximum order of function to compute
936 x : float
937 Argument at which to evaluate
939 Returns
940 -------
941 yn : ndarray
942 Value of y0(x), ..., yn(x)
943 ynp : ndarray
944 First derivative y0'(x), ..., yn'(x)
946 Notes
947 -----
948 The computation is carried out via ascending recurrence, using the
949 relation DLMF 10.51.1 [2]_.
951 Wrapper for a Fortran routine created by Shanjie Zhang and Jianming
952 Jin [1]_.
954 References
955 ----------
956 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
957 Functions", John Wiley and Sons, 1996.
958 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html
959 .. [2] NIST Digital Library of Mathematical Functions.
960 https://dlmf.nist.gov/10.51.E1
962 """
963 if not (isscalar(n) and isscalar(x)):
964 raise ValueError("arguments must be scalars.")
965 n = _nonneg_int_or_fail(n, 'n', strict=False)
966 if (n == 0):
967 n1 = 1
968 else:
969 n1 = n
970 nm, jn, jnp = specfun.rcty(n1, x)
971 return jn[:(n+1)], jnp[:(n+1)]
974def erf_zeros(nt):
975 """Compute the first nt zero in the first quadrant, ordered by absolute value.
977 Zeros in the other quadrants can be obtained by using the symmetries erf(-z) = erf(z) and
978 erf(conj(z)) = conj(erf(z)).
981 Parameters
982 ----------
983 nt : int
984 The number of zeros to compute
986 Returns
987 -------
988 The locations of the zeros of erf : ndarray (complex)
989 Complex values at which zeros of erf(z)
991 Examples
992 --------
993 >>> from scipy import special
994 >>> special.erf_zeros(1)
995 array([1.45061616+1.880943j])
997 Check that erf is (close to) zero for the value returned by erf_zeros
999 >>> special.erf(special.erf_zeros(1))
1000 array([4.95159469e-14-1.16407394e-16j])
1002 References
1003 ----------
1004 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
1005 Functions", John Wiley and Sons, 1996.
1006 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html
1008 """
1009 if (floor(nt) != nt) or (nt <= 0) or not isscalar(nt):
1010 raise ValueError("Argument must be positive scalar integer.")
1011 return specfun.cerzo(nt)
1014def fresnelc_zeros(nt):
1015 """Compute nt complex zeros of cosine Fresnel integral C(z).
1017 References
1018 ----------
1019 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
1020 Functions", John Wiley and Sons, 1996.
1021 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html
1023 """
1024 if (floor(nt) != nt) or (nt <= 0) or not isscalar(nt):
1025 raise ValueError("Argument must be positive scalar integer.")
1026 return specfun.fcszo(1, nt)
1029def fresnels_zeros(nt):
1030 """Compute nt complex zeros of sine Fresnel integral S(z).
1032 References
1033 ----------
1034 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
1035 Functions", John Wiley and Sons, 1996.
1036 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html
1038 """
1039 if (floor(nt) != nt) or (nt <= 0) or not isscalar(nt):
1040 raise ValueError("Argument must be positive scalar integer.")
1041 return specfun.fcszo(2, nt)
1044def fresnel_zeros(nt):
1045 """Compute nt complex zeros of sine and cosine Fresnel integrals S(z) and C(z).
1047 References
1048 ----------
1049 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
1050 Functions", John Wiley and Sons, 1996.
1051 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html
1053 """
1054 if (floor(nt) != nt) or (nt <= 0) or not isscalar(nt):
1055 raise ValueError("Argument must be positive scalar integer.")
1056 return specfun.fcszo(2, nt), specfun.fcszo(1, nt)
1059def assoc_laguerre(x, n, k=0.0):
1060 """Compute the generalized (associated) Laguerre polynomial of degree n and order k.
1062 The polynomial :math:`L^{(k)}_n(x)` is orthogonal over ``[0, inf)``,
1063 with weighting function ``exp(-x) * x**k`` with ``k > -1``.
1065 Notes
1066 -----
1067 `assoc_laguerre` is a simple wrapper around `eval_genlaguerre`, with
1068 reversed argument order ``(x, n, k=0.0) --> (n, k, x)``.
1070 """
1071 return orthogonal.eval_genlaguerre(n, k, x)
1074digamma = psi
1077def polygamma(n, x):
1078 r"""Polygamma functions.
1080 Defined as :math:`\psi^{(n)}(x)` where :math:`\psi` is the
1081 `digamma` function. See [dlmf]_ for details.
1083 Parameters
1084 ----------
1085 n : array_like
1086 The order of the derivative of the digamma function; must be
1087 integral
1088 x : array_like
1089 Real valued input
1091 Returns
1092 -------
1093 ndarray
1094 Function results
1096 See Also
1097 --------
1098 digamma
1100 References
1101 ----------
1102 .. [dlmf] NIST, Digital Library of Mathematical Functions,
1103 https://dlmf.nist.gov/5.15
1105 Examples
1106 --------
1107 >>> from scipy import special
1108 >>> x = [2, 3, 25.5]
1109 >>> special.polygamma(1, x)
1110 array([ 0.64493407, 0.39493407, 0.03999467])
1111 >>> special.polygamma(0, x) == special.psi(x)
1112 array([ True, True, True], dtype=bool)
1114 """
1115 n, x = asarray(n), asarray(x)
1116 fac2 = (-1.0)**(n+1) * gamma(n+1.0) * zeta(n+1, x)
1117 return where(n == 0, psi(x), fac2)
1120def mathieu_even_coef(m, q):
1121 r"""Fourier coefficients for even Mathieu and modified Mathieu functions.
1123 The Fourier series of the even solutions of the Mathieu differential
1124 equation are of the form
1126 .. math:: \mathrm{ce}_{2n}(z, q) = \sum_{k=0}^{\infty} A_{(2n)}^{(2k)} \cos 2kz
1128 .. math:: \mathrm{ce}_{2n+1}(z, q) = \sum_{k=0}^{\infty} A_{(2n+1)}^{(2k+1)} \cos (2k+1)z
1130 This function returns the coefficients :math:`A_{(2n)}^{(2k)}` for even
1131 input m=2n, and the coefficients :math:`A_{(2n+1)}^{(2k+1)}` for odd input
1132 m=2n+1.
1134 Parameters
1135 ----------
1136 m : int
1137 Order of Mathieu functions. Must be non-negative.
1138 q : float (>=0)
1139 Parameter of Mathieu functions. Must be non-negative.
1141 Returns
1142 -------
1143 Ak : ndarray
1144 Even or odd Fourier coefficients, corresponding to even or odd m.
1146 References
1147 ----------
1148 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
1149 Functions", John Wiley and Sons, 1996.
1150 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html
1151 .. [2] NIST Digital Library of Mathematical Functions
1152 https://dlmf.nist.gov/28.4#i
1154 """
1155 if not (isscalar(m) and isscalar(q)):
1156 raise ValueError("m and q must be scalars.")
1157 if (q < 0):
1158 raise ValueError("q >=0")
1159 if (m != floor(m)) or (m < 0):
1160 raise ValueError("m must be an integer >=0.")
1162 if (q <= 1):
1163 qm = 7.5 + 56.1*sqrt(q) - 134.7*q + 90.7*sqrt(q)*q
1164 else:
1165 qm = 17.0 + 3.1*sqrt(q) - .126*q + .0037*sqrt(q)*q
1166 km = int(qm + 0.5*m)
1167 if km > 251:
1168 print("Warning, too many predicted coefficients.")
1169 kd = 1
1170 m = int(floor(m))
1171 if m % 2:
1172 kd = 2
1174 a = mathieu_a(m, q)
1175 fc = specfun.fcoef(kd, m, q, a)
1176 return fc[:km]
1179def mathieu_odd_coef(m, q):
1180 r"""Fourier coefficients for even Mathieu and modified Mathieu functions.
1182 The Fourier series of the odd solutions of the Mathieu differential
1183 equation are of the form
1185 .. math:: \mathrm{se}_{2n+1}(z, q) = \sum_{k=0}^{\infty} B_{(2n+1)}^{(2k+1)} \sin (2k+1)z
1187 .. math:: \mathrm{se}_{2n+2}(z, q) = \sum_{k=0}^{\infty} B_{(2n+2)}^{(2k+2)} \sin (2k+2)z
1189 This function returns the coefficients :math:`B_{(2n+2)}^{(2k+2)}` for even
1190 input m=2n+2, and the coefficients :math:`B_{(2n+1)}^{(2k+1)}` for odd
1191 input m=2n+1.
1193 Parameters
1194 ----------
1195 m : int
1196 Order of Mathieu functions. Must be non-negative.
1197 q : float (>=0)
1198 Parameter of Mathieu functions. Must be non-negative.
1200 Returns
1201 -------
1202 Bk : ndarray
1203 Even or odd Fourier coefficients, corresponding to even or odd m.
1205 References
1206 ----------
1207 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
1208 Functions", John Wiley and Sons, 1996.
1209 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html
1211 """
1212 if not (isscalar(m) and isscalar(q)):
1213 raise ValueError("m and q must be scalars.")
1214 if (q < 0):
1215 raise ValueError("q >=0")
1216 if (m != floor(m)) or (m <= 0):
1217 raise ValueError("m must be an integer > 0")
1219 if (q <= 1):
1220 qm = 7.5 + 56.1*sqrt(q) - 134.7*q + 90.7*sqrt(q)*q
1221 else:
1222 qm = 17.0 + 3.1*sqrt(q) - .126*q + .0037*sqrt(q)*q
1223 km = int(qm + 0.5*m)
1224 if km > 251:
1225 print("Warning, too many predicted coefficients.")
1226 kd = 4
1227 m = int(floor(m))
1228 if m % 2:
1229 kd = 3
1231 b = mathieu_b(m, q)
1232 fc = specfun.fcoef(kd, m, q, b)
1233 return fc[:km]
1236def lpmn(m, n, z):
1237 """Sequence of associated Legendre functions of the first kind.
1239 Computes the associated Legendre function of the first kind of order m and
1240 degree n, ``Pmn(z)`` = :math:`P_n^m(z)`, and its derivative, ``Pmn'(z)``.
1241 Returns two arrays of size ``(m+1, n+1)`` containing ``Pmn(z)`` and
1242 ``Pmn'(z)`` for all orders from ``0..m`` and degrees from ``0..n``.
1244 This function takes a real argument ``z``. For complex arguments ``z``
1245 use clpmn instead.
1247 Parameters
1248 ----------
1249 m : int
1250 ``|m| <= n``; the order of the Legendre function.
1251 n : int
1252 where ``n >= 0``; the degree of the Legendre function. Often
1253 called ``l`` (lower case L) in descriptions of the associated
1254 Legendre function
1255 z : float
1256 Input value.
1258 Returns
1259 -------
1260 Pmn_z : (m+1, n+1) array
1261 Values for all orders 0..m and degrees 0..n
1262 Pmn_d_z : (m+1, n+1) array
1263 Derivatives for all orders 0..m and degrees 0..n
1265 See Also
1266 --------
1267 clpmn: associated Legendre functions of the first kind for complex z
1269 Notes
1270 -----
1271 In the interval (-1, 1), Ferrer's function of the first kind is
1272 returned. The phase convention used for the intervals (1, inf)
1273 and (-inf, -1) is such that the result is always real.
1275 References
1276 ----------
1277 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
1278 Functions", John Wiley and Sons, 1996.
1279 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html
1280 .. [2] NIST Digital Library of Mathematical Functions
1281 https://dlmf.nist.gov/14.3
1283 """
1284 if not isscalar(m) or (abs(m) > n):
1285 raise ValueError("m must be <= n.")
1286 if not isscalar(n) or (n < 0):
1287 raise ValueError("n must be a non-negative integer.")
1288 if not isscalar(z):
1289 raise ValueError("z must be scalar.")
1290 if iscomplex(z):
1291 raise ValueError("Argument must be real. Use clpmn instead.")
1292 if (m < 0):
1293 mp = -m
1294 mf, nf = mgrid[0:mp+1, 0:n+1]
1295 with ufuncs.errstate(all='ignore'):
1296 if abs(z) < 1:
1297 # Ferrer function; DLMF 14.9.3
1298 fixarr = where(mf > nf, 0.0,
1299 (-1)**mf * gamma(nf-mf+1) / gamma(nf+mf+1))
1300 else:
1301 # Match to clpmn; DLMF 14.9.13
1302 fixarr = where(mf > nf, 0.0, gamma(nf-mf+1) / gamma(nf+mf+1))
1303 else:
1304 mp = m
1305 p, pd = specfun.lpmn(mp, n, z)
1306 if (m < 0):
1307 p = p * fixarr
1308 pd = pd * fixarr
1309 return p, pd
1312def clpmn(m, n, z, type=3):
1313 """Associated Legendre function of the first kind for complex arguments.
1315 Computes the associated Legendre function of the first kind of order m and
1316 degree n, ``Pmn(z)`` = :math:`P_n^m(z)`, and its derivative, ``Pmn'(z)``.
1317 Returns two arrays of size ``(m+1, n+1)`` containing ``Pmn(z)`` and
1318 ``Pmn'(z)`` for all orders from ``0..m`` and degrees from ``0..n``.
1320 Parameters
1321 ----------
1322 m : int
1323 ``|m| <= n``; the order of the Legendre function.
1324 n : int
1325 where ``n >= 0``; the degree of the Legendre function. Often
1326 called ``l`` (lower case L) in descriptions of the associated
1327 Legendre function
1328 z : float or complex
1329 Input value.
1330 type : int, optional
1331 takes values 2 or 3
1332 2: cut on the real axis ``|x| > 1``
1333 3: cut on the real axis ``-1 < x < 1`` (default)
1335 Returns
1336 -------
1337 Pmn_z : (m+1, n+1) array
1338 Values for all orders ``0..m`` and degrees ``0..n``
1339 Pmn_d_z : (m+1, n+1) array
1340 Derivatives for all orders ``0..m`` and degrees ``0..n``
1342 See Also
1343 --------
1344 lpmn: associated Legendre functions of the first kind for real z
1346 Notes
1347 -----
1348 By default, i.e. for ``type=3``, phase conventions are chosen according
1349 to [1]_ such that the function is analytic. The cut lies on the interval
1350 (-1, 1). Approaching the cut from above or below in general yields a phase
1351 factor with respect to Ferrer's function of the first kind
1352 (cf. `lpmn`).
1354 For ``type=2`` a cut at ``|x| > 1`` is chosen. Approaching the real values
1355 on the interval (-1, 1) in the complex plane yields Ferrer's function
1356 of the first kind.
1358 References
1359 ----------
1360 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
1361 Functions", John Wiley and Sons, 1996.
1362 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html
1363 .. [2] NIST Digital Library of Mathematical Functions
1364 https://dlmf.nist.gov/14.21
1366 """
1367 if not isscalar(m) or (abs(m) > n):
1368 raise ValueError("m must be <= n.")
1369 if not isscalar(n) or (n < 0):
1370 raise ValueError("n must be a non-negative integer.")
1371 if not isscalar(z):
1372 raise ValueError("z must be scalar.")
1373 if not(type == 2 or type == 3):
1374 raise ValueError("type must be either 2 or 3.")
1375 if (m < 0):
1376 mp = -m
1377 mf, nf = mgrid[0:mp+1, 0:n+1]
1378 with ufuncs.errstate(all='ignore'):
1379 if type == 2:
1380 fixarr = where(mf > nf, 0.0,
1381 (-1)**mf * gamma(nf-mf+1) / gamma(nf+mf+1))
1382 else:
1383 fixarr = where(mf > nf, 0.0, gamma(nf-mf+1) / gamma(nf+mf+1))
1384 else:
1385 mp = m
1386 p, pd = specfun.clpmn(mp, n, real(z), imag(z), type)
1387 if (m < 0):
1388 p = p * fixarr
1389 pd = pd * fixarr
1390 return p, pd
1393def lqmn(m, n, z):
1394 """Sequence of associated Legendre functions of the second kind.
1396 Computes the associated Legendre function of the second kind of order m and
1397 degree n, ``Qmn(z)`` = :math:`Q_n^m(z)`, and its derivative, ``Qmn'(z)``.
1398 Returns two arrays of size ``(m+1, n+1)`` containing ``Qmn(z)`` and
1399 ``Qmn'(z)`` for all orders from ``0..m`` and degrees from ``0..n``.
1401 Parameters
1402 ----------
1403 m : int
1404 ``|m| <= n``; the order of the Legendre function.
1405 n : int
1406 where ``n >= 0``; the degree of the Legendre function. Often
1407 called ``l`` (lower case L) in descriptions of the associated
1408 Legendre function
1409 z : complex
1410 Input value.
1412 Returns
1413 -------
1414 Qmn_z : (m+1, n+1) array
1415 Values for all orders 0..m and degrees 0..n
1416 Qmn_d_z : (m+1, n+1) array
1417 Derivatives for all orders 0..m and degrees 0..n
1419 References
1420 ----------
1421 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
1422 Functions", John Wiley and Sons, 1996.
1423 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html
1425 """
1426 if not isscalar(m) or (m < 0):
1427 raise ValueError("m must be a non-negative integer.")
1428 if not isscalar(n) or (n < 0):
1429 raise ValueError("n must be a non-negative integer.")
1430 if not isscalar(z):
1431 raise ValueError("z must be scalar.")
1432 m = int(m)
1433 n = int(n)
1435 # Ensure neither m nor n == 0
1436 mm = max(1, m)
1437 nn = max(1, n)
1439 if iscomplex(z):
1440 q, qd = specfun.clqmn(mm, nn, z)
1441 else:
1442 q, qd = specfun.lqmn(mm, nn, z)
1443 return q[:(m+1), :(n+1)], qd[:(m+1), :(n+1)]
1446def bernoulli(n):
1447 """Bernoulli numbers B0..Bn (inclusive).
1449 References
1450 ----------
1451 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
1452 Functions", John Wiley and Sons, 1996.
1453 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html
1455 """
1456 if not isscalar(n) or (n < 0):
1457 raise ValueError("n must be a non-negative integer.")
1458 n = int(n)
1459 if (n < 2):
1460 n1 = 2
1461 else:
1462 n1 = n
1463 return specfun.bernob(int(n1))[:(n+1)]
1466def euler(n):
1467 """Euler numbers E(0), E(1), ..., E(n).
1469 The Euler numbers [1]_ are also known as the secant numbers.
1471 Because ``euler(n)`` returns floating point values, it does not give
1472 exact values for large `n`. The first inexact value is E(22).
1474 Parameters
1475 ----------
1476 n : int
1477 The highest index of the Euler number to be returned.
1479 Returns
1480 -------
1481 ndarray
1482 The Euler numbers [E(0), E(1), ..., E(n)].
1483 The odd Euler numbers, which are all zero, are included.
1485 References
1486 ----------
1487 .. [1] Sequence A122045, The On-Line Encyclopedia of Integer Sequences,
1488 https://oeis.org/A122045
1489 .. [2] Zhang, Shanjie and Jin, Jianming. "Computation of Special
1490 Functions", John Wiley and Sons, 1996.
1491 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html
1493 Examples
1494 --------
1495 >>> from scipy.special import euler
1496 >>> euler(6)
1497 array([ 1., 0., -1., 0., 5., 0., -61.])
1499 >>> euler(13).astype(np.int64)
1500 array([ 1, 0, -1, 0, 5, 0, -61,
1501 0, 1385, 0, -50521, 0, 2702765, 0])
1503 >>> euler(22)[-1] # Exact value of E(22) is -69348874393137901.
1504 -69348874393137976.0
1506 """
1507 if not isscalar(n) or (n < 0):
1508 raise ValueError("n must be a non-negative integer.")
1509 n = int(n)
1510 if (n < 2):
1511 n1 = 2
1512 else:
1513 n1 = n
1514 return specfun.eulerb(n1)[:(n+1)]
1517def lpn(n, z):
1518 """Legendre function of the first kind.
1520 Compute sequence of Legendre functions of the first kind (polynomials),
1521 Pn(z) and derivatives for all degrees from 0 to n (inclusive).
1523 See also special.legendre for polynomial class.
1525 References
1526 ----------
1527 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
1528 Functions", John Wiley and Sons, 1996.
1529 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html
1531 """
1532 if not (isscalar(n) and isscalar(z)):
1533 raise ValueError("arguments must be scalars.")
1534 n = _nonneg_int_or_fail(n, 'n', strict=False)
1535 if (n < 1):
1536 n1 = 1
1537 else:
1538 n1 = n
1539 if iscomplex(z):
1540 pn, pd = specfun.clpn(n1, z)
1541 else:
1542 pn, pd = specfun.lpn(n1, z)
1543 return pn[:(n+1)], pd[:(n+1)]
1546def lqn(n, z):
1547 """Legendre function of the second kind.
1549 Compute sequence of Legendre functions of the second kind, Qn(z) and
1550 derivatives for all degrees from 0 to n (inclusive).
1552 References
1553 ----------
1554 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
1555 Functions", John Wiley and Sons, 1996.
1556 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html
1558 """
1559 if not (isscalar(n) and isscalar(z)):
1560 raise ValueError("arguments must be scalars.")
1561 n = _nonneg_int_or_fail(n, 'n', strict=False)
1562 if (n < 1):
1563 n1 = 1
1564 else:
1565 n1 = n
1566 if iscomplex(z):
1567 qn, qd = specfun.clqn(n1, z)
1568 else:
1569 qn, qd = specfun.lqnb(n1, z)
1570 return qn[:(n+1)], qd[:(n+1)]
1573def ai_zeros(nt):
1574 """
1575 Compute `nt` zeros and values of the Airy function Ai and its derivative.
1577 Computes the first `nt` zeros, `a`, of the Airy function Ai(x);
1578 first `nt` zeros, `ap`, of the derivative of the Airy function Ai'(x);
1579 the corresponding values Ai(a');
1580 and the corresponding values Ai'(a).
1582 Parameters
1583 ----------
1584 nt : int
1585 Number of zeros to compute
1587 Returns
1588 -------
1589 a : ndarray
1590 First `nt` zeros of Ai(x)
1591 ap : ndarray
1592 First `nt` zeros of Ai'(x)
1593 ai : ndarray
1594 Values of Ai(x) evaluated at first `nt` zeros of Ai'(x)
1595 aip : ndarray
1596 Values of Ai'(x) evaluated at first `nt` zeros of Ai(x)
1598 Examples
1599 --------
1600 >>> from scipy import special
1601 >>> a, ap, ai, aip = special.ai_zeros(3)
1602 >>> a
1603 array([-2.33810741, -4.08794944, -5.52055983])
1604 >>> ap
1605 array([-1.01879297, -3.24819758, -4.82009921])
1606 >>> ai
1607 array([ 0.53565666, -0.41901548, 0.38040647])
1608 >>> aip
1609 array([ 0.70121082, -0.80311137, 0.86520403])
1611 References
1612 ----------
1613 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
1614 Functions", John Wiley and Sons, 1996.
1615 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html
1617 """
1618 kf = 1
1619 if not isscalar(nt) or (floor(nt) != nt) or (nt <= 0):
1620 raise ValueError("nt must be a positive integer scalar.")
1621 return specfun.airyzo(nt, kf)
1624def bi_zeros(nt):
1625 """
1626 Compute `nt` zeros and values of the Airy function Bi and its derivative.
1628 Computes the first `nt` zeros, b, of the Airy function Bi(x);
1629 first `nt` zeros, b', of the derivative of the Airy function Bi'(x);
1630 the corresponding values Bi(b');
1631 and the corresponding values Bi'(b).
1633 Parameters
1634 ----------
1635 nt : int
1636 Number of zeros to compute
1638 Returns
1639 -------
1640 b : ndarray
1641 First `nt` zeros of Bi(x)
1642 bp : ndarray
1643 First `nt` zeros of Bi'(x)
1644 bi : ndarray
1645 Values of Bi(x) evaluated at first `nt` zeros of Bi'(x)
1646 bip : ndarray
1647 Values of Bi'(x) evaluated at first `nt` zeros of Bi(x)
1649 Examples
1650 --------
1651 >>> from scipy import special
1652 >>> b, bp, bi, bip = special.bi_zeros(3)
1653 >>> b
1654 array([-1.17371322, -3.2710933 , -4.83073784])
1655 >>> bp
1656 array([-2.29443968, -4.07315509, -5.51239573])
1657 >>> bi
1658 array([-0.45494438, 0.39652284, -0.36796916])
1659 >>> bip
1660 array([ 0.60195789, -0.76031014, 0.83699101])
1662 References
1663 ----------
1664 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
1665 Functions", John Wiley and Sons, 1996.
1666 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html
1668 """
1669 kf = 2
1670 if not isscalar(nt) or (floor(nt) != nt) or (nt <= 0):
1671 raise ValueError("nt must be a positive integer scalar.")
1672 return specfun.airyzo(nt, kf)
1675def lmbda(v, x):
1676 r"""Jahnke-Emden Lambda function, Lambdav(x).
1678 This function is defined as [2]_,
1680 .. math:: \Lambda_v(x) = \Gamma(v+1) \frac{J_v(x)}{(x/2)^v},
1682 where :math:`\Gamma` is the gamma function and :math:`J_v` is the
1683 Bessel function of the first kind.
1685 Parameters
1686 ----------
1687 v : float
1688 Order of the Lambda function
1689 x : float
1690 Value at which to evaluate the function and derivatives
1692 Returns
1693 -------
1694 vl : ndarray
1695 Values of Lambda_vi(x), for vi=v-int(v), vi=1+v-int(v), ..., vi=v.
1696 dl : ndarray
1697 Derivatives Lambda_vi'(x), for vi=v-int(v), vi=1+v-int(v), ..., vi=v.
1699 References
1700 ----------
1701 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
1702 Functions", John Wiley and Sons, 1996.
1703 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html
1704 .. [2] Jahnke, E. and Emde, F. "Tables of Functions with Formulae and
1705 Curves" (4th ed.), Dover, 1945
1706 """
1707 if not (isscalar(v) and isscalar(x)):
1708 raise ValueError("arguments must be scalars.")
1709 if (v < 0):
1710 raise ValueError("argument must be > 0.")
1711 n = int(v)
1712 v0 = v - n
1713 if (n < 1):
1714 n1 = 1
1715 else:
1716 n1 = n
1717 v1 = n1 + v0
1718 if (v != floor(v)):
1719 vm, vl, dl = specfun.lamv(v1, x)
1720 else:
1721 vm, vl, dl = specfun.lamn(v1, x)
1722 return vl[:(n+1)], dl[:(n+1)]
1725def pbdv_seq(v, x):
1726 """Parabolic cylinder functions Dv(x) and derivatives.
1728 Parameters
1729 ----------
1730 v : float
1731 Order of the parabolic cylinder function
1732 x : float
1733 Value at which to evaluate the function and derivatives
1735 Returns
1736 -------
1737 dv : ndarray
1738 Values of D_vi(x), for vi=v-int(v), vi=1+v-int(v), ..., vi=v.
1739 dp : ndarray
1740 Derivatives D_vi'(x), for vi=v-int(v), vi=1+v-int(v), ..., vi=v.
1742 References
1743 ----------
1744 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
1745 Functions", John Wiley and Sons, 1996, chapter 13.
1746 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html
1748 """
1749 if not (isscalar(v) and isscalar(x)):
1750 raise ValueError("arguments must be scalars.")
1751 n = int(v)
1752 v0 = v-n
1753 if (n < 1):
1754 n1 = 1
1755 else:
1756 n1 = n
1757 v1 = n1 + v0
1758 dv, dp, pdf, pdd = specfun.pbdv(v1, x)
1759 return dv[:n1+1], dp[:n1+1]
1762def pbvv_seq(v, x):
1763 """Parabolic cylinder functions Vv(x) and derivatives.
1765 Parameters
1766 ----------
1767 v : float
1768 Order of the parabolic cylinder function
1769 x : float
1770 Value at which to evaluate the function and derivatives
1772 Returns
1773 -------
1774 dv : ndarray
1775 Values of V_vi(x), for vi=v-int(v), vi=1+v-int(v), ..., vi=v.
1776 dp : ndarray
1777 Derivatives V_vi'(x), for vi=v-int(v), vi=1+v-int(v), ..., vi=v.
1779 References
1780 ----------
1781 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
1782 Functions", John Wiley and Sons, 1996, chapter 13.
1783 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html
1785 """
1786 if not (isscalar(v) and isscalar(x)):
1787 raise ValueError("arguments must be scalars.")
1788 n = int(v)
1789 v0 = v-n
1790 if (n <= 1):
1791 n1 = 1
1792 else:
1793 n1 = n
1794 v1 = n1 + v0
1795 dv, dp, pdf, pdd = specfun.pbvv(v1, x)
1796 return dv[:n1+1], dp[:n1+1]
1799def pbdn_seq(n, z):
1800 """Parabolic cylinder functions Dn(z) and derivatives.
1802 Parameters
1803 ----------
1804 n : int
1805 Order of the parabolic cylinder function
1806 z : complex
1807 Value at which to evaluate the function and derivatives
1809 Returns
1810 -------
1811 dv : ndarray
1812 Values of D_i(z), for i=0, ..., i=n.
1813 dp : ndarray
1814 Derivatives D_i'(z), for i=0, ..., i=n.
1816 References
1817 ----------
1818 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
1819 Functions", John Wiley and Sons, 1996, chapter 13.
1820 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html
1822 """
1823 if not (isscalar(n) and isscalar(z)):
1824 raise ValueError("arguments must be scalars.")
1825 if (floor(n) != n):
1826 raise ValueError("n must be an integer.")
1827 if (abs(n) <= 1):
1828 n1 = 1
1829 else:
1830 n1 = n
1831 cpb, cpd = specfun.cpbdn(n1, z)
1832 return cpb[:n1+1], cpd[:n1+1]
1835def ber_zeros(nt):
1836 """Compute nt zeros of the Kelvin function ber.
1838 Parameters
1839 ----------
1840 nt : int
1841 Number of zeros to compute. Must be positive.
1843 Returns
1844 -------
1845 ndarray
1846 First `nt` zeros of the Kelvin function.
1848 See Also
1849 --------
1850 ber
1852 References
1853 ----------
1854 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
1855 Functions", John Wiley and Sons, 1996.
1856 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html
1858 """
1859 if not isscalar(nt) or (floor(nt) != nt) or (nt <= 0):
1860 raise ValueError("nt must be positive integer scalar.")
1861 return specfun.klvnzo(nt, 1)
1864def bei_zeros(nt):
1865 """Compute nt zeros of the Kelvin function bei.
1867 Parameters
1868 ----------
1869 nt : int
1870 Number of zeros to compute. Must be positive.
1872 Returns
1873 -------
1874 ndarray
1875 First `nt` zeros of the Kelvin function.
1877 See Also
1878 --------
1879 bei
1881 References
1882 ----------
1883 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
1884 Functions", John Wiley and Sons, 1996.
1885 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html
1887 """
1888 if not isscalar(nt) or (floor(nt) != nt) or (nt <= 0):
1889 raise ValueError("nt must be positive integer scalar.")
1890 return specfun.klvnzo(nt, 2)
1893def ker_zeros(nt):
1894 """Compute nt zeros of the Kelvin function ker.
1896 Parameters
1897 ----------
1898 nt : int
1899 Number of zeros to compute. Must be positive.
1901 Returns
1902 -------
1903 ndarray
1904 First `nt` zeros of the Kelvin function.
1906 See Also
1907 --------
1908 ker
1910 References
1911 ----------
1912 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
1913 Functions", John Wiley and Sons, 1996.
1914 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html
1916 """
1917 if not isscalar(nt) or (floor(nt) != nt) or (nt <= 0):
1918 raise ValueError("nt must be positive integer scalar.")
1919 return specfun.klvnzo(nt, 3)
1922def kei_zeros(nt):
1923 """Compute nt zeros of the Kelvin function kei.
1925 Parameters
1926 ----------
1927 nt : int
1928 Number of zeros to compute. Must be positive.
1930 Returns
1931 -------
1932 ndarray
1933 First `nt` zeros of the Kelvin function.
1935 See Also
1936 --------
1937 kei
1939 References
1940 ----------
1941 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
1942 Functions", John Wiley and Sons, 1996.
1943 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html
1945 """
1946 if not isscalar(nt) or (floor(nt) != nt) or (nt <= 0):
1947 raise ValueError("nt must be positive integer scalar.")
1948 return specfun.klvnzo(nt, 4)
1951def berp_zeros(nt):
1952 """Compute nt zeros of the derivative of the Kelvin function ber.
1954 Parameters
1955 ----------
1956 nt : int
1957 Number of zeros to compute. Must be positive.
1959 Returns
1960 -------
1961 ndarray
1962 First `nt` zeros of the derivative of the Kelvin function.
1964 See Also
1965 --------
1966 ber, berp
1968 References
1969 ----------
1970 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
1971 Functions", John Wiley and Sons, 1996.
1972 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html
1974 """
1975 if not isscalar(nt) or (floor(nt) != nt) or (nt <= 0):
1976 raise ValueError("nt must be positive integer scalar.")
1977 return specfun.klvnzo(nt, 5)
1980def beip_zeros(nt):
1981 """Compute nt zeros of the derivative of the Kelvin function bei.
1983 Parameters
1984 ----------
1985 nt : int
1986 Number of zeros to compute. Must be positive.
1988 Returns
1989 -------
1990 ndarray
1991 First `nt` zeros of the derivative of the Kelvin function.
1993 See Also
1994 --------
1995 bei, beip
1997 References
1998 ----------
1999 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
2000 Functions", John Wiley and Sons, 1996.
2001 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html
2003 """
2004 if not isscalar(nt) or (floor(nt) != nt) or (nt <= 0):
2005 raise ValueError("nt must be positive integer scalar.")
2006 return specfun.klvnzo(nt, 6)
2009def kerp_zeros(nt):
2010 """Compute nt zeros of the derivative of the Kelvin function ker.
2012 Parameters
2013 ----------
2014 nt : int
2015 Number of zeros to compute. Must be positive.
2017 Returns
2018 -------
2019 ndarray
2020 First `nt` zeros of the derivative of the Kelvin function.
2022 See Also
2023 --------
2024 ker, kerp
2026 References
2027 ----------
2028 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
2029 Functions", John Wiley and Sons, 1996.
2030 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html
2032 """
2033 if not isscalar(nt) or (floor(nt) != nt) or (nt <= 0):
2034 raise ValueError("nt must be positive integer scalar.")
2035 return specfun.klvnzo(nt, 7)
2038def keip_zeros(nt):
2039 """Compute nt zeros of the derivative of the Kelvin function kei.
2041 Parameters
2042 ----------
2043 nt : int
2044 Number of zeros to compute. Must be positive.
2046 Returns
2047 -------
2048 ndarray
2049 First `nt` zeros of the derivative of the Kelvin function.
2051 See Also
2052 --------
2053 kei, keip
2055 References
2056 ----------
2057 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
2058 Functions", John Wiley and Sons, 1996.
2059 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html
2061 """
2062 if not isscalar(nt) or (floor(nt) != nt) or (nt <= 0):
2063 raise ValueError("nt must be positive integer scalar.")
2064 return specfun.klvnzo(nt, 8)
2067def kelvin_zeros(nt):
2068 """Compute nt zeros of all Kelvin functions.
2070 Returned in a length-8 tuple of arrays of length nt. The tuple contains
2071 the arrays of zeros of (ber, bei, ker, kei, ber', bei', ker', kei').
2073 References
2074 ----------
2075 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
2076 Functions", John Wiley and Sons, 1996.
2077 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html
2079 """
2080 if not isscalar(nt) or (floor(nt) != nt) or (nt <= 0):
2081 raise ValueError("nt must be positive integer scalar.")
2082 return (specfun.klvnzo(nt, 1),
2083 specfun.klvnzo(nt, 2),
2084 specfun.klvnzo(nt, 3),
2085 specfun.klvnzo(nt, 4),
2086 specfun.klvnzo(nt, 5),
2087 specfun.klvnzo(nt, 6),
2088 specfun.klvnzo(nt, 7),
2089 specfun.klvnzo(nt, 8))
2092def pro_cv_seq(m, n, c):
2093 """Characteristic values for prolate spheroidal wave functions.
2095 Compute a sequence of characteristic values for the prolate
2096 spheroidal wave functions for mode m and n'=m..n and spheroidal
2097 parameter c.
2099 References
2100 ----------
2101 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
2102 Functions", John Wiley and Sons, 1996.
2103 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html
2105 """
2106 if not (isscalar(m) and isscalar(n) and isscalar(c)):
2107 raise ValueError("Arguments must be scalars.")
2108 if (n != floor(n)) or (m != floor(m)):
2109 raise ValueError("Modes must be integers.")
2110 if (n-m > 199):
2111 raise ValueError("Difference between n and m is too large.")
2112 maxL = n-m+1
2113 return specfun.segv(m, n, c, 1)[1][:maxL]
2116def obl_cv_seq(m, n, c):
2117 """Characteristic values for oblate spheroidal wave functions.
2119 Compute a sequence of characteristic values for the oblate
2120 spheroidal wave functions for mode m and n'=m..n and spheroidal
2121 parameter c.
2123 References
2124 ----------
2125 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
2126 Functions", John Wiley and Sons, 1996.
2127 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html
2129 """
2130 if not (isscalar(m) and isscalar(n) and isscalar(c)):
2131 raise ValueError("Arguments must be scalars.")
2132 if (n != floor(n)) or (m != floor(m)):
2133 raise ValueError("Modes must be integers.")
2134 if (n-m > 199):
2135 raise ValueError("Difference between n and m is too large.")
2136 maxL = n-m+1
2137 return specfun.segv(m, n, c, -1)[1][:maxL]
2140def comb(N, k, exact=False, repetition=False):
2141 """The number of combinations of N things taken k at a time.
2143 This is often expressed as "N choose k".
2145 Parameters
2146 ----------
2147 N : int, ndarray
2148 Number of things.
2149 k : int, ndarray
2150 Number of elements taken.
2151 exact : bool, optional
2152 If `exact` is False, then floating point precision is used, otherwise
2153 exact long integer is computed.
2154 repetition : bool, optional
2155 If `repetition` is True, then the number of combinations with
2156 repetition is computed.
2158 Returns
2159 -------
2160 val : int, float, ndarray
2161 The total number of combinations.
2163 See Also
2164 --------
2165 binom : Binomial coefficient ufunc
2167 Notes
2168 -----
2169 - Array arguments accepted only for exact=False case.
2170 - If N < 0, or k < 0, then 0 is returned.
2171 - If k > N and repetition=False, then 0 is returned.
2173 Examples
2174 --------
2175 >>> from scipy.special import comb
2176 >>> k = np.array([3, 4])
2177 >>> n = np.array([10, 10])
2178 >>> comb(n, k, exact=False)
2179 array([ 120., 210.])
2180 >>> comb(10, 3, exact=True)
2181 120
2182 >>> comb(10, 3, exact=True, repetition=True)
2183 220
2185 """
2186 if repetition:
2187 return comb(N + k - 1, k, exact)
2188 if exact:
2189 return _comb_int(N, k)
2190 else:
2191 k, N = asarray(k), asarray(N)
2192 cond = (k <= N) & (N >= 0) & (k >= 0)
2193 vals = binom(N, k)
2194 if isinstance(vals, np.ndarray):
2195 vals[~cond] = 0
2196 elif not cond:
2197 vals = np.float64(0)
2198 return vals
2201def perm(N, k, exact=False):
2202 """Permutations of N things taken k at a time, i.e., k-permutations of N.
2204 It's also known as "partial permutations".
2206 Parameters
2207 ----------
2208 N : int, ndarray
2209 Number of things.
2210 k : int, ndarray
2211 Number of elements taken.
2212 exact : bool, optional
2213 If `exact` is False, then floating point precision is used, otherwise
2214 exact long integer is computed.
2216 Returns
2217 -------
2218 val : int, ndarray
2219 The number of k-permutations of N.
2221 Notes
2222 -----
2223 - Array arguments accepted only for exact=False case.
2224 - If k > N, N < 0, or k < 0, then a 0 is returned.
2226 Examples
2227 --------
2228 >>> from scipy.special import perm
2229 >>> k = np.array([3, 4])
2230 >>> n = np.array([10, 10])
2231 >>> perm(n, k)
2232 array([ 720., 5040.])
2233 >>> perm(10, 3, exact=True)
2234 720
2236 """
2237 if exact:
2238 if (k > N) or (N < 0) or (k < 0):
2239 return 0
2240 val = 1
2241 for i in range(N - k + 1, N + 1):
2242 val *= i
2243 return val
2244 else:
2245 k, N = asarray(k), asarray(N)
2246 cond = (k <= N) & (N >= 0) & (k >= 0)
2247 vals = poch(N - k + 1, k)
2248 if isinstance(vals, np.ndarray):
2249 vals[~cond] = 0
2250 elif not cond:
2251 vals = np.float64(0)
2252 return vals
2255# https://stackoverflow.com/a/16327037
2256def _range_prod(lo, hi):
2257 """
2258 Product of a range of numbers.
2260 Returns the product of
2261 lo * (lo+1) * (lo+2) * ... * (hi-2) * (hi-1) * hi
2262 = hi! / (lo-1)!
2264 Breaks into smaller products first for speed:
2265 _range_prod(2, 9) = ((2*3)*(4*5))*((6*7)*(8*9))
2266 """
2267 if lo + 1 < hi:
2268 mid = (hi + lo) // 2
2269 return _range_prod(lo, mid) * _range_prod(mid + 1, hi)
2270 if lo == hi:
2271 return lo
2272 return lo * hi
2275def factorial(n, exact=False):
2276 """
2277 The factorial of a number or array of numbers.
2279 The factorial of non-negative integer `n` is the product of all
2280 positive integers less than or equal to `n`::
2282 n! = n * (n - 1) * (n - 2) * ... * 1
2284 Parameters
2285 ----------
2286 n : int or array_like of ints
2287 Input values. If ``n < 0``, the return value is 0.
2288 exact : bool, optional
2289 If True, calculate the answer exactly using long integer arithmetic.
2290 If False, result is approximated in floating point rapidly using the
2291 `gamma` function.
2292 Default is False.
2294 Returns
2295 -------
2296 nf : float or int or ndarray
2297 Factorial of `n`, as integer or float depending on `exact`.
2299 Notes
2300 -----
2301 For arrays with ``exact=True``, the factorial is computed only once, for
2302 the largest input, with each other result computed in the process.
2303 The output dtype is increased to ``int64`` or ``object`` if necessary.
2305 With ``exact=False`` the factorial is approximated using the gamma
2306 function:
2308 .. math:: n! = \\Gamma(n+1)
2310 Examples
2311 --------
2312 >>> from scipy.special import factorial
2313 >>> arr = np.array([3, 4, 5])
2314 >>> factorial(arr, exact=False)
2315 array([ 6., 24., 120.])
2316 >>> factorial(arr, exact=True)
2317 array([ 6, 24, 120])
2318 >>> factorial(5, exact=True)
2319 120
2321 """
2322 if exact:
2323 if np.ndim(n) == 0:
2324 if np.isnan(n):
2325 return n
2326 return 0 if n < 0 else math.factorial(n)
2327 else:
2328 n = asarray(n)
2329 un = np.unique(n).astype(object)
2331 # Convert to object array of long ints if np.int can't handle size
2332 if np.isnan(n).any():
2333 dt = float
2334 elif un[-1] > 20:
2335 dt = object
2336 elif un[-1] > 12:
2337 dt = np.int64
2338 else:
2339 dt = np.int
2341 out = np.empty_like(n, dtype=dt)
2343 # Handle invalid/trivial values
2344 # Ignore runtime warning when less operator used w/np.nan
2345 with np.errstate(all='ignore'):
2346 un = un[un > 1]
2347 out[n < 2] = 1
2348 out[n < 0] = 0
2350 # Calculate products of each range of numbers
2351 if un.size:
2352 val = math.factorial(un[0])
2353 out[n == un[0]] = val
2354 for i in range(len(un) - 1):
2355 prev = un[i] + 1
2356 current = un[i + 1]
2357 val *= _range_prod(prev, current)
2358 out[n == current] = val
2360 if np.isnan(n).any():
2361 out = out.astype(np.float64)
2362 out[np.isnan(n)] = n[np.isnan(n)]
2363 return out
2364 else:
2365 out = ufuncs._factorial(n)
2366 return out
2369def factorial2(n, exact=False):
2370 """Double factorial.
2372 This is the factorial with every second value skipped. E.g., ``7!! = 7 * 5
2373 * 3 * 1``. It can be approximated numerically as::
2375 n!! = special.gamma(n/2+1)*2**((m+1)/2)/sqrt(pi) n odd
2376 = 2**(n/2) * (n/2)! n even
2378 Parameters
2379 ----------
2380 n : int or array_like
2381 Calculate ``n!!``. Arrays are only supported with `exact` set
2382 to False. If ``n < 0``, the return value is 0.
2383 exact : bool, optional
2384 The result can be approximated rapidly using the gamma-formula
2385 above (default). If `exact` is set to True, calculate the
2386 answer exactly using integer arithmetic.
2388 Returns
2389 -------
2390 nff : float or int
2391 Double factorial of `n`, as an int or a float depending on
2392 `exact`.
2394 Examples
2395 --------
2396 >>> from scipy.special import factorial2
2397 >>> factorial2(7, exact=False)
2398 array(105.00000000000001)
2399 >>> factorial2(7, exact=True)
2400 105
2402 """
2403 if exact:
2404 if n < -1:
2405 return 0
2406 if n <= 0:
2407 return 1
2408 val = 1
2409 for k in range(n, 0, -2):
2410 val *= k
2411 return val
2412 else:
2413 n = asarray(n)
2414 vals = zeros(n.shape, 'd')
2415 cond1 = (n % 2) & (n >= -1)
2416 cond2 = (1-(n % 2)) & (n >= -1)
2417 oddn = extract(cond1, n)
2418 evenn = extract(cond2, n)
2419 nd2o = oddn / 2.0
2420 nd2e = evenn / 2.0
2421 place(vals, cond1, gamma(nd2o + 1) / sqrt(pi) * pow(2.0, nd2o + 0.5))
2422 place(vals, cond2, gamma(nd2e + 1) * pow(2.0, nd2e))
2423 return vals
2426def factorialk(n, k, exact=True):
2427 """Multifactorial of n of order k, n(!!...!).
2429 This is the multifactorial of n skipping k values. For example,
2431 factorialk(17, 4) = 17!!!! = 17 * 13 * 9 * 5 * 1
2433 In particular, for any integer ``n``, we have
2435 factorialk(n, 1) = factorial(n)
2437 factorialk(n, 2) = factorial2(n)
2439 Parameters
2440 ----------
2441 n : int
2442 Calculate multifactorial. If `n` < 0, the return value is 0.
2443 k : int
2444 Order of multifactorial.
2445 exact : bool, optional
2446 If exact is set to True, calculate the answer exactly using
2447 integer arithmetic.
2449 Returns
2450 -------
2451 val : int
2452 Multifactorial of `n`.
2454 Raises
2455 ------
2456 NotImplementedError
2457 Raises when exact is False
2459 Examples
2460 --------
2461 >>> from scipy.special import factorialk
2462 >>> factorialk(5, 1, exact=True)
2463 120
2464 >>> factorialk(5, 3, exact=True)
2465 10
2467 """
2468 if exact:
2469 if n < 1-k:
2470 return 0
2471 if n <= 0:
2472 return 1
2473 val = 1
2474 for j in range(n, 0, -k):
2475 val = val*j
2476 return val
2477 else:
2478 raise NotImplementedError
2481def zeta(x, q=None, out=None):
2482 r"""
2483 Riemann or Hurwitz zeta function.
2485 Parameters
2486 ----------
2487 x : array_like of float
2488 Input data, must be real
2489 q : array_like of float, optional
2490 Input data, must be real. Defaults to Riemann zeta.
2491 out : ndarray, optional
2492 Output array for the computed values.
2494 Returns
2495 -------
2496 out : array_like
2497 Values of zeta(x).
2499 Notes
2500 -----
2501 The two-argument version is the Hurwitz zeta function
2503 .. math::
2505 \zeta(x, q) = \sum_{k=0}^{\infty} \frac{1}{(k + q)^x};
2507 see [dlmf]_ for details. The Riemann zeta function corresponds to
2508 the case when ``q = 1``.
2510 See Also
2511 --------
2512 zetac
2514 References
2515 ----------
2516 .. [dlmf] NIST, Digital Library of Mathematical Functions,
2517 https://dlmf.nist.gov/25.11#i
2519 Examples
2520 --------
2521 >>> from scipy.special import zeta, polygamma, factorial
2523 Some specific values:
2525 >>> zeta(2), np.pi**2/6
2526 (1.6449340668482266, 1.6449340668482264)
2528 >>> zeta(4), np.pi**4/90
2529 (1.0823232337111381, 1.082323233711138)
2531 Relation to the `polygamma` function:
2533 >>> m = 3
2534 >>> x = 1.25
2535 >>> polygamma(m, x)
2536 array(2.782144009188397)
2537 >>> (-1)**(m+1) * factorial(m) * zeta(m+1, x)
2538 2.7821440091883969
2540 """
2541 if q is None:
2542 return ufuncs._riemann_zeta(x, out)
2543 else:
2544 return ufuncs._zeta(x, q, out)