Hide keyboard shortcuts

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# 

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)) 

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 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 

1454 

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)] 

1464 

1465 

1466def euler(n): 

1467 """Euler numbers E(0), E(1), ..., E(n). 

1468 

1469 The Euler numbers [1]_ are also known as the secant numbers. 

1470 

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). 

1473 

1474 Parameters 

1475 ---------- 

1476 n : int 

1477 The highest index of the Euler number to be returned. 

1478 

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. 

1484 

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 

1492 

1493 Examples 

1494 -------- 

1495 >>> from scipy.special import euler 

1496 >>> euler(6) 

1497 array([ 1., 0., -1., 0., 5., 0., -61.]) 

1498 

1499 >>> euler(13).astype(np.int64) 

1500 array([ 1, 0, -1, 0, 5, 0, -61, 

1501 0, 1385, 0, -50521, 0, 2702765, 0]) 

1502 

1503 >>> euler(22)[-1] # Exact value of E(22) is -69348874393137901. 

1504 -69348874393137976.0 

1505 

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)] 

1515 

1516 

1517def lpn(n, z): 

1518 """Legendre function of the first kind. 

1519 

1520 Compute sequence of Legendre functions of the first kind (polynomials), 

1521 Pn(z) and derivatives for all degrees from 0 to n (inclusive). 

1522 

1523 See also special.legendre for polynomial class. 

1524 

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 

1530 

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)] 

1544 

1545 

1546def lqn(n, z): 

1547 """Legendre function of the second kind. 

1548 

1549 Compute sequence of Legendre functions of the second kind, Qn(z) and 

1550 derivatives for all degrees from 0 to n (inclusive). 

1551 

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 

1557 

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)] 

1571 

1572 

1573def ai_zeros(nt): 

1574 """ 

1575 Compute `nt` zeros and values of the Airy function Ai and its derivative. 

1576 

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). 

1581 

1582 Parameters 

1583 ---------- 

1584 nt : int 

1585 Number of zeros to compute 

1586 

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) 

1597 

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]) 

1610 

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 

1616 

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) 

1622 

1623 

1624def bi_zeros(nt): 

1625 """ 

1626 Compute `nt` zeros and values of the Airy function Bi and its derivative. 

1627 

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). 

1632 

1633 Parameters 

1634 ---------- 

1635 nt : int 

1636 Number of zeros to compute 

1637 

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) 

1648 

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]) 

1661 

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 

1667 

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) 

1673 

1674 

1675def lmbda(v, x): 

1676 r"""Jahnke-Emden Lambda function, Lambdav(x). 

1677 

1678 This function is defined as [2]_, 

1679 

1680 .. math:: \Lambda_v(x) = \Gamma(v+1) \frac{J_v(x)}{(x/2)^v}, 

1681 

1682 where :math:`\Gamma` is the gamma function and :math:`J_v` is the 

1683 Bessel function of the first kind. 

1684 

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 

1691 

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. 

1698 

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)] 

1723 

1724 

1725def pbdv_seq(v, x): 

1726 """Parabolic cylinder functions Dv(x) and derivatives. 

1727 

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 

1734 

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. 

1741 

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 

1747 

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] 

1760 

1761 

1762def pbvv_seq(v, x): 

1763 """Parabolic cylinder functions Vv(x) and derivatives. 

1764 

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 

1771 

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. 

1778 

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 

1784 

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] 

1797 

1798 

1799def pbdn_seq(n, z): 

1800 """Parabolic cylinder functions Dn(z) and derivatives. 

1801 

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 

1808 

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. 

1815 

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 

1821 

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] 

1833 

1834 

1835def ber_zeros(nt): 

1836 """Compute nt zeros of the Kelvin function ber. 

1837 

1838 Parameters 

1839 ---------- 

1840 nt : int 

1841 Number of zeros to compute. Must be positive. 

1842 

1843 Returns 

1844 ------- 

1845 ndarray 

1846 First `nt` zeros of the Kelvin function. 

1847 

1848 See Also 

1849 -------- 

1850 ber 

1851 

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 

1857 

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) 

1862 

1863 

1864def bei_zeros(nt): 

1865 """Compute nt zeros of the Kelvin function bei. 

1866 

1867 Parameters 

1868 ---------- 

1869 nt : int 

1870 Number of zeros to compute. Must be positive. 

1871 

1872 Returns 

1873 ------- 

1874 ndarray 

1875 First `nt` zeros of the Kelvin function. 

1876 

1877 See Also 

1878 -------- 

1879 bei 

1880 

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 

1886 

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) 

1891 

1892 

1893def ker_zeros(nt): 

1894 """Compute nt zeros of the Kelvin function ker. 

1895 

1896 Parameters 

1897 ---------- 

1898 nt : int 

1899 Number of zeros to compute. Must be positive. 

1900 

1901 Returns 

1902 ------- 

1903 ndarray 

1904 First `nt` zeros of the Kelvin function. 

1905 

1906 See Also 

1907 -------- 

1908 ker 

1909 

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 

1915 

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) 

1920 

1921 

1922def kei_zeros(nt): 

1923 """Compute nt zeros of the Kelvin function kei. 

1924 

1925 Parameters 

1926 ---------- 

1927 nt : int 

1928 Number of zeros to compute. Must be positive. 

1929 

1930 Returns 

1931 ------- 

1932 ndarray 

1933 First `nt` zeros of the Kelvin function. 

1934 

1935 See Also 

1936 -------- 

1937 kei 

1938 

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 

1944 

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) 

1949 

1950 

1951def berp_zeros(nt): 

1952 """Compute nt zeros of the derivative of the Kelvin function ber. 

1953 

1954 Parameters 

1955 ---------- 

1956 nt : int 

1957 Number of zeros to compute. Must be positive. 

1958 

1959 Returns 

1960 ------- 

1961 ndarray 

1962 First `nt` zeros of the derivative of the Kelvin function. 

1963 

1964 See Also 

1965 -------- 

1966 ber, berp 

1967 

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 

1973 

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) 

1978 

1979 

1980def beip_zeros(nt): 

1981 """Compute nt zeros of the derivative of the Kelvin function bei. 

1982 

1983 Parameters 

1984 ---------- 

1985 nt : int 

1986 Number of zeros to compute. Must be positive. 

1987 

1988 Returns 

1989 ------- 

1990 ndarray 

1991 First `nt` zeros of the derivative of the Kelvin function. 

1992 

1993 See Also 

1994 -------- 

1995 bei, beip 

1996 

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 

2002 

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) 

2007 

2008 

2009def kerp_zeros(nt): 

2010 """Compute nt zeros of the derivative of the Kelvin function ker. 

2011 

2012 Parameters 

2013 ---------- 

2014 nt : int 

2015 Number of zeros to compute. Must be positive. 

2016 

2017 Returns 

2018 ------- 

2019 ndarray 

2020 First `nt` zeros of the derivative of the Kelvin function. 

2021 

2022 See Also 

2023 -------- 

2024 ker, kerp 

2025 

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 

2031 

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) 

2036 

2037 

2038def keip_zeros(nt): 

2039 """Compute nt zeros of the derivative of the Kelvin function kei. 

2040 

2041 Parameters 

2042 ---------- 

2043 nt : int 

2044 Number of zeros to compute. Must be positive. 

2045 

2046 Returns 

2047 ------- 

2048 ndarray 

2049 First `nt` zeros of the derivative of the Kelvin function. 

2050 

2051 See Also 

2052 -------- 

2053 kei, keip 

2054 

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 

2060 

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) 

2065 

2066 

2067def kelvin_zeros(nt): 

2068 """Compute nt zeros of all Kelvin functions. 

2069 

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'). 

2072 

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 

2078 

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)) 

2090 

2091 

2092def pro_cv_seq(m, n, c): 

2093 """Characteristic values for prolate spheroidal wave functions. 

2094 

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. 

2098 

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 

2104 

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] 

2114 

2115 

2116def obl_cv_seq(m, n, c): 

2117 """Characteristic values for oblate spheroidal wave functions. 

2118 

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. 

2122 

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 

2128 

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] 

2138 

2139 

2140def comb(N, k, exact=False, repetition=False): 

2141 """The number of combinations of N things taken k at a time. 

2142 

2143 This is often expressed as "N choose k". 

2144 

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. 

2157 

2158 Returns 

2159 ------- 

2160 val : int, float, ndarray 

2161 The total number of combinations. 

2162 

2163 See Also 

2164 -------- 

2165 binom : Binomial coefficient ufunc 

2166 

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. 

2172 

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 

2184 

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 

2199 

2200 

2201def perm(N, k, exact=False): 

2202 """Permutations of N things taken k at a time, i.e., k-permutations of N. 

2203 

2204 It's also known as "partial permutations". 

2205 

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. 

2215 

2216 Returns 

2217 ------- 

2218 val : int, ndarray 

2219 The number of k-permutations of N. 

2220 

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. 

2225 

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 

2235 

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 

2253 

2254 

2255# https://stackoverflow.com/a/16327037 

2256def _range_prod(lo, hi): 

2257 """ 

2258 Product of a range of numbers. 

2259 

2260 Returns the product of 

2261 lo * (lo+1) * (lo+2) * ... * (hi-2) * (hi-1) * hi 

2262 = hi! / (lo-1)! 

2263 

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 

2273 

2274 

2275def factorial(n, exact=False): 

2276 """ 

2277 The factorial of a number or array of numbers. 

2278 

2279 The factorial of non-negative integer `n` is the product of all 

2280 positive integers less than or equal to `n`:: 

2281 

2282 n! = n * (n - 1) * (n - 2) * ... * 1 

2283 

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. 

2293 

2294 Returns 

2295 ------- 

2296 nf : float or int or ndarray 

2297 Factorial of `n`, as integer or float depending on `exact`. 

2298 

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. 

2304 

2305 With ``exact=False`` the factorial is approximated using the gamma 

2306 function: 

2307 

2308 .. math:: n! = \\Gamma(n+1) 

2309 

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 

2320 

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) 

2330 

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 

2340 

2341 out = np.empty_like(n, dtype=dt) 

2342 

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 

2349 

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 

2359 

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 

2367 

2368 

2369def factorial2(n, exact=False): 

2370 """Double factorial. 

2371 

2372 This is the factorial with every second value skipped. E.g., ``7!! = 7 * 5 

2373 * 3 * 1``. It can be approximated numerically as:: 

2374 

2375 n!! = special.gamma(n/2+1)*2**((m+1)/2)/sqrt(pi) n odd 

2376 = 2**(n/2) * (n/2)! n even 

2377 

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. 

2387 

2388 Returns 

2389 ------- 

2390 nff : float or int 

2391 Double factorial of `n`, as an int or a float depending on 

2392 `exact`. 

2393 

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 

2401 

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 

2424 

2425 

2426def factorialk(n, k, exact=True): 

2427 """Multifactorial of n of order k, n(!!...!). 

2428 

2429 This is the multifactorial of n skipping k values. For example, 

2430 

2431 factorialk(17, 4) = 17!!!! = 17 * 13 * 9 * 5 * 1 

2432 

2433 In particular, for any integer ``n``, we have 

2434 

2435 factorialk(n, 1) = factorial(n) 

2436 

2437 factorialk(n, 2) = factorial2(n) 

2438 

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. 

2448 

2449 Returns 

2450 ------- 

2451 val : int 

2452 Multifactorial of `n`. 

2453 

2454 Raises 

2455 ------ 

2456 NotImplementedError 

2457 Raises when exact is False 

2458 

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 

2466 

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 

2479 

2480 

2481def zeta(x, q=None, out=None): 

2482 r""" 

2483 Riemann or Hurwitz zeta function. 

2484 

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. 

2493 

2494 Returns 

2495 ------- 

2496 out : array_like 

2497 Values of zeta(x). 

2498 

2499 Notes 

2500 ----- 

2501 The two-argument version is the Hurwitz zeta function 

2502 

2503 .. math:: 

2504 

2505 \zeta(x, q) = \sum_{k=0}^{\infty} \frac{1}{(k + q)^x}; 

2506 

2507 see [dlmf]_ for details. The Riemann zeta function corresponds to 

2508 the case when ``q = 1``. 

2509 

2510 See Also 

2511 -------- 

2512 zetac 

2513 

2514 References 

2515 ---------- 

2516 .. [dlmf] NIST, Digital Library of Mathematical Functions, 

2517 https://dlmf.nist.gov/25.11#i 

2518 

2519 Examples 

2520 -------- 

2521 >>> from scipy.special import zeta, polygamma, factorial 

2522 

2523 Some specific values: 

2524 

2525 >>> zeta(2), np.pi**2/6 

2526 (1.6449340668482266, 1.6449340668482264) 

2527 

2528 >>> zeta(4), np.pi**4/90 

2529 (1.0823232337111381, 1.082323233711138) 

2530 

2531 Relation to the `polygamma` function: 

2532 

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 

2539 

2540 """ 

2541 if q is None: 

2542 return ufuncs._riemann_zeta(x, out) 

2543 else: 

2544 return ufuncs._zeta(x, q, out)