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