1# 2# Author: Travis Oliphant, 2002 3# 4 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 18 19 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] 92 93 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)) from err 107 return n 108 109 110def diric(x, n): 111 """Periodic sinc function, also called the Dirichlet function. 112 113 The Dirichlet function is defined as:: 114 115 diric(x, n) = sin(x * n/2) / (n * sin(x / 2)), 116 117 where `n` is a positive integer. 118 119 Parameters 120 ---------- 121 x : array_like 122 Input data 123 n : int 124 Integer defining the periodicity. 125 126 Returns 127 ------- 128 diric : ndarray 129 130 Examples 131 -------- 132 >>> from scipy import special 133 >>> import matplotlib.pyplot as plt 134 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() 142 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. 146 147 Suppress output of values that are effectively 0: 148 149 >>> np.set_printoptions(suppress=True) 150 151 Create a signal `x` of length `m` with `k` ones: 152 153 >>> m = 8 154 >>> k = 3 155 >>> x = np.zeros(m) 156 >>> x[:k] = 1 157 158 Use the FFT to compute the Fourier transform of `x`, and 159 inspect the magnitudes of the coefficients: 160 161 >>> np.abs(np.fft.fft(x)) 162 array([ 3. , 2.41421356, 1. , 0.41421356, 1. , 163 0.41421356, 1. , 2.41421356]) 164 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`: 168 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) 182 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 191 192 mask1 = (n <= 0) | (n != floor(n)) 193 place(y, mask1, nan) 194 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))) 202 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 209 210 211def jnjnp_zeros(nt): 212 """Compute zeros of integer-order Bessel functions Jn and Jn'. 213 214 Results are arranged in order of the magnitudes of the zeros. 215 216 Parameters 217 ---------- 218 nt : int 219 Number (<=1200) of zeros to compute 220 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`. 233 234 See Also 235 -------- 236 jn_zeros, jnp_zeros : to get separated arrays of zeros. 237 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 243 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] 250 251 252def jnyn_zeros(n, nt): 253 """Compute nt zeros of Bessel functions Jn(x), Jn'(x), Yn(x), and Yn'(x). 254 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. 258 259 Parameters 260 ---------- 261 n : int 262 Order of the Bessel functions 263 nt : int 264 Number (<=1200) of zeros to compute 265 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' 276 277 See Also 278 -------- 279 jn_zeros, jnp_zeros, yn_zeros, ynp_zeros 280 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 286 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) 295 296 297def jn_zeros(n, nt): 298 r"""Compute zeros of integer-order Bessel functions Jn. 299 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`. 304 305 Parameters 306 ---------- 307 n : int 308 Order of Bessel function 309 nt : int 310 Number of zeros to return 311 312 Returns 313 ------- 314 ndarray 315 First `n` zeros of the Bessel function. 316 317 See Also 318 -------- 319 jv 320 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 326 327 Examples 328 -------- 329 >>> import scipy.special as sc 330 331 We can check that we are getting approximations of the zeros by 332 evaluating them with `jv`. 333 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]) 340 341 Note that the zero at ``x = 0`` for ``n > 0`` is not included. 342 343 >>> sc.jv(1, 0) 344 0.0 345 346 """ 347 return jnyn_zeros(n, nt)[0] 348 349 350def jnp_zeros(n, nt): 351 r"""Compute zeros of integer-order Bessel function derivatives Jn'. 352 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`. 357 358 Parameters 359 ---------- 360 n : int 361 Order of Bessel function 362 nt : int 363 Number of zeros to return 364 365 Returns 366 ------- 367 ndarray 368 First `n` zeros of the Bessel function. 369 370 See Also 371 -------- 372 jvp, jv 373 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 379 380 Examples 381 -------- 382 >>> import scipy.special as sc 383 384 We can check that we are getting approximations of the zeros by 385 evaluating them with `jvp`. 386 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]) 393 394 Note that the zero at ``x = 0`` for ``n > 1`` is not included. 395 396 >>> sc.jvp(n, 0) 397 0.0 398 399 """ 400 return jnyn_zeros(n, nt)[1] 401 402 403def yn_zeros(n, nt): 404 r"""Compute zeros of integer-order Bessel function Yn(x). 405 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. 408 409 Parameters 410 ---------- 411 n : int 412 Order of Bessel function 413 nt : int 414 Number of zeros to return 415 416 Returns 417 ------- 418 ndarray 419 First `n` zeros of the Bessel function. 420 421 See Also 422 -------- 423 yn, yv 424 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 430 431 Examples 432 -------- 433 >>> import scipy.special as sc 434 435 We can check that we are getting approximations of the zeros by 436 evaluating them with `yn`. 437 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]) 444 445 """ 446 return jnyn_zeros(n, nt)[2] 447 448 449def ynp_zeros(n, nt): 450 r"""Compute zeros of integer-order Bessel function derivatives Yn'(x). 451 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. 455 456 Parameters 457 ---------- 458 n : int 459 Order of Bessel function 460 nt : int 461 Number of zeros to return 462 463 See Also 464 -------- 465 yvp 466 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 472 473 Examples 474 -------- 475 >>> import scipy.special as sc 476 477 We can check that we are getting approximations of the zeros by 478 evaluating them with `yvp`. 479 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]) 486 487 """ 488 return jnyn_zeros(n, nt)[3] 489 490 491def y0_zeros(nt, complex=False): 492 """Compute nt zeros of Bessel function Y0(z), and derivative at each zero. 493 494 The derivatives are given by Y0'(z0) = -Y1(z0) at each zero z0. 495 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. 505 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 512 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 518 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) 525 526 527def y1_zeros(nt, complex=False): 528 """Compute nt zeros of Bessel function Y1(z), and derivative at each zero. 529 530 The derivatives are given by Y1'(z1) = Y0(z1) at each zero z1. 531 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. 541 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 548 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 554 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) 561 562 563def y1p_zeros(nt, complex=False): 564 """Compute nt zeros of Bessel derivative Y1'(z), and value at each zero. 565 566 The values are given by Y1(z1) at each z1 where Y1'(z1)=0. 567 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. 577 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 584 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 590 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) 597 598 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) 611 612 613def jvp(v, z, n=1): 614 """Compute derivatives of Bessel functions of the first kind. 615 616 Compute the nth derivative of the Bessel function `Jv` with 617 respect to `z`. 618 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 628 629 Returns 630 ------- 631 scalar or ndarray 632 Values of the derivative of the Bessel function. 633 634 Notes 635 ----- 636 The derivative is computed using the relation DLFM 10.6.7 [2]_. 637 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 645 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) 652 653 654def yvp(v, z, n=1): 655 """Compute derivatives of Bessel functions of the second kind. 656 657 Compute the nth derivative of the Bessel function `Yv` with 658 respect to `z`. 659 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 668 669 Returns 670 ------- 671 scalar or ndarray 672 nth derivative of the Bessel function. 673 674 Notes 675 ----- 676 The derivative is computed using the relation DLFM 10.6.7 [2]_. 677 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 685 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) 692 693 694def kvp(v, z, n=1): 695 """Compute nth derivative of real-order modified Bessel function Kv(z) 696 697 Kv(z) is the modified Bessel function of the second kind. 698 Derivative is calculated with respect to `z`. 699 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. 708 709 Returns 710 ------- 711 out : ndarray 712 The results 713 714 Examples 715 -------- 716 Calculate multiple values at order 5: 717 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]) 722 723 724 Calculate for a single value at multiple orders: 725 726 >>> kvp((4, 4.5, 5), 1) 727 array([ -184.0309, -568.9585, -1849.0354]) 728 729 Notes 730 ----- 731 The derivative is computed using the relation DLFM 10.29.5 [2]_. 732 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 740 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) 747 748 749def ivp(v, z, n=1): 750 """Compute derivatives of modified Bessel functions of the first kind. 751 752 Compute the nth derivative of the modified Bessel function `Iv` 753 with respect to `z`. 754 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 764 765 Returns 766 ------- 767 scalar or ndarray 768 nth derivative of the modified Bessel function. 769 770 See Also 771 -------- 772 iv 773 774 Notes 775 ----- 776 The derivative is computed using the relation DLFM 10.29.5 [2]_. 777 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 785 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) 792 793 794def h1vp(v, z, n=1): 795 """Compute nth derivative of Hankel function H1v(z) with respect to `z`. 796 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 806 807 Returns 808 ------- 809 scalar or ndarray 810 Values of the derivative of the Hankel function. 811 812 Notes 813 ----- 814 The derivative is computed using the relation DLFM 10.6.7 [2]_. 815 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 823 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) 830 831 832def h2vp(v, z, n=1): 833 """Compute nth derivative of Hankel function H2v(z) with respect to `z`. 834 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 844 845 Returns 846 ------- 847 scalar or ndarray 848 Values of the derivative of the Hankel function. 849 850 Notes 851 ----- 852 The derivative is computed using the relation DLFM 10.6.7 [2]_. 853 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 861 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) 868 869 870def riccati_jn(n, x): 871 r"""Compute Ricatti-Bessel function of the first kind and its derivative. 872 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`. 876 877 This function computes the value and first derivative of the 878 Ricatti-Bessel function for all orders up to and including `n`. 879 880 Parameters 881 ---------- 882 n : int 883 Maximum order of function to compute 884 x : float 885 Argument at which to evaluate 886 887 Returns 888 ------- 889 jn : ndarray 890 Value of j0(x), ..., jn(x) 891 jnp : ndarray 892 First derivative j0'(x), ..., jn'(x) 893 894 Notes 895 ----- 896 The computation is carried out via backward recurrence, using the 897 relation DLMF 10.51.1 [2]_. 898 899 Wrapper for a Fortran routine created by Shanjie Zhang and Jianming 900 Jin [1]_. 901 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 909 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)] 920 921 922def riccati_yn(n, x): 923 """Compute Ricatti-Bessel function of the second kind and its derivative. 924 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`. 928 929 This function computes the value and first derivative of the function for 930 all orders up to and including `n`. 931 932 Parameters 933 ---------- 934 n : int 935 Maximum order of function to compute 936 x : float 937 Argument at which to evaluate 938 939 Returns 940 ------- 941 yn : ndarray 942 Value of y0(x), ..., yn(x) 943 ynp : ndarray 944 First derivative y0'(x), ..., yn'(x) 945 946 Notes 947 ----- 948 The computation is carried out via ascending recurrence, using the 949 relation DLMF 10.51.1 [2]_. 950 951 Wrapper for a Fortran routine created by Shanjie Zhang and Jianming 952 Jin [1]_. 953 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 961 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)] 972 973 974def erf_zeros(nt): 975 """Compute the first nt zero in the first quadrant, ordered by absolute value. 976 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)). 979 980 981 Parameters 982 ---------- 983 nt : int 984 The number of zeros to compute 985 986 Returns 987 ------- 988 The locations of the zeros of erf : ndarray (complex) 989 Complex values at which zeros of erf(z) 990 991 Examples 992 -------- 993 >>> from scipy import special 994 >>> special.erf_zeros(1) 995 array([1.45061616+1.880943j]) 996 997 Check that erf is (close to) zero for the value returned by erf_zeros 998 999 >>> special.erf(special.erf_zeros(1)) 1000 array([4.95159469e-14-1.16407394e-16j]) 1001 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 1007 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) 1012 1013 1014def fresnelc_zeros(nt): 1015 """Compute nt complex zeros of cosine Fresnel integral C(z). 1016 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 1022 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) 1027 1028 1029def fresnels_zeros(nt): 1030 """Compute nt complex zeros of sine Fresnel integral S(z). 1031 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 1037 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) 1042 1043 1044def fresnel_zeros(nt): 1045 """Compute nt complex zeros of sine and cosine Fresnel integrals S(z) and C(z). 1046 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 1052 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) 1057 1058 1059def assoc_laguerre(x, n, k=0.0): 1060 """Compute the generalized (associated) Laguerre polynomial of degree n and order k. 1061 1062 The polynomial :math:`L^{(k)}_n(x)` is orthogonal over ``[0, inf)``, 1063 with weighting function ``exp(-x) * x**k`` with ``k > -1``. 1064 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)``. 1069 1070 """ 1071 return orthogonal.eval_genlaguerre(n, k, x) 1072 1073 1074digamma = psi 1075 1076 1077def polygamma(n, x): 1078 r"""Polygamma functions. 1079 1080 Defined as :math:`\psi^{(n)}(x)` where :math:`\psi` is the 1081 `digamma` function. See [dlmf]_ for details. 1082 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 1090 1091 Returns 1092 ------- 1093 ndarray 1094 Function results 1095 1096 See Also 1097 -------- 1098 digamma 1099 1100 References 1101 ---------- 1102 .. [dlmf] NIST, Digital Library of Mathematical Functions, 1103 https://dlmf.nist.gov/5.15 1104 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) 1113 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) 1118 1119 1120def mathieu_even_coef(m, q): 1121 r"""Fourier coefficients for even Mathieu and modified Mathieu functions. 1122 1123 The Fourier series of the even solutions of the Mathieu differential 1124 equation are of the form 1125 1126 .. math:: \mathrm{ce}_{2n}(z, q) = \sum_{k=0}^{\infty} A_{(2n)}^{(2k)} \cos 2kz 1127 1128 .. math:: \mathrm{ce}_{2n+1}(z, q) = \sum_{k=0}^{\infty} A_{(2n+1)}^{(2k+1)} \cos (2k+1)z 1129 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. 1133 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. 1140 1141 Returns 1142 ------- 1143 Ak : ndarray 1144 Even or odd Fourier coefficients, corresponding to even or odd m. 1145 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 1153 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.") 1161 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 1173 1174 a = mathieu_a(m, q) 1175 fc = specfun.fcoef(kd, m, q, a) 1176 return fc[:km] 1177 1178 1179def mathieu_odd_coef(m, q): 1180 r"""Fourier coefficients for even Mathieu and modified Mathieu functions. 1181 1182 The Fourier series of the odd solutions of the Mathieu differential 1183 equation are of the form 1184 1185 .. math:: \mathrm{se}_{2n+1}(z, q) = \sum_{k=0}^{\infty} B_{(2n+1)}^{(2k+1)} \sin (2k+1)z 1186 1187 .. math:: \mathrm{se}_{2n+2}(z, q) = \sum_{k=0}^{\infty} B_{(2n+2)}^{(2k+2)} \sin (2k+2)z 1188 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. 1192 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. 1199 1200 Returns 1201 ------- 1202 Bk : ndarray 1203 Even or odd Fourier coefficients, corresponding to even or odd m. 1204 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 1210 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") 1218 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 1230 1231 b = mathieu_b(m, q) 1232 fc = specfun.fcoef(kd, m, q, b) 1233 return fc[:km] 1234 1235 1236def lpmn(m, n, z): 1237 """Sequence of associated Legendre functions of the first kind. 1238 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``. 1243 1244 This function takes a real argument ``z``. For complex arguments ``z`` 1245 use clpmn instead. 1246 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. 1257 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 1264 1265 See Also 1266 -------- 1267 clpmn: associated Legendre functions of the first kind for complex z 1268 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. 1274 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 1282 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 1310 1311 1312def clpmn(m, n, z, type=3): 1313 """Associated Legendre function of the first kind for complex arguments. 1314 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``. 1319 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) 1334 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`` 1341 1342 See Also 1343 -------- 1344 lpmn: associated Legendre functions of the first kind for real z 1345 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`). 1353 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. 1357 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 1365 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 1391 1392 1393def lqmn(m, n, z): 1394 """Sequence of associated Legendre functions of the second kind. 1395 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``. 1400 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. 1411 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 1418 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 1424 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) 1434 1435 # Ensure neither m nor n == 0 1436 mm = max(1, m) 1437 nn = max(1, n) 1438 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)] 1444 1445 1446def bernoulli(n): 1447 """Bernoulli numbers B0..Bn (inclusive). 1448 1449 Parameters 1450 ---------- 1451 n : int 1452 Indicated the number of terms in the Bernoulli series to generate. 1453 1454 Returns 1455 ------- 1456 ndarray 1457 The Bernoulli numbers ``[B(0), B(1), ..., B(n)]``. 1458 1459 References 1460 ---------- 1461 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special 1462 Functions", John Wiley and Sons, 1996. 1463 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html 1464 .. [2] "Bernoulli number", Wikipedia, https://en.wikipedia.org/wiki/Bernoulli_number 1465 1466 Examples 1467 -------- 1468 >>> from scipy.special import bernoulli, zeta 1469 >>> bernoulli(4) 1470 array([ 1. , -0.5 , 0.16666667, 0. , -0.03333333]) 1471 1472 The Wikipedia article ([2]_) points out the relationship between the 1473 Bernoulli numbers and the zeta function, ``B_n^+ = -n * zeta(1 - n)`` 1474 for ``n > 0``: 1475 1476 >>> n = np.arange(1, 5) 1477 >>> -n * zeta(1 - n) 1478 array([ 0.5 , 0.16666667, -0. , -0.03333333]) 1479 1480 Note that, in the notation used in the wikipedia article, 1481 `bernoulli` computes ``B_n^-`` (i.e. it used the convention that 1482 ``B_1`` is -1/2). The relation given above is for ``B_n^+``, so the 1483 sign of 0.5 does not match the output of ``bernoulli(4)``. 1484 1485 """ 1486 if not isscalar(n) or (n < 0): 1487 raise ValueError("n must be a non-negative integer.") 1488 n = int(n) 1489 if (n < 2): 1490 n1 = 2 1491 else: 1492 n1 = n 1493 return specfun.bernob(int(n1))[:(n+1)] 1494 1495 1496def euler(n): 1497 """Euler numbers E(0), E(1), ..., E(n). 1498 1499 The Euler numbers [1]_ are also known as the secant numbers. 1500 1501 Because ``euler(n)`` returns floating point values, it does not give 1502 exact values for large `n`. The first inexact value is E(22). 1503 1504 Parameters 1505 ---------- 1506 n : int 1507 The highest index of the Euler number to be returned. 1508 1509 Returns 1510 ------- 1511 ndarray 1512 The Euler numbers [E(0), E(1), ..., E(n)]. 1513 The odd Euler numbers, which are all zero, are included. 1514 1515 References 1516 ---------- 1517 .. [1] Sequence A122045, The On-Line Encyclopedia of Integer Sequences, 1518 https://oeis.org/A122045 1519 .. [2] Zhang, Shanjie and Jin, Jianming. "Computation of Special 1520 Functions", John Wiley and Sons, 1996. 1521 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html 1522 1523 Examples 1524 -------- 1525 >>> from scipy.special import euler 1526 >>> euler(6) 1527 array([ 1., 0., -1., 0., 5., 0., -61.]) 1528 1529 >>> euler(13).astype(np.int64) 1530 array([ 1, 0, -1, 0, 5, 0, -61, 1531 0, 1385, 0, -50521, 0, 2702765, 0]) 1532 1533 >>> euler(22)[-1] # Exact value of E(22) is -69348874393137901. 1534 -69348874393137976.0 1535 1536 """ 1537 if not isscalar(n) or (n < 0): 1538 raise ValueError("n must be a non-negative integer.") 1539 n = int(n) 1540 if (n < 2): 1541 n1 = 2 1542 else: 1543 n1 = n 1544 return specfun.eulerb(n1)[:(n+1)] 1545 1546 1547def lpn(n, z): 1548 """Legendre function of the first kind. 1549 1550 Compute sequence of Legendre functions of the first kind (polynomials), 1551 Pn(z) and derivatives for all degrees from 0 to n (inclusive). 1552 1553 See also special.legendre for polynomial class. 1554 1555 References 1556 ---------- 1557 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special 1558 Functions", John Wiley and Sons, 1996. 1559 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html 1560 1561 """ 1562 if not (isscalar(n) and isscalar(z)): 1563 raise ValueError("arguments must be scalars.") 1564 n = _nonneg_int_or_fail(n, 'n', strict=False) 1565 if (n < 1): 1566 n1 = 1 1567 else: 1568 n1 = n 1569 if iscomplex(z): 1570 pn, pd = specfun.clpn(n1, z) 1571 else: 1572 pn, pd = specfun.lpn(n1, z) 1573 return pn[:(n+1)], pd[:(n+1)] 1574 1575 1576def lqn(n, z): 1577 """Legendre function of the second kind. 1578 1579 Compute sequence of Legendre functions of the second kind, Qn(z) and 1580 derivatives for all degrees from 0 to n (inclusive). 1581 1582 References 1583 ---------- 1584 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special 1585 Functions", John Wiley and Sons, 1996. 1586 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html 1587 1588 """ 1589 if not (isscalar(n) and isscalar(z)): 1590 raise ValueError("arguments must be scalars.") 1591 n = _nonneg_int_or_fail(n, 'n', strict=False) 1592 if (n < 1): 1593 n1 = 1 1594 else: 1595 n1 = n 1596 if iscomplex(z): 1597 qn, qd = specfun.clqn(n1, z) 1598 else: 1599 qn, qd = specfun.lqnb(n1, z) 1600 return qn[:(n+1)], qd[:(n+1)] 1601 1602 1603def ai_zeros(nt): 1604 """ 1605 Compute `nt` zeros and values of the Airy function Ai and its derivative. 1606 1607 Computes the first `nt` zeros, `a`, of the Airy function Ai(x); 1608 first `nt` zeros, `ap`, of the derivative of the Airy function Ai'(x); 1609 the corresponding values Ai(a'); 1610 and the corresponding values Ai'(a). 1611 1612 Parameters 1613 ---------- 1614 nt : int 1615 Number of zeros to compute 1616 1617 Returns 1618 ------- 1619 a : ndarray 1620 First `nt` zeros of Ai(x) 1621 ap : ndarray 1622 First `nt` zeros of Ai'(x) 1623 ai : ndarray 1624 Values of Ai(x) evaluated at first `nt` zeros of Ai'(x) 1625 aip : ndarray 1626 Values of Ai'(x) evaluated at first `nt` zeros of Ai(x) 1627 1628 Examples 1629 -------- 1630 >>> from scipy import special 1631 >>> a, ap, ai, aip = special.ai_zeros(3) 1632 >>> a 1633 array([-2.33810741, -4.08794944, -5.52055983]) 1634 >>> ap 1635 array([-1.01879297, -3.24819758, -4.82009921]) 1636 >>> ai 1637 array([ 0.53565666, -0.41901548, 0.38040647]) 1638 >>> aip 1639 array([ 0.70121082, -0.80311137, 0.86520403]) 1640 1641 References 1642 ---------- 1643 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special 1644 Functions", John Wiley and Sons, 1996. 1645 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html 1646 1647 """ 1648 kf = 1 1649 if not isscalar(nt) or (floor(nt) != nt) or (nt <= 0): 1650 raise ValueError("nt must be a positive integer scalar.") 1651 return specfun.airyzo(nt, kf) 1652 1653 1654def bi_zeros(nt): 1655 """ 1656 Compute `nt` zeros and values of the Airy function Bi and its derivative. 1657 1658 Computes the first `nt` zeros, b, of the Airy function Bi(x); 1659 first `nt` zeros, b', of the derivative of the Airy function Bi'(x); 1660 the corresponding values Bi(b'); 1661 and the corresponding values Bi'(b). 1662 1663 Parameters 1664 ---------- 1665 nt : int 1666 Number of zeros to compute 1667 1668 Returns 1669 ------- 1670 b : ndarray 1671 First `nt` zeros of Bi(x) 1672 bp : ndarray 1673 First `nt` zeros of Bi'(x) 1674 bi : ndarray 1675 Values of Bi(x) evaluated at first `nt` zeros of Bi'(x) 1676 bip : ndarray 1677 Values of Bi'(x) evaluated at first `nt` zeros of Bi(x) 1678 1679 Examples 1680 -------- 1681 >>> from scipy import special 1682 >>> b, bp, bi, bip = special.bi_zeros(3) 1683 >>> b 1684 array([-1.17371322, -3.2710933 , -4.83073784]) 1685 >>> bp 1686 array([-2.29443968, -4.07315509, -5.51239573]) 1687 >>> bi 1688 array([-0.45494438, 0.39652284, -0.36796916]) 1689 >>> bip 1690 array([ 0.60195789, -0.76031014, 0.83699101]) 1691 1692 References 1693 ---------- 1694 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special 1695 Functions", John Wiley and Sons, 1996. 1696 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html 1697 1698 """ 1699 kf = 2 1700 if not isscalar(nt) or (floor(nt) != nt) or (nt <= 0): 1701 raise ValueError("nt must be a positive integer scalar.") 1702 return specfun.airyzo(nt, kf) 1703 1704 1705def lmbda(v, x): 1706 r"""Jahnke-Emden Lambda function, Lambdav(x). 1707 1708 This function is defined as [2]_, 1709 1710 .. math:: \Lambda_v(x) = \Gamma(v+1) \frac{J_v(x)}{(x/2)^v}, 1711 1712 where :math:`\Gamma` is the gamma function and :math:`J_v` is the 1713 Bessel function of the first kind. 1714 1715 Parameters 1716 ---------- 1717 v : float 1718 Order of the Lambda function 1719 x : float 1720 Value at which to evaluate the function and derivatives 1721 1722 Returns 1723 ------- 1724 vl : ndarray 1725 Values of Lambda_vi(x), for vi=v-int(v), vi=1+v-int(v), ..., vi=v. 1726 dl : ndarray 1727 Derivatives Lambda_vi'(x), for vi=v-int(v), vi=1+v-int(v), ..., vi=v. 1728 1729 References 1730 ---------- 1731 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special 1732 Functions", John Wiley and Sons, 1996. 1733 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html 1734 .. [2] Jahnke, E. and Emde, F. "Tables of Functions with Formulae and 1735 Curves" (4th ed.), Dover, 1945 1736 """ 1737 if not (isscalar(v) and isscalar(x)): 1738 raise ValueError("arguments must be scalars.") 1739 if (v < 0): 1740 raise ValueError("argument must be > 0.") 1741 n = int(v) 1742 v0 = v - n 1743 if (n < 1): 1744 n1 = 1 1745 else: 1746 n1 = n 1747 v1 = n1 + v0 1748 if (v != floor(v)): 1749 vm, vl, dl = specfun.lamv(v1, x) 1750 else: 1751 vm, vl, dl = specfun.lamn(v1, x) 1752 return vl[:(n+1)], dl[:(n+1)] 1753 1754 1755def pbdv_seq(v, x): 1756 """Parabolic cylinder functions Dv(x) and derivatives. 1757 1758 Parameters 1759 ---------- 1760 v : float 1761 Order of the parabolic cylinder function 1762 x : float 1763 Value at which to evaluate the function and derivatives 1764 1765 Returns 1766 ------- 1767 dv : ndarray 1768 Values of D_vi(x), for vi=v-int(v), vi=1+v-int(v), ..., vi=v. 1769 dp : ndarray 1770 Derivatives D_vi'(x), for vi=v-int(v), vi=1+v-int(v), ..., vi=v. 1771 1772 References 1773 ---------- 1774 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special 1775 Functions", John Wiley and Sons, 1996, chapter 13. 1776 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html 1777 1778 """ 1779 if not (isscalar(v) and isscalar(x)): 1780 raise ValueError("arguments must be scalars.") 1781 n = int(v) 1782 v0 = v-n 1783 if (n < 1): 1784 n1 = 1 1785 else: 1786 n1 = n 1787 v1 = n1 + v0 1788 dv, dp, pdf, pdd = specfun.pbdv(v1, x) 1789 return dv[:n1+1], dp[:n1+1] 1790 1791 1792def pbvv_seq(v, x): 1793 """Parabolic cylinder functions Vv(x) and derivatives. 1794 1795 Parameters 1796 ---------- 1797 v : float 1798 Order of the parabolic cylinder function 1799 x : float 1800 Value at which to evaluate the function and derivatives 1801 1802 Returns 1803 ------- 1804 dv : ndarray 1805 Values of V_vi(x), for vi=v-int(v), vi=1+v-int(v), ..., vi=v. 1806 dp : ndarray 1807 Derivatives V_vi'(x), for vi=v-int(v), vi=1+v-int(v), ..., vi=v. 1808 1809 References 1810 ---------- 1811 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special 1812 Functions", John Wiley and Sons, 1996, chapter 13. 1813 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html 1814 1815 """ 1816 if not (isscalar(v) and isscalar(x)): 1817 raise ValueError("arguments must be scalars.") 1818 n = int(v) 1819 v0 = v-n 1820 if (n <= 1): 1821 n1 = 1 1822 else: 1823 n1 = n 1824 v1 = n1 + v0 1825 dv, dp, pdf, pdd = specfun.pbvv(v1, x) 1826 return dv[:n1+1], dp[:n1+1] 1827 1828 1829def pbdn_seq(n, z): 1830 """Parabolic cylinder functions Dn(z) and derivatives. 1831 1832 Parameters 1833 ---------- 1834 n : int 1835 Order of the parabolic cylinder function 1836 z : complex 1837 Value at which to evaluate the function and derivatives 1838 1839 Returns 1840 ------- 1841 dv : ndarray 1842 Values of D_i(z), for i=0, ..., i=n. 1843 dp : ndarray 1844 Derivatives D_i'(z), for i=0, ..., i=n. 1845 1846 References 1847 ---------- 1848 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special 1849 Functions", John Wiley and Sons, 1996, chapter 13. 1850 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html 1851 1852 """ 1853 if not (isscalar(n) and isscalar(z)): 1854 raise ValueError("arguments must be scalars.") 1855 if (floor(n) != n): 1856 raise ValueError("n must be an integer.") 1857 if (abs(n) <= 1): 1858 n1 = 1 1859 else: 1860 n1 = n 1861 cpb, cpd = specfun.cpbdn(n1, z) 1862 return cpb[:n1+1], cpd[:n1+1] 1863 1864 1865def ber_zeros(nt): 1866 """Compute nt zeros of the Kelvin function ber. 1867 1868 Parameters 1869 ---------- 1870 nt : int 1871 Number of zeros to compute. Must be positive. 1872 1873 Returns 1874 ------- 1875 ndarray 1876 First `nt` zeros of the Kelvin function. 1877 1878 See Also 1879 -------- 1880 ber 1881 1882 References 1883 ---------- 1884 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special 1885 Functions", John Wiley and Sons, 1996. 1886 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html 1887 1888 """ 1889 if not isscalar(nt) or (floor(nt) != nt) or (nt <= 0): 1890 raise ValueError("nt must be positive integer scalar.") 1891 return specfun.klvnzo(nt, 1) 1892 1893 1894def bei_zeros(nt): 1895 """Compute nt zeros of the Kelvin function bei. 1896 1897 Parameters 1898 ---------- 1899 nt : int 1900 Number of zeros to compute. Must be positive. 1901 1902 Returns 1903 ------- 1904 ndarray 1905 First `nt` zeros of the Kelvin function. 1906 1907 See Also 1908 -------- 1909 bei 1910 1911 References 1912 ---------- 1913 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special 1914 Functions", John Wiley and Sons, 1996. 1915 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html 1916 1917 """ 1918 if not isscalar(nt) or (floor(nt) != nt) or (nt <= 0): 1919 raise ValueError("nt must be positive integer scalar.") 1920 return specfun.klvnzo(nt, 2) 1921 1922 1923def ker_zeros(nt): 1924 """Compute nt zeros of the Kelvin function ker. 1925 1926 Parameters 1927 ---------- 1928 nt : int 1929 Number of zeros to compute. Must be positive. 1930 1931 Returns 1932 ------- 1933 ndarray 1934 First `nt` zeros of the Kelvin function. 1935 1936 See Also 1937 -------- 1938 ker 1939 1940 References 1941 ---------- 1942 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special 1943 Functions", John Wiley and Sons, 1996. 1944 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html 1945 1946 """ 1947 if not isscalar(nt) or (floor(nt) != nt) or (nt <= 0): 1948 raise ValueError("nt must be positive integer scalar.") 1949 return specfun.klvnzo(nt, 3) 1950 1951 1952def kei_zeros(nt): 1953 """Compute nt zeros of the Kelvin function kei. 1954 1955 Parameters 1956 ---------- 1957 nt : int 1958 Number of zeros to compute. Must be positive. 1959 1960 Returns 1961 ------- 1962 ndarray 1963 First `nt` zeros of the Kelvin function. 1964 1965 See Also 1966 -------- 1967 kei 1968 1969 References 1970 ---------- 1971 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special 1972 Functions", John Wiley and Sons, 1996. 1973 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html 1974 1975 """ 1976 if not isscalar(nt) or (floor(nt) != nt) or (nt <= 0): 1977 raise ValueError("nt must be positive integer scalar.") 1978 return specfun.klvnzo(nt, 4) 1979 1980 1981def berp_zeros(nt): 1982 """Compute nt zeros of the derivative of the Kelvin function ber. 1983 1984 Parameters 1985 ---------- 1986 nt : int 1987 Number of zeros to compute. Must be positive. 1988 1989 Returns 1990 ------- 1991 ndarray 1992 First `nt` zeros of the derivative of the Kelvin function. 1993 1994 See Also 1995 -------- 1996 ber, berp 1997 1998 References 1999 ---------- 2000 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special 2001 Functions", John Wiley and Sons, 1996. 2002 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html 2003 2004 """ 2005 if not isscalar(nt) or (floor(nt) != nt) or (nt <= 0): 2006 raise ValueError("nt must be positive integer scalar.") 2007 return specfun.klvnzo(nt, 5) 2008 2009 2010def beip_zeros(nt): 2011 """Compute nt zeros of the derivative of the Kelvin function bei. 2012 2013 Parameters 2014 ---------- 2015 nt : int 2016 Number of zeros to compute. Must be positive. 2017 2018 Returns 2019 ------- 2020 ndarray 2021 First `nt` zeros of the derivative of the Kelvin function. 2022 2023 See Also 2024 -------- 2025 bei, beip 2026 2027 References 2028 ---------- 2029 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special 2030 Functions", John Wiley and Sons, 1996. 2031 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html 2032 2033 """ 2034 if not isscalar(nt) or (floor(nt) != nt) or (nt <= 0): 2035 raise ValueError("nt must be positive integer scalar.") 2036 return specfun.klvnzo(nt, 6) 2037 2038 2039def kerp_zeros(nt): 2040 """Compute nt zeros of the derivative of the Kelvin function ker. 2041 2042 Parameters 2043 ---------- 2044 nt : int 2045 Number of zeros to compute. Must be positive. 2046 2047 Returns 2048 ------- 2049 ndarray 2050 First `nt` zeros of the derivative of the Kelvin function. 2051 2052 See Also 2053 -------- 2054 ker, kerp 2055 2056 References 2057 ---------- 2058 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special 2059 Functions", John Wiley and Sons, 1996. 2060 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html 2061 2062 """ 2063 if not isscalar(nt) or (floor(nt) != nt) or (nt <= 0): 2064 raise ValueError("nt must be positive integer scalar.") 2065 return specfun.klvnzo(nt, 7) 2066 2067 2068def keip_zeros(nt): 2069 """Compute nt zeros of the derivative of the Kelvin function kei. 2070 2071 Parameters 2072 ---------- 2073 nt : int 2074 Number of zeros to compute. Must be positive. 2075 2076 Returns 2077 ------- 2078 ndarray 2079 First `nt` zeros of the derivative of the Kelvin function. 2080 2081 See Also 2082 -------- 2083 kei, keip 2084 2085 References 2086 ---------- 2087 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special 2088 Functions", John Wiley and Sons, 1996. 2089 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html 2090 2091 """ 2092 if not isscalar(nt) or (floor(nt) != nt) or (nt <= 0): 2093 raise ValueError("nt must be positive integer scalar.") 2094 return specfun.klvnzo(nt, 8) 2095 2096 2097def kelvin_zeros(nt): 2098 """Compute nt zeros of all Kelvin functions. 2099 2100 Returned in a length-8 tuple of arrays of length nt. The tuple contains 2101 the arrays of zeros of (ber, bei, ker, kei, ber', bei', ker', kei'). 2102 2103 References 2104 ---------- 2105 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special 2106 Functions", John Wiley and Sons, 1996. 2107 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html 2108 2109 """ 2110 if not isscalar(nt) or (floor(nt) != nt) or (nt <= 0): 2111 raise ValueError("nt must be positive integer scalar.") 2112 return (specfun.klvnzo(nt, 1), 2113 specfun.klvnzo(nt, 2), 2114 specfun.klvnzo(nt, 3), 2115 specfun.klvnzo(nt, 4), 2116 specfun.klvnzo(nt, 5), 2117 specfun.klvnzo(nt, 6), 2118 specfun.klvnzo(nt, 7), 2119 specfun.klvnzo(nt, 8)) 2120 2121 2122def pro_cv_seq(m, n, c): 2123 """Characteristic values for prolate spheroidal wave functions. 2124 2125 Compute a sequence of characteristic values for the prolate 2126 spheroidal wave functions for mode m and n'=m..n and spheroidal 2127 parameter c. 2128 2129 References 2130 ---------- 2131 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special 2132 Functions", John Wiley and Sons, 1996. 2133 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html 2134 2135 """ 2136 if not (isscalar(m) and isscalar(n) and isscalar(c)): 2137 raise ValueError("Arguments must be scalars.") 2138 if (n != floor(n)) or (m != floor(m)): 2139 raise ValueError("Modes must be integers.") 2140 if (n-m > 199): 2141 raise ValueError("Difference between n and m is too large.") 2142 maxL = n-m+1 2143 return specfun.segv(m, n, c, 1)[1][:maxL] 2144 2145 2146def obl_cv_seq(m, n, c): 2147 """Characteristic values for oblate spheroidal wave functions. 2148 2149 Compute a sequence of characteristic values for the oblate 2150 spheroidal wave functions for mode m and n'=m..n and spheroidal 2151 parameter c. 2152 2153 References 2154 ---------- 2155 .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special 2156 Functions", John Wiley and Sons, 1996. 2157 https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html 2158 2159 """ 2160 if not (isscalar(m) and isscalar(n) and isscalar(c)): 2161 raise ValueError("Arguments must be scalars.") 2162 if (n != floor(n)) or (m != floor(m)): 2163 raise ValueError("Modes must be integers.") 2164 if (n-m > 199): 2165 raise ValueError("Difference between n and m is too large.") 2166 maxL = n-m+1 2167 return specfun.segv(m, n, c, -1)[1][:maxL] 2168 2169 2170def comb(N, k, exact=False, repetition=False): 2171 """The number of combinations of N things taken k at a time. 2172 2173 This is often expressed as "N choose k". 2174 2175 Parameters 2176 ---------- 2177 N : int, ndarray 2178 Number of things. 2179 k : int, ndarray 2180 Number of elements taken. 2181 exact : bool, optional 2182 If `exact` is False, then floating point precision is used, otherwise 2183 exact long integer is computed. 2184 repetition : bool, optional 2185 If `repetition` is True, then the number of combinations with 2186 repetition is computed. 2187 2188 Returns 2189 ------- 2190 val : int, float, ndarray 2191 The total number of combinations. 2192 2193 See Also 2194 -------- 2195 binom : Binomial coefficient ufunc 2196 2197 Notes 2198 ----- 2199 - Array arguments accepted only for exact=False case. 2200 - If N < 0, or k < 0, then 0 is returned. 2201 - If k > N and repetition=False, then 0 is returned. 2202 2203 Examples 2204 -------- 2205 >>> from scipy.special import comb 2206 >>> k = np.array([3, 4]) 2207 >>> n = np.array([10, 10]) 2208 >>> comb(n, k, exact=False) 2209 array([ 120., 210.]) 2210 >>> comb(10, 3, exact=True) 2211 120 2212 >>> comb(10, 3, exact=True, repetition=True) 2213 220 2214 2215 """ 2216 if repetition: 2217 return comb(N + k - 1, k, exact) 2218 if exact: 2219 return _comb_int(N, k) 2220 else: 2221 k, N = asarray(k), asarray(N) 2222 cond = (k <= N) & (N >= 0) & (k >= 0) 2223 vals = binom(N, k) 2224 if isinstance(vals, np.ndarray): 2225 vals[~cond] = 0 2226 elif not cond: 2227 vals = np.float64(0) 2228 return vals 2229 2230 2231def perm(N, k, exact=False): 2232 """Permutations of N things taken k at a time, i.e., k-permutations of N. 2233 2234 It's also known as "partial permutations". 2235 2236 Parameters 2237 ---------- 2238 N : int, ndarray 2239 Number of things. 2240 k : int, ndarray 2241 Number of elements taken. 2242 exact : bool, optional 2243 If `exact` is False, then floating point precision is used, otherwise 2244 exact long integer is computed. 2245 2246 Returns 2247 ------- 2248 val : int, ndarray 2249 The number of k-permutations of N. 2250 2251 Notes 2252 ----- 2253 - Array arguments accepted only for exact=False case. 2254 - If k > N, N < 0, or k < 0, then a 0 is returned. 2255 2256 Examples 2257 -------- 2258 >>> from scipy.special import perm 2259 >>> k = np.array([3, 4]) 2260 >>> n = np.array([10, 10]) 2261 >>> perm(n, k) 2262 array([ 720., 5040.]) 2263 >>> perm(10, 3, exact=True) 2264 720 2265 2266 """ 2267 if exact: 2268 if (k > N) or (N < 0) or (k < 0): 2269 return 0 2270 val = 1 2271 for i in range(N - k + 1, N + 1): 2272 val *= i 2273 return val 2274 else: 2275 k, N = asarray(k), asarray(N) 2276 cond = (k <= N) & (N >= 0) & (k >= 0) 2277 vals = poch(N - k + 1, k) 2278 if isinstance(vals, np.ndarray): 2279 vals[~cond] = 0 2280 elif not cond: 2281 vals = np.float64(0) 2282 return vals 2283 2284 2285# https://stackoverflow.com/a/16327037 2286def _range_prod(lo, hi): 2287 """ 2288 Product of a range of numbers. 2289 2290 Returns the product of 2291 lo * (lo+1) * (lo+2) * ... * (hi-2) * (hi-1) * hi 2292 = hi! / (lo-1)! 2293 2294 Breaks into smaller products first for speed: 2295 _range_prod(2, 9) = ((2*3)*(4*5))*((6*7)*(8*9)) 2296 """ 2297 if lo + 1 < hi: 2298 mid = (hi + lo) // 2 2299 return _range_prod(lo, mid) * _range_prod(mid + 1, hi) 2300 if lo == hi: 2301 return lo 2302 return lo * hi 2303 2304 2305def factorial(n, exact=False): 2306 """ 2307 The factorial of a number or array of numbers. 2308 2309 The factorial of non-negative integer `n` is the product of all 2310 positive integers less than or equal to `n`:: 2311 2312 n! = n * (n - 1) * (n - 2) * ... * 1 2313 2314 Parameters 2315 ---------- 2316 n : int or array_like of ints 2317 Input values. If ``n < 0``, the return value is 0. 2318 exact : bool, optional 2319 If True, calculate the answer exactly using long integer arithmetic. 2320 If False, result is approximated in floating point rapidly using the 2321 `gamma` function. 2322 Default is False. 2323 2324 Returns 2325 ------- 2326 nf : float or int or ndarray 2327 Factorial of `n`, as integer or float depending on `exact`. 2328 2329 Notes 2330 ----- 2331 For arrays with ``exact=True``, the factorial is computed only once, for 2332 the largest input, with each other result computed in the process. 2333 The output dtype is increased to ``int64`` or ``object`` if necessary. 2334 2335 With ``exact=False`` the factorial is approximated using the gamma 2336 function: 2337 2338 .. math:: n! = \\Gamma(n+1) 2339 2340 Examples 2341 -------- 2342 >>> from scipy.special import factorial 2343 >>> arr = np.array([3, 4, 5]) 2344 >>> factorial(arr, exact=False) 2345 array([ 6., 24., 120.]) 2346 >>> factorial(arr, exact=True) 2347 array([ 6, 24, 120]) 2348 >>> factorial(5, exact=True) 2349 120 2350 2351 """ 2352 if exact: 2353 if np.ndim(n) == 0: 2354 if np.isnan(n): 2355 return n 2356 return 0 if n < 0 else math.factorial(n) 2357 else: 2358 n = asarray(n) 2359 un = np.unique(n).astype(object) 2360 2361 # Convert to object array of long ints if np.int_ can't handle size 2362 if np.isnan(n).any(): 2363 dt = float 2364 elif un[-1] > 20: 2365 dt = object 2366 elif un[-1] > 12: 2367 dt = np.int64 2368 else: 2369 dt = np.int_ 2370 2371 out = np.empty_like(n, dtype=dt) 2372 2373 # Handle invalid/trivial values 2374 # Ignore runtime warning when less operator used w/np.nan 2375 with np.errstate(all='ignore'): 2376 un = un[un > 1] 2377 out[n < 2] = 1 2378 out[n < 0] = 0 2379 2380 # Calculate products of each range of numbers 2381 if un.size: 2382 val = math.factorial(un[0]) 2383 out[n == un[0]] = val 2384 for i in range(len(un) - 1): 2385 prev = un[i] + 1 2386 current = un[i + 1] 2387 val *= _range_prod(prev, current) 2388 out[n == current] = val 2389 2390 if np.isnan(n).any(): 2391 out = out.astype(np.float64) 2392 out[np.isnan(n)] = n[np.isnan(n)] 2393 return out 2394 else: 2395 out = ufuncs._factorial(n) 2396 return out 2397 2398 2399def factorial2(n, exact=False): 2400 """Double factorial. 2401 2402 This is the factorial with every second value skipped. E.g., ``7!! = 7 * 5 2403 * 3 * 1``. It can be approximated numerically as:: 2404 2405 n!! = special.gamma(n/2+1)*2**((m+1)/2)/sqrt(pi) n odd 2406 = 2**(n/2) * (n/2)! n even 2407 2408 Parameters 2409 ---------- 2410 n : int or array_like 2411 Calculate ``n!!``. Arrays are only supported with `exact` set 2412 to False. If ``n < 0``, the return value is 0. 2413 exact : bool, optional 2414 The result can be approximated rapidly using the gamma-formula 2415 above (default). If `exact` is set to True, calculate the 2416 answer exactly using integer arithmetic. 2417 2418 Returns 2419 ------- 2420 nff : float or int 2421 Double factorial of `n`, as an int or a float depending on 2422 `exact`. 2423 2424 Examples 2425 -------- 2426 >>> from scipy.special import factorial2 2427 >>> factorial2(7, exact=False) 2428 array(105.00000000000001) 2429 >>> factorial2(7, exact=True) 2430 105 2431 2432 """ 2433 if exact: 2434 if n < -1: 2435 return 0 2436 if n <= 0: 2437 return 1 2438 val = 1 2439 for k in range(n, 0, -2): 2440 val *= k 2441 return val 2442 else: 2443 n = asarray(n) 2444 vals = zeros(n.shape, 'd') 2445 cond1 = (n % 2) & (n >= -1) 2446 cond2 = (1-(n % 2)) & (n >= -1) 2447 oddn = extract(cond1, n) 2448 evenn = extract(cond2, n) 2449 nd2o = oddn / 2.0 2450 nd2e = evenn / 2.0 2451 place(vals, cond1, gamma(nd2o + 1) / sqrt(pi) * pow(2.0, nd2o + 0.5)) 2452 place(vals, cond2, gamma(nd2e + 1) * pow(2.0, nd2e)) 2453 return vals 2454 2455 2456def factorialk(n, k, exact=True): 2457 """Multifactorial of n of order k, n(!!...!). 2458 2459 This is the multifactorial of n skipping k values. For example, 2460 2461 factorialk(17, 4) = 17!!!! = 17 * 13 * 9 * 5 * 1 2462 2463 In particular, for any integer ``n``, we have 2464 2465 factorialk(n, 1) = factorial(n) 2466 2467 factorialk(n, 2) = factorial2(n) 2468 2469 Parameters 2470 ---------- 2471 n : int 2472 Calculate multifactorial. If `n` < 0, the return value is 0. 2473 k : int 2474 Order of multifactorial. 2475 exact : bool, optional 2476 If exact is set to True, calculate the answer exactly using 2477 integer arithmetic. 2478 2479 Returns 2480 ------- 2481 val : int 2482 Multifactorial of `n`. 2483 2484 Raises 2485 ------ 2486 NotImplementedError 2487 Raises when exact is False 2488 2489 Examples 2490 -------- 2491 >>> from scipy.special import factorialk 2492 >>> factorialk(5, 1, exact=True) 2493 120 2494 >>> factorialk(5, 3, exact=True) 2495 10 2496 2497 """ 2498 if exact: 2499 if n < 1-k: 2500 return 0 2501 if n <= 0: 2502 return 1 2503 val = 1 2504 for j in range(n, 0, -k): 2505 val = val*j 2506 return val 2507 else: 2508 raise NotImplementedError 2509 2510 2511def zeta(x, q=None, out=None): 2512 r""" 2513 Riemann or Hurwitz zeta function. 2514 2515 Parameters 2516 ---------- 2517 x : array_like of float 2518 Input data, must be real 2519 q : array_like of float, optional 2520 Input data, must be real. Defaults to Riemann zeta. 2521 out : ndarray, optional 2522 Output array for the computed values. 2523 2524 Returns 2525 ------- 2526 out : array_like 2527 Values of zeta(x). 2528 2529 Notes 2530 ----- 2531 The two-argument version is the Hurwitz zeta function 2532 2533 .. math:: 2534 2535 \zeta(x, q) = \sum_{k=0}^{\infty} \frac{1}{(k + q)^x}; 2536 2537 see [dlmf]_ for details. The Riemann zeta function corresponds to 2538 the case when ``q = 1``. 2539 2540 See Also 2541 -------- 2542 zetac 2543 2544 References 2545 ---------- 2546 .. [dlmf] NIST, Digital Library of Mathematical Functions, 2547 https://dlmf.nist.gov/25.11#i 2548 2549 Examples 2550 -------- 2551 >>> from scipy.special import zeta, polygamma, factorial 2552 2553 Some specific values: 2554 2555 >>> zeta(2), np.pi**2/6 2556 (1.6449340668482266, 1.6449340668482264) 2557 2558 >>> zeta(4), np.pi**4/90 2559 (1.0823232337111381, 1.082323233711138) 2560 2561 Relation to the `polygamma` function: 2562 2563 >>> m = 3 2564 >>> x = 1.25 2565 >>> polygamma(m, x) 2566 array(2.782144009188397) 2567 >>> (-1)**(m+1) * factorial(m) * zeta(m+1, x) 2568 2.7821440091883969 2569 2570 """ 2571 if q is None: 2572 return ufuncs._riemann_zeta(x, out) 2573 else: 2574 return ufuncs._zeta(x, q, out) 2575