1# this program corresponds to special.py 2 3### Means test is not done yet 4# E Means test is giving error (E) 5# F Means test is failing (F) 6# EF Means test is giving error and Failing 7#! Means test is segfaulting 8# 8 Means test runs forever 9 10### test_besselpoly 11### test_mathieu_a 12### test_mathieu_even_coef 13### test_mathieu_odd_coef 14### test_modfresnelp 15### test_modfresnelm 16# test_pbdv_seq 17### test_pbvv_seq 18### test_sph_harm 19 20import itertools 21import platform 22import sys 23 24import numpy as np 25from numpy import (array, isnan, r_, arange, finfo, pi, sin, cos, tan, exp, 26 log, zeros, sqrt, asarray, inf, nan_to_num, real, arctan, float_) 27 28import pytest 29from pytest import raises as assert_raises 30from numpy.testing import (assert_equal, assert_almost_equal, 31 assert_array_equal, assert_array_almost_equal, assert_approx_equal, 32 assert_, assert_allclose, assert_array_almost_equal_nulp, 33 suppress_warnings) 34 35from scipy import special 36import scipy.special._ufuncs as cephes 37from scipy.special import ellipk 38 39from scipy.special._testutils import with_special_errors, \ 40 assert_func_equal, FuncData 41 42import math 43 44 45class TestCephes: 46 def test_airy(self): 47 cephes.airy(0) 48 49 def test_airye(self): 50 cephes.airye(0) 51 52 def test_binom(self): 53 n = np.array([0.264, 4, 5.2, 17]) 54 k = np.array([2, 0.4, 7, 3.3]) 55 nk = np.array(np.broadcast_arrays(n[:,None], k[None,:]) 56 ).reshape(2, -1).T 57 rknown = np.array([[-0.097152, 0.9263051596159367, 0.01858423645695389, 58 -0.007581020651518199],[6, 2.0214389119675666, 0, 2.9827344527963846], 59 [10.92, 2.22993515861399, -0.00585728, 10.468891352063146], 60 [136, 3.5252179590758828, 19448, 1024.5526916174495]]) 61 assert_func_equal(cephes.binom, rknown.ravel(), nk, rtol=1e-13) 62 63 # Test branches in implementation 64 np.random.seed(1234) 65 n = np.r_[np.arange(-7, 30), 1000*np.random.rand(30) - 500] 66 k = np.arange(0, 102) 67 nk = np.array(np.broadcast_arrays(n[:,None], k[None,:]) 68 ).reshape(2, -1).T 69 70 assert_func_equal(cephes.binom, 71 cephes.binom(nk[:,0], nk[:,1] * (1 + 1e-15)), 72 nk, 73 atol=1e-10, rtol=1e-10) 74 75 def test_binom_2(self): 76 # Test branches in implementation 77 np.random.seed(1234) 78 n = np.r_[np.logspace(1, 300, 20)] 79 k = np.arange(0, 102) 80 nk = np.array(np.broadcast_arrays(n[:,None], k[None,:]) 81 ).reshape(2, -1).T 82 83 assert_func_equal(cephes.binom, 84 cephes.binom(nk[:,0], nk[:,1] * (1 + 1e-15)), 85 nk, 86 atol=1e-10, rtol=1e-10) 87 88 def test_binom_exact(self): 89 @np.vectorize 90 def binom_int(n, k): 91 n = int(n) 92 k = int(k) 93 num = int(1) 94 den = int(1) 95 for i in range(1, k+1): 96 num *= i + n - k 97 den *= i 98 return float(num/den) 99 100 np.random.seed(1234) 101 n = np.arange(1, 15) 102 k = np.arange(0, 15) 103 nk = np.array(np.broadcast_arrays(n[:,None], k[None,:]) 104 ).reshape(2, -1).T 105 nk = nk[nk[:,0] >= nk[:,1]] 106 assert_func_equal(cephes.binom, 107 binom_int(nk[:,0], nk[:,1]), 108 nk, 109 atol=0, rtol=0) 110 111 def test_binom_nooverflow_8346(self): 112 # Test (binom(n, k) doesn't overflow prematurely */ 113 dataset = [ 114 (1000, 500, 2.70288240945436551e+299), 115 (1002, 501, 1.08007396880791225e+300), 116 (1004, 502, 4.31599279169058121e+300), 117 (1006, 503, 1.72468101616263781e+301), 118 (1008, 504, 6.89188009236419153e+301), 119 (1010, 505, 2.75402257948335448e+302), 120 (1012, 506, 1.10052048531923757e+303), 121 (1014, 507, 4.39774063758732849e+303), 122 (1016, 508, 1.75736486108312519e+304), 123 (1018, 509, 7.02255427788423734e+304), 124 (1020, 510, 2.80626776829962255e+305), 125 (1022, 511, 1.12140876377061240e+306), 126 (1024, 512, 4.48125455209897109e+306), 127 (1026, 513, 1.79075474304149900e+307), 128 (1028, 514, 7.15605105487789676e+307) 129 ] 130 dataset = np.asarray(dataset) 131 FuncData(cephes.binom, dataset, (0, 1), 2, rtol=1e-12).check() 132 133 def test_bdtr(self): 134 assert_equal(cephes.bdtr(1,1,0.5),1.0) 135 136 def test_bdtri(self): 137 assert_equal(cephes.bdtri(1,3,0.5),0.5) 138 139 def test_bdtrc(self): 140 assert_equal(cephes.bdtrc(1,3,0.5),0.5) 141 142 def test_bdtrin(self): 143 assert_equal(cephes.bdtrin(1,0,1),5.0) 144 145 def test_bdtrik(self): 146 cephes.bdtrik(1,3,0.5) 147 148 def test_bei(self): 149 assert_equal(cephes.bei(0),0.0) 150 151 def test_beip(self): 152 assert_equal(cephes.beip(0),0.0) 153 154 def test_ber(self): 155 assert_equal(cephes.ber(0),1.0) 156 157 def test_berp(self): 158 assert_equal(cephes.berp(0),0.0) 159 160 def test_besselpoly(self): 161 assert_equal(cephes.besselpoly(0,0,0),1.0) 162 163 def test_beta(self): 164 assert_equal(cephes.beta(1,1),1.0) 165 assert_allclose(cephes.beta(-100.3, 1e-200), cephes.gamma(1e-200)) 166 assert_allclose(cephes.beta(0.0342, 171), 24.070498359873497, 167 rtol=1e-13, atol=0) 168 169 def test_betainc(self): 170 assert_equal(cephes.betainc(1,1,1),1.0) 171 assert_allclose(cephes.betainc(0.0342, 171, 1e-10), 0.55269916901806648) 172 173 def test_betaln(self): 174 assert_equal(cephes.betaln(1,1),0.0) 175 assert_allclose(cephes.betaln(-100.3, 1e-200), cephes.gammaln(1e-200)) 176 assert_allclose(cephes.betaln(0.0342, 170), 3.1811881124242447, 177 rtol=1e-14, atol=0) 178 179 def test_betaincinv(self): 180 assert_equal(cephes.betaincinv(1,1,1),1.0) 181 assert_allclose(cephes.betaincinv(0.0342, 171, 0.25), 182 8.4231316935498957e-21, rtol=3e-12, atol=0) 183 184 def test_beta_inf(self): 185 assert_(np.isinf(special.beta(-1, 2))) 186 187 def test_btdtr(self): 188 assert_equal(cephes.btdtr(1,1,1),1.0) 189 190 def test_btdtri(self): 191 assert_equal(cephes.btdtri(1,1,1),1.0) 192 193 def test_btdtria(self): 194 assert_equal(cephes.btdtria(1,1,1),5.0) 195 196 def test_btdtrib(self): 197 assert_equal(cephes.btdtrib(1,1,1),5.0) 198 199 def test_cbrt(self): 200 assert_approx_equal(cephes.cbrt(1),1.0) 201 202 def test_chdtr(self): 203 assert_equal(cephes.chdtr(1,0),0.0) 204 205 def test_chdtrc(self): 206 assert_equal(cephes.chdtrc(1,0),1.0) 207 208 def test_chdtri(self): 209 assert_equal(cephes.chdtri(1,1),0.0) 210 211 def test_chdtriv(self): 212 assert_equal(cephes.chdtriv(0,0),5.0) 213 214 def test_chndtr(self): 215 assert_equal(cephes.chndtr(0,1,0),0.0) 216 217 # Each row holds (x, nu, lam, expected_value) 218 # These values were computed using Wolfram Alpha with 219 # CDF[NoncentralChiSquareDistribution[nu, lam], x] 220 values = np.array([ 221 [25.00, 20.0, 400, 4.1210655112396197139e-57], 222 [25.00, 8.00, 250, 2.3988026526832425878e-29], 223 [0.001, 8.00, 40., 5.3761806201366039084e-24], 224 [0.010, 8.00, 40., 5.45396231055999457039e-20], 225 [20.00, 2.00, 107, 1.39390743555819597802e-9], 226 [22.50, 2.00, 107, 7.11803307138105870671e-9], 227 [25.00, 2.00, 107, 3.11041244829864897313e-8], 228 [3.000, 2.00, 1.0, 0.62064365321954362734], 229 [350.0, 300., 10., 0.93880128006276407710], 230 [100.0, 13.5, 10., 0.99999999650104210949], 231 [700.0, 20.0, 400, 0.99999999925680650105], 232 [150.0, 13.5, 10., 0.99999999999999983046], 233 [160.0, 13.5, 10., 0.99999999999999999518], # 1.0 234 ]) 235 cdf = cephes.chndtr(values[:, 0], values[:, 1], values[:, 2]) 236 assert_allclose(cdf, values[:, 3], rtol=1e-12) 237 238 assert_almost_equal(cephes.chndtr(np.inf, np.inf, 0), 2.0) 239 assert_almost_equal(cephes.chndtr(2, 1, np.inf), 0.0) 240 assert_(np.isnan(cephes.chndtr(np.nan, 1, 2))) 241 assert_(np.isnan(cephes.chndtr(5, np.nan, 2))) 242 assert_(np.isnan(cephes.chndtr(5, 1, np.nan))) 243 244 def test_chndtridf(self): 245 assert_equal(cephes.chndtridf(0,0,1),5.0) 246 247 def test_chndtrinc(self): 248 assert_equal(cephes.chndtrinc(0,1,0),5.0) 249 250 def test_chndtrix(self): 251 assert_equal(cephes.chndtrix(0,1,0),0.0) 252 253 def test_cosdg(self): 254 assert_equal(cephes.cosdg(0),1.0) 255 256 def test_cosm1(self): 257 assert_equal(cephes.cosm1(0),0.0) 258 259 def test_cotdg(self): 260 assert_almost_equal(cephes.cotdg(45),1.0) 261 262 def test_dawsn(self): 263 assert_equal(cephes.dawsn(0),0.0) 264 assert_allclose(cephes.dawsn(1.23), 0.50053727749081767) 265 266 def test_diric(self): 267 # Test behavior near multiples of 2pi. Regression test for issue 268 # described in gh-4001. 269 n_odd = [1, 5, 25] 270 x = np.array(2*np.pi + 5e-5).astype(np.float32) 271 assert_almost_equal(special.diric(x, n_odd), 1.0, decimal=7) 272 x = np.array(2*np.pi + 1e-9).astype(np.float64) 273 assert_almost_equal(special.diric(x, n_odd), 1.0, decimal=15) 274 x = np.array(2*np.pi + 1e-15).astype(np.float64) 275 assert_almost_equal(special.diric(x, n_odd), 1.0, decimal=15) 276 if hasattr(np, 'float128'): 277 # No float128 available in 32-bit numpy 278 x = np.array(2*np.pi + 1e-12).astype(np.float128) 279 assert_almost_equal(special.diric(x, n_odd), 1.0, decimal=19) 280 281 n_even = [2, 4, 24] 282 x = np.array(2*np.pi + 1e-9).astype(np.float64) 283 assert_almost_equal(special.diric(x, n_even), -1.0, decimal=15) 284 285 # Test at some values not near a multiple of pi 286 x = np.arange(0.2*np.pi, 1.0*np.pi, 0.2*np.pi) 287 octave_result = [0.872677996249965, 0.539344662916632, 288 0.127322003750035, -0.206011329583298] 289 assert_almost_equal(special.diric(x, 3), octave_result, decimal=15) 290 291 def test_diric_broadcasting(self): 292 x = np.arange(5) 293 n = np.array([1, 3, 7]) 294 assert_(special.diric(x[:, np.newaxis], n).shape == (x.size, n.size)) 295 296 def test_ellipe(self): 297 assert_equal(cephes.ellipe(1),1.0) 298 299 def test_ellipeinc(self): 300 assert_equal(cephes.ellipeinc(0,1),0.0) 301 302 def test_ellipj(self): 303 cephes.ellipj(0,1) 304 305 def test_ellipk(self): 306 assert_allclose(ellipk(0), pi/2) 307 308 def test_ellipkinc(self): 309 assert_equal(cephes.ellipkinc(0,0),0.0) 310 311 def test_erf(self): 312 assert_equal(cephes.erf(0), 0.0) 313 314 def test_erf_symmetry(self): 315 x = 5.905732037710919 316 assert_equal(cephes.erf(x) + cephes.erf(-x), 0.0) 317 318 def test_erfc(self): 319 assert_equal(cephes.erfc(0), 1.0) 320 321 def test_exp10(self): 322 assert_approx_equal(cephes.exp10(2),100.0) 323 324 def test_exp2(self): 325 assert_equal(cephes.exp2(2),4.0) 326 327 def test_expm1(self): 328 assert_equal(cephes.expm1(0),0.0) 329 assert_equal(cephes.expm1(np.inf), np.inf) 330 assert_equal(cephes.expm1(-np.inf), -1) 331 assert_equal(cephes.expm1(np.nan), np.nan) 332 333 def test_expm1_complex(self): 334 expm1 = cephes.expm1 335 assert_equal(expm1(0 + 0j), 0 + 0j) 336 assert_equal(expm1(complex(np.inf, 0)), complex(np.inf, 0)) 337 assert_equal(expm1(complex(np.inf, 1)), complex(np.inf, np.inf)) 338 assert_equal(expm1(complex(np.inf, 2)), complex(-np.inf, np.inf)) 339 assert_equal(expm1(complex(np.inf, 4)), complex(-np.inf, -np.inf)) 340 assert_equal(expm1(complex(np.inf, 5)), complex(np.inf, -np.inf)) 341 assert_equal(expm1(complex(1, np.inf)), complex(np.nan, np.nan)) 342 assert_equal(expm1(complex(0, np.inf)), complex(np.nan, np.nan)) 343 assert_equal(expm1(complex(np.inf, np.inf)), complex(np.inf, np.nan)) 344 assert_equal(expm1(complex(-np.inf, np.inf)), complex(-1, 0)) 345 assert_equal(expm1(complex(-np.inf, np.nan)), complex(-1, 0)) 346 assert_equal(expm1(complex(np.inf, np.nan)), complex(np.inf, np.nan)) 347 assert_equal(expm1(complex(0, np.nan)), complex(np.nan, np.nan)) 348 assert_equal(expm1(complex(1, np.nan)), complex(np.nan, np.nan)) 349 assert_equal(expm1(complex(np.nan, 1)), complex(np.nan, np.nan)) 350 assert_equal(expm1(complex(np.nan, np.nan)), complex(np.nan, np.nan)) 351 352 @pytest.mark.xfail(reason='The real part of expm1(z) bad at these points') 353 def test_expm1_complex_hard(self): 354 # The real part of this function is difficult to evaluate when 355 # z.real = -log(cos(z.imag)). 356 y = np.array([0.1, 0.2, 0.3, 5, 11, 20]) 357 x = -np.log(np.cos(y)) 358 z = x + 1j*y 359 360 # evaluate using mpmath.expm1 with dps=1000 361 expected = np.array([-5.5507901846769623e-17+0.10033467208545054j, 362 2.4289354732893695e-18+0.20271003550867248j, 363 4.5235500262585768e-17+0.30933624960962319j, 364 7.8234305217489006e-17-3.3805150062465863j, 365 -1.3685191953697676e-16-225.95084645419513j, 366 8.7175620481291045e-17+2.2371609442247422j]) 367 found = cephes.expm1(z) 368 # this passes. 369 assert_array_almost_equal_nulp(found.imag, expected.imag, 3) 370 # this fails. 371 assert_array_almost_equal_nulp(found.real, expected.real, 20) 372 373 def test_fdtr(self): 374 assert_equal(cephes.fdtr(1, 1, 0), 0.0) 375 # Computed using Wolfram Alpha: CDF[FRatioDistribution[1e-6, 5], 10] 376 assert_allclose(cephes.fdtr(1e-6, 5, 10), 0.9999940790193488, 377 rtol=1e-12) 378 379 def test_fdtrc(self): 380 assert_equal(cephes.fdtrc(1, 1, 0), 1.0) 381 # Computed using Wolfram Alpha: 382 # 1 - CDF[FRatioDistribution[2, 1/10], 1e10] 383 assert_allclose(cephes.fdtrc(2, 0.1, 1e10), 0.27223784621293512, 384 rtol=1e-12) 385 386 def test_fdtri(self): 387 assert_allclose(cephes.fdtri(1, 1, [0.499, 0.501]), 388 array([0.9937365, 1.00630298]), rtol=1e-6) 389 # From Wolfram Alpha: 390 # CDF[FRatioDistribution[1/10, 1], 3] = 0.8756751669632105666874... 391 p = 0.8756751669632105666874 392 assert_allclose(cephes.fdtri(0.1, 1, p), 3, rtol=1e-12) 393 394 @pytest.mark.xfail(reason='Returns nan on i686.') 395 def test_fdtri_mysterious_failure(self): 396 assert_allclose(cephes.fdtri(1, 1, 0.5), 1) 397 398 def test_fdtridfd(self): 399 assert_equal(cephes.fdtridfd(1,0,0),5.0) 400 401 def test_fresnel(self): 402 assert_equal(cephes.fresnel(0),(0.0,0.0)) 403 404 def test_gamma(self): 405 assert_equal(cephes.gamma(5),24.0) 406 407 def test_gammainccinv(self): 408 assert_equal(cephes.gammainccinv(5,1),0.0) 409 410 def test_gammaln(self): 411 cephes.gammaln(10) 412 413 def test_gammasgn(self): 414 vals = np.array([-4, -3.5, -2.3, 1, 4.2], np.float64) 415 assert_array_equal(cephes.gammasgn(vals), np.sign(cephes.rgamma(vals))) 416 417 def test_gdtr(self): 418 assert_equal(cephes.gdtr(1,1,0),0.0) 419 420 def test_gdtr_inf(self): 421 assert_equal(cephes.gdtr(1,1,np.inf),1.0) 422 423 def test_gdtrc(self): 424 assert_equal(cephes.gdtrc(1,1,0),1.0) 425 426 def test_gdtria(self): 427 assert_equal(cephes.gdtria(0,1,1),0.0) 428 429 def test_gdtrib(self): 430 cephes.gdtrib(1,0,1) 431 # assert_equal(cephes.gdtrib(1,0,1),5.0) 432 433 def test_gdtrix(self): 434 cephes.gdtrix(1,1,.1) 435 436 def test_hankel1(self): 437 cephes.hankel1(1,1) 438 439 def test_hankel1e(self): 440 cephes.hankel1e(1,1) 441 442 def test_hankel2(self): 443 cephes.hankel2(1,1) 444 445 def test_hankel2e(self): 446 cephes.hankel2e(1,1) 447 448 def test_hyp1f1(self): 449 assert_approx_equal(cephes.hyp1f1(1,1,1), exp(1.0)) 450 assert_approx_equal(cephes.hyp1f1(3,4,-6), 0.026056422099537251095) 451 cephes.hyp1f1(1,1,1) 452 453 def test_hyp2f1(self): 454 assert_equal(cephes.hyp2f1(1,1,1,0),1.0) 455 456 def test_i0(self): 457 assert_equal(cephes.i0(0),1.0) 458 459 def test_i0e(self): 460 assert_equal(cephes.i0e(0),1.0) 461 462 def test_i1(self): 463 assert_equal(cephes.i1(0),0.0) 464 465 def test_i1e(self): 466 assert_equal(cephes.i1e(0),0.0) 467 468 def test_it2i0k0(self): 469 cephes.it2i0k0(1) 470 471 def test_it2j0y0(self): 472 cephes.it2j0y0(1) 473 474 def test_it2struve0(self): 475 cephes.it2struve0(1) 476 477 def test_itairy(self): 478 cephes.itairy(1) 479 480 def test_iti0k0(self): 481 assert_equal(cephes.iti0k0(0),(0.0,0.0)) 482 483 def test_itj0y0(self): 484 assert_equal(cephes.itj0y0(0),(0.0,0.0)) 485 486 def test_itmodstruve0(self): 487 assert_equal(cephes.itmodstruve0(0),0.0) 488 489 def test_itstruve0(self): 490 assert_equal(cephes.itstruve0(0),0.0) 491 492 def test_iv(self): 493 assert_equal(cephes.iv(1,0),0.0) 494 495 def _check_ive(self): 496 assert_equal(cephes.ive(1,0),0.0) 497 498 def test_j0(self): 499 assert_equal(cephes.j0(0),1.0) 500 501 def test_j1(self): 502 assert_equal(cephes.j1(0),0.0) 503 504 def test_jn(self): 505 assert_equal(cephes.jn(0,0),1.0) 506 507 def test_jv(self): 508 assert_equal(cephes.jv(0,0),1.0) 509 510 def _check_jve(self): 511 assert_equal(cephes.jve(0,0),1.0) 512 513 def test_k0(self): 514 cephes.k0(2) 515 516 def test_k0e(self): 517 cephes.k0e(2) 518 519 def test_k1(self): 520 cephes.k1(2) 521 522 def test_k1e(self): 523 cephes.k1e(2) 524 525 def test_kei(self): 526 cephes.kei(2) 527 528 def test_keip(self): 529 assert_equal(cephes.keip(0),0.0) 530 531 def test_ker(self): 532 cephes.ker(2) 533 534 def test_kerp(self): 535 cephes.kerp(2) 536 537 def _check_kelvin(self): 538 cephes.kelvin(2) 539 540 def test_kn(self): 541 cephes.kn(1,1) 542 543 def test_kolmogi(self): 544 assert_equal(cephes.kolmogi(1),0.0) 545 assert_(np.isnan(cephes.kolmogi(np.nan))) 546 547 def test_kolmogorov(self): 548 assert_equal(cephes.kolmogorov(0), 1.0) 549 550 def test_kolmogp(self): 551 assert_equal(cephes._kolmogp(0), -0.0) 552 553 def test_kolmogc(self): 554 assert_equal(cephes._kolmogc(0), 0.0) 555 556 def test_kolmogci(self): 557 assert_equal(cephes._kolmogci(0), 0.0) 558 assert_(np.isnan(cephes._kolmogci(np.nan))) 559 560 def _check_kv(self): 561 cephes.kv(1,1) 562 563 def _check_kve(self): 564 cephes.kve(1,1) 565 566 def test_log1p(self): 567 log1p = cephes.log1p 568 assert_equal(log1p(0), 0.0) 569 assert_equal(log1p(-1), -np.inf) 570 assert_equal(log1p(-2), np.nan) 571 assert_equal(log1p(np.inf), np.inf) 572 573 def test_log1p_complex(self): 574 log1p = cephes.log1p 575 c = complex 576 assert_equal(log1p(0 + 0j), 0 + 0j) 577 assert_equal(log1p(c(-1, 0)), c(-np.inf, 0)) 578 with suppress_warnings() as sup: 579 sup.filter(RuntimeWarning, "invalid value encountered in multiply") 580 assert_allclose(log1p(c(1, np.inf)), c(np.inf, np.pi/2)) 581 assert_equal(log1p(c(1, np.nan)), c(np.nan, np.nan)) 582 assert_allclose(log1p(c(-np.inf, 1)), c(np.inf, np.pi)) 583 assert_equal(log1p(c(np.inf, 1)), c(np.inf, 0)) 584 assert_allclose(log1p(c(-np.inf, np.inf)), c(np.inf, 3*np.pi/4)) 585 assert_allclose(log1p(c(np.inf, np.inf)), c(np.inf, np.pi/4)) 586 assert_equal(log1p(c(np.inf, np.nan)), c(np.inf, np.nan)) 587 assert_equal(log1p(c(-np.inf, np.nan)), c(np.inf, np.nan)) 588 assert_equal(log1p(c(np.nan, np.inf)), c(np.inf, np.nan)) 589 assert_equal(log1p(c(np.nan, 1)), c(np.nan, np.nan)) 590 assert_equal(log1p(c(np.nan, np.nan)), c(np.nan, np.nan)) 591 592 def test_lpmv(self): 593 assert_equal(cephes.lpmv(0,0,1),1.0) 594 595 def test_mathieu_a(self): 596 assert_equal(cephes.mathieu_a(1,0),1.0) 597 598 def test_mathieu_b(self): 599 assert_equal(cephes.mathieu_b(1,0),1.0) 600 601 def test_mathieu_cem(self): 602 assert_equal(cephes.mathieu_cem(1,0,0),(1.0,0.0)) 603 604 # Test AMS 20.2.27 605 @np.vectorize 606 def ce_smallq(m, q, z): 607 z *= np.pi/180 608 if m == 0: 609 return 2**(-0.5) * (1 - .5*q*cos(2*z)) # + O(q^2) 610 elif m == 1: 611 return cos(z) - q/8 * cos(3*z) # + O(q^2) 612 elif m == 2: 613 return cos(2*z) - q*(cos(4*z)/12 - 1/4) # + O(q^2) 614 else: 615 return cos(m*z) - q*(cos((m+2)*z)/(4*(m+1)) - cos((m-2)*z)/(4*(m-1))) # + O(q^2) 616 m = np.arange(0, 100) 617 q = np.r_[0, np.logspace(-30, -9, 10)] 618 assert_allclose(cephes.mathieu_cem(m[:,None], q[None,:], 0.123)[0], 619 ce_smallq(m[:,None], q[None,:], 0.123), 620 rtol=1e-14, atol=0) 621 622 def test_mathieu_sem(self): 623 assert_equal(cephes.mathieu_sem(1,0,0),(0.0,1.0)) 624 625 # Test AMS 20.2.27 626 @np.vectorize 627 def se_smallq(m, q, z): 628 z *= np.pi/180 629 if m == 1: 630 return sin(z) - q/8 * sin(3*z) # + O(q^2) 631 elif m == 2: 632 return sin(2*z) - q*sin(4*z)/12 # + O(q^2) 633 else: 634 return sin(m*z) - q*(sin((m+2)*z)/(4*(m+1)) - sin((m-2)*z)/(4*(m-1))) # + O(q^2) 635 m = np.arange(1, 100) 636 q = np.r_[0, np.logspace(-30, -9, 10)] 637 assert_allclose(cephes.mathieu_sem(m[:,None], q[None,:], 0.123)[0], 638 se_smallq(m[:,None], q[None,:], 0.123), 639 rtol=1e-14, atol=0) 640 641 def test_mathieu_modcem1(self): 642 assert_equal(cephes.mathieu_modcem1(1,0,0),(0.0,0.0)) 643 644 def test_mathieu_modcem2(self): 645 cephes.mathieu_modcem2(1,1,1) 646 647 # Test reflection relation AMS 20.6.19 648 m = np.arange(0, 4)[:,None,None] 649 q = np.r_[np.logspace(-2, 2, 10)][None,:,None] 650 z = np.linspace(0, 1, 7)[None,None,:] 651 652 y1 = cephes.mathieu_modcem2(m, q, -z)[0] 653 654 fr = -cephes.mathieu_modcem2(m, q, 0)[0] / cephes.mathieu_modcem1(m, q, 0)[0] 655 y2 = -cephes.mathieu_modcem2(m, q, z)[0] - 2*fr*cephes.mathieu_modcem1(m, q, z)[0] 656 657 assert_allclose(y1, y2, rtol=1e-10) 658 659 def test_mathieu_modsem1(self): 660 assert_equal(cephes.mathieu_modsem1(1,0,0),(0.0,0.0)) 661 662 def test_mathieu_modsem2(self): 663 cephes.mathieu_modsem2(1,1,1) 664 665 # Test reflection relation AMS 20.6.20 666 m = np.arange(1, 4)[:,None,None] 667 q = np.r_[np.logspace(-2, 2, 10)][None,:,None] 668 z = np.linspace(0, 1, 7)[None,None,:] 669 670 y1 = cephes.mathieu_modsem2(m, q, -z)[0] 671 fr = cephes.mathieu_modsem2(m, q, 0)[1] / cephes.mathieu_modsem1(m, q, 0)[1] 672 y2 = cephes.mathieu_modsem2(m, q, z)[0] - 2*fr*cephes.mathieu_modsem1(m, q, z)[0] 673 assert_allclose(y1, y2, rtol=1e-10) 674 675 def test_mathieu_overflow(self): 676 # Check that these return NaNs instead of causing a SEGV 677 assert_equal(cephes.mathieu_cem(10000, 0, 1.3), (np.nan, np.nan)) 678 assert_equal(cephes.mathieu_sem(10000, 0, 1.3), (np.nan, np.nan)) 679 assert_equal(cephes.mathieu_cem(10000, 1.5, 1.3), (np.nan, np.nan)) 680 assert_equal(cephes.mathieu_sem(10000, 1.5, 1.3), (np.nan, np.nan)) 681 assert_equal(cephes.mathieu_modcem1(10000, 1.5, 1.3), (np.nan, np.nan)) 682 assert_equal(cephes.mathieu_modsem1(10000, 1.5, 1.3), (np.nan, np.nan)) 683 assert_equal(cephes.mathieu_modcem2(10000, 1.5, 1.3), (np.nan, np.nan)) 684 assert_equal(cephes.mathieu_modsem2(10000, 1.5, 1.3), (np.nan, np.nan)) 685 686 def test_mathieu_ticket_1847(self): 687 # Regression test --- this call had some out-of-bounds access 688 # and could return nan occasionally 689 for k in range(60): 690 v = cephes.mathieu_modsem2(2, 100, -1) 691 # Values from ACM TOMS 804 (derivate by numerical differentiation) 692 assert_allclose(v[0], 0.1431742913063671074347, rtol=1e-10) 693 assert_allclose(v[1], 0.9017807375832909144719, rtol=1e-4) 694 695 def test_modfresnelm(self): 696 cephes.modfresnelm(0) 697 698 def test_modfresnelp(self): 699 cephes.modfresnelp(0) 700 701 def _check_modstruve(self): 702 assert_equal(cephes.modstruve(1,0),0.0) 703 704 def test_nbdtr(self): 705 assert_equal(cephes.nbdtr(1,1,1),1.0) 706 707 def test_nbdtrc(self): 708 assert_equal(cephes.nbdtrc(1,1,1),0.0) 709 710 def test_nbdtri(self): 711 assert_equal(cephes.nbdtri(1,1,1),1.0) 712 713 def __check_nbdtrik(self): 714 cephes.nbdtrik(1,.4,.5) 715 716 def test_nbdtrin(self): 717 assert_equal(cephes.nbdtrin(1,0,0),5.0) 718 719 def test_ncfdtr(self): 720 assert_equal(cephes.ncfdtr(1,1,1,0),0.0) 721 722 def test_ncfdtri(self): 723 assert_equal(cephes.ncfdtri(1, 1, 1, 0), 0.0) 724 f = [0.5, 1, 1.5] 725 p = cephes.ncfdtr(2, 3, 1.5, f) 726 assert_allclose(cephes.ncfdtri(2, 3, 1.5, p), f) 727 728 def test_ncfdtridfd(self): 729 dfd = [1, 2, 3] 730 p = cephes.ncfdtr(2, dfd, 0.25, 15) 731 assert_allclose(cephes.ncfdtridfd(2, p, 0.25, 15), dfd) 732 733 def test_ncfdtridfn(self): 734 dfn = [0.1, 1, 2, 3, 1e4] 735 p = cephes.ncfdtr(dfn, 2, 0.25, 15) 736 assert_allclose(cephes.ncfdtridfn(p, 2, 0.25, 15), dfn, rtol=1e-5) 737 738 def test_ncfdtrinc(self): 739 nc = [0.5, 1.5, 2.0] 740 p = cephes.ncfdtr(2, 3, nc, 15) 741 assert_allclose(cephes.ncfdtrinc(2, 3, p, 15), nc) 742 743 def test_nctdtr(self): 744 assert_equal(cephes.nctdtr(1,0,0),0.5) 745 assert_equal(cephes.nctdtr(9, 65536, 45), 0.0) 746 747 assert_approx_equal(cephes.nctdtr(np.inf, 1., 1.), 0.5, 5) 748 assert_(np.isnan(cephes.nctdtr(2., np.inf, 10.))) 749 assert_approx_equal(cephes.nctdtr(2., 1., np.inf), 1.) 750 751 assert_(np.isnan(cephes.nctdtr(np.nan, 1., 1.))) 752 assert_(np.isnan(cephes.nctdtr(2., np.nan, 1.))) 753 assert_(np.isnan(cephes.nctdtr(2., 1., np.nan))) 754 755 def __check_nctdtridf(self): 756 cephes.nctdtridf(1,0.5,0) 757 758 def test_nctdtrinc(self): 759 cephes.nctdtrinc(1,0,0) 760 761 def test_nctdtrit(self): 762 cephes.nctdtrit(.1,0.2,.5) 763 764 def test_nrdtrimn(self): 765 assert_approx_equal(cephes.nrdtrimn(0.5,1,1),1.0) 766 767 def test_nrdtrisd(self): 768 assert_allclose(cephes.nrdtrisd(0.5,0.5,0.5), 0.0, 769 atol=0, rtol=0) 770 771 def test_obl_ang1(self): 772 cephes.obl_ang1(1,1,1,0) 773 774 def test_obl_ang1_cv(self): 775 result = cephes.obl_ang1_cv(1,1,1,1,0) 776 assert_almost_equal(result[0],1.0) 777 assert_almost_equal(result[1],0.0) 778 779 def _check_obl_cv(self): 780 assert_equal(cephes.obl_cv(1,1,0),2.0) 781 782 def test_obl_rad1(self): 783 cephes.obl_rad1(1,1,1,0) 784 785 def test_obl_rad1_cv(self): 786 cephes.obl_rad1_cv(1,1,1,1,0) 787 788 def test_obl_rad2(self): 789 cephes.obl_rad2(1,1,1,0) 790 791 def test_obl_rad2_cv(self): 792 cephes.obl_rad2_cv(1,1,1,1,0) 793 794 def test_pbdv(self): 795 assert_equal(cephes.pbdv(1,0),(0.0,1.0)) 796 797 def test_pbvv(self): 798 cephes.pbvv(1,0) 799 800 def test_pbwa(self): 801 cephes.pbwa(1,0) 802 803 def test_pdtr(self): 804 val = cephes.pdtr(0, 1) 805 assert_almost_equal(val, np.exp(-1)) 806 # Edge case: m = 0. 807 val = cephes.pdtr([0, 1, 2], 0) 808 assert_array_equal(val, [1, 1, 1]) 809 810 def test_pdtrc(self): 811 val = cephes.pdtrc(0, 1) 812 assert_almost_equal(val, 1 - np.exp(-1)) 813 # Edge case: m = 0. 814 val = cephes.pdtrc([0, 1, 2], 0.0) 815 assert_array_equal(val, [0, 0, 0]) 816 817 def test_pdtri(self): 818 with suppress_warnings() as sup: 819 sup.filter(RuntimeWarning, "floating point number truncated to an integer") 820 cephes.pdtri(0.5,0.5) 821 822 def test_pdtrik(self): 823 k = cephes.pdtrik(0.5, 1) 824 assert_almost_equal(cephes.gammaincc(k + 1, 1), 0.5) 825 # Edge case: m = 0 or very small. 826 k = cephes.pdtrik([[0], [0.25], [0.95]], [0, 1e-20, 1e-6]) 827 assert_array_equal(k, np.zeros((3, 3))) 828 829 def test_pro_ang1(self): 830 cephes.pro_ang1(1,1,1,0) 831 832 def test_pro_ang1_cv(self): 833 assert_array_almost_equal(cephes.pro_ang1_cv(1,1,1,1,0), 834 array((1.0,0.0))) 835 836 def _check_pro_cv(self): 837 assert_equal(cephes.pro_cv(1,1,0),2.0) 838 839 def test_pro_rad1(self): 840 cephes.pro_rad1(1,1,1,0.1) 841 842 def test_pro_rad1_cv(self): 843 cephes.pro_rad1_cv(1,1,1,1,0) 844 845 def test_pro_rad2(self): 846 cephes.pro_rad2(1,1,1,0) 847 848 def test_pro_rad2_cv(self): 849 cephes.pro_rad2_cv(1,1,1,1,0) 850 851 def test_psi(self): 852 cephes.psi(1) 853 854 def test_radian(self): 855 assert_equal(cephes.radian(0,0,0),0) 856 857 def test_rgamma(self): 858 assert_equal(cephes.rgamma(1),1.0) 859 860 def test_round(self): 861 assert_equal(cephes.round(3.4),3.0) 862 assert_equal(cephes.round(-3.4),-3.0) 863 assert_equal(cephes.round(3.6),4.0) 864 assert_equal(cephes.round(-3.6),-4.0) 865 assert_equal(cephes.round(3.5),4.0) 866 assert_equal(cephes.round(-3.5),-4.0) 867 868 def test_shichi(self): 869 cephes.shichi(1) 870 871 def test_sici(self): 872 cephes.sici(1) 873 874 s, c = cephes.sici(np.inf) 875 assert_almost_equal(s, np.pi * 0.5) 876 assert_almost_equal(c, 0) 877 878 s, c = cephes.sici(-np.inf) 879 assert_almost_equal(s, -np.pi * 0.5) 880 assert_(np.isnan(c), "cosine integral(-inf) is not nan") 881 882 def test_sindg(self): 883 assert_equal(cephes.sindg(90),1.0) 884 885 def test_smirnov(self): 886 assert_equal(cephes.smirnov(1,.1),0.9) 887 assert_(np.isnan(cephes.smirnov(1,np.nan))) 888 889 def test_smirnovp(self): 890 assert_equal(cephes._smirnovp(1, .1), -1) 891 assert_equal(cephes._smirnovp(2, 0.75), -2*(0.25)**(2-1)) 892 assert_equal(cephes._smirnovp(3, 0.75), -3*(0.25)**(3-1)) 893 assert_(np.isnan(cephes._smirnovp(1, np.nan))) 894 895 def test_smirnovc(self): 896 assert_equal(cephes._smirnovc(1,.1),0.1) 897 assert_(np.isnan(cephes._smirnovc(1,np.nan))) 898 x10 = np.linspace(0, 1, 11, endpoint=True) 899 assert_almost_equal(cephes._smirnovc(3, x10), 1-cephes.smirnov(3, x10)) 900 x4 = np.linspace(0, 1, 5, endpoint=True) 901 assert_almost_equal(cephes._smirnovc(4, x4), 1-cephes.smirnov(4, x4)) 902 903 def test_smirnovi(self): 904 assert_almost_equal(cephes.smirnov(1,cephes.smirnovi(1,0.4)),0.4) 905 assert_almost_equal(cephes.smirnov(1,cephes.smirnovi(1,0.6)),0.6) 906 assert_(np.isnan(cephes.smirnovi(1,np.nan))) 907 908 def test_smirnovci(self): 909 assert_almost_equal(cephes._smirnovc(1,cephes._smirnovci(1,0.4)),0.4) 910 assert_almost_equal(cephes._smirnovc(1,cephes._smirnovci(1,0.6)),0.6) 911 assert_(np.isnan(cephes._smirnovci(1,np.nan))) 912 913 def test_spence(self): 914 assert_equal(cephes.spence(1),0.0) 915 916 def test_stdtr(self): 917 assert_equal(cephes.stdtr(1,0),0.5) 918 assert_almost_equal(cephes.stdtr(1,1), 0.75) 919 assert_almost_equal(cephes.stdtr(1,2), 0.852416382349) 920 921 def test_stdtridf(self): 922 cephes.stdtridf(0.7,1) 923 924 def test_stdtrit(self): 925 cephes.stdtrit(1,0.7) 926 927 def test_struve(self): 928 assert_equal(cephes.struve(0,0),0.0) 929 930 def test_tandg(self): 931 assert_equal(cephes.tandg(45),1.0) 932 933 def test_tklmbda(self): 934 assert_almost_equal(cephes.tklmbda(1,1),1.0) 935 936 def test_y0(self): 937 cephes.y0(1) 938 939 def test_y1(self): 940 cephes.y1(1) 941 942 def test_yn(self): 943 cephes.yn(1,1) 944 945 def test_yv(self): 946 cephes.yv(1,1) 947 948 def _check_yve(self): 949 cephes.yve(1,1) 950 951 def test_wofz(self): 952 z = [complex(624.2,-0.26123), complex(-0.4,3.), complex(0.6,2.), 953 complex(-1.,1.), complex(-1.,-9.), complex(-1.,9.), 954 complex(-0.0000000234545,1.1234), complex(-3.,5.1), 955 complex(-53,30.1), complex(0.0,0.12345), 956 complex(11,1), complex(-22,-2), complex(9,-28), 957 complex(21,-33), complex(1e5,1e5), complex(1e14,1e14) 958 ] 959 w = [ 960 complex(-3.78270245518980507452677445620103199303131110e-7, 961 0.000903861276433172057331093754199933411710053155), 962 complex(0.1764906227004816847297495349730234591778719532788, 963 -0.02146550539468457616788719893991501311573031095617), 964 complex(0.2410250715772692146133539023007113781272362309451, 965 0.06087579663428089745895459735240964093522265589350), 966 complex(0.30474420525691259245713884106959496013413834051768, 967 -0.20821893820283162728743734725471561394145872072738), 968 complex(7.317131068972378096865595229600561710140617977e34, 969 8.321873499714402777186848353320412813066170427e34), 970 complex(0.0615698507236323685519612934241429530190806818395, 971 -0.00676005783716575013073036218018565206070072304635), 972 complex(0.3960793007699874918961319170187598400134746631, 973 -5.593152259116644920546186222529802777409274656e-9), 974 complex(0.08217199226739447943295069917990417630675021771804, 975 -0.04701291087643609891018366143118110965272615832184), 976 complex(0.00457246000350281640952328010227885008541748668738, 977 -0.00804900791411691821818731763401840373998654987934), 978 complex(0.8746342859608052666092782112565360755791467973338452, 979 0.), 980 complex(0.00468190164965444174367477874864366058339647648741, 981 0.0510735563901306197993676329845149741675029197050), 982 complex(-0.0023193175200187620902125853834909543869428763219, 983 -0.025460054739731556004902057663500272721780776336), 984 complex(9.11463368405637174660562096516414499772662584e304, 985 3.97101807145263333769664875189354358563218932e305), 986 complex(-4.4927207857715598976165541011143706155432296e281, 987 -2.8019591213423077494444700357168707775769028e281), 988 complex(2.820947917809305132678577516325951485807107151e-6, 989 2.820947917668257736791638444590253942253354058e-6), 990 complex(2.82094791773878143474039725787438662716372268e-15, 991 2.82094791773878143474039725773333923127678361e-15) 992 ] 993 assert_func_equal(cephes.wofz, w, z, rtol=1e-13) 994 995 996class TestAiry: 997 def test_airy(self): 998 # This tests the airy function to ensure 8 place accuracy in computation 999 1000 x = special.airy(.99) 1001 assert_array_almost_equal(x,array([0.13689066,-0.16050153,1.19815925,0.92046818]),8) 1002 x = special.airy(.41) 1003 assert_array_almost_equal(x,array([0.25238916,-.23480512,0.80686202,0.51053919]),8) 1004 x = special.airy(-.36) 1005 assert_array_almost_equal(x,array([0.44508477,-0.23186773,0.44939534,0.48105354]),8) 1006 1007 def test_airye(self): 1008 a = special.airye(0.01) 1009 b = special.airy(0.01) 1010 b1 = [None]*4 1011 for n in range(2): 1012 b1[n] = b[n]*exp(2.0/3.0*0.01*sqrt(0.01)) 1013 for n in range(2,4): 1014 b1[n] = b[n]*exp(-abs(real(2.0/3.0*0.01*sqrt(0.01)))) 1015 assert_array_almost_equal(a,b1,6) 1016 1017 def test_bi_zeros(self): 1018 bi = special.bi_zeros(2) 1019 bia = (array([-1.17371322, -3.2710930]), 1020 array([-2.29443968, -4.07315509]), 1021 array([-0.45494438, 0.39652284]), 1022 array([0.60195789, -0.76031014])) 1023 assert_array_almost_equal(bi,bia,4) 1024 1025 bi = special.bi_zeros(5) 1026 assert_array_almost_equal(bi[0],array([-1.173713222709127, 1027 -3.271093302836352, 1028 -4.830737841662016, 1029 -6.169852128310251, 1030 -7.376762079367764]),11) 1031 1032 assert_array_almost_equal(bi[1],array([-2.294439682614122, 1033 -4.073155089071828, 1034 -5.512395729663599, 1035 -6.781294445990305, 1036 -7.940178689168587]),10) 1037 1038 assert_array_almost_equal(bi[2],array([-0.454944383639657, 1039 0.396522836094465, 1040 -0.367969161486959, 1041 0.349499116831805, 1042 -0.336026240133662]),11) 1043 1044 assert_array_almost_equal(bi[3],array([0.601957887976239, 1045 -0.760310141492801, 1046 0.836991012619261, 1047 -0.88947990142654, 1048 0.929983638568022]),10) 1049 1050 def test_ai_zeros(self): 1051 ai = special.ai_zeros(1) 1052 assert_array_almost_equal(ai,(array([-2.33810741]), 1053 array([-1.01879297]), 1054 array([0.5357]), 1055 array([0.7012])),4) 1056 1057 def test_ai_zeros_big(self): 1058 z, zp, ai_zpx, aip_zx = special.ai_zeros(50000) 1059 ai_z, aip_z, _, _ = special.airy(z) 1060 ai_zp, aip_zp, _, _ = special.airy(zp) 1061 1062 ai_envelope = 1/abs(z)**(1./4) 1063 aip_envelope = abs(zp)**(1./4) 1064 1065 # Check values 1066 assert_allclose(ai_zpx, ai_zp, rtol=1e-10) 1067 assert_allclose(aip_zx, aip_z, rtol=1e-10) 1068 1069 # Check they are zeros 1070 assert_allclose(ai_z/ai_envelope, 0, atol=1e-10, rtol=0) 1071 assert_allclose(aip_zp/aip_envelope, 0, atol=1e-10, rtol=0) 1072 1073 # Check first zeros, DLMF 9.9.1 1074 assert_allclose(z[:6], 1075 [-2.3381074105, -4.0879494441, -5.5205598281, 1076 -6.7867080901, -7.9441335871, -9.0226508533], rtol=1e-10) 1077 assert_allclose(zp[:6], 1078 [-1.0187929716, -3.2481975822, -4.8200992112, 1079 -6.1633073556, -7.3721772550, -8.4884867340], rtol=1e-10) 1080 1081 def test_bi_zeros_big(self): 1082 z, zp, bi_zpx, bip_zx = special.bi_zeros(50000) 1083 _, _, bi_z, bip_z = special.airy(z) 1084 _, _, bi_zp, bip_zp = special.airy(zp) 1085 1086 bi_envelope = 1/abs(z)**(1./4) 1087 bip_envelope = abs(zp)**(1./4) 1088 1089 # Check values 1090 assert_allclose(bi_zpx, bi_zp, rtol=1e-10) 1091 assert_allclose(bip_zx, bip_z, rtol=1e-10) 1092 1093 # Check they are zeros 1094 assert_allclose(bi_z/bi_envelope, 0, atol=1e-10, rtol=0) 1095 assert_allclose(bip_zp/bip_envelope, 0, atol=1e-10, rtol=0) 1096 1097 # Check first zeros, DLMF 9.9.2 1098 assert_allclose(z[:6], 1099 [-1.1737132227, -3.2710933028, -4.8307378417, 1100 -6.1698521283, -7.3767620794, -8.4919488465], rtol=1e-10) 1101 assert_allclose(zp[:6], 1102 [-2.2944396826, -4.0731550891, -5.5123957297, 1103 -6.7812944460, -7.9401786892, -9.0195833588], rtol=1e-10) 1104 1105 1106class TestAssocLaguerre: 1107 def test_assoc_laguerre(self): 1108 a1 = special.genlaguerre(11,1) 1109 a2 = special.assoc_laguerre(.2,11,1) 1110 assert_array_almost_equal(a2,a1(.2),8) 1111 a2 = special.assoc_laguerre(1,11,1) 1112 assert_array_almost_equal(a2,a1(1),8) 1113 1114 1115class TestBesselpoly: 1116 def test_besselpoly(self): 1117 pass 1118 1119 1120class TestKelvin: 1121 def test_bei(self): 1122 mbei = special.bei(2) 1123 assert_almost_equal(mbei, 0.9722916273066613,5) # this may not be exact 1124 1125 def test_beip(self): 1126 mbeip = special.beip(2) 1127 assert_almost_equal(mbeip,0.91701361338403631,5) # this may not be exact 1128 1129 def test_ber(self): 1130 mber = special.ber(2) 1131 assert_almost_equal(mber,0.75173418271380821,5) # this may not be exact 1132 1133 def test_berp(self): 1134 mberp = special.berp(2) 1135 assert_almost_equal(mberp,-0.49306712470943909,5) # this may not be exact 1136 1137 def test_bei_zeros(self): 1138 # Abramowitz & Stegun, Table 9.12 1139 bi = special.bei_zeros(5) 1140 assert_array_almost_equal(bi,array([5.02622, 1141 9.45541, 1142 13.89349, 1143 18.33398, 1144 22.77544]),4) 1145 1146 def test_beip_zeros(self): 1147 bip = special.beip_zeros(5) 1148 assert_array_almost_equal(bip,array([3.772673304934953, 1149 8.280987849760042, 1150 12.742147523633703, 1151 17.193431752512542, 1152 21.641143941167325]),8) 1153 1154 def test_ber_zeros(self): 1155 ber = special.ber_zeros(5) 1156 assert_array_almost_equal(ber,array([2.84892, 1157 7.23883, 1158 11.67396, 1159 16.11356, 1160 20.55463]),4) 1161 1162 def test_berp_zeros(self): 1163 brp = special.berp_zeros(5) 1164 assert_array_almost_equal(brp,array([6.03871, 1165 10.51364, 1166 14.96844, 1167 19.41758, 1168 23.86430]),4) 1169 1170 def test_kelvin(self): 1171 mkelv = special.kelvin(2) 1172 assert_array_almost_equal(mkelv,(special.ber(2) + special.bei(2)*1j, 1173 special.ker(2) + special.kei(2)*1j, 1174 special.berp(2) + special.beip(2)*1j, 1175 special.kerp(2) + special.keip(2)*1j),8) 1176 1177 def test_kei(self): 1178 mkei = special.kei(2) 1179 assert_almost_equal(mkei,-0.20240006776470432,5) 1180 1181 def test_keip(self): 1182 mkeip = special.keip(2) 1183 assert_almost_equal(mkeip,0.21980790991960536,5) 1184 1185 def test_ker(self): 1186 mker = special.ker(2) 1187 assert_almost_equal(mker,-0.041664513991509472,5) 1188 1189 def test_kerp(self): 1190 mkerp = special.kerp(2) 1191 assert_almost_equal(mkerp,-0.10660096588105264,5) 1192 1193 def test_kei_zeros(self): 1194 kei = special.kei_zeros(5) 1195 assert_array_almost_equal(kei,array([3.91467, 1196 8.34422, 1197 12.78256, 1198 17.22314, 1199 21.66464]),4) 1200 1201 def test_keip_zeros(self): 1202 keip = special.keip_zeros(5) 1203 assert_array_almost_equal(keip,array([4.93181, 1204 9.40405, 1205 13.85827, 1206 18.30717, 1207 22.75379]),4) 1208 1209 # numbers come from 9.9 of A&S pg. 381 1210 def test_kelvin_zeros(self): 1211 tmp = special.kelvin_zeros(5) 1212 berz,beiz,kerz,keiz,berpz,beipz,kerpz,keipz = tmp 1213 assert_array_almost_equal(berz,array([2.84892, 1214 7.23883, 1215 11.67396, 1216 16.11356, 1217 20.55463]),4) 1218 assert_array_almost_equal(beiz,array([5.02622, 1219 9.45541, 1220 13.89349, 1221 18.33398, 1222 22.77544]),4) 1223 assert_array_almost_equal(kerz,array([1.71854, 1224 6.12728, 1225 10.56294, 1226 15.00269, 1227 19.44382]),4) 1228 assert_array_almost_equal(keiz,array([3.91467, 1229 8.34422, 1230 12.78256, 1231 17.22314, 1232 21.66464]),4) 1233 assert_array_almost_equal(berpz,array([6.03871, 1234 10.51364, 1235 14.96844, 1236 19.41758, 1237 23.86430]),4) 1238 assert_array_almost_equal(beipz,array([3.77267, 1239 # table from 1927 had 3.77320 1240 # but this is more accurate 1241 8.28099, 1242 12.74215, 1243 17.19343, 1244 21.64114]),4) 1245 assert_array_almost_equal(kerpz,array([2.66584, 1246 7.17212, 1247 11.63218, 1248 16.08312, 1249 20.53068]),4) 1250 assert_array_almost_equal(keipz,array([4.93181, 1251 9.40405, 1252 13.85827, 1253 18.30717, 1254 22.75379]),4) 1255 1256 def test_ker_zeros(self): 1257 ker = special.ker_zeros(5) 1258 assert_array_almost_equal(ker,array([1.71854, 1259 6.12728, 1260 10.56294, 1261 15.00269, 1262 19.44381]),4) 1263 1264 def test_kerp_zeros(self): 1265 kerp = special.kerp_zeros(5) 1266 assert_array_almost_equal(kerp,array([2.66584, 1267 7.17212, 1268 11.63218, 1269 16.08312, 1270 20.53068]),4) 1271 1272 1273class TestBernoulli: 1274 def test_bernoulli(self): 1275 brn = special.bernoulli(5) 1276 assert_array_almost_equal(brn,array([1.0000, 1277 -0.5000, 1278 0.1667, 1279 0.0000, 1280 -0.0333, 1281 0.0000]),4) 1282 1283 1284class TestBeta: 1285 def test_beta(self): 1286 bet = special.beta(2,4) 1287 betg = (special.gamma(2)*special.gamma(4))/special.gamma(6) 1288 assert_almost_equal(bet,betg,8) 1289 1290 def test_betaln(self): 1291 betln = special.betaln(2,4) 1292 bet = log(abs(special.beta(2,4))) 1293 assert_almost_equal(betln,bet,8) 1294 1295 def test_betainc(self): 1296 btinc = special.betainc(1,1,.2) 1297 assert_almost_equal(btinc,0.2,8) 1298 1299 def test_betaincinv(self): 1300 y = special.betaincinv(2,4,.5) 1301 comp = special.betainc(2,4,y) 1302 assert_almost_equal(comp,.5,5) 1303 1304 1305class TestCombinatorics: 1306 def test_comb(self): 1307 assert_array_almost_equal(special.comb([10, 10], [3, 4]), [120., 210.]) 1308 assert_almost_equal(special.comb(10, 3), 120.) 1309 assert_equal(special.comb(10, 3, exact=True), 120) 1310 assert_equal(special.comb(10, 3, exact=True, repetition=True), 220) 1311 1312 assert_allclose([special.comb(20, k, exact=True) for k in range(21)], 1313 special.comb(20, list(range(21))), atol=1e-15) 1314 1315 ii = np.iinfo(int).max + 1 1316 assert_equal(special.comb(ii, ii-1, exact=True), ii) 1317 1318 expected = 100891344545564193334812497256 1319 assert_equal(special.comb(100, 50, exact=True), expected) 1320 1321 def test_comb_with_np_int64(self): 1322 n = 70 1323 k = 30 1324 np_n = np.int64(n) 1325 np_k = np.int64(k) 1326 assert_equal(special.comb(np_n, np_k, exact=True), 1327 special.comb(n, k, exact=True)) 1328 1329 def test_comb_zeros(self): 1330 assert_equal(special.comb(2, 3, exact=True), 0) 1331 assert_equal(special.comb(-1, 3, exact=True), 0) 1332 assert_equal(special.comb(2, -1, exact=True), 0) 1333 assert_equal(special.comb(2, -1, exact=False), 0) 1334 assert_array_almost_equal(special.comb([2, -1, 2, 10], [3, 3, -1, 3]), 1335 [0., 0., 0., 120.]) 1336 1337 def test_perm(self): 1338 assert_array_almost_equal(special.perm([10, 10], [3, 4]), [720., 5040.]) 1339 assert_almost_equal(special.perm(10, 3), 720.) 1340 assert_equal(special.perm(10, 3, exact=True), 720) 1341 1342 def test_perm_zeros(self): 1343 assert_equal(special.perm(2, 3, exact=True), 0) 1344 assert_equal(special.perm(-1, 3, exact=True), 0) 1345 assert_equal(special.perm(2, -1, exact=True), 0) 1346 assert_equal(special.perm(2, -1, exact=False), 0) 1347 assert_array_almost_equal(special.perm([2, -1, 2, 10], [3, 3, -1, 3]), 1348 [0., 0., 0., 720.]) 1349 1350 1351class TestTrigonometric: 1352 def test_cbrt(self): 1353 cb = special.cbrt(27) 1354 cbrl = 27**(1.0/3.0) 1355 assert_approx_equal(cb,cbrl) 1356 1357 def test_cbrtmore(self): 1358 cb1 = special.cbrt(27.9) 1359 cbrl1 = 27.9**(1.0/3.0) 1360 assert_almost_equal(cb1,cbrl1,8) 1361 1362 def test_cosdg(self): 1363 cdg = special.cosdg(90) 1364 cdgrl = cos(pi/2.0) 1365 assert_almost_equal(cdg,cdgrl,8) 1366 1367 def test_cosdgmore(self): 1368 cdgm = special.cosdg(30) 1369 cdgmrl = cos(pi/6.0) 1370 assert_almost_equal(cdgm,cdgmrl,8) 1371 1372 def test_cosm1(self): 1373 cs = (special.cosm1(0),special.cosm1(.3),special.cosm1(pi/10)) 1374 csrl = (cos(0)-1,cos(.3)-1,cos(pi/10)-1) 1375 assert_array_almost_equal(cs,csrl,8) 1376 1377 def test_cotdg(self): 1378 ct = special.cotdg(30) 1379 ctrl = tan(pi/6.0)**(-1) 1380 assert_almost_equal(ct,ctrl,8) 1381 1382 def test_cotdgmore(self): 1383 ct1 = special.cotdg(45) 1384 ctrl1 = tan(pi/4.0)**(-1) 1385 assert_almost_equal(ct1,ctrl1,8) 1386 1387 def test_specialpoints(self): 1388 assert_almost_equal(special.cotdg(45), 1.0, 14) 1389 assert_almost_equal(special.cotdg(-45), -1.0, 14) 1390 assert_almost_equal(special.cotdg(90), 0.0, 14) 1391 assert_almost_equal(special.cotdg(-90), 0.0, 14) 1392 assert_almost_equal(special.cotdg(135), -1.0, 14) 1393 assert_almost_equal(special.cotdg(-135), 1.0, 14) 1394 assert_almost_equal(special.cotdg(225), 1.0, 14) 1395 assert_almost_equal(special.cotdg(-225), -1.0, 14) 1396 assert_almost_equal(special.cotdg(270), 0.0, 14) 1397 assert_almost_equal(special.cotdg(-270), 0.0, 14) 1398 assert_almost_equal(special.cotdg(315), -1.0, 14) 1399 assert_almost_equal(special.cotdg(-315), 1.0, 14) 1400 assert_almost_equal(special.cotdg(765), 1.0, 14) 1401 1402 def test_sinc(self): 1403 # the sinc implementation and more extensive sinc tests are in numpy 1404 assert_array_equal(special.sinc([0]), 1) 1405 assert_equal(special.sinc(0.0), 1.0) 1406 1407 def test_sindg(self): 1408 sn = special.sindg(90) 1409 assert_equal(sn,1.0) 1410 1411 def test_sindgmore(self): 1412 snm = special.sindg(30) 1413 snmrl = sin(pi/6.0) 1414 assert_almost_equal(snm,snmrl,8) 1415 snm1 = special.sindg(45) 1416 snmrl1 = sin(pi/4.0) 1417 assert_almost_equal(snm1,snmrl1,8) 1418 1419 1420class TestTandg: 1421 1422 def test_tandg(self): 1423 tn = special.tandg(30) 1424 tnrl = tan(pi/6.0) 1425 assert_almost_equal(tn,tnrl,8) 1426 1427 def test_tandgmore(self): 1428 tnm = special.tandg(45) 1429 tnmrl = tan(pi/4.0) 1430 assert_almost_equal(tnm,tnmrl,8) 1431 tnm1 = special.tandg(60) 1432 tnmrl1 = tan(pi/3.0) 1433 assert_almost_equal(tnm1,tnmrl1,8) 1434 1435 def test_specialpoints(self): 1436 assert_almost_equal(special.tandg(0), 0.0, 14) 1437 assert_almost_equal(special.tandg(45), 1.0, 14) 1438 assert_almost_equal(special.tandg(-45), -1.0, 14) 1439 assert_almost_equal(special.tandg(135), -1.0, 14) 1440 assert_almost_equal(special.tandg(-135), 1.0, 14) 1441 assert_almost_equal(special.tandg(180), 0.0, 14) 1442 assert_almost_equal(special.tandg(-180), 0.0, 14) 1443 assert_almost_equal(special.tandg(225), 1.0, 14) 1444 assert_almost_equal(special.tandg(-225), -1.0, 14) 1445 assert_almost_equal(special.tandg(315), -1.0, 14) 1446 assert_almost_equal(special.tandg(-315), 1.0, 14) 1447 1448 1449class TestEllip: 1450 def test_ellipj_nan(self): 1451 """Regression test for #912.""" 1452 special.ellipj(0.5, np.nan) 1453 1454 def test_ellipj(self): 1455 el = special.ellipj(0.2,0) 1456 rel = [sin(0.2),cos(0.2),1.0,0.20] 1457 assert_array_almost_equal(el,rel,13) 1458 1459 def test_ellipk(self): 1460 elk = special.ellipk(.2) 1461 assert_almost_equal(elk,1.659623598610528,11) 1462 1463 assert_equal(special.ellipkm1(0.0), np.inf) 1464 assert_equal(special.ellipkm1(1.0), pi/2) 1465 assert_equal(special.ellipkm1(np.inf), 0.0) 1466 assert_equal(special.ellipkm1(np.nan), np.nan) 1467 assert_equal(special.ellipkm1(-1), np.nan) 1468 assert_allclose(special.ellipk(-10), 0.7908718902387385) 1469 1470 def test_ellipkinc(self): 1471 elkinc = special.ellipkinc(pi/2,.2) 1472 elk = special.ellipk(0.2) 1473 assert_almost_equal(elkinc,elk,15) 1474 alpha = 20*pi/180 1475 phi = 45*pi/180 1476 m = sin(alpha)**2 1477 elkinc = special.ellipkinc(phi,m) 1478 assert_almost_equal(elkinc,0.79398143,8) 1479 # From pg. 614 of A & S 1480 1481 assert_equal(special.ellipkinc(pi/2, 0.0), pi/2) 1482 assert_equal(special.ellipkinc(pi/2, 1.0), np.inf) 1483 assert_equal(special.ellipkinc(pi/2, -np.inf), 0.0) 1484 assert_equal(special.ellipkinc(pi/2, np.nan), np.nan) 1485 assert_equal(special.ellipkinc(pi/2, 2), np.nan) 1486 assert_equal(special.ellipkinc(0, 0.5), 0.0) 1487 assert_equal(special.ellipkinc(np.inf, 0.5), np.inf) 1488 assert_equal(special.ellipkinc(-np.inf, 0.5), -np.inf) 1489 assert_equal(special.ellipkinc(np.inf, np.inf), np.nan) 1490 assert_equal(special.ellipkinc(np.inf, -np.inf), np.nan) 1491 assert_equal(special.ellipkinc(-np.inf, -np.inf), np.nan) 1492 assert_equal(special.ellipkinc(-np.inf, np.inf), np.nan) 1493 assert_equal(special.ellipkinc(np.nan, 0.5), np.nan) 1494 assert_equal(special.ellipkinc(np.nan, np.nan), np.nan) 1495 1496 assert_allclose(special.ellipkinc(0.38974112035318718, 1), 0.4, rtol=1e-14) 1497 assert_allclose(special.ellipkinc(1.5707, -10), 0.79084284661724946) 1498 1499 def test_ellipkinc_2(self): 1500 # Regression test for gh-3550 1501 # ellipkinc(phi, mbad) was NaN and mvals[2:6] were twice the correct value 1502 mbad = 0.68359375000000011 1503 phi = 0.9272952180016123 1504 m = np.nextafter(mbad, 0) 1505 mvals = [] 1506 for j in range(10): 1507 mvals.append(m) 1508 m = np.nextafter(m, 1) 1509 f = special.ellipkinc(phi, mvals) 1510 assert_array_almost_equal_nulp(f, np.full_like(f, 1.0259330100195334), 1) 1511 # this bug also appears at phi + n * pi for at least small n 1512 f1 = special.ellipkinc(phi + pi, mvals) 1513 assert_array_almost_equal_nulp(f1, np.full_like(f1, 5.1296650500976675), 2) 1514 1515 def test_ellipkinc_singular(self): 1516 # ellipkinc(phi, 1) has closed form and is finite only for phi in (-pi/2, pi/2) 1517 xlog = np.logspace(-300, -17, 25) 1518 xlin = np.linspace(1e-17, 0.1, 25) 1519 xlin2 = np.linspace(0.1, pi/2, 25, endpoint=False) 1520 1521 assert_allclose(special.ellipkinc(xlog, 1), np.arcsinh(np.tan(xlog)), rtol=1e14) 1522 assert_allclose(special.ellipkinc(xlin, 1), np.arcsinh(np.tan(xlin)), rtol=1e14) 1523 assert_allclose(special.ellipkinc(xlin2, 1), np.arcsinh(np.tan(xlin2)), rtol=1e14) 1524 assert_equal(special.ellipkinc(np.pi/2, 1), np.inf) 1525 assert_allclose(special.ellipkinc(-xlog, 1), np.arcsinh(np.tan(-xlog)), rtol=1e14) 1526 assert_allclose(special.ellipkinc(-xlin, 1), np.arcsinh(np.tan(-xlin)), rtol=1e14) 1527 assert_allclose(special.ellipkinc(-xlin2, 1), np.arcsinh(np.tan(-xlin2)), rtol=1e14) 1528 assert_equal(special.ellipkinc(-np.pi/2, 1), np.inf) 1529 1530 def test_ellipe(self): 1531 ele = special.ellipe(.2) 1532 assert_almost_equal(ele,1.4890350580958529,8) 1533 1534 assert_equal(special.ellipe(0.0), pi/2) 1535 assert_equal(special.ellipe(1.0), 1.0) 1536 assert_equal(special.ellipe(-np.inf), np.inf) 1537 assert_equal(special.ellipe(np.nan), np.nan) 1538 assert_equal(special.ellipe(2), np.nan) 1539 assert_allclose(special.ellipe(-10), 3.6391380384177689) 1540 1541 def test_ellipeinc(self): 1542 eleinc = special.ellipeinc(pi/2,.2) 1543 ele = special.ellipe(0.2) 1544 assert_almost_equal(eleinc,ele,14) 1545 # pg 617 of A & S 1546 alpha, phi = 52*pi/180,35*pi/180 1547 m = sin(alpha)**2 1548 eleinc = special.ellipeinc(phi,m) 1549 assert_almost_equal(eleinc, 0.58823065, 8) 1550 1551 assert_equal(special.ellipeinc(pi/2, 0.0), pi/2) 1552 assert_equal(special.ellipeinc(pi/2, 1.0), 1.0) 1553 assert_equal(special.ellipeinc(pi/2, -np.inf), np.inf) 1554 assert_equal(special.ellipeinc(pi/2, np.nan), np.nan) 1555 assert_equal(special.ellipeinc(pi/2, 2), np.nan) 1556 assert_equal(special.ellipeinc(0, 0.5), 0.0) 1557 assert_equal(special.ellipeinc(np.inf, 0.5), np.inf) 1558 assert_equal(special.ellipeinc(-np.inf, 0.5), -np.inf) 1559 assert_equal(special.ellipeinc(np.inf, -np.inf), np.inf) 1560 assert_equal(special.ellipeinc(-np.inf, -np.inf), -np.inf) 1561 assert_equal(special.ellipeinc(np.inf, np.inf), np.nan) 1562 assert_equal(special.ellipeinc(-np.inf, np.inf), np.nan) 1563 assert_equal(special.ellipeinc(np.nan, 0.5), np.nan) 1564 assert_equal(special.ellipeinc(np.nan, np.nan), np.nan) 1565 assert_allclose(special.ellipeinc(1.5707, -10), 3.6388185585822876) 1566 1567 def test_ellipeinc_2(self): 1568 # Regression test for gh-3550 1569 # ellipeinc(phi, mbad) was NaN and mvals[2:6] were twice the correct value 1570 mbad = 0.68359375000000011 1571 phi = 0.9272952180016123 1572 m = np.nextafter(mbad, 0) 1573 mvals = [] 1574 for j in range(10): 1575 mvals.append(m) 1576 m = np.nextafter(m, 1) 1577 f = special.ellipeinc(phi, mvals) 1578 assert_array_almost_equal_nulp(f, np.full_like(f, 0.84442884574781019), 2) 1579 # this bug also appears at phi + n * pi for at least small n 1580 f1 = special.ellipeinc(phi + pi, mvals) 1581 assert_array_almost_equal_nulp(f1, np.full_like(f1, 3.3471442287390509), 4) 1582 1583 1584class TestErf: 1585 1586 def test_erf(self): 1587 er = special.erf(.25) 1588 assert_almost_equal(er,0.2763263902,8) 1589 1590 def test_erf_zeros(self): 1591 erz = special.erf_zeros(5) 1592 erzr = array([1.45061616+1.88094300j, 1593 2.24465928+2.61657514j, 1594 2.83974105+3.17562810j, 1595 3.33546074+3.64617438j, 1596 3.76900557+4.06069723j]) 1597 assert_array_almost_equal(erz,erzr,4) 1598 1599 def _check_variant_func(self, func, other_func, rtol, atol=0): 1600 np.random.seed(1234) 1601 n = 10000 1602 x = np.random.pareto(0.02, n) * (2*np.random.randint(0, 2, n) - 1) 1603 y = np.random.pareto(0.02, n) * (2*np.random.randint(0, 2, n) - 1) 1604 z = x + 1j*y 1605 1606 with np.errstate(all='ignore'): 1607 w = other_func(z) 1608 w_real = other_func(x).real 1609 1610 mask = np.isfinite(w) 1611 w = w[mask] 1612 z = z[mask] 1613 1614 mask = np.isfinite(w_real) 1615 w_real = w_real[mask] 1616 x = x[mask] 1617 1618 # test both real and complex variants 1619 assert_func_equal(func, w, z, rtol=rtol, atol=atol) 1620 assert_func_equal(func, w_real, x, rtol=rtol, atol=atol) 1621 1622 def test_erfc_consistent(self): 1623 self._check_variant_func( 1624 cephes.erfc, 1625 lambda z: 1 - cephes.erf(z), 1626 rtol=1e-12, 1627 atol=1e-14 # <- the test function loses precision 1628 ) 1629 1630 def test_erfcx_consistent(self): 1631 self._check_variant_func( 1632 cephes.erfcx, 1633 lambda z: np.exp(z*z) * cephes.erfc(z), 1634 rtol=1e-12 1635 ) 1636 1637 def test_erfi_consistent(self): 1638 self._check_variant_func( 1639 cephes.erfi, 1640 lambda z: -1j * cephes.erf(1j*z), 1641 rtol=1e-12 1642 ) 1643 1644 def test_dawsn_consistent(self): 1645 self._check_variant_func( 1646 cephes.dawsn, 1647 lambda z: sqrt(pi)/2 * np.exp(-z*z) * cephes.erfi(z), 1648 rtol=1e-12 1649 ) 1650 1651 def test_erf_nan_inf(self): 1652 vals = [np.nan, -np.inf, np.inf] 1653 expected = [np.nan, -1, 1] 1654 assert_allclose(special.erf(vals), expected, rtol=1e-15) 1655 1656 def test_erfc_nan_inf(self): 1657 vals = [np.nan, -np.inf, np.inf] 1658 expected = [np.nan, 2, 0] 1659 assert_allclose(special.erfc(vals), expected, rtol=1e-15) 1660 1661 def test_erfcx_nan_inf(self): 1662 vals = [np.nan, -np.inf, np.inf] 1663 expected = [np.nan, np.inf, 0] 1664 assert_allclose(special.erfcx(vals), expected, rtol=1e-15) 1665 1666 def test_erfi_nan_inf(self): 1667 vals = [np.nan, -np.inf, np.inf] 1668 expected = [np.nan, -np.inf, np.inf] 1669 assert_allclose(special.erfi(vals), expected, rtol=1e-15) 1670 1671 def test_dawsn_nan_inf(self): 1672 vals = [np.nan, -np.inf, np.inf] 1673 expected = [np.nan, -0.0, 0.0] 1674 assert_allclose(special.dawsn(vals), expected, rtol=1e-15) 1675 1676 def test_wofz_nan_inf(self): 1677 vals = [np.nan, -np.inf, np.inf] 1678 expected = [np.nan + np.nan * 1.j, 0.-0.j, 0.+0.j] 1679 assert_allclose(special.wofz(vals), expected, rtol=1e-15) 1680 1681 1682class TestEuler: 1683 def test_euler(self): 1684 eu0 = special.euler(0) 1685 eu1 = special.euler(1) 1686 eu2 = special.euler(2) # just checking segfaults 1687 assert_allclose(eu0, [1], rtol=1e-15) 1688 assert_allclose(eu1, [1, 0], rtol=1e-15) 1689 assert_allclose(eu2, [1, 0, -1], rtol=1e-15) 1690 eu24 = special.euler(24) 1691 mathworld = [1,1,5,61,1385,50521,2702765,199360981, 1692 19391512145,2404879675441, 1693 370371188237525,69348874393137901, 1694 15514534163557086905] 1695 correct = zeros((25,),'d') 1696 for k in range(0,13): 1697 if (k % 2): 1698 correct[2*k] = -float(mathworld[k]) 1699 else: 1700 correct[2*k] = float(mathworld[k]) 1701 with np.errstate(all='ignore'): 1702 err = nan_to_num((eu24-correct)/correct) 1703 errmax = max(err) 1704 assert_almost_equal(errmax, 0.0, 14) 1705 1706 1707class TestExp: 1708 def test_exp2(self): 1709 ex = special.exp2(2) 1710 exrl = 2**2 1711 assert_equal(ex,exrl) 1712 1713 def test_exp2more(self): 1714 exm = special.exp2(2.5) 1715 exmrl = 2**(2.5) 1716 assert_almost_equal(exm,exmrl,8) 1717 1718 def test_exp10(self): 1719 ex = special.exp10(2) 1720 exrl = 10**2 1721 assert_approx_equal(ex,exrl) 1722 1723 def test_exp10more(self): 1724 exm = special.exp10(2.5) 1725 exmrl = 10**(2.5) 1726 assert_almost_equal(exm,exmrl,8) 1727 1728 def test_expm1(self): 1729 ex = (special.expm1(2),special.expm1(3),special.expm1(4)) 1730 exrl = (exp(2)-1,exp(3)-1,exp(4)-1) 1731 assert_array_almost_equal(ex,exrl,8) 1732 1733 def test_expm1more(self): 1734 ex1 = (special.expm1(2),special.expm1(2.1),special.expm1(2.2)) 1735 exrl1 = (exp(2)-1,exp(2.1)-1,exp(2.2)-1) 1736 assert_array_almost_equal(ex1,exrl1,8) 1737 1738 1739class TestFactorialFunctions: 1740 def test_factorial(self): 1741 # Some known values, float math 1742 assert_array_almost_equal(special.factorial(0), 1) 1743 assert_array_almost_equal(special.factorial(1), 1) 1744 assert_array_almost_equal(special.factorial(2), 2) 1745 assert_array_almost_equal([6., 24., 120.], 1746 special.factorial([3, 4, 5], exact=False)) 1747 assert_array_almost_equal(special.factorial([[5, 3], [4, 3]]), 1748 [[120, 6], [24, 6]]) 1749 1750 # Some known values, integer math 1751 assert_equal(special.factorial(0, exact=True), 1) 1752 assert_equal(special.factorial(1, exact=True), 1) 1753 assert_equal(special.factorial(2, exact=True), 2) 1754 assert_equal(special.factorial(5, exact=True), 120) 1755 assert_equal(special.factorial(15, exact=True), 1307674368000) 1756 1757 # ndarray shape is maintained 1758 assert_equal(special.factorial([7, 4, 15, 10], exact=True), 1759 [5040, 24, 1307674368000, 3628800]) 1760 1761 assert_equal(special.factorial([[5, 3], [4, 3]], True), 1762 [[120, 6], [24, 6]]) 1763 1764 # object arrays 1765 assert_equal(special.factorial(np.arange(-3, 22), True), 1766 special.factorial(np.arange(-3, 22), False)) 1767 1768 # int64 array 1769 assert_equal(special.factorial(np.arange(-3, 15), True), 1770 special.factorial(np.arange(-3, 15), False)) 1771 1772 # int32 array 1773 assert_equal(special.factorial(np.arange(-3, 5), True), 1774 special.factorial(np.arange(-3, 5), False)) 1775 1776 # Consistent output for n < 0 1777 for exact in (True, False): 1778 assert_array_equal(0, special.factorial(-3, exact)) 1779 assert_array_equal([1, 2, 0, 0], 1780 special.factorial([1, 2, -5, -4], exact)) 1781 1782 for n in range(0, 22): 1783 # Compare all with math.factorial 1784 correct = math.factorial(n) 1785 assert_array_equal(correct, special.factorial(n, True)) 1786 assert_array_equal(correct, special.factorial([n], True)[0]) 1787 1788 assert_allclose(float(correct), special.factorial(n, False)) 1789 assert_allclose(float(correct), special.factorial([n], False)[0]) 1790 1791 # Compare exact=True vs False, scalar vs array 1792 assert_array_equal(special.factorial(n, True), 1793 special.factorial(n, False)) 1794 1795 assert_array_equal(special.factorial([n], True), 1796 special.factorial([n], False)) 1797 1798 @pytest.mark.parametrize('x, exact', [ 1799 (1, True), 1800 (1, False), 1801 (np.array(1), True), 1802 (np.array(1), False), 1803 ]) 1804 def test_factorial_0d_return_type(self, x, exact): 1805 assert np.isscalar(special.factorial(x, exact=exact)) 1806 1807 def test_factorial2(self): 1808 assert_array_almost_equal([105., 384., 945.], 1809 special.factorial2([7, 8, 9], exact=False)) 1810 assert_equal(special.factorial2(7, exact=True), 105) 1811 1812 def test_factorialk(self): 1813 assert_equal(special.factorialk(5, 1, exact=True), 120) 1814 assert_equal(special.factorialk(5, 3, exact=True), 10) 1815 1816 @pytest.mark.parametrize('x, exact', [ 1817 (np.nan, True), 1818 (np.nan, False), 1819 (np.array([np.nan]), True), 1820 (np.array([np.nan]), False), 1821 ]) 1822 def test_nan_inputs(self, x, exact): 1823 result = special.factorial(x, exact=exact) 1824 assert_(np.isnan(result)) 1825 1826 # GH-13122: special.factorial() argument should be an array of integers. 1827 # On Python 3.10, math.factorial() reject float. 1828 # On Python 3.9, a DeprecationWarning is emitted. 1829 # A numpy array casts all integers to float if the array contains a 1830 # single NaN. 1831 @pytest.mark.skipif(sys.version_info >= (3, 10), 1832 reason="Python 3.10+ math.factorial() requires int") 1833 def test_mixed_nan_inputs(self): 1834 x = np.array([np.nan, 1, 2, 3, np.nan]) 1835 with suppress_warnings() as sup: 1836 sup.filter(DeprecationWarning, "Using factorial\\(\\) with floats is deprecated") 1837 result = special.factorial(x, exact=True) 1838 assert_equal(np.array([np.nan, 1, 2, 6, np.nan]), result) 1839 result = special.factorial(x, exact=False) 1840 assert_equal(np.array([np.nan, 1, 2, 6, np.nan]), result) 1841 1842 1843class TestFresnel: 1844 def test_fresnel(self): 1845 frs = array(special.fresnel(.5)) 1846 assert_array_almost_equal(frs,array([0.064732432859999287, 0.49234422587144644]),8) 1847 1848 def test_fresnel_inf1(self): 1849 frs = special.fresnel(np.inf) 1850 assert_equal(frs, (0.5, 0.5)) 1851 1852 def test_fresnel_inf2(self): 1853 frs = special.fresnel(-np.inf) 1854 assert_equal(frs, (-0.5, -0.5)) 1855 1856 # values from pg 329 Table 7.11 of A & S 1857 # slightly corrected in 4th decimal place 1858 def test_fresnel_zeros(self): 1859 szo, czo = special.fresnel_zeros(5) 1860 assert_array_almost_equal(szo, 1861 array([2.0093+0.2885j, 1862 2.8335+0.2443j, 1863 3.4675+0.2185j, 1864 4.0026+0.2009j, 1865 4.4742+0.1877j]),3) 1866 assert_array_almost_equal(czo, 1867 array([1.7437+0.3057j, 1868 2.6515+0.2529j, 1869 3.3204+0.2240j, 1870 3.8757+0.2047j, 1871 4.3611+0.1907j]),3) 1872 vals1 = special.fresnel(szo)[0] 1873 vals2 = special.fresnel(czo)[1] 1874 assert_array_almost_equal(vals1,0,14) 1875 assert_array_almost_equal(vals2,0,14) 1876 1877 def test_fresnelc_zeros(self): 1878 szo, czo = special.fresnel_zeros(6) 1879 frc = special.fresnelc_zeros(6) 1880 assert_array_almost_equal(frc,czo,12) 1881 1882 def test_fresnels_zeros(self): 1883 szo, czo = special.fresnel_zeros(5) 1884 frs = special.fresnels_zeros(5) 1885 assert_array_almost_equal(frs,szo,12) 1886 1887 1888class TestGamma: 1889 def test_gamma(self): 1890 gam = special.gamma(5) 1891 assert_equal(gam,24.0) 1892 1893 def test_gammaln(self): 1894 gamln = special.gammaln(3) 1895 lngam = log(special.gamma(3)) 1896 assert_almost_equal(gamln,lngam,8) 1897 1898 def test_gammainccinv(self): 1899 gccinv = special.gammainccinv(.5,.5) 1900 gcinv = special.gammaincinv(.5,.5) 1901 assert_almost_equal(gccinv,gcinv,8) 1902 1903 @with_special_errors 1904 def test_gammaincinv(self): 1905 y = special.gammaincinv(.4,.4) 1906 x = special.gammainc(.4,y) 1907 assert_almost_equal(x,0.4,1) 1908 y = special.gammainc(10, 0.05) 1909 x = special.gammaincinv(10, 2.5715803516000736e-20) 1910 assert_almost_equal(0.05, x, decimal=10) 1911 assert_almost_equal(y, 2.5715803516000736e-20, decimal=10) 1912 x = special.gammaincinv(50, 8.20754777388471303050299243573393e-18) 1913 assert_almost_equal(11.0, x, decimal=10) 1914 1915 @with_special_errors 1916 def test_975(self): 1917 # Regression test for ticket #975 -- switch point in algorithm 1918 # check that things work OK at the point, immediately next floats 1919 # around it, and a bit further away 1920 pts = [0.25, 1921 np.nextafter(0.25, 0), 0.25 - 1e-12, 1922 np.nextafter(0.25, 1), 0.25 + 1e-12] 1923 for xp in pts: 1924 y = special.gammaincinv(.4, xp) 1925 x = special.gammainc(0.4, y) 1926 assert_allclose(x, xp, rtol=1e-12) 1927 1928 def test_rgamma(self): 1929 rgam = special.rgamma(8) 1930 rlgam = 1/special.gamma(8) 1931 assert_almost_equal(rgam,rlgam,8) 1932 1933 def test_infinity(self): 1934 assert_(np.isinf(special.gamma(-1))) 1935 assert_equal(special.rgamma(-1), 0) 1936 1937 1938class TestHankel: 1939 1940 def test_negv1(self): 1941 assert_almost_equal(special.hankel1(-3,2), -special.hankel1(3,2), 14) 1942 1943 def test_hankel1(self): 1944 hank1 = special.hankel1(1,.1) 1945 hankrl = (special.jv(1,.1) + special.yv(1,.1)*1j) 1946 assert_almost_equal(hank1,hankrl,8) 1947 1948 def test_negv1e(self): 1949 assert_almost_equal(special.hankel1e(-3,2), -special.hankel1e(3,2), 14) 1950 1951 def test_hankel1e(self): 1952 hank1e = special.hankel1e(1,.1) 1953 hankrle = special.hankel1(1,.1)*exp(-.1j) 1954 assert_almost_equal(hank1e,hankrle,8) 1955 1956 def test_negv2(self): 1957 assert_almost_equal(special.hankel2(-3,2), -special.hankel2(3,2), 14) 1958 1959 def test_hankel2(self): 1960 hank2 = special.hankel2(1,.1) 1961 hankrl2 = (special.jv(1,.1) - special.yv(1,.1)*1j) 1962 assert_almost_equal(hank2,hankrl2,8) 1963 1964 def test_neg2e(self): 1965 assert_almost_equal(special.hankel2e(-3,2), -special.hankel2e(3,2), 14) 1966 1967 def test_hankl2e(self): 1968 hank2e = special.hankel2e(1,.1) 1969 hankrl2e = special.hankel2e(1,.1) 1970 assert_almost_equal(hank2e,hankrl2e,8) 1971 1972 1973class TestHyper: 1974 def test_h1vp(self): 1975 h1 = special.h1vp(1,.1) 1976 h1real = (special.jvp(1,.1) + special.yvp(1,.1)*1j) 1977 assert_almost_equal(h1,h1real,8) 1978 1979 def test_h2vp(self): 1980 h2 = special.h2vp(1,.1) 1981 h2real = (special.jvp(1,.1) - special.yvp(1,.1)*1j) 1982 assert_almost_equal(h2,h2real,8) 1983 1984 def test_hyp0f1(self): 1985 # scalar input 1986 assert_allclose(special.hyp0f1(2.5, 0.5), 1.21482702689997, rtol=1e-12) 1987 assert_allclose(special.hyp0f1(2.5, 0), 1.0, rtol=1e-15) 1988 1989 # float input, expected values match mpmath 1990 x = special.hyp0f1(3.0, [-1.5, -1, 0, 1, 1.5]) 1991 expected = np.array([0.58493659229143, 0.70566805723127, 1.0, 1992 1.37789689539747, 1.60373685288480]) 1993 assert_allclose(x, expected, rtol=1e-12) 1994 1995 # complex input 1996 x = special.hyp0f1(3.0, np.array([-1.5, -1, 0, 1, 1.5]) + 0.j) 1997 assert_allclose(x, expected.astype(complex), rtol=1e-12) 1998 1999 # test broadcasting 2000 x1 = [0.5, 1.5, 2.5] 2001 x2 = [0, 1, 0.5] 2002 x = special.hyp0f1(x1, x2) 2003 expected = [1.0, 1.8134302039235093, 1.21482702689997] 2004 assert_allclose(x, expected, rtol=1e-12) 2005 x = special.hyp0f1(np.row_stack([x1] * 2), x2) 2006 assert_allclose(x, np.row_stack([expected] * 2), rtol=1e-12) 2007 assert_raises(ValueError, special.hyp0f1, 2008 np.row_stack([x1] * 3), [0, 1]) 2009 2010 def test_hyp0f1_gh5764(self): 2011 # Just checks the point that failed; there's a more systematic 2012 # test in test_mpmath 2013 res = special.hyp0f1(0.8, 0.5 + 0.5*1J) 2014 # The expected value was generated using mpmath 2015 assert_almost_equal(res, 1.6139719776441115 + 1J*0.80893054061790665) 2016 2017 def test_hyp1f1(self): 2018 hyp1 = special.hyp1f1(.1,.1,.3) 2019 assert_almost_equal(hyp1, 1.3498588075760032,7) 2020 2021 # test contributed by Moritz Deger (2008-05-29) 2022 # https://github.com/scipy/scipy/issues/1186 (Trac #659) 2023 2024 # reference data obtained from mathematica [ a, b, x, m(a,b,x)]: 2025 # produced with test_hyp1f1.nb 2026 ref_data = array([[-8.38132975e+00, -1.28436461e+01, -2.91081397e+01, 1.04178330e+04], 2027 [2.91076882e+00, -6.35234333e+00, -1.27083993e+01, 6.68132725e+00], 2028 [-1.42938258e+01, 1.80869131e-01, 1.90038728e+01, 1.01385897e+05], 2029 [5.84069088e+00, 1.33187908e+01, 2.91290106e+01, 1.59469411e+08], 2030 [-2.70433202e+01, -1.16274873e+01, -2.89582384e+01, 1.39900152e+24], 2031 [4.26344966e+00, -2.32701773e+01, 1.91635759e+01, 6.13816915e+21], 2032 [1.20514340e+01, -3.40260240e+00, 7.26832235e+00, 1.17696112e+13], 2033 [2.77372955e+01, -1.99424687e+00, 3.61332246e+00, 3.07419615e+13], 2034 [1.50310939e+01, -2.91198675e+01, -1.53581080e+01, -3.79166033e+02], 2035 [1.43995827e+01, 9.84311196e+00, 1.93204553e+01, 2.55836264e+10], 2036 [-4.08759686e+00, 1.34437025e+01, -1.42072843e+01, 1.70778449e+01], 2037 [8.05595738e+00, -1.31019838e+01, 1.52180721e+01, 3.06233294e+21], 2038 [1.81815804e+01, -1.42908793e+01, 9.57868793e+00, -2.84771348e+20], 2039 [-2.49671396e+01, 1.25082843e+01, -1.71562286e+01, 2.36290426e+07], 2040 [2.67277673e+01, 1.70315414e+01, 6.12701450e+00, 7.77917232e+03], 2041 [2.49565476e+01, 2.91694684e+01, 6.29622660e+00, 2.35300027e+02], 2042 [6.11924542e+00, -1.59943768e+00, 9.57009289e+00, 1.32906326e+11], 2043 [-1.47863653e+01, 2.41691301e+01, -1.89981821e+01, 2.73064953e+03], 2044 [2.24070483e+01, -2.93647433e+00, 8.19281432e+00, -6.42000372e+17], 2045 [8.04042600e-01, 1.82710085e+01, -1.97814534e+01, 5.48372441e-01], 2046 [1.39590390e+01, 1.97318686e+01, 2.37606635e+00, 5.51923681e+00], 2047 [-4.66640483e+00, -2.00237930e+01, 7.40365095e+00, 4.50310752e+00], 2048 [2.76821999e+01, -6.36563968e+00, 1.11533984e+01, -9.28725179e+23], 2049 [-2.56764457e+01, 1.24544906e+00, 1.06407572e+01, 1.25922076e+01], 2050 [3.20447808e+00, 1.30874383e+01, 2.26098014e+01, 2.03202059e+04], 2051 [-1.24809647e+01, 4.15137113e+00, -2.92265700e+01, 2.39621411e+08], 2052 [2.14778108e+01, -2.35162960e+00, -1.13758664e+01, 4.46882152e-01], 2053 [-9.85469168e+00, -3.28157680e+00, 1.67447548e+01, -1.07342390e+07], 2054 [1.08122310e+01, -2.47353236e+01, -1.15622349e+01, -2.91733796e+03], 2055 [-2.67933347e+01, -3.39100709e+00, 2.56006986e+01, -5.29275382e+09], 2056 [-8.60066776e+00, -8.02200924e+00, 1.07231926e+01, 1.33548320e+06], 2057 [-1.01724238e-01, -1.18479709e+01, -2.55407104e+01, 1.55436570e+00], 2058 [-3.93356771e+00, 2.11106818e+01, -2.57598485e+01, 2.13467840e+01], 2059 [3.74750503e+00, 1.55687633e+01, -2.92841720e+01, 1.43873509e-02], 2060 [6.99726781e+00, 2.69855571e+01, -1.63707771e+01, 3.08098673e-02], 2061 [-2.31996011e+01, 3.47631054e+00, 9.75119815e-01, 1.79971073e-02], 2062 [2.38951044e+01, -2.91460190e+01, -2.50774708e+00, 9.56934814e+00], 2063 [1.52730825e+01, 5.77062507e+00, 1.21922003e+01, 1.32345307e+09], 2064 [1.74673917e+01, 1.89723426e+01, 4.94903250e+00, 9.90859484e+01], 2065 [1.88971241e+01, 2.86255413e+01, 5.52360109e-01, 1.44165360e+00], 2066 [1.02002319e+01, -1.66855152e+01, -2.55426235e+01, 6.56481554e+02], 2067 [-1.79474153e+01, 1.22210200e+01, -1.84058212e+01, 8.24041812e+05], 2068 [-1.36147103e+01, 1.32365492e+00, -7.22375200e+00, 9.92446491e+05], 2069 [7.57407832e+00, 2.59738234e+01, -1.34139168e+01, 3.64037761e-02], 2070 [2.21110169e+00, 1.28012666e+01, 1.62529102e+01, 1.33433085e+02], 2071 [-2.64297569e+01, -1.63176658e+01, -1.11642006e+01, -2.44797251e+13], 2072 [-2.46622944e+01, -3.02147372e+00, 8.29159315e+00, -3.21799070e+05], 2073 [-1.37215095e+01, -1.96680183e+01, 2.91940118e+01, 3.21457520e+12], 2074 [-5.45566105e+00, 2.81292086e+01, 1.72548215e-01, 9.66973000e-01], 2075 [-1.55751298e+00, -8.65703373e+00, 2.68622026e+01, -3.17190834e+16], 2076 [2.45393609e+01, -2.70571903e+01, 1.96815505e+01, 1.80708004e+37], 2077 [5.77482829e+00, 1.53203143e+01, 2.50534322e+01, 1.14304242e+06], 2078 [-1.02626819e+01, 2.36887658e+01, -2.32152102e+01, 7.28965646e+02], 2079 [-1.30833446e+00, -1.28310210e+01, 1.87275544e+01, -9.33487904e+12], 2080 [5.83024676e+00, -1.49279672e+01, 2.44957538e+01, -7.61083070e+27], 2081 [-2.03130747e+01, 2.59641715e+01, -2.06174328e+01, 4.54744859e+04], 2082 [1.97684551e+01, -2.21410519e+01, -2.26728740e+01, 3.53113026e+06], 2083 [2.73673444e+01, 2.64491725e+01, 1.57599882e+01, 1.07385118e+07], 2084 [5.73287971e+00, 1.21111904e+01, 1.33080171e+01, 2.63220467e+03], 2085 [-2.82751072e+01, 2.08605881e+01, 9.09838900e+00, -6.60957033e-07], 2086 [1.87270691e+01, -1.74437016e+01, 1.52413599e+01, 6.59572851e+27], 2087 [6.60681457e+00, -2.69449855e+00, 9.78972047e+00, -2.38587870e+12], 2088 [1.20895561e+01, -2.51355765e+01, 2.30096101e+01, 7.58739886e+32], 2089 [-2.44682278e+01, 2.10673441e+01, -1.36705538e+01, 4.54213550e+04], 2090 [-4.50665152e+00, 3.72292059e+00, -4.83403707e+00, 2.68938214e+01], 2091 [-7.46540049e+00, -1.08422222e+01, -1.72203805e+01, -2.09402162e+02], 2092 [-2.00307551e+01, -7.50604431e+00, -2.78640020e+01, 4.15985444e+19], 2093 [1.99890876e+01, 2.20677419e+01, -2.51301778e+01, 1.23840297e-09], 2094 [2.03183823e+01, -7.66942559e+00, 2.10340070e+01, 1.46285095e+31], 2095 [-2.90315825e+00, -2.55785967e+01, -9.58779316e+00, 2.65714264e-01], 2096 [2.73960829e+01, -1.80097203e+01, -2.03070131e+00, 2.52908999e+02], 2097 [-2.11708058e+01, -2.70304032e+01, 2.48257944e+01, 3.09027527e+08], 2098 [2.21959758e+01, 4.00258675e+00, -1.62853977e+01, -9.16280090e-09], 2099 [1.61661840e+01, -2.26845150e+01, 2.17226940e+01, -8.24774394e+33], 2100 [-3.35030306e+00, 1.32670581e+00, 9.39711214e+00, -1.47303163e+01], 2101 [7.23720726e+00, -2.29763909e+01, 2.34709682e+01, -9.20711735e+29], 2102 [2.71013568e+01, 1.61951087e+01, -7.11388906e-01, 2.98750911e-01], 2103 [8.40057933e+00, -7.49665220e+00, 2.95587388e+01, 6.59465635e+29], 2104 [-1.51603423e+01, 1.94032322e+01, -7.60044357e+00, 1.05186941e+02], 2105 [-8.83788031e+00, -2.72018313e+01, 1.88269907e+00, 1.81687019e+00], 2106 [-1.87283712e+01, 5.87479570e+00, -1.91210203e+01, 2.52235612e+08], 2107 [-5.61338513e-01, 2.69490237e+01, 1.16660111e-01, 9.97567783e-01], 2108 [-5.44354025e+00, -1.26721408e+01, -4.66831036e+00, 1.06660735e-01], 2109 [-2.18846497e+00, 2.33299566e+01, 9.62564397e+00, 3.03842061e-01], 2110 [6.65661299e+00, -2.39048713e+01, 1.04191807e+01, 4.73700451e+13], 2111 [-2.57298921e+01, -2.60811296e+01, 2.74398110e+01, -5.32566307e+11], 2112 [-1.11431826e+01, -1.59420160e+01, -1.84880553e+01, -1.01514747e+02], 2113 [6.50301931e+00, 2.59859051e+01, -2.33270137e+01, 1.22760500e-02], 2114 [-1.94987891e+01, -2.62123262e+01, 3.90323225e+00, 1.71658894e+01], 2115 [7.26164601e+00, -1.41469402e+01, 2.81499763e+01, -2.50068329e+31], 2116 [-1.52424040e+01, 2.99719005e+01, -2.85753678e+01, 1.31906693e+04], 2117 [5.24149291e+00, -1.72807223e+01, 2.22129493e+01, 2.50748475e+25], 2118 [3.63207230e-01, -9.54120862e-02, -2.83874044e+01, 9.43854939e-01], 2119 [-2.11326457e+00, -1.25707023e+01, 1.17172130e+00, 1.20812698e+00], 2120 [2.48513582e+00, 1.03652647e+01, -1.84625148e+01, 6.47910997e-02], 2121 [2.65395942e+01, 2.74794672e+01, 1.29413428e+01, 2.89306132e+05], 2122 [-9.49445460e+00, 1.59930921e+01, -1.49596331e+01, 3.27574841e+02], 2123 [-5.89173945e+00, 9.96742426e+00, 2.60318889e+01, -3.15842908e-01], 2124 [-1.15387239e+01, -2.21433107e+01, -2.17686413e+01, 1.56724718e-01], 2125 [-5.30592244e+00, -2.42752190e+01, 1.29734035e+00, 1.31985534e+00]]) 2126 2127 for a,b,c,expected in ref_data: 2128 result = special.hyp1f1(a,b,c) 2129 assert_(abs(expected - result)/expected < 1e-4) 2130 2131 def test_hyp1f1_gh2957(self): 2132 hyp1 = special.hyp1f1(0.5, 1.5, -709.7827128933) 2133 hyp2 = special.hyp1f1(0.5, 1.5, -709.7827128934) 2134 assert_almost_equal(hyp1, hyp2, 12) 2135 2136 def test_hyp1f1_gh2282(self): 2137 hyp = special.hyp1f1(0.5, 1.5, -1000) 2138 assert_almost_equal(hyp, 0.028024956081989643, 12) 2139 2140 def test_hyp2f1(self): 2141 # a collection of special cases taken from AMS 55 2142 values = [[0.5, 1, 1.5, 0.2**2, 0.5/0.2*log((1+0.2)/(1-0.2))], 2143 [0.5, 1, 1.5, -0.2**2, 1./0.2*arctan(0.2)], 2144 [1, 1, 2, 0.2, -1/0.2*log(1-0.2)], 2145 [3, 3.5, 1.5, 0.2**2, 2146 0.5/0.2/(-5)*((1+0.2)**(-5)-(1-0.2)**(-5))], 2147 [-3, 3, 0.5, sin(0.2)**2, cos(2*3*0.2)], 2148 [3, 4, 8, 1, special.gamma(8)*special.gamma(8-4-3)/special.gamma(8-3)/special.gamma(8-4)], 2149 [3, 2, 3-2+1, -1, 1./2**3*sqrt(pi) * 2150 special.gamma(1+3-2)/special.gamma(1+0.5*3-2)/special.gamma(0.5+0.5*3)], 2151 [5, 2, 5-2+1, -1, 1./2**5*sqrt(pi) * 2152 special.gamma(1+5-2)/special.gamma(1+0.5*5-2)/special.gamma(0.5+0.5*5)], 2153 [4, 0.5+4, 1.5-2*4, -1./3, (8./9)**(-2*4)*special.gamma(4./3) * 2154 special.gamma(1.5-2*4)/special.gamma(3./2)/special.gamma(4./3-2*4)], 2155 # and some others 2156 # ticket #424 2157 [1.5, -0.5, 1.0, -10.0, 4.1300097765277476484], 2158 # negative integer a or b, with c-a-b integer and x > 0.9 2159 [-2,3,1,0.95,0.715], 2160 [2,-3,1,0.95,-0.007], 2161 [-6,3,1,0.95,0.0000810625], 2162 [2,-5,1,0.95,-0.000029375], 2163 # huge negative integers 2164 (10, -900, 10.5, 0.99, 1.91853705796607664803709475658e-24), 2165 (10, -900, -10.5, 0.99, 3.54279200040355710199058559155e-18), 2166 ] 2167 for i, (a, b, c, x, v) in enumerate(values): 2168 cv = special.hyp2f1(a, b, c, x) 2169 assert_almost_equal(cv, v, 8, err_msg='test #%d' % i) 2170 2171 def test_hyperu(self): 2172 val1 = special.hyperu(1,0.1,100) 2173 assert_almost_equal(val1,0.0098153,7) 2174 a,b = [0.3,0.6,1.2,-2.7],[1.5,3.2,-0.4,-3.2] 2175 a,b = asarray(a), asarray(b) 2176 z = 0.5 2177 hypu = special.hyperu(a,b,z) 2178 hprl = (pi/sin(pi*b))*(special.hyp1f1(a,b,z) / 2179 (special.gamma(1+a-b)*special.gamma(b)) - 2180 z**(1-b)*special.hyp1f1(1+a-b,2-b,z) 2181 / (special.gamma(a)*special.gamma(2-b))) 2182 assert_array_almost_equal(hypu,hprl,12) 2183 2184 def test_hyperu_gh2287(self): 2185 assert_almost_equal(special.hyperu(1, 1.5, 20.2), 2186 0.048360918656699191, 12) 2187 2188 2189class TestBessel: 2190 def test_itj0y0(self): 2191 it0 = array(special.itj0y0(.2)) 2192 assert_array_almost_equal(it0,array([0.19933433254006822, -0.34570883800412566]),8) 2193 2194 def test_it2j0y0(self): 2195 it2 = array(special.it2j0y0(.2)) 2196 assert_array_almost_equal(it2,array([0.0049937546274601858, -0.43423067011231614]),8) 2197 2198 def test_negv_iv(self): 2199 assert_equal(special.iv(3,2), special.iv(-3,2)) 2200 2201 def test_j0(self): 2202 oz = special.j0(.1) 2203 ozr = special.jn(0,.1) 2204 assert_almost_equal(oz,ozr,8) 2205 2206 def test_j1(self): 2207 o1 = special.j1(.1) 2208 o1r = special.jn(1,.1) 2209 assert_almost_equal(o1,o1r,8) 2210 2211 def test_jn(self): 2212 jnnr = special.jn(1,.2) 2213 assert_almost_equal(jnnr,0.099500832639235995,8) 2214 2215 def test_negv_jv(self): 2216 assert_almost_equal(special.jv(-3,2), -special.jv(3,2), 14) 2217 2218 def test_jv(self): 2219 values = [[0, 0.1, 0.99750156206604002], 2220 [2./3, 1e-8, 0.3239028506761532e-5], 2221 [2./3, 1e-10, 0.1503423854873779e-6], 2222 [3.1, 1e-10, 0.1711956265409013e-32], 2223 [2./3, 4.0, -0.2325440850267039], 2224 ] 2225 for i, (v, x, y) in enumerate(values): 2226 yc = special.jv(v, x) 2227 assert_almost_equal(yc, y, 8, err_msg='test #%d' % i) 2228 2229 def test_negv_jve(self): 2230 assert_almost_equal(special.jve(-3,2), -special.jve(3,2), 14) 2231 2232 def test_jve(self): 2233 jvexp = special.jve(1,.2) 2234 assert_almost_equal(jvexp,0.099500832639235995,8) 2235 jvexp1 = special.jve(1,.2+1j) 2236 z = .2+1j 2237 jvexpr = special.jv(1,z)*exp(-abs(z.imag)) 2238 assert_almost_equal(jvexp1,jvexpr,8) 2239 2240 def test_jn_zeros(self): 2241 jn0 = special.jn_zeros(0,5) 2242 jn1 = special.jn_zeros(1,5) 2243 assert_array_almost_equal(jn0,array([2.4048255577, 2244 5.5200781103, 2245 8.6537279129, 2246 11.7915344391, 2247 14.9309177086]),4) 2248 assert_array_almost_equal(jn1,array([3.83171, 2249 7.01559, 2250 10.17347, 2251 13.32369, 2252 16.47063]),4) 2253 2254 jn102 = special.jn_zeros(102,5) 2255 assert_allclose(jn102, array([110.89174935992040343, 2256 117.83464175788308398, 2257 123.70194191713507279, 2258 129.02417238949092824, 2259 134.00114761868422559]), rtol=1e-13) 2260 2261 jn301 = special.jn_zeros(301,5) 2262 assert_allclose(jn301, array([313.59097866698830153, 2263 323.21549776096288280, 2264 331.22338738656748796, 2265 338.39676338872084500, 2266 345.03284233056064157]), rtol=1e-13) 2267 2268 def test_jn_zeros_slow(self): 2269 jn0 = special.jn_zeros(0, 300) 2270 assert_allclose(jn0[260-1], 816.02884495068867280, rtol=1e-13) 2271 assert_allclose(jn0[280-1], 878.86068707124422606, rtol=1e-13) 2272 assert_allclose(jn0[300-1], 941.69253065317954064, rtol=1e-13) 2273 2274 jn10 = special.jn_zeros(10, 300) 2275 assert_allclose(jn10[260-1], 831.67668514305631151, rtol=1e-13) 2276 assert_allclose(jn10[280-1], 894.51275095371316931, rtol=1e-13) 2277 assert_allclose(jn10[300-1], 957.34826370866539775, rtol=1e-13) 2278 2279 jn3010 = special.jn_zeros(3010,5) 2280 assert_allclose(jn3010, array([3036.86590780927, 2281 3057.06598526482, 2282 3073.66360690272, 2283 3088.37736494778, 2284 3101.86438139042]), rtol=1e-8) 2285 2286 def test_jnjnp_zeros(self): 2287 jn = special.jn 2288 2289 def jnp(n, x): 2290 return (jn(n-1,x) - jn(n+1,x))/2 2291 for nt in range(1, 30): 2292 z, n, m, t = special.jnjnp_zeros(nt) 2293 for zz, nn, tt in zip(z, n, t): 2294 if tt == 0: 2295 assert_allclose(jn(nn, zz), 0, atol=1e-6) 2296 elif tt == 1: 2297 assert_allclose(jnp(nn, zz), 0, atol=1e-6) 2298 else: 2299 raise AssertionError("Invalid t return for nt=%d" % nt) 2300 2301 def test_jnp_zeros(self): 2302 jnp = special.jnp_zeros(1,5) 2303 assert_array_almost_equal(jnp, array([1.84118, 2304 5.33144, 2305 8.53632, 2306 11.70600, 2307 14.86359]),4) 2308 jnp = special.jnp_zeros(443,5) 2309 assert_allclose(special.jvp(443, jnp), 0, atol=1e-15) 2310 2311 def test_jnyn_zeros(self): 2312 jnz = special.jnyn_zeros(1,5) 2313 assert_array_almost_equal(jnz,(array([3.83171, 2314 7.01559, 2315 10.17347, 2316 13.32369, 2317 16.47063]), 2318 array([1.84118, 2319 5.33144, 2320 8.53632, 2321 11.70600, 2322 14.86359]), 2323 array([2.19714, 2324 5.42968, 2325 8.59601, 2326 11.74915, 2327 14.89744]), 2328 array([3.68302, 2329 6.94150, 2330 10.12340, 2331 13.28576, 2332 16.44006])),5) 2333 2334 def test_jvp(self): 2335 jvprim = special.jvp(2,2) 2336 jv0 = (special.jv(1,2)-special.jv(3,2))/2 2337 assert_almost_equal(jvprim,jv0,10) 2338 2339 def test_k0(self): 2340 ozk = special.k0(.1) 2341 ozkr = special.kv(0,.1) 2342 assert_almost_equal(ozk,ozkr,8) 2343 2344 def test_k0e(self): 2345 ozke = special.k0e(.1) 2346 ozker = special.kve(0,.1) 2347 assert_almost_equal(ozke,ozker,8) 2348 2349 def test_k1(self): 2350 o1k = special.k1(.1) 2351 o1kr = special.kv(1,.1) 2352 assert_almost_equal(o1k,o1kr,8) 2353 2354 def test_k1e(self): 2355 o1ke = special.k1e(.1) 2356 o1ker = special.kve(1,.1) 2357 assert_almost_equal(o1ke,o1ker,8) 2358 2359 def test_jacobi(self): 2360 a = 5*np.random.random() - 1 2361 b = 5*np.random.random() - 1 2362 P0 = special.jacobi(0,a,b) 2363 P1 = special.jacobi(1,a,b) 2364 P2 = special.jacobi(2,a,b) 2365 P3 = special.jacobi(3,a,b) 2366 2367 assert_array_almost_equal(P0.c,[1],13) 2368 assert_array_almost_equal(P1.c,array([a+b+2,a-b])/2.0,13) 2369 cp = [(a+b+3)*(a+b+4), 4*(a+b+3)*(a+2), 4*(a+1)*(a+2)] 2370 p2c = [cp[0],cp[1]-2*cp[0],cp[2]-cp[1]+cp[0]] 2371 assert_array_almost_equal(P2.c,array(p2c)/8.0,13) 2372 cp = [(a+b+4)*(a+b+5)*(a+b+6),6*(a+b+4)*(a+b+5)*(a+3), 2373 12*(a+b+4)*(a+2)*(a+3),8*(a+1)*(a+2)*(a+3)] 2374 p3c = [cp[0],cp[1]-3*cp[0],cp[2]-2*cp[1]+3*cp[0],cp[3]-cp[2]+cp[1]-cp[0]] 2375 assert_array_almost_equal(P3.c,array(p3c)/48.0,13) 2376 2377 def test_kn(self): 2378 kn1 = special.kn(0,.2) 2379 assert_almost_equal(kn1,1.7527038555281462,8) 2380 2381 def test_negv_kv(self): 2382 assert_equal(special.kv(3.0, 2.2), special.kv(-3.0, 2.2)) 2383 2384 def test_kv0(self): 2385 kv0 = special.kv(0,.2) 2386 assert_almost_equal(kv0, 1.7527038555281462, 10) 2387 2388 def test_kv1(self): 2389 kv1 = special.kv(1,0.2) 2390 assert_almost_equal(kv1, 4.775972543220472, 10) 2391 2392 def test_kv2(self): 2393 kv2 = special.kv(2,0.2) 2394 assert_almost_equal(kv2, 49.51242928773287, 10) 2395 2396 def test_kn_largeorder(self): 2397 assert_allclose(special.kn(32, 1), 1.7516596664574289e+43) 2398 2399 def test_kv_largearg(self): 2400 assert_equal(special.kv(0, 1e19), 0) 2401 2402 def test_negv_kve(self): 2403 assert_equal(special.kve(3.0, 2.2), special.kve(-3.0, 2.2)) 2404 2405 def test_kve(self): 2406 kve1 = special.kve(0,.2) 2407 kv1 = special.kv(0,.2)*exp(.2) 2408 assert_almost_equal(kve1,kv1,8) 2409 z = .2+1j 2410 kve2 = special.kve(0,z) 2411 kv2 = special.kv(0,z)*exp(z) 2412 assert_almost_equal(kve2,kv2,8) 2413 2414 def test_kvp_v0n1(self): 2415 z = 2.2 2416 assert_almost_equal(-special.kv(1,z), special.kvp(0,z, n=1), 10) 2417 2418 def test_kvp_n1(self): 2419 v = 3. 2420 z = 2.2 2421 xc = -special.kv(v+1,z) + v/z*special.kv(v,z) 2422 x = special.kvp(v,z, n=1) 2423 assert_almost_equal(xc, x, 10) # this function (kvp) is broken 2424 2425 def test_kvp_n2(self): 2426 v = 3. 2427 z = 2.2 2428 xc = (z**2+v**2-v)/z**2 * special.kv(v,z) + special.kv(v+1,z)/z 2429 x = special.kvp(v, z, n=2) 2430 assert_almost_equal(xc, x, 10) 2431 2432 def test_y0(self): 2433 oz = special.y0(.1) 2434 ozr = special.yn(0,.1) 2435 assert_almost_equal(oz,ozr,8) 2436 2437 def test_y1(self): 2438 o1 = special.y1(.1) 2439 o1r = special.yn(1,.1) 2440 assert_almost_equal(o1,o1r,8) 2441 2442 def test_y0_zeros(self): 2443 yo,ypo = special.y0_zeros(2) 2444 zo,zpo = special.y0_zeros(2,complex=1) 2445 all = r_[yo,zo] 2446 allval = r_[ypo,zpo] 2447 assert_array_almost_equal(abs(special.yv(0.0,all)),0.0,11) 2448 assert_array_almost_equal(abs(special.yv(1,all)-allval),0.0,11) 2449 2450 def test_y1_zeros(self): 2451 y1 = special.y1_zeros(1) 2452 assert_array_almost_equal(y1,(array([2.19714]),array([0.52079])),5) 2453 2454 def test_y1p_zeros(self): 2455 y1p = special.y1p_zeros(1,complex=1) 2456 assert_array_almost_equal(y1p,(array([0.5768+0.904j]), array([-0.7635+0.5892j])),3) 2457 2458 def test_yn_zeros(self): 2459 an = special.yn_zeros(4,2) 2460 assert_array_almost_equal(an,array([5.64515, 9.36162]),5) 2461 an = special.yn_zeros(443,5) 2462 assert_allclose(an, [450.13573091578090314, 463.05692376675001542, 2463 472.80651546418663566, 481.27353184725625838, 2464 488.98055964441374646], rtol=1e-15) 2465 2466 def test_ynp_zeros(self): 2467 ao = special.ynp_zeros(0,2) 2468 assert_array_almost_equal(ao,array([2.19714133, 5.42968104]),6) 2469 ao = special.ynp_zeros(43,5) 2470 assert_allclose(special.yvp(43, ao), 0, atol=1e-15) 2471 ao = special.ynp_zeros(443,5) 2472 assert_allclose(special.yvp(443, ao), 0, atol=1e-9) 2473 2474 def test_ynp_zeros_large_order(self): 2475 ao = special.ynp_zeros(443,5) 2476 assert_allclose(special.yvp(443, ao), 0, atol=1e-14) 2477 2478 def test_yn(self): 2479 yn2n = special.yn(1,.2) 2480 assert_almost_equal(yn2n,-3.3238249881118471,8) 2481 2482 def test_negv_yv(self): 2483 assert_almost_equal(special.yv(-3,2), -special.yv(3,2), 14) 2484 2485 def test_yv(self): 2486 yv2 = special.yv(1,.2) 2487 assert_almost_equal(yv2,-3.3238249881118471,8) 2488 2489 def test_negv_yve(self): 2490 assert_almost_equal(special.yve(-3,2), -special.yve(3,2), 14) 2491 2492 def test_yve(self): 2493 yve2 = special.yve(1,.2) 2494 assert_almost_equal(yve2,-3.3238249881118471,8) 2495 yve2r = special.yv(1,.2+1j)*exp(-1) 2496 yve22 = special.yve(1,.2+1j) 2497 assert_almost_equal(yve22,yve2r,8) 2498 2499 def test_yvp(self): 2500 yvpr = (special.yv(1,.2) - special.yv(3,.2))/2.0 2501 yvp1 = special.yvp(2,.2) 2502 assert_array_almost_equal(yvp1,yvpr,10) 2503 2504 def _cephes_vs_amos_points(self): 2505 """Yield points at which to compare Cephes implementation to AMOS""" 2506 # check several points, including large-amplitude ones 2507 v = [-120, -100.3, -20., -10., -1., -.5, 0., 1., 12.49, 120., 301] 2508 z = [-1300, -11, -10, -1, 1., 10., 200.5, 401., 600.5, 700.6, 1300, 2509 10003] 2510 yield from itertools.product(v, z) 2511 2512 # check half-integers; these are problematic points at least 2513 # for cephes/iv 2514 yield from itertools.product(0.5 + arange(-60, 60), [3.5]) 2515 2516 def check_cephes_vs_amos(self, f1, f2, rtol=1e-11, atol=0, skip=None): 2517 for v, z in self._cephes_vs_amos_points(): 2518 if skip is not None and skip(v, z): 2519 continue 2520 c1, c2, c3 = f1(v, z), f1(v,z+0j), f2(int(v), z) 2521 if np.isinf(c1): 2522 assert_(np.abs(c2) >= 1e300, (v, z)) 2523 elif np.isnan(c1): 2524 assert_(c2.imag != 0, (v, z)) 2525 else: 2526 assert_allclose(c1, c2, err_msg=(v, z), rtol=rtol, atol=atol) 2527 if v == int(v): 2528 assert_allclose(c3, c2, err_msg=(v, z), 2529 rtol=rtol, atol=atol) 2530 2531 @pytest.mark.xfail(platform.machine() == 'ppc64le', 2532 reason="fails on ppc64le") 2533 def test_jv_cephes_vs_amos(self): 2534 self.check_cephes_vs_amos(special.jv, special.jn, rtol=1e-10, atol=1e-305) 2535 2536 @pytest.mark.xfail(platform.machine() == 'ppc64le', 2537 reason="fails on ppc64le") 2538 def test_yv_cephes_vs_amos(self): 2539 self.check_cephes_vs_amos(special.yv, special.yn, rtol=1e-11, atol=1e-305) 2540 2541 def test_yv_cephes_vs_amos_only_small_orders(self): 2542 skipper = lambda v, z: (abs(v) > 50) 2543 self.check_cephes_vs_amos(special.yv, special.yn, rtol=1e-11, atol=1e-305, skip=skipper) 2544 2545 def test_iv_cephes_vs_amos(self): 2546 with np.errstate(all='ignore'): 2547 self.check_cephes_vs_amos(special.iv, special.iv, rtol=5e-9, atol=1e-305) 2548 2549 @pytest.mark.slow 2550 def test_iv_cephes_vs_amos_mass_test(self): 2551 N = 1000000 2552 np.random.seed(1) 2553 v = np.random.pareto(0.5, N) * (-1)**np.random.randint(2, size=N) 2554 x = np.random.pareto(0.2, N) * (-1)**np.random.randint(2, size=N) 2555 2556 imsk = (np.random.randint(8, size=N) == 0) 2557 v[imsk] = v[imsk].astype(int) 2558 2559 with np.errstate(all='ignore'): 2560 c1 = special.iv(v, x) 2561 c2 = special.iv(v, x+0j) 2562 2563 # deal with differences in the inf and zero cutoffs 2564 c1[abs(c1) > 1e300] = np.inf 2565 c2[abs(c2) > 1e300] = np.inf 2566 c1[abs(c1) < 1e-300] = 0 2567 c2[abs(c2) < 1e-300] = 0 2568 2569 dc = abs(c1/c2 - 1) 2570 dc[np.isnan(dc)] = 0 2571 2572 k = np.argmax(dc) 2573 2574 # Most error apparently comes from AMOS and not our implementation; 2575 # there are some problems near integer orders there 2576 assert_(dc[k] < 2e-7, (v[k], x[k], special.iv(v[k], x[k]), special.iv(v[k], x[k]+0j))) 2577 2578 def test_kv_cephes_vs_amos(self): 2579 self.check_cephes_vs_amos(special.kv, special.kn, rtol=1e-9, atol=1e-305) 2580 self.check_cephes_vs_amos(special.kv, special.kv, rtol=1e-9, atol=1e-305) 2581 2582 def test_ticket_623(self): 2583 assert_allclose(special.jv(3, 4), 0.43017147387562193) 2584 assert_allclose(special.jv(301, 1300), 0.0183487151115275) 2585 assert_allclose(special.jv(301, 1296.0682), -0.0224174325312048) 2586 2587 def test_ticket_853(self): 2588 """Negative-order Bessels""" 2589 # cephes 2590 assert_allclose(special.jv(-1, 1), -0.4400505857449335) 2591 assert_allclose(special.jv(-2, 1), 0.1149034849319005) 2592 assert_allclose(special.yv(-1, 1), 0.7812128213002887) 2593 assert_allclose(special.yv(-2, 1), -1.650682606816255) 2594 assert_allclose(special.iv(-1, 1), 0.5651591039924851) 2595 assert_allclose(special.iv(-2, 1), 0.1357476697670383) 2596 assert_allclose(special.kv(-1, 1), 0.6019072301972347) 2597 assert_allclose(special.kv(-2, 1), 1.624838898635178) 2598 assert_allclose(special.jv(-0.5, 1), 0.43109886801837607952) 2599 assert_allclose(special.yv(-0.5, 1), 0.6713967071418031) 2600 assert_allclose(special.iv(-0.5, 1), 1.231200214592967) 2601 assert_allclose(special.kv(-0.5, 1), 0.4610685044478945) 2602 # amos 2603 assert_allclose(special.jv(-1, 1+0j), -0.4400505857449335) 2604 assert_allclose(special.jv(-2, 1+0j), 0.1149034849319005) 2605 assert_allclose(special.yv(-1, 1+0j), 0.7812128213002887) 2606 assert_allclose(special.yv(-2, 1+0j), -1.650682606816255) 2607 2608 assert_allclose(special.iv(-1, 1+0j), 0.5651591039924851) 2609 assert_allclose(special.iv(-2, 1+0j), 0.1357476697670383) 2610 assert_allclose(special.kv(-1, 1+0j), 0.6019072301972347) 2611 assert_allclose(special.kv(-2, 1+0j), 1.624838898635178) 2612 2613 assert_allclose(special.jv(-0.5, 1+0j), 0.43109886801837607952) 2614 assert_allclose(special.jv(-0.5, 1+1j), 0.2628946385649065-0.827050182040562j) 2615 assert_allclose(special.yv(-0.5, 1+0j), 0.6713967071418031) 2616 assert_allclose(special.yv(-0.5, 1+1j), 0.967901282890131+0.0602046062142816j) 2617 2618 assert_allclose(special.iv(-0.5, 1+0j), 1.231200214592967) 2619 assert_allclose(special.iv(-0.5, 1+1j), 0.77070737376928+0.39891821043561j) 2620 assert_allclose(special.kv(-0.5, 1+0j), 0.4610685044478945) 2621 assert_allclose(special.kv(-0.5, 1+1j), 0.06868578341999-0.38157825981268j) 2622 2623 assert_allclose(special.jve(-0.5,1+0.3j), special.jv(-0.5, 1+0.3j)*exp(-0.3)) 2624 assert_allclose(special.yve(-0.5,1+0.3j), special.yv(-0.5, 1+0.3j)*exp(-0.3)) 2625 assert_allclose(special.ive(-0.5,0.3+1j), special.iv(-0.5, 0.3+1j)*exp(-0.3)) 2626 assert_allclose(special.kve(-0.5,0.3+1j), special.kv(-0.5, 0.3+1j)*exp(0.3+1j)) 2627 2628 assert_allclose(special.hankel1(-0.5, 1+1j), special.jv(-0.5, 1+1j) + 1j*special.yv(-0.5,1+1j)) 2629 assert_allclose(special.hankel2(-0.5, 1+1j), special.jv(-0.5, 1+1j) - 1j*special.yv(-0.5,1+1j)) 2630 2631 def test_ticket_854(self): 2632 """Real-valued Bessel domains""" 2633 assert_(isnan(special.jv(0.5, -1))) 2634 assert_(isnan(special.iv(0.5, -1))) 2635 assert_(isnan(special.yv(0.5, -1))) 2636 assert_(isnan(special.yv(1, -1))) 2637 assert_(isnan(special.kv(0.5, -1))) 2638 assert_(isnan(special.kv(1, -1))) 2639 assert_(isnan(special.jve(0.5, -1))) 2640 assert_(isnan(special.ive(0.5, -1))) 2641 assert_(isnan(special.yve(0.5, -1))) 2642 assert_(isnan(special.yve(1, -1))) 2643 assert_(isnan(special.kve(0.5, -1))) 2644 assert_(isnan(special.kve(1, -1))) 2645 assert_(isnan(special.airye(-1)[0:2]).all(), special.airye(-1)) 2646 assert_(not isnan(special.airye(-1)[2:4]).any(), special.airye(-1)) 2647 2648 def test_gh_7909(self): 2649 assert_(special.kv(1.5, 0) == np.inf) 2650 assert_(special.kve(1.5, 0) == np.inf) 2651 2652 def test_ticket_503(self): 2653 """Real-valued Bessel I overflow""" 2654 assert_allclose(special.iv(1, 700), 1.528500390233901e302) 2655 assert_allclose(special.iv(1000, 1120), 1.301564549405821e301) 2656 2657 def test_iv_hyperg_poles(self): 2658 assert_allclose(special.iv(-0.5, 1), 1.231200214592967) 2659 2660 def iv_series(self, v, z, n=200): 2661 k = arange(0, n).astype(float_) 2662 r = (v+2*k)*log(.5*z) - special.gammaln(k+1) - special.gammaln(v+k+1) 2663 r[isnan(r)] = inf 2664 r = exp(r) 2665 err = abs(r).max() * finfo(float_).eps * n + abs(r[-1])*10 2666 return r.sum(), err 2667 2668 def test_i0_series(self): 2669 for z in [1., 10., 200.5]: 2670 value, err = self.iv_series(0, z) 2671 assert_allclose(special.i0(z), value, atol=err, err_msg=z) 2672 2673 def test_i1_series(self): 2674 for z in [1., 10., 200.5]: 2675 value, err = self.iv_series(1, z) 2676 assert_allclose(special.i1(z), value, atol=err, err_msg=z) 2677 2678 def test_iv_series(self): 2679 for v in [-20., -10., -1., 0., 1., 12.49, 120.]: 2680 for z in [1., 10., 200.5, -1+2j]: 2681 value, err = self.iv_series(v, z) 2682 assert_allclose(special.iv(v, z), value, atol=err, err_msg=(v, z)) 2683 2684 def test_i0(self): 2685 values = [[0.0, 1.0], 2686 [1e-10, 1.0], 2687 [0.1, 0.9071009258], 2688 [0.5, 0.6450352706], 2689 [1.0, 0.4657596077], 2690 [2.5, 0.2700464416], 2691 [5.0, 0.1835408126], 2692 [20.0, 0.0897803119], 2693 ] 2694 for i, (x, v) in enumerate(values): 2695 cv = special.i0(x) * exp(-x) 2696 assert_almost_equal(cv, v, 8, err_msg='test #%d' % i) 2697 2698 def test_i0e(self): 2699 oize = special.i0e(.1) 2700 oizer = special.ive(0,.1) 2701 assert_almost_equal(oize,oizer,8) 2702 2703 def test_i1(self): 2704 values = [[0.0, 0.0], 2705 [1e-10, 0.4999999999500000e-10], 2706 [0.1, 0.0452984468], 2707 [0.5, 0.1564208032], 2708 [1.0, 0.2079104154], 2709 [5.0, 0.1639722669], 2710 [20.0, 0.0875062222], 2711 ] 2712 for i, (x, v) in enumerate(values): 2713 cv = special.i1(x) * exp(-x) 2714 assert_almost_equal(cv, v, 8, err_msg='test #%d' % i) 2715 2716 def test_i1e(self): 2717 oi1e = special.i1e(.1) 2718 oi1er = special.ive(1,.1) 2719 assert_almost_equal(oi1e,oi1er,8) 2720 2721 def test_iti0k0(self): 2722 iti0 = array(special.iti0k0(5)) 2723 assert_array_almost_equal(iti0,array([31.848667776169801, 1.5673873907283657]),5) 2724 2725 def test_it2i0k0(self): 2726 it2k = special.it2i0k0(.1) 2727 assert_array_almost_equal(it2k,array([0.0012503906973464409, 3.3309450354686687]),6) 2728 2729 def test_iv(self): 2730 iv1 = special.iv(0,.1)*exp(-.1) 2731 assert_almost_equal(iv1,0.90710092578230106,10) 2732 2733 def test_negv_ive(self): 2734 assert_equal(special.ive(3,2), special.ive(-3,2)) 2735 2736 def test_ive(self): 2737 ive1 = special.ive(0,.1) 2738 iv1 = special.iv(0,.1)*exp(-.1) 2739 assert_almost_equal(ive1,iv1,10) 2740 2741 def test_ivp0(self): 2742 assert_almost_equal(special.iv(1,2), special.ivp(0,2), 10) 2743 2744 def test_ivp(self): 2745 y = (special.iv(0,2) + special.iv(2,2))/2 2746 x = special.ivp(1,2) 2747 assert_almost_equal(x,y,10) 2748 2749 2750class TestLaguerre: 2751 def test_laguerre(self): 2752 lag0 = special.laguerre(0) 2753 lag1 = special.laguerre(1) 2754 lag2 = special.laguerre(2) 2755 lag3 = special.laguerre(3) 2756 lag4 = special.laguerre(4) 2757 lag5 = special.laguerre(5) 2758 assert_array_almost_equal(lag0.c,[1],13) 2759 assert_array_almost_equal(lag1.c,[-1,1],13) 2760 assert_array_almost_equal(lag2.c,array([1,-4,2])/2.0,13) 2761 assert_array_almost_equal(lag3.c,array([-1,9,-18,6])/6.0,13) 2762 assert_array_almost_equal(lag4.c,array([1,-16,72,-96,24])/24.0,13) 2763 assert_array_almost_equal(lag5.c,array([-1,25,-200,600,-600,120])/120.0,13) 2764 2765 def test_genlaguerre(self): 2766 k = 5*np.random.random() - 0.9 2767 lag0 = special.genlaguerre(0,k) 2768 lag1 = special.genlaguerre(1,k) 2769 lag2 = special.genlaguerre(2,k) 2770 lag3 = special.genlaguerre(3,k) 2771 assert_equal(lag0.c,[1]) 2772 assert_equal(lag1.c,[-1,k+1]) 2773 assert_almost_equal(lag2.c,array([1,-2*(k+2),(k+1.)*(k+2.)])/2.0) 2774 assert_almost_equal(lag3.c,array([-1,3*(k+3),-3*(k+2)*(k+3),(k+1)*(k+2)*(k+3)])/6.0) 2775 2776 2777# Base polynomials come from Abrahmowitz and Stegan 2778class TestLegendre: 2779 def test_legendre(self): 2780 leg0 = special.legendre(0) 2781 leg1 = special.legendre(1) 2782 leg2 = special.legendre(2) 2783 leg3 = special.legendre(3) 2784 leg4 = special.legendre(4) 2785 leg5 = special.legendre(5) 2786 assert_equal(leg0.c, [1]) 2787 assert_equal(leg1.c, [1,0]) 2788 assert_almost_equal(leg2.c, array([3,0,-1])/2.0, decimal=13) 2789 assert_almost_equal(leg3.c, array([5,0,-3,0])/2.0) 2790 assert_almost_equal(leg4.c, array([35,0,-30,0,3])/8.0) 2791 assert_almost_equal(leg5.c, array([63,0,-70,0,15,0])/8.0) 2792 2793 2794class TestLambda: 2795 def test_lmbda(self): 2796 lam = special.lmbda(1,.1) 2797 lamr = (array([special.jn(0,.1), 2*special.jn(1,.1)/.1]), 2798 array([special.jvp(0,.1), -2*special.jv(1,.1)/.01 + 2*special.jvp(1,.1)/.1])) 2799 assert_array_almost_equal(lam,lamr,8) 2800 2801 2802class TestLog1p: 2803 def test_log1p(self): 2804 l1p = (special.log1p(10), special.log1p(11), special.log1p(12)) 2805 l1prl = (log(11), log(12), log(13)) 2806 assert_array_almost_equal(l1p,l1prl,8) 2807 2808 def test_log1pmore(self): 2809 l1pm = (special.log1p(1), special.log1p(1.1), special.log1p(1.2)) 2810 l1pmrl = (log(2),log(2.1),log(2.2)) 2811 assert_array_almost_equal(l1pm,l1pmrl,8) 2812 2813 2814class TestLegendreFunctions: 2815 def test_clpmn(self): 2816 z = 0.5+0.3j 2817 clp = special.clpmn(2, 2, z, 3) 2818 assert_array_almost_equal(clp, 2819 (array([[1.0000, z, 0.5*(3*z*z-1)], 2820 [0.0000, sqrt(z*z-1), 3*z*sqrt(z*z-1)], 2821 [0.0000, 0.0000, 3*(z*z-1)]]), 2822 array([[0.0000, 1.0000, 3*z], 2823 [0.0000, z/sqrt(z*z-1), 3*(2*z*z-1)/sqrt(z*z-1)], 2824 [0.0000, 0.0000, 6*z]])), 2825 7) 2826 2827 def test_clpmn_close_to_real_2(self): 2828 eps = 1e-10 2829 m = 1 2830 n = 3 2831 x = 0.5 2832 clp_plus = special.clpmn(m, n, x+1j*eps, 2)[0][m, n] 2833 clp_minus = special.clpmn(m, n, x-1j*eps, 2)[0][m, n] 2834 assert_array_almost_equal(array([clp_plus, clp_minus]), 2835 array([special.lpmv(m, n, x), 2836 special.lpmv(m, n, x)]), 2837 7) 2838 2839 def test_clpmn_close_to_real_3(self): 2840 eps = 1e-10 2841 m = 1 2842 n = 3 2843 x = 0.5 2844 clp_plus = special.clpmn(m, n, x+1j*eps, 3)[0][m, n] 2845 clp_minus = special.clpmn(m, n, x-1j*eps, 3)[0][m, n] 2846 assert_array_almost_equal(array([clp_plus, clp_minus]), 2847 array([special.lpmv(m, n, x)*np.exp(-0.5j*m*np.pi), 2848 special.lpmv(m, n, x)*np.exp(0.5j*m*np.pi)]), 2849 7) 2850 2851 def test_clpmn_across_unit_circle(self): 2852 eps = 1e-7 2853 m = 1 2854 n = 1 2855 x = 1j 2856 for type in [2, 3]: 2857 assert_almost_equal(special.clpmn(m, n, x+1j*eps, type)[0][m, n], 2858 special.clpmn(m, n, x-1j*eps, type)[0][m, n], 6) 2859 2860 def test_inf(self): 2861 for z in (1, -1): 2862 for n in range(4): 2863 for m in range(1, n): 2864 lp = special.clpmn(m, n, z) 2865 assert_(np.isinf(lp[1][1,1:]).all()) 2866 lp = special.lpmn(m, n, z) 2867 assert_(np.isinf(lp[1][1,1:]).all()) 2868 2869 def test_deriv_clpmn(self): 2870 # data inside and outside of the unit circle 2871 zvals = [0.5+0.5j, -0.5+0.5j, -0.5-0.5j, 0.5-0.5j, 2872 1+1j, -1+1j, -1-1j, 1-1j] 2873 m = 2 2874 n = 3 2875 for type in [2, 3]: 2876 for z in zvals: 2877 for h in [1e-3, 1e-3j]: 2878 approx_derivative = (special.clpmn(m, n, z+0.5*h, type)[0] 2879 - special.clpmn(m, n, z-0.5*h, type)[0])/h 2880 assert_allclose(special.clpmn(m, n, z, type)[1], 2881 approx_derivative, 2882 rtol=1e-4) 2883 2884 def test_lpmn(self): 2885 lp = special.lpmn(0,2,.5) 2886 assert_array_almost_equal(lp,(array([[1.00000, 2887 0.50000, 2888 -0.12500]]), 2889 array([[0.00000, 2890 1.00000, 2891 1.50000]])),4) 2892 2893 def test_lpn(self): 2894 lpnf = special.lpn(2,.5) 2895 assert_array_almost_equal(lpnf,(array([1.00000, 2896 0.50000, 2897 -0.12500]), 2898 array([0.00000, 2899 1.00000, 2900 1.50000])),4) 2901 2902 def test_lpmv(self): 2903 lp = special.lpmv(0,2,.5) 2904 assert_almost_equal(lp,-0.125,7) 2905 lp = special.lpmv(0,40,.001) 2906 assert_almost_equal(lp,0.1252678976534484,7) 2907 2908 # XXX: this is outside the domain of the current implementation, 2909 # so ensure it returns a NaN rather than a wrong answer. 2910 with np.errstate(all='ignore'): 2911 lp = special.lpmv(-1,-1,.001) 2912 assert_(lp != 0 or np.isnan(lp)) 2913 2914 def test_lqmn(self): 2915 lqmnf = special.lqmn(0,2,.5) 2916 lqf = special.lqn(2,.5) 2917 assert_array_almost_equal(lqmnf[0][0],lqf[0],4) 2918 assert_array_almost_equal(lqmnf[1][0],lqf[1],4) 2919 2920 def test_lqmn_gt1(self): 2921 """algorithm for real arguments changes at 1.0001 2922 test against analytical result for m=2, n=1 2923 """ 2924 x0 = 1.0001 2925 delta = 0.00002 2926 for x in (x0-delta, x0+delta): 2927 lq = special.lqmn(2, 1, x)[0][-1, -1] 2928 expected = 2/(x*x-1) 2929 assert_almost_equal(lq, expected) 2930 2931 def test_lqmn_shape(self): 2932 a, b = special.lqmn(4, 4, 1.1) 2933 assert_equal(a.shape, (5, 5)) 2934 assert_equal(b.shape, (5, 5)) 2935 2936 a, b = special.lqmn(4, 0, 1.1) 2937 assert_equal(a.shape, (5, 1)) 2938 assert_equal(b.shape, (5, 1)) 2939 2940 def test_lqn(self): 2941 lqf = special.lqn(2,.5) 2942 assert_array_almost_equal(lqf,(array([0.5493, -0.7253, -0.8187]), 2943 array([1.3333, 1.216, -0.8427])),4) 2944 2945 2946class TestMathieu: 2947 2948 def test_mathieu_a(self): 2949 pass 2950 2951 def test_mathieu_even_coef(self): 2952 special.mathieu_even_coef(2,5) 2953 # Q not defined broken and cannot figure out proper reporting order 2954 2955 def test_mathieu_odd_coef(self): 2956 # same problem as above 2957 pass 2958 2959 2960class TestFresnelIntegral: 2961 2962 def test_modfresnelp(self): 2963 pass 2964 2965 def test_modfresnelm(self): 2966 pass 2967 2968 2969class TestOblCvSeq: 2970 def test_obl_cv_seq(self): 2971 obl = special.obl_cv_seq(0,3,1) 2972 assert_array_almost_equal(obl,array([-0.348602, 2973 1.393206, 2974 5.486800, 2975 11.492120]),5) 2976 2977 2978class TestParabolicCylinder: 2979 def test_pbdn_seq(self): 2980 pb = special.pbdn_seq(1,.1) 2981 assert_array_almost_equal(pb,(array([0.9975, 2982 0.0998]), 2983 array([-0.0499, 2984 0.9925])),4) 2985 2986 def test_pbdv(self): 2987 special.pbdv(1,.2) 2988 1/2*(.2)*special.pbdv(1,.2)[0] - special.pbdv(0,.2)[0] 2989 2990 def test_pbdv_seq(self): 2991 pbn = special.pbdn_seq(1,.1) 2992 pbv = special.pbdv_seq(1,.1) 2993 assert_array_almost_equal(pbv,(real(pbn[0]),real(pbn[1])),4) 2994 2995 def test_pbdv_points(self): 2996 # simple case 2997 eta = np.linspace(-10, 10, 5) 2998 z = 2**(eta/2)*np.sqrt(np.pi)/special.gamma(.5-.5*eta) 2999 assert_allclose(special.pbdv(eta, 0.)[0], z, rtol=1e-14, atol=1e-14) 3000 3001 # some points 3002 assert_allclose(special.pbdv(10.34, 20.44)[0], 1.3731383034455e-32, rtol=1e-12) 3003 assert_allclose(special.pbdv(-9.53, 3.44)[0], 3.166735001119246e-8, rtol=1e-12) 3004 3005 def test_pbdv_gradient(self): 3006 x = np.linspace(-4, 4, 8)[:,None] 3007 eta = np.linspace(-10, 10, 5)[None,:] 3008 3009 p = special.pbdv(eta, x) 3010 eps = 1e-7 + 1e-7*abs(x) 3011 dp = (special.pbdv(eta, x + eps)[0] - special.pbdv(eta, x - eps)[0]) / eps / 2. 3012 assert_allclose(p[1], dp, rtol=1e-6, atol=1e-6) 3013 3014 def test_pbvv_gradient(self): 3015 x = np.linspace(-4, 4, 8)[:,None] 3016 eta = np.linspace(-10, 10, 5)[None,:] 3017 3018 p = special.pbvv(eta, x) 3019 eps = 1e-7 + 1e-7*abs(x) 3020 dp = (special.pbvv(eta, x + eps)[0] - special.pbvv(eta, x - eps)[0]) / eps / 2. 3021 assert_allclose(p[1], dp, rtol=1e-6, atol=1e-6) 3022 3023 3024class TestPolygamma: 3025 # from Table 6.2 (pg. 271) of A&S 3026 def test_polygamma(self): 3027 poly2 = special.polygamma(2,1) 3028 poly3 = special.polygamma(3,1) 3029 assert_almost_equal(poly2,-2.4041138063,10) 3030 assert_almost_equal(poly3,6.4939394023,10) 3031 3032 # Test polygamma(0, x) == psi(x) 3033 x = [2, 3, 1.1e14] 3034 assert_almost_equal(special.polygamma(0, x), special.psi(x)) 3035 3036 # Test broadcasting 3037 n = [0, 1, 2] 3038 x = [0.5, 1.5, 2.5] 3039 expected = [-1.9635100260214238, 0.93480220054467933, 3040 -0.23620405164172739] 3041 assert_almost_equal(special.polygamma(n, x), expected) 3042 expected = np.row_stack([expected]*2) 3043 assert_almost_equal(special.polygamma(n, np.row_stack([x]*2)), 3044 expected) 3045 assert_almost_equal(special.polygamma(np.row_stack([n]*2), x), 3046 expected) 3047 3048 3049class TestProCvSeq: 3050 def test_pro_cv_seq(self): 3051 prol = special.pro_cv_seq(0,3,1) 3052 assert_array_almost_equal(prol,array([0.319000, 3053 2.593084, 3054 6.533471, 3055 12.514462]),5) 3056 3057 3058class TestPsi: 3059 def test_psi(self): 3060 ps = special.psi(1) 3061 assert_almost_equal(ps,-0.57721566490153287,8) 3062 3063 3064class TestRadian: 3065 def test_radian(self): 3066 rad = special.radian(90,0,0) 3067 assert_almost_equal(rad,pi/2.0,5) 3068 3069 def test_radianmore(self): 3070 rad1 = special.radian(90,1,60) 3071 assert_almost_equal(rad1,pi/2+0.0005816135199345904,5) 3072 3073 3074class TestRiccati: 3075 def test_riccati_jn(self): 3076 N, x = 2, 0.2 3077 S = np.empty((N, N)) 3078 for n in range(N): 3079 j = special.spherical_jn(n, x) 3080 jp = special.spherical_jn(n, x, derivative=True) 3081 S[0,n] = x*j 3082 S[1,n] = x*jp + j 3083 assert_array_almost_equal(S, special.riccati_jn(n, x), 8) 3084 3085 def test_riccati_yn(self): 3086 N, x = 2, 0.2 3087 C = np.empty((N, N)) 3088 for n in range(N): 3089 y = special.spherical_yn(n, x) 3090 yp = special.spherical_yn(n, x, derivative=True) 3091 C[0,n] = x*y 3092 C[1,n] = x*yp + y 3093 assert_array_almost_equal(C, special.riccati_yn(n, x), 8) 3094 3095 3096class TestRound: 3097 def test_round(self): 3098 rnd = list(map(int,(special.round(10.1),special.round(10.4),special.round(10.5),special.round(10.6)))) 3099 3100 # Note: According to the documentation, scipy.special.round is 3101 # supposed to round to the nearest even number if the fractional 3102 # part is exactly 0.5. On some platforms, this does not appear 3103 # to work and thus this test may fail. However, this unit test is 3104 # correctly written. 3105 rndrl = (10,10,10,11) 3106 assert_array_equal(rnd,rndrl) 3107 3108 3109def test_sph_harm(): 3110 # Tests derived from tables in 3111 # https://en.wikipedia.org/wiki/Table_of_spherical_harmonics 3112 sh = special.sph_harm 3113 pi = np.pi 3114 exp = np.exp 3115 sqrt = np.sqrt 3116 sin = np.sin 3117 cos = np.cos 3118 assert_array_almost_equal(sh(0,0,0,0), 3119 0.5/sqrt(pi)) 3120 assert_array_almost_equal(sh(-2,2,0.,pi/4), 3121 0.25*sqrt(15./(2.*pi)) * 3122 (sin(pi/4))**2.) 3123 assert_array_almost_equal(sh(-2,2,0.,pi/2), 3124 0.25*sqrt(15./(2.*pi))) 3125 assert_array_almost_equal(sh(2,2,pi,pi/2), 3126 0.25*sqrt(15/(2.*pi)) * 3127 exp(0+2.*pi*1j)*sin(pi/2.)**2.) 3128 assert_array_almost_equal(sh(2,4,pi/4.,pi/3.), 3129 (3./8.)*sqrt(5./(2.*pi)) * 3130 exp(0+2.*pi/4.*1j) * 3131 sin(pi/3.)**2. * 3132 (7.*cos(pi/3.)**2.-1)) 3133 assert_array_almost_equal(sh(4,4,pi/8.,pi/6.), 3134 (3./16.)*sqrt(35./(2.*pi)) * 3135 exp(0+4.*pi/8.*1j)*sin(pi/6.)**4.) 3136 3137 3138def test_sph_harm_ufunc_loop_selection(): 3139 # see https://github.com/scipy/scipy/issues/4895 3140 dt = np.dtype(np.complex128) 3141 assert_equal(special.sph_harm(0, 0, 0, 0).dtype, dt) 3142 assert_equal(special.sph_harm([0], 0, 0, 0).dtype, dt) 3143 assert_equal(special.sph_harm(0, [0], 0, 0).dtype, dt) 3144 assert_equal(special.sph_harm(0, 0, [0], 0).dtype, dt) 3145 assert_equal(special.sph_harm(0, 0, 0, [0]).dtype, dt) 3146 assert_equal(special.sph_harm([0], [0], [0], [0]).dtype, dt) 3147 3148 3149class TestStruve: 3150 def _series(self, v, z, n=100): 3151 """Compute Struve function & error estimate from its power series.""" 3152 k = arange(0, n) 3153 r = (-1)**k * (.5*z)**(2*k+v+1)/special.gamma(k+1.5)/special.gamma(k+v+1.5) 3154 err = abs(r).max() * finfo(float_).eps * n 3155 return r.sum(), err 3156 3157 def test_vs_series(self): 3158 """Check Struve function versus its power series""" 3159 for v in [-20, -10, -7.99, -3.4, -1, 0, 1, 3.4, 12.49, 16]: 3160 for z in [1, 10, 19, 21, 30]: 3161 value, err = self._series(v, z) 3162 assert_allclose(special.struve(v, z), value, rtol=0, atol=err), (v, z) 3163 3164 def test_some_values(self): 3165 assert_allclose(special.struve(-7.99, 21), 0.0467547614113, rtol=1e-7) 3166 assert_allclose(special.struve(-8.01, 21), 0.0398716951023, rtol=1e-8) 3167 assert_allclose(special.struve(-3.0, 200), 0.0142134427432, rtol=1e-12) 3168 assert_allclose(special.struve(-8.0, -41), 0.0192469727846, rtol=1e-11) 3169 assert_equal(special.struve(-12, -41), -special.struve(-12, 41)) 3170 assert_equal(special.struve(+12, -41), -special.struve(+12, 41)) 3171 assert_equal(special.struve(-11, -41), +special.struve(-11, 41)) 3172 assert_equal(special.struve(+11, -41), +special.struve(+11, 41)) 3173 3174 assert_(isnan(special.struve(-7.1, -1))) 3175 assert_(isnan(special.struve(-10.1, -1))) 3176 3177 def test_regression_679(self): 3178 """Regression test for #679""" 3179 assert_allclose(special.struve(-1.0, 20 - 1e-8), special.struve(-1.0, 20 + 1e-8)) 3180 assert_allclose(special.struve(-2.0, 20 - 1e-8), special.struve(-2.0, 20 + 1e-8)) 3181 assert_allclose(special.struve(-4.3, 20 - 1e-8), special.struve(-4.3, 20 + 1e-8)) 3182 3183 3184def test_chi2_smalldf(): 3185 assert_almost_equal(special.chdtr(0.6,3), 0.957890536704110) 3186 3187 3188def test_ch2_inf(): 3189 assert_equal(special.chdtr(0.7,np.inf), 1.0) 3190 3191 3192def test_chi2c_smalldf(): 3193 assert_almost_equal(special.chdtrc(0.6,3), 1-0.957890536704110) 3194 3195 3196def test_chi2_inv_smalldf(): 3197 assert_almost_equal(special.chdtri(0.6,1-0.957890536704110), 3) 3198 3199 3200def test_agm_simple(): 3201 rtol = 1e-13 3202 3203 # Gauss's constant 3204 assert_allclose(1/special.agm(1, np.sqrt(2)), 0.834626841674073186, 3205 rtol=rtol) 3206 3207 # These values were computed using Wolfram Alpha, with the 3208 # function ArithmeticGeometricMean[a, b]. 3209 agm13 = 1.863616783244897 3210 agm15 = 2.604008190530940 3211 agm35 = 3.936235503649555 3212 assert_allclose(special.agm([[1], [3]], [1, 3, 5]), 3213 [[1, agm13, agm15], 3214 [agm13, 3, agm35]], rtol=rtol) 3215 3216 # Computed by the iteration formula using mpmath, 3217 # with mpmath.mp.prec = 1000: 3218 agm12 = 1.4567910310469068 3219 assert_allclose(special.agm(1, 2), agm12, rtol=rtol) 3220 assert_allclose(special.agm(2, 1), agm12, rtol=rtol) 3221 assert_allclose(special.agm(-1, -2), -agm12, rtol=rtol) 3222 assert_allclose(special.agm(24, 6), 13.458171481725614, rtol=rtol) 3223 assert_allclose(special.agm(13, 123456789.5), 11111458.498599306, 3224 rtol=rtol) 3225 assert_allclose(special.agm(1e30, 1), 2.229223055945383e+28, rtol=rtol) 3226 assert_allclose(special.agm(1e-22, 1), 0.030182566420169886, rtol=rtol) 3227 assert_allclose(special.agm(1e150, 1e180), 2.229223055945383e+178, 3228 rtol=rtol) 3229 assert_allclose(special.agm(1e180, 1e-150), 2.0634722510162677e+177, 3230 rtol=rtol) 3231 assert_allclose(special.agm(1e-150, 1e-170), 3.3112619670463756e-152, 3232 rtol=rtol) 3233 fi = np.finfo(1.0) 3234 assert_allclose(special.agm(fi.tiny, fi.max), 1.9892072050015473e+305, 3235 rtol=rtol) 3236 assert_allclose(special.agm(0.75*fi.max, fi.max), 1.564904312298045e+308, 3237 rtol=rtol) 3238 assert_allclose(special.agm(fi.tiny, 3*fi.tiny), 4.1466849866735005e-308, 3239 rtol=rtol) 3240 3241 # zero, nan and inf cases. 3242 assert_equal(special.agm(0, 0), 0) 3243 assert_equal(special.agm(99, 0), 0) 3244 3245 assert_equal(special.agm(-1, 10), np.nan) 3246 assert_equal(special.agm(0, np.inf), np.nan) 3247 assert_equal(special.agm(np.inf, 0), np.nan) 3248 assert_equal(special.agm(0, -np.inf), np.nan) 3249 assert_equal(special.agm(-np.inf, 0), np.nan) 3250 assert_equal(special.agm(np.inf, -np.inf), np.nan) 3251 assert_equal(special.agm(-np.inf, np.inf), np.nan) 3252 assert_equal(special.agm(1, np.nan), np.nan) 3253 assert_equal(special.agm(np.nan, -1), np.nan) 3254 3255 assert_equal(special.agm(1, np.inf), np.inf) 3256 assert_equal(special.agm(np.inf, 1), np.inf) 3257 assert_equal(special.agm(-1, -np.inf), -np.inf) 3258 assert_equal(special.agm(-np.inf, -1), -np.inf) 3259 3260 3261def test_legacy(): 3262 # Legacy behavior: truncating arguments to integers 3263 with suppress_warnings() as sup: 3264 sup.filter(RuntimeWarning, "floating point number truncated to an integer") 3265 assert_equal(special.expn(1, 0.3), special.expn(1.8, 0.3)) 3266 assert_equal(special.nbdtrc(1, 2, 0.3), special.nbdtrc(1.8, 2.8, 0.3)) 3267 assert_equal(special.nbdtr(1, 2, 0.3), special.nbdtr(1.8, 2.8, 0.3)) 3268 assert_equal(special.nbdtri(1, 2, 0.3), special.nbdtri(1.8, 2.8, 0.3)) 3269 assert_equal(special.pdtri(1, 0.3), special.pdtri(1.8, 0.3)) 3270 assert_equal(special.kn(1, 0.3), special.kn(1.8, 0.3)) 3271 assert_equal(special.yn(1, 0.3), special.yn(1.8, 0.3)) 3272 assert_equal(special.smirnov(1, 0.3), special.smirnov(1.8, 0.3)) 3273 assert_equal(special.smirnovi(1, 0.3), special.smirnovi(1.8, 0.3)) 3274 3275 3276@with_special_errors 3277def test_error_raising(): 3278 assert_raises(special.SpecialFunctionError, special.iv, 1, 1e99j) 3279 3280 3281def test_xlogy(): 3282 def xfunc(x, y): 3283 with np.errstate(invalid='ignore'): 3284 if x == 0 and not np.isnan(y): 3285 return x 3286 else: 3287 return x*np.log(y) 3288 3289 z1 = np.asarray([(0,0), (0, np.nan), (0, np.inf), (1.0, 2.0)], dtype=float) 3290 z2 = np.r_[z1, [(0, 1j), (1, 1j)]] 3291 3292 w1 = np.vectorize(xfunc)(z1[:,0], z1[:,1]) 3293 assert_func_equal(special.xlogy, w1, z1, rtol=1e-13, atol=1e-13) 3294 w2 = np.vectorize(xfunc)(z2[:,0], z2[:,1]) 3295 assert_func_equal(special.xlogy, w2, z2, rtol=1e-13, atol=1e-13) 3296 3297 3298def test_xlog1py(): 3299 def xfunc(x, y): 3300 with np.errstate(invalid='ignore'): 3301 if x == 0 and not np.isnan(y): 3302 return x 3303 else: 3304 return x * np.log1p(y) 3305 3306 z1 = np.asarray([(0,0), (0, np.nan), (0, np.inf), (1.0, 2.0), 3307 (1, 1e-30)], dtype=float) 3308 w1 = np.vectorize(xfunc)(z1[:,0], z1[:,1]) 3309 assert_func_equal(special.xlog1py, w1, z1, rtol=1e-13, atol=1e-13) 3310 3311 3312def test_entr(): 3313 def xfunc(x): 3314 if x < 0: 3315 return -np.inf 3316 else: 3317 return -special.xlogy(x, x) 3318 values = (0, 0.5, 1.0, np.inf) 3319 signs = [-1, 1] 3320 arr = [] 3321 for sgn, v in itertools.product(signs, values): 3322 arr.append(sgn * v) 3323 z = np.array(arr, dtype=float) 3324 w = np.vectorize(xfunc, otypes=[np.float64])(z) 3325 assert_func_equal(special.entr, w, z, rtol=1e-13, atol=1e-13) 3326 3327 3328def test_kl_div(): 3329 def xfunc(x, y): 3330 if x < 0 or y < 0 or (y == 0 and x != 0): 3331 # extension of natural domain to preserve convexity 3332 return np.inf 3333 elif np.isposinf(x) or np.isposinf(y): 3334 # limits within the natural domain 3335 return np.inf 3336 elif x == 0: 3337 return y 3338 else: 3339 return special.xlogy(x, x/y) - x + y 3340 values = (0, 0.5, 1.0) 3341 signs = [-1, 1] 3342 arr = [] 3343 for sgna, va, sgnb, vb in itertools.product(signs, values, signs, values): 3344 arr.append((sgna*va, sgnb*vb)) 3345 z = np.array(arr, dtype=float) 3346 w = np.vectorize(xfunc, otypes=[np.float64])(z[:,0], z[:,1]) 3347 assert_func_equal(special.kl_div, w, z, rtol=1e-13, atol=1e-13) 3348 3349 3350def test_rel_entr(): 3351 def xfunc(x, y): 3352 if x > 0 and y > 0: 3353 return special.xlogy(x, x/y) 3354 elif x == 0 and y >= 0: 3355 return 0 3356 else: 3357 return np.inf 3358 values = (0, 0.5, 1.0) 3359 signs = [-1, 1] 3360 arr = [] 3361 for sgna, va, sgnb, vb in itertools.product(signs, values, signs, values): 3362 arr.append((sgna*va, sgnb*vb)) 3363 z = np.array(arr, dtype=float) 3364 w = np.vectorize(xfunc, otypes=[np.float64])(z[:,0], z[:,1]) 3365 assert_func_equal(special.rel_entr, w, z, rtol=1e-13, atol=1e-13) 3366 3367 3368def test_huber(): 3369 assert_equal(special.huber(-1, 1.5), np.inf) 3370 assert_allclose(special.huber(2, 1.5), 0.5 * np.square(1.5)) 3371 assert_allclose(special.huber(2, 2.5), 2 * (2.5 - 0.5 * 2)) 3372 3373 def xfunc(delta, r): 3374 if delta < 0: 3375 return np.inf 3376 elif np.abs(r) < delta: 3377 return 0.5 * np.square(r) 3378 else: 3379 return delta * (np.abs(r) - 0.5 * delta) 3380 3381 z = np.random.randn(10, 2) 3382 w = np.vectorize(xfunc, otypes=[np.float64])(z[:,0], z[:,1]) 3383 assert_func_equal(special.huber, w, z, rtol=1e-13, atol=1e-13) 3384 3385 3386def test_pseudo_huber(): 3387 def xfunc(delta, r): 3388 if delta < 0: 3389 return np.inf 3390 elif (not delta) or (not r): 3391 return 0 3392 else: 3393 return delta**2 * (np.sqrt(1 + (r/delta)**2) - 1) 3394 3395 z = np.array(np.random.randn(10, 2).tolist() + [[0, 0.5], [0.5, 0]]) 3396 w = np.vectorize(xfunc, otypes=[np.float64])(z[:,0], z[:,1]) 3397 assert_func_equal(special.pseudo_huber, w, z, rtol=1e-13, atol=1e-13) 3398