specialfunctions.h

Go to the documentation of this file.
1 /*************************************************************************
2 ALGLIB 3.11.0 (source code generated 2017-05-11)
3 Copyright (c) Sergey Bochkanov (ALGLIB project).
4 
5 >>> SOURCE LICENSE >>>
6 This program is free software; you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation (www.fsf.org); either version 2 of the
9 License, or (at your option) any later version.
10 
11 This program is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15 
16 A copy of the GNU General Public License is available at
17 http://www.fsf.org/licensing/licenses
18 >>> END OF LICENSE >>>
19 *************************************************************************/
20 #ifndef _specialfunctions_pkg_h
21 #define _specialfunctions_pkg_h
22 #include "ap.h"
23 #include "alglibinternal.h"
24 
26 //
27 // THIS SECTION CONTAINS COMPUTATIONAL CORE DECLARATIONS (DATATYPES)
28 //
30 namespace alglib_impl
31 {
32 
33 }
34 
36 //
37 // THIS SECTION CONTAINS C++ INTERFACE
38 //
40 namespace alglib
41 {
42 
43 
44 /*************************************************************************
45 Gamma function
46 
47 Input parameters:
48  X - argument
49 
50 Domain:
51  0 < X < 171.6
52  -170 < X < 0, X is not an integer.
53 
54 Relative error:
55  arithmetic domain # trials peak rms
56  IEEE -170,-33 20000 2.3e-15 3.3e-16
57  IEEE -33, 33 20000 9.4e-16 2.2e-16
58  IEEE 33, 171.6 20000 2.3e-15 3.2e-16
59 
60 Cephes Math Library Release 2.8: June, 2000
61 Original copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier
62 Translated to AlgoPascal by Bochkanov Sergey (2005, 2006, 2007).
63 *************************************************************************/
64 double gammafunction(const double x);
65 
66 
67 /*************************************************************************
68 Natural logarithm of gamma function
69 
70 Input parameters:
71  X - argument
72 
73 Result:
74  logarithm of the absolute value of the Gamma(X).
75 
76 Output parameters:
77  SgnGam - sign(Gamma(X))
78 
79 Domain:
80  0 < X < 2.55e305
81  -2.55e305 < X < 0, X is not an integer.
82 
83 ACCURACY:
84 arithmetic domain # trials peak rms
85  IEEE 0, 3 28000 5.4e-16 1.1e-16
86  IEEE 2.718, 2.556e305 40000 3.5e-16 8.3e-17
87 The error criterion was relative when the function magnitude
88 was greater than one but absolute when it was less than one.
89 
90 The following test used the relative error criterion, though
91 at certain points the relative error could be much higher than
92 indicated.
93  IEEE -200, -4 10000 4.8e-16 1.3e-16
94 
95 Cephes Math Library Release 2.8: June, 2000
96 Copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier
97 Translated to AlgoPascal by Bochkanov Sergey (2005, 2006, 2007).
98 *************************************************************************/
99 double lngamma(const double x, double &sgngam);
100 
101 /*************************************************************************
102 Error function
103 
104 The integral is
105 
106  x
107  -
108  2 | | 2
109  erf(x) = -------- | exp( - t ) dt.
110  sqrt(pi) | |
111  -
112  0
113 
114 For 0 <= |x| < 1, erf(x) = x * P4(x**2)/Q5(x**2); otherwise
115 erf(x) = 1 - erfc(x).
116 
117 
118 ACCURACY:
119 
120  Relative error:
121 arithmetic domain # trials peak rms
122  IEEE 0,1 30000 3.7e-16 1.0e-16
123 
124 Cephes Math Library Release 2.8: June, 2000
125 Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier
126 *************************************************************************/
127 double errorfunction(const double x);
128 
129 
130 /*************************************************************************
131 Complementary error function
132 
133  1 - erf(x) =
134 
135  inf.
136  -
137  2 | | 2
138  erfc(x) = -------- | exp( - t ) dt
139  sqrt(pi) | |
140  -
141  x
142 
143 
144 For small x, erfc(x) = 1 - erf(x); otherwise rational
145 approximations are computed.
146 
147 
148 ACCURACY:
149 
150  Relative error:
151 arithmetic domain # trials peak rms
152  IEEE 0,26.6417 30000 5.7e-14 1.5e-14
153 
154 Cephes Math Library Release 2.8: June, 2000
155 Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier
156 *************************************************************************/
157 double errorfunctionc(const double x);
158 
159 
160 /*************************************************************************
161 Normal distribution function
162 
163 Returns the area under the Gaussian probability density
164 function, integrated from minus infinity to x:
165 
166  x
167  -
168  1 | | 2
169  ndtr(x) = --------- | exp( - t /2 ) dt
170  sqrt(2pi) | |
171  -
172  -inf.
173 
174  = ( 1 + erf(z) ) / 2
175  = erfc(z) / 2
176 
177 where z = x/sqrt(2). Computation is via the functions
178 erf and erfc.
179 
180 
181 ACCURACY:
182 
183  Relative error:
184 arithmetic domain # trials peak rms
185  IEEE -13,0 30000 3.4e-14 6.7e-15
186 
187 Cephes Math Library Release 2.8: June, 2000
188 Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier
189 *************************************************************************/
190 double normaldistribution(const double x);
191 
192 
193 /*************************************************************************
194 Inverse of the error function
195 
196 Cephes Math Library Release 2.8: June, 2000
197 Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier
198 *************************************************************************/
199 double inverf(const double e);
200 
201 
202 /*************************************************************************
203 Inverse of Normal distribution function
204 
205 Returns the argument, x, for which the area under the
206 Gaussian probability density function (integrated from
207 minus infinity to x) is equal to y.
208 
209 
210 For small arguments 0 < y < exp(-2), the program computes
211 z = sqrt( -2.0 * log(y) ); then the approximation is
212 x = z - log(z)/z - (1/z) P(1/z) / Q(1/z).
213 There are two rational functions P/Q, one for 0 < y < exp(-32)
214 and the other for y up to exp(-2). For larger arguments,
215 w = y - 0.5, and x/sqrt(2pi) = w + w**3 R(w**2)/S(w**2)).
216 
217 ACCURACY:
218 
219  Relative error:
220 arithmetic domain # trials peak rms
221  IEEE 0.125, 1 20000 7.2e-16 1.3e-16
222  IEEE 3e-308, 0.135 50000 4.6e-16 9.8e-17
223 
224 Cephes Math Library Release 2.8: June, 2000
225 Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier
226 *************************************************************************/
227 double invnormaldistribution(const double y0);
228 
229 /*************************************************************************
230 Incomplete gamma integral
231 
232 The function is defined by
233 
234  x
235  -
236  1 | | -t a-1
237  igam(a,x) = ----- | e t dt.
238  - | |
239  | (a) -
240  0
241 
242 
243 In this implementation both arguments must be positive.
244 The integral is evaluated by either a power series or
245 continued fraction expansion, depending on the relative
246 values of a and x.
247 
248 ACCURACY:
249 
250  Relative error:
251 arithmetic domain # trials peak rms
252  IEEE 0,30 200000 3.6e-14 2.9e-15
253  IEEE 0,100 300000 9.9e-14 1.5e-14
254 
255 Cephes Math Library Release 2.8: June, 2000
256 Copyright 1985, 1987, 2000 by Stephen L. Moshier
257 *************************************************************************/
258 double incompletegamma(const double a, const double x);
259 
260 
261 /*************************************************************************
262 Complemented incomplete gamma integral
263 
264 The function is defined by
265 
266 
267  igamc(a,x) = 1 - igam(a,x)
268 
269  inf.
270  -
271  1 | | -t a-1
272  = ----- | e t dt.
273  - | |
274  | (a) -
275  x
276 
277 
278 In this implementation both arguments must be positive.
279 The integral is evaluated by either a power series or
280 continued fraction expansion, depending on the relative
281 values of a and x.
282 
283 ACCURACY:
284 
285 Tested at random a, x.
286  a x Relative error:
287 arithmetic domain domain # trials peak rms
288  IEEE 0.5,100 0,100 200000 1.9e-14 1.7e-15
289  IEEE 0.01,0.5 0,100 200000 1.4e-13 1.6e-15
290 
291 Cephes Math Library Release 2.8: June, 2000
292 Copyright 1985, 1987, 2000 by Stephen L. Moshier
293 *************************************************************************/
294 double incompletegammac(const double a, const double x);
295 
296 
297 /*************************************************************************
298 Inverse of complemented imcomplete gamma integral
299 
300 Given p, the function finds x such that
301 
302  igamc( a, x ) = p.
303 
304 Starting with the approximate value
305 
306  3
307  x = a t
308 
309  where
310 
311  t = 1 - d - ndtri(p) sqrt(d)
312 
313 and
314 
315  d = 1/9a,
316 
317 the routine performs up to 10 Newton iterations to find the
318 root of igamc(a,x) - p = 0.
319 
320 ACCURACY:
321 
322 Tested at random a, p in the intervals indicated.
323 
324  a p Relative error:
325 arithmetic domain domain # trials peak rms
326  IEEE 0.5,100 0,0.5 100000 1.0e-14 1.7e-15
327  IEEE 0.01,0.5 0,0.5 100000 9.0e-14 3.4e-15
328  IEEE 0.5,10000 0,0.5 20000 2.3e-13 3.8e-14
329 
330 Cephes Math Library Release 2.8: June, 2000
331 Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier
332 *************************************************************************/
333 double invincompletegammac(const double a, const double y0);
334 
335 /*************************************************************************
336 Complete elliptic integral of the first kind
337 
338 Approximates the integral
339 
340 
341 
342  pi/2
343  -
344  | |
345  | dt
346 K(m) = | ------------------
347  | 2
348  | | sqrt( 1 - m sin t )
349  -
350  0
351 
352 using the approximation
353 
354  P(x) - log x Q(x).
355 
356 ACCURACY:
357 
358  Relative error:
359 arithmetic domain # trials peak rms
360  IEEE 0,1 30000 2.5e-16 6.8e-17
361 
362 Cephes Math Library, Release 2.8: June, 2000
363 Copyright 1984, 1987, 2000 by Stephen L. Moshier
364 *************************************************************************/
365 double ellipticintegralk(const double m);
366 
367 
368 /*************************************************************************
369 Complete elliptic integral of the first kind
370 
371 Approximates the integral
372 
373 
374 
375  pi/2
376  -
377  | |
378  | dt
379 K(m) = | ------------------
380  | 2
381  | | sqrt( 1 - m sin t )
382  -
383  0
384 
385 where m = 1 - m1, using the approximation
386 
387  P(x) - log x Q(x).
388 
389 The argument m1 is used rather than m so that the logarithmic
390 singularity at m = 1 will be shifted to the origin; this
391 preserves maximum accuracy.
392 
393 K(0) = pi/2.
394 
395 ACCURACY:
396 
397  Relative error:
398 arithmetic domain # trials peak rms
399  IEEE 0,1 30000 2.5e-16 6.8e-17
400 
401 Cephes Math Library, Release 2.8: June, 2000
402 Copyright 1984, 1987, 2000 by Stephen L. Moshier
403 *************************************************************************/
404 double ellipticintegralkhighprecision(const double m1);
405 
406 
407 /*************************************************************************
408 Incomplete elliptic integral of the first kind F(phi|m)
409 
410 Approximates the integral
411 
412 
413 
414  phi
415  -
416  | |
417  | dt
418 F(phi_\m) = | ------------------
419  | 2
420  | | sqrt( 1 - m sin t )
421  -
422  0
423 
424 of amplitude phi and modulus m, using the arithmetic -
425 geometric mean algorithm.
426 
427 
428 
429 
430 ACCURACY:
431 
432 Tested at random points with m in [0, 1] and phi as indicated.
433 
434  Relative error:
435 arithmetic domain # trials peak rms
436  IEEE -10,10 200000 7.4e-16 1.0e-16
437 
438 Cephes Math Library Release 2.8: June, 2000
439 Copyright 1984, 1987, 2000 by Stephen L. Moshier
440 *************************************************************************/
441 double incompleteellipticintegralk(const double phi, const double m);
442 
443 
444 /*************************************************************************
445 Complete elliptic integral of the second kind
446 
447 Approximates the integral
448 
449 
450  pi/2
451  -
452  | | 2
453 E(m) = | sqrt( 1 - m sin t ) dt
454  | |
455  -
456  0
457 
458 using the approximation
459 
460  P(x) - x log x Q(x).
461 
462 ACCURACY:
463 
464  Relative error:
465 arithmetic domain # trials peak rms
466  IEEE 0, 1 10000 2.1e-16 7.3e-17
467 
468 Cephes Math Library, Release 2.8: June, 2000
469 Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
470 *************************************************************************/
471 double ellipticintegrale(const double m);
472 
473 
474 /*************************************************************************
475 Incomplete elliptic integral of the second kind
476 
477 Approximates the integral
478 
479 
480  phi
481  -
482  | |
483  | 2
484 E(phi_\m) = | sqrt( 1 - m sin t ) dt
485  |
486  | |
487  -
488  0
489 
490 of amplitude phi and modulus m, using the arithmetic -
491 geometric mean algorithm.
492 
493 ACCURACY:
494 
495 Tested at random arguments with phi in [-10, 10] and m in
496 [0, 1].
497  Relative error:
498 arithmetic domain # trials peak rms
499  IEEE -10,10 150000 3.3e-15 1.4e-16
500 
501 Cephes Math Library Release 2.8: June, 2000
502 Copyright 1984, 1987, 1993, 2000 by Stephen L. Moshier
503 *************************************************************************/
504 double incompleteellipticintegrale(const double phi, const double m);
505 
506 /*************************************************************************
507 Calculation of the value of the Hermite polynomial.
508 
509 Parameters:
510  n - degree, n>=0
511  x - argument
512 
513 Result:
514  the value of the Hermite polynomial Hn at x
515 *************************************************************************/
516 double hermitecalculate(const ae_int_t n, const double x);
517 
518 
519 /*************************************************************************
520 Summation of Hermite polynomials using Clenshaw's recurrence formula.
521 
522 This routine calculates
523  c[0]*H0(x) + c[1]*H1(x) + ... + c[N]*HN(x)
524 
525 Parameters:
526  n - degree, n>=0
527  x - argument
528 
529 Result:
530  the value of the Hermite polynomial at x
531 *************************************************************************/
532 double hermitesum(const real_1d_array &c, const ae_int_t n, const double x);
533 
534 
535 /*************************************************************************
536 Representation of Hn as C[0] + C[1]*X + ... + C[N]*X^N
537 
538 Input parameters:
539  N - polynomial degree, n>=0
540 
541 Output parameters:
542  C - coefficients
543 *************************************************************************/
545 
546 /*************************************************************************
547 Dawson's Integral
548 
549 Approximates the integral
550 
551  x
552  -
553  2 | | 2
554  dawsn(x) = exp( -x ) | exp( t ) dt
555  | |
556  -
557  0
558 
559 Three different rational approximations are employed, for
560 the intervals 0 to 3.25; 3.25 to 6.25; and 6.25 up.
561 
562 ACCURACY:
563 
564  Relative error:
565 arithmetic domain # trials peak rms
566  IEEE 0,10 10000 6.9e-16 1.0e-16
567 
568 Cephes Math Library Release 2.8: June, 2000
569 Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
570 *************************************************************************/
571 double dawsonintegral(const double x);
572 
573 /*************************************************************************
574 Sine and cosine integrals
575 
576 Evaluates the integrals
577 
578  x
579  -
580  | cos t - 1
581  Ci(x) = eul + ln x + | --------- dt,
582  | t
583  -
584  0
585  x
586  -
587  | sin t
588  Si(x) = | ----- dt
589  | t
590  -
591  0
592 
593 where eul = 0.57721566490153286061 is Euler's constant.
594 The integrals are approximated by rational functions.
595 For x > 8 auxiliary functions f(x) and g(x) are employed
596 such that
597 
598 Ci(x) = f(x) sin(x) - g(x) cos(x)
599 Si(x) = pi/2 - f(x) cos(x) - g(x) sin(x)
600 
601 
602 ACCURACY:
603  Test interval = [0,50].
604 Absolute error, except relative when > 1:
605 arithmetic function # trials peak rms
606  IEEE Si 30000 4.4e-16 7.3e-17
607  IEEE Ci 30000 6.9e-16 5.1e-17
608 
609 Cephes Math Library Release 2.1: January, 1989
610 Copyright 1984, 1987, 1989 by Stephen L. Moshier
611 *************************************************************************/
612 void sinecosineintegrals(const double x, double &si, double &ci);
613 
614 
615 /*************************************************************************
616 Hyperbolic sine and cosine integrals
617 
618 Approximates the integrals
619 
620  x
621  -
622  | | cosh t - 1
623  Chi(x) = eul + ln x + | ----------- dt,
624  | | t
625  -
626  0
627 
628  x
629  -
630  | | sinh t
631  Shi(x) = | ------ dt
632  | | t
633  -
634  0
635 
636 where eul = 0.57721566490153286061 is Euler's constant.
637 The integrals are evaluated by power series for x < 8
638 and by Chebyshev expansions for x between 8 and 88.
639 For large x, both functions approach exp(x)/2x.
640 Arguments greater than 88 in magnitude return MAXNUM.
641 
642 
643 ACCURACY:
644 
645 Test interval 0 to 88.
646  Relative error:
647 arithmetic function # trials peak rms
648  IEEE Shi 30000 6.9e-16 1.6e-16
649  Absolute error, except relative when |Chi| > 1:
650  IEEE Chi 30000 8.4e-16 1.4e-16
651 
652 Cephes Math Library Release 2.8: June, 2000
653 Copyright 1984, 1987, 2000 by Stephen L. Moshier
654 *************************************************************************/
655 void hyperbolicsinecosineintegrals(const double x, double &shi, double &chi);
656 
657 /*************************************************************************
658 Poisson distribution
659 
660 Returns the sum of the first k+1 terms of the Poisson
661 distribution:
662 
663  k j
664  -- -m m
665  > e --
666  -- j!
667  j=0
668 
669 The terms are not summed directly; instead the incomplete
670 gamma integral is employed, according to the relation
671 
672 y = pdtr( k, m ) = igamc( k+1, m ).
673 
674 The arguments must both be positive.
675 ACCURACY:
676 
677 See incomplete gamma function
678 
679 Cephes Math Library Release 2.8: June, 2000
680 Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier
681 *************************************************************************/
682 double poissondistribution(const ae_int_t k, const double m);
683 
684 
685 /*************************************************************************
686 Complemented Poisson distribution
687 
688 Returns the sum of the terms k+1 to infinity of the Poisson
689 distribution:
690 
691  inf. j
692  -- -m m
693  > e --
694  -- j!
695  j=k+1
696 
697 The terms are not summed directly; instead the incomplete
698 gamma integral is employed, according to the formula
699 
700 y = pdtrc( k, m ) = igam( k+1, m ).
701 
702 The arguments must both be positive.
703 
704 ACCURACY:
705 
706 See incomplete gamma function
707 
708 Cephes Math Library Release 2.8: June, 2000
709 Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier
710 *************************************************************************/
711 double poissoncdistribution(const ae_int_t k, const double m);
712 
713 
714 /*************************************************************************
715 Inverse Poisson distribution
716 
717 Finds the Poisson variable x such that the integral
718 from 0 to x of the Poisson density is equal to the
719 given probability y.
720 
721 This is accomplished using the inverse gamma integral
722 function and the relation
723 
724  m = igami( k+1, y ).
725 
726 ACCURACY:
727 
728 See inverse incomplete gamma function
729 
730 Cephes Math Library Release 2.8: June, 2000
731 Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier
732 *************************************************************************/
733 double invpoissondistribution(const ae_int_t k, const double y);
734 
735 /*************************************************************************
736 Bessel function of order zero
737 
738 Returns Bessel function of order zero of the argument.
739 
740 The domain is divided into the intervals [0, 5] and
741 (5, infinity). In the first interval the following rational
742 approximation is used:
743 
744 
745  2 2
746 (w - r ) (w - r ) P (w) / Q (w)
747  1 2 3 8
748 
749  2
750 where w = x and the two r's are zeros of the function.
751 
752 In the second interval, the Hankel asymptotic expansion
753 is employed with two rational functions of degree 6/6
754 and 7/7.
755 
756 ACCURACY:
757 
758  Absolute error:
759 arithmetic domain # trials peak rms
760  IEEE 0, 30 60000 4.2e-16 1.1e-16
761 
762 Cephes Math Library Release 2.8: June, 2000
763 Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
764 *************************************************************************/
765 double besselj0(const double x);
766 
767 
768 /*************************************************************************
769 Bessel function of order one
770 
771 Returns Bessel function of order one of the argument.
772 
773 The domain is divided into the intervals [0, 8] and
774 (8, infinity). In the first interval a 24 term Chebyshev
775 expansion is used. In the second, the asymptotic
776 trigonometric representation is employed using two
777 rational functions of degree 5/5.
778 
779 ACCURACY:
780 
781  Absolute error:
782 arithmetic domain # trials peak rms
783  IEEE 0, 30 30000 2.6e-16 1.1e-16
784 
785 Cephes Math Library Release 2.8: June, 2000
786 Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
787 *************************************************************************/
788 double besselj1(const double x);
789 
790 
791 /*************************************************************************
792 Bessel function of integer order
793 
794 Returns Bessel function of order n, where n is a
795 (possibly negative) integer.
796 
797 The ratio of jn(x) to j0(x) is computed by backward
798 recurrence. First the ratio jn/jn-1 is found by a
799 continued fraction expansion. Then the recurrence
800 relating successive orders is applied until j0 or j1 is
801 reached.
802 
803 If n = 0 or 1 the routine for j0 or j1 is called
804 directly.
805 
806 ACCURACY:
807 
808  Absolute error:
809 arithmetic range # trials peak rms
810  IEEE 0, 30 5000 4.4e-16 7.9e-17
811 
812 
813 Not suitable for large n or x. Use jv() (fractional order) instead.
814 
815 Cephes Math Library Release 2.8: June, 2000
816 Copyright 1984, 1987, 2000 by Stephen L. Moshier
817 *************************************************************************/
818 double besseljn(const ae_int_t n, const double x);
819 
820 
821 /*************************************************************************
822 Bessel function of the second kind, order zero
823 
824 Returns Bessel function of the second kind, of order
825 zero, of the argument.
826 
827 The domain is divided into the intervals [0, 5] and
828 (5, infinity). In the first interval a rational approximation
829 R(x) is employed to compute
830  y0(x) = R(x) + 2 * log(x) * j0(x) / PI.
831 Thus a call to j0() is required.
832 
833 In the second interval, the Hankel asymptotic expansion
834 is employed with two rational functions of degree 6/6
835 and 7/7.
836 
837 
838 
839 ACCURACY:
840 
841  Absolute error, when y0(x) < 1; else relative error:
842 
843 arithmetic domain # trials peak rms
844  IEEE 0, 30 30000 1.3e-15 1.6e-16
845 
846 Cephes Math Library Release 2.8: June, 2000
847 Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
848 *************************************************************************/
849 double bessely0(const double x);
850 
851 
852 /*************************************************************************
853 Bessel function of second kind of order one
854 
855 Returns Bessel function of the second kind of order one
856 of the argument.
857 
858 The domain is divided into the intervals [0, 8] and
859 (8, infinity). In the first interval a 25 term Chebyshev
860 expansion is used, and a call to j1() is required.
861 In the second, the asymptotic trigonometric representation
862 is employed using two rational functions of degree 5/5.
863 
864 ACCURACY:
865 
866  Absolute error:
867 arithmetic domain # trials peak rms
868  IEEE 0, 30 30000 1.0e-15 1.3e-16
869 
870 Cephes Math Library Release 2.8: June, 2000
871 Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
872 *************************************************************************/
873 double bessely1(const double x);
874 
875 
876 /*************************************************************************
877 Bessel function of second kind of integer order
878 
879 Returns Bessel function of order n, where n is a
880 (possibly negative) integer.
881 
882 The function is evaluated by forward recurrence on
883 n, starting with values computed by the routines
884 y0() and y1().
885 
886 If n = 0 or 1 the routine for y0 or y1 is called
887 directly.
888 
889 ACCURACY:
890  Absolute error, except relative
891  when y > 1:
892 arithmetic domain # trials peak rms
893  IEEE 0, 30 30000 3.4e-15 4.3e-16
894 
895 Cephes Math Library Release 2.8: June, 2000
896 Copyright 1984, 1987, 2000 by Stephen L. Moshier
897 *************************************************************************/
898 double besselyn(const ae_int_t n, const double x);
899 
900 
901 /*************************************************************************
902 Modified Bessel function of order zero
903 
904 Returns modified Bessel function of order zero of the
905 argument.
906 
907 The function is defined as i0(x) = j0( ix ).
908 
909 The range is partitioned into the two intervals [0,8] and
910 (8, infinity). Chebyshev polynomial expansions are employed
911 in each interval.
912 
913 ACCURACY:
914 
915  Relative error:
916 arithmetic domain # trials peak rms
917  IEEE 0,30 30000 5.8e-16 1.4e-16
918 
919 Cephes Math Library Release 2.8: June, 2000
920 Copyright 1984, 1987, 2000 by Stephen L. Moshier
921 *************************************************************************/
922 double besseli0(const double x);
923 
924 
925 /*************************************************************************
926 Modified Bessel function of order one
927 
928 Returns modified Bessel function of order one of the
929 argument.
930 
931 The function is defined as i1(x) = -i j1( ix ).
932 
933 The range is partitioned into the two intervals [0,8] and
934 (8, infinity). Chebyshev polynomial expansions are employed
935 in each interval.
936 
937 ACCURACY:
938 
939  Relative error:
940 arithmetic domain # trials peak rms
941  IEEE 0, 30 30000 1.9e-15 2.1e-16
942 
943 Cephes Math Library Release 2.8: June, 2000
944 Copyright 1985, 1987, 2000 by Stephen L. Moshier
945 *************************************************************************/
946 double besseli1(const double x);
947 
948 
949 /*************************************************************************
950 Modified Bessel function, second kind, order zero
951 
952 Returns modified Bessel function of the second kind
953 of order zero of the argument.
954 
955 The range is partitioned into the two intervals [0,8] and
956 (8, infinity). Chebyshev polynomial expansions are employed
957 in each interval.
958 
959 ACCURACY:
960 
961 Tested at 2000 random points between 0 and 8. Peak absolute
962 error (relative when K0 > 1) was 1.46e-14; rms, 4.26e-15.
963  Relative error:
964 arithmetic domain # trials peak rms
965  IEEE 0, 30 30000 1.2e-15 1.6e-16
966 
967 Cephes Math Library Release 2.8: June, 2000
968 Copyright 1984, 1987, 2000 by Stephen L. Moshier
969 *************************************************************************/
970 double besselk0(const double x);
971 
972 
973 /*************************************************************************
974 Modified Bessel function, second kind, order one
975 
976 Computes the modified Bessel function of the second kind
977 of order one of the argument.
978 
979 The range is partitioned into the two intervals [0,2] and
980 (2, infinity). Chebyshev polynomial expansions are employed
981 in each interval.
982 
983 ACCURACY:
984 
985  Relative error:
986 arithmetic domain # trials peak rms
987  IEEE 0, 30 30000 1.2e-15 1.6e-16
988 
989 Cephes Math Library Release 2.8: June, 2000
990 Copyright 1984, 1987, 2000 by Stephen L. Moshier
991 *************************************************************************/
992 double besselk1(const double x);
993 
994 
995 /*************************************************************************
996 Modified Bessel function, second kind, integer order
997 
998 Returns modified Bessel function of the second kind
999 of order n of the argument.
1000 
1001 The range is partitioned into the two intervals [0,9.55] and
1002 (9.55, infinity). An ascending power series is used in the
1003 low range, and an asymptotic expansion in the high range.
1004 
1005 ACCURACY:
1006 
1007  Relative error:
1008 arithmetic domain # trials peak rms
1009  IEEE 0,30 90000 1.8e-8 3.0e-10
1010 
1011 Error is high only near the crossover point x = 9.55
1012 between the two expansions used.
1013 
1014 Cephes Math Library Release 2.8: June, 2000
1015 Copyright 1984, 1987, 1988, 2000 by Stephen L. Moshier
1016 *************************************************************************/
1017 double besselkn(const ae_int_t nn, const double x);
1018 
1019 /*************************************************************************
1020 Incomplete beta integral
1021 
1022 Returns incomplete beta integral of the arguments, evaluated
1023 from zero to x. The function is defined as
1024 
1025  x
1026  - -
1027  | (a+b) | | a-1 b-1
1028  ----------- | t (1-t) dt.
1029  - - | |
1030  | (a) | (b) -
1031  0
1032 
1033 The domain of definition is 0 <= x <= 1. In this
1034 implementation a and b are restricted to positive values.
1035 The integral from x to 1 may be obtained by the symmetry
1036 relation
1037 
1038  1 - incbet( a, b, x ) = incbet( b, a, 1-x ).
1039 
1040 The integral is evaluated by a continued fraction expansion
1041 or, when b*x is small, by a power series.
1042 
1043 ACCURACY:
1044 
1045 Tested at uniformly distributed random points (a,b,x) with a and b
1046 in "domain" and x between 0 and 1.
1047  Relative error
1048 arithmetic domain # trials peak rms
1049  IEEE 0,5 10000 6.9e-15 4.5e-16
1050  IEEE 0,85 250000 2.2e-13 1.7e-14
1051  IEEE 0,1000 30000 5.3e-12 6.3e-13
1052  IEEE 0,10000 250000 9.3e-11 7.1e-12
1053  IEEE 0,100000 10000 8.7e-10 4.8e-11
1054 Outputs smaller than the IEEE gradual underflow threshold
1055 were excluded from these statistics.
1056 
1057 Cephes Math Library, Release 2.8: June, 2000
1058 Copyright 1984, 1995, 2000 by Stephen L. Moshier
1059 *************************************************************************/
1060 double incompletebeta(const double a, const double b, const double x);
1061 
1062 
1063 /*************************************************************************
1064 Inverse of imcomplete beta integral
1065 
1066 Given y, the function finds x such that
1067 
1068  incbet( a, b, x ) = y .
1069 
1070 The routine performs interval halving or Newton iterations to find the
1071 root of incbet(a,b,x) - y = 0.
1072 
1073 
1074 ACCURACY:
1075 
1076  Relative error:
1077  x a,b
1078 arithmetic domain domain # trials peak rms
1079  IEEE 0,1 .5,10000 50000 5.8e-12 1.3e-13
1080  IEEE 0,1 .25,100 100000 1.8e-13 3.9e-15
1081  IEEE 0,1 0,5 50000 1.1e-12 5.5e-15
1082 With a and b constrained to half-integer or integer values:
1083  IEEE 0,1 .5,10000 50000 5.8e-12 1.1e-13
1084  IEEE 0,1 .5,100 100000 1.7e-14 7.9e-16
1085 With a = .5, b constrained to half-integer or integer values:
1086  IEEE 0,1 .5,10000 10000 8.3e-11 1.0e-11
1087 
1088 Cephes Math Library Release 2.8: June, 2000
1089 Copyright 1984, 1996, 2000 by Stephen L. Moshier
1090 *************************************************************************/
1091 double invincompletebeta(const double a, const double b, const double y);
1092 
1093 /*************************************************************************
1094 F distribution
1095 
1096 Returns the area from zero to x under the F density
1097 function (also known as Snedcor's density or the
1098 variance ratio density). This is the density
1099 of x = (u1/df1)/(u2/df2), where u1 and u2 are random
1100 variables having Chi square distributions with df1
1101 and df2 degrees of freedom, respectively.
1102 The incomplete beta integral is used, according to the
1103 formula
1104 
1105 P(x) = incbet( df1/2, df2/2, (df1*x/(df2 + df1*x) ).
1106 
1107 
1108 The arguments a and b are greater than zero, and x is
1109 nonnegative.
1110 
1111 ACCURACY:
1112 
1113 Tested at random points (a,b,x).
1114 
1115  x a,b Relative error:
1116 arithmetic domain domain # trials peak rms
1117  IEEE 0,1 0,100 100000 9.8e-15 1.7e-15
1118  IEEE 1,5 0,100 100000 6.5e-15 3.5e-16
1119  IEEE 0,1 1,10000 100000 2.2e-11 3.3e-12
1120  IEEE 1,5 1,10000 100000 1.1e-11 1.7e-13
1121 
1122 Cephes Math Library Release 2.8: June, 2000
1123 Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier
1124 *************************************************************************/
1125 double fdistribution(const ae_int_t a, const ae_int_t b, const double x);
1126 
1127 
1128 /*************************************************************************
1129 Complemented F distribution
1130 
1131 Returns the area from x to infinity under the F density
1132 function (also known as Snedcor's density or the
1133 variance ratio density).
1134 
1135 
1136  inf.
1137  -
1138  1 | | a-1 b-1
1139 1-P(x) = ------ | t (1-t) dt
1140  B(a,b) | |
1141  -
1142  x
1143 
1144 
1145 The incomplete beta integral is used, according to the
1146 formula
1147 
1148 P(x) = incbet( df2/2, df1/2, (df2/(df2 + df1*x) ).
1149 
1150 
1151 ACCURACY:
1152 
1153 Tested at random points (a,b,x) in the indicated intervals.
1154  x a,b Relative error:
1155 arithmetic domain domain # trials peak rms
1156  IEEE 0,1 1,100 100000 3.7e-14 5.9e-16
1157  IEEE 1,5 1,100 100000 8.0e-15 1.6e-15
1158  IEEE 0,1 1,10000 100000 1.8e-11 3.5e-13
1159  IEEE 1,5 1,10000 100000 2.0e-11 3.0e-12
1160 
1161 Cephes Math Library Release 2.8: June, 2000
1162 Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier
1163 *************************************************************************/
1164 double fcdistribution(const ae_int_t a, const ae_int_t b, const double x);
1165 
1166 
1167 /*************************************************************************
1168 Inverse of complemented F distribution
1169 
1170 Finds the F density argument x such that the integral
1171 from x to infinity of the F density is equal to the
1172 given probability p.
1173 
1174 This is accomplished using the inverse beta integral
1175 function and the relations
1176 
1177  z = incbi( df2/2, df1/2, p )
1178  x = df2 (1-z) / (df1 z).
1179 
1180 Note: the following relations hold for the inverse of
1181 the uncomplemented F distribution:
1182 
1183  z = incbi( df1/2, df2/2, p )
1184  x = df2 z / (df1 (1-z)).
1185 
1186 ACCURACY:
1187 
1188 Tested at random points (a,b,p).
1189 
1190  a,b Relative error:
1191 arithmetic domain # trials peak rms
1192  For p between .001 and 1:
1193  IEEE 1,100 100000 8.3e-15 4.7e-16
1194  IEEE 1,10000 100000 2.1e-11 1.4e-13
1195  For p between 10^-6 and 10^-3:
1196  IEEE 1,100 50000 1.3e-12 8.4e-15
1197  IEEE 1,10000 50000 3.0e-12 4.8e-14
1198 
1199 Cephes Math Library Release 2.8: June, 2000
1200 Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier
1201 *************************************************************************/
1202 double invfdistribution(const ae_int_t a, const ae_int_t b, const double y);
1203 
1204 /*************************************************************************
1205 Fresnel integral
1206 
1207 Evaluates the Fresnel integrals
1208 
1209  x
1210  -
1211  | |
1212 C(x) = | cos(pi/2 t**2) dt,
1213  | |
1214  -
1215  0
1216 
1217  x
1218  -
1219  | |
1220 S(x) = | sin(pi/2 t**2) dt.
1221  | |
1222  -
1223  0
1224 
1225 
1226 The integrals are evaluated by a power series for x < 1.
1227 For x >= 1 auxiliary functions f(x) and g(x) are employed
1228 such that
1229 
1230 C(x) = 0.5 + f(x) sin( pi/2 x**2 ) - g(x) cos( pi/2 x**2 )
1231 S(x) = 0.5 - f(x) cos( pi/2 x**2 ) - g(x) sin( pi/2 x**2 )
1232 
1233 
1234 
1235 ACCURACY:
1236 
1237  Relative error.
1238 
1239 Arithmetic function domain # trials peak rms
1240  IEEE S(x) 0, 10 10000 2.0e-15 3.2e-16
1241  IEEE C(x) 0, 10 10000 1.8e-15 3.3e-16
1242 
1243 Cephes Math Library Release 2.8: June, 2000
1244 Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
1245 *************************************************************************/
1246 void fresnelintegral(const double x, double &c, double &s);
1247 
1248 /*************************************************************************
1249 Jacobian Elliptic Functions
1250 
1251 Evaluates the Jacobian elliptic functions sn(u|m), cn(u|m),
1252 and dn(u|m) of parameter m between 0 and 1, and real
1253 argument u.
1254 
1255 These functions are periodic, with quarter-period on the
1256 real axis equal to the complete elliptic integral
1257 ellpk(1.0-m).
1258 
1259 Relation to incomplete elliptic integral:
1260 If u = ellik(phi,m), then sn(u|m) = sin(phi),
1261 and cn(u|m) = cos(phi). Phi is called the amplitude of u.
1262 
1263 Computation is by means of the arithmetic-geometric mean
1264 algorithm, except when m is within 1e-9 of 0 or 1. In the
1265 latter case with m close to 1, the approximation applies
1266 only for phi < pi/2.
1267 
1268 ACCURACY:
1269 
1270 Tested at random points with u between 0 and 10, m between
1271 0 and 1.
1272 
1273  Absolute error (* = relative error):
1274 arithmetic function # trials peak rms
1275  IEEE phi 10000 9.2e-16* 1.4e-16*
1276  IEEE sn 50000 4.1e-15 4.6e-16
1277  IEEE cn 40000 3.6e-15 4.4e-16
1278  IEEE dn 10000 1.3e-12 1.8e-14
1279 
1280  Peak error observed in consistency check using addition
1281 theorem for sn(u+v) was 4e-16 (absolute). Also tested by
1282 the above relation to the incomplete elliptic integral.
1283 Accuracy deteriorates when u is large.
1284 
1285 Cephes Math Library Release 2.8: June, 2000
1286 Copyright 1984, 1987, 2000 by Stephen L. Moshier
1287 *************************************************************************/
1288 void jacobianellipticfunctions(const double u, const double m, double &sn, double &cn, double &dn, double &ph);
1289 
1290 /*************************************************************************
1291 Psi (digamma) function
1292 
1293  d -
1294  psi(x) = -- ln | (x)
1295  dx
1296 
1297 is the logarithmic derivative of the gamma function.
1298 For integer x,
1299  n-1
1300  -
1301 psi(n) = -EUL + > 1/k.
1302  -
1303  k=1
1304 
1305 This formula is used for 0 < n <= 10. If x is negative, it
1306 is transformed to a positive argument by the reflection
1307 formula psi(1-x) = psi(x) + pi cot(pi x).
1308 For general positive x, the argument is made greater than 10
1309 using the recurrence psi(x+1) = psi(x) + 1/x.
1310 Then the following asymptotic expansion is applied:
1311 
1312  inf. B
1313  - 2k
1314 psi(x) = log(x) - 1/2x - > -------
1315  - 2k
1316  k=1 2k x
1317 
1318 where the B2k are Bernoulli numbers.
1319 
1320 ACCURACY:
1321  Relative error (except absolute when |psi| < 1):
1322 arithmetic domain # trials peak rms
1323  IEEE 0,30 30000 1.3e-15 1.4e-16
1324  IEEE -30,0 40000 1.5e-15 2.2e-16
1325 
1326 Cephes Math Library Release 2.8: June, 2000
1327 Copyright 1984, 1987, 1992, 2000 by Stephen L. Moshier
1328 *************************************************************************/
1329 double psi(const double x);
1330 
1331 /*************************************************************************
1332 Exponential integral Ei(x)
1333 
1334  x
1335  - t
1336  | | e
1337  Ei(x) = -|- --- dt .
1338  | | t
1339  -
1340  -inf
1341 
1342 Not defined for x <= 0.
1343 See also expn.c.
1344 
1345 
1346 
1347 ACCURACY:
1348 
1349  Relative error:
1350 arithmetic domain # trials peak rms
1351  IEEE 0,100 50000 8.6e-16 1.3e-16
1352 
1353 Cephes Math Library Release 2.8: May, 1999
1354 Copyright 1999 by Stephen L. Moshier
1355 *************************************************************************/
1356 double exponentialintegralei(const double x);
1357 
1358 
1359 /*************************************************************************
1360 Exponential integral En(x)
1361 
1362 Evaluates the exponential integral
1363 
1364  inf.
1365  -
1366  | | -xt
1367  | e
1368  E (x) = | ---- dt.
1369  n | n
1370  | | t
1371  -
1372  1
1373 
1374 
1375 Both n and x must be nonnegative.
1376 
1377 The routine employs either a power series, a continued
1378 fraction, or an asymptotic formula depending on the
1379 relative values of n and x.
1380 
1381 ACCURACY:
1382 
1383  Relative error:
1384 arithmetic domain # trials peak rms
1385  IEEE 0, 30 10000 1.7e-15 3.6e-16
1386 
1387 Cephes Math Library Release 2.8: June, 2000
1388 Copyright 1985, 2000 by Stephen L. Moshier
1389 *************************************************************************/
1390 double exponentialintegralen(const double x, const ae_int_t n);
1391 
1392 /*************************************************************************
1393 Calculation of the value of the Laguerre polynomial.
1394 
1395 Parameters:
1396  n - degree, n>=0
1397  x - argument
1398 
1399 Result:
1400  the value of the Laguerre polynomial Ln at x
1401 *************************************************************************/
1402 double laguerrecalculate(const ae_int_t n, const double x);
1403 
1404 
1405 /*************************************************************************
1406 Summation of Laguerre polynomials using Clenshaw's recurrence formula.
1407 
1408 This routine calculates c[0]*L0(x) + c[1]*L1(x) + ... + c[N]*LN(x)
1409 
1410 Parameters:
1411  n - degree, n>=0
1412  x - argument
1413 
1414 Result:
1415  the value of the Laguerre polynomial at x
1416 *************************************************************************/
1417 double laguerresum(const real_1d_array &c, const ae_int_t n, const double x);
1418 
1419 
1420 /*************************************************************************
1421 Representation of Ln as C[0] + C[1]*X + ... + C[N]*X^N
1422 
1423 Input parameters:
1424  N - polynomial degree, n>=0
1425 
1426 Output parameters:
1427  C - coefficients
1428 *************************************************************************/
1430 
1431 /*************************************************************************
1432 Chi-square distribution
1433 
1434 Returns the area under the left hand tail (from 0 to x)
1435 of the Chi square probability density function with
1436 v degrees of freedom.
1437 
1438 
1439  x
1440  -
1441  1 | | v/2-1 -t/2
1442  P( x | v ) = ----------- | t e dt
1443  v/2 - | |
1444  2 | (v/2) -
1445  0
1446 
1447 where x is the Chi-square variable.
1448 
1449 The incomplete gamma integral is used, according to the
1450 formula
1451 
1452 y = chdtr( v, x ) = igam( v/2.0, x/2.0 ).
1453 
1454 The arguments must both be positive.
1455 
1456 ACCURACY:
1457 
1458 See incomplete gamma function
1459 
1460 
1461 Cephes Math Library Release 2.8: June, 2000
1462 Copyright 1984, 1987, 2000 by Stephen L. Moshier
1463 *************************************************************************/
1464 double chisquaredistribution(const double v, const double x);
1465 
1466 
1467 /*************************************************************************
1468 Complemented Chi-square distribution
1469 
1470 Returns the area under the right hand tail (from x to
1471 infinity) of the Chi square probability density function
1472 with v degrees of freedom:
1473 
1474  inf.
1475  -
1476  1 | | v/2-1 -t/2
1477  P( x | v ) = ----------- | t e dt
1478  v/2 - | |
1479  2 | (v/2) -
1480  x
1481 
1482 where x is the Chi-square variable.
1483 
1484 The incomplete gamma integral is used, according to the
1485 formula
1486 
1487 y = chdtr( v, x ) = igamc( v/2.0, x/2.0 ).
1488 
1489 The arguments must both be positive.
1490 
1491 ACCURACY:
1492 
1493 See incomplete gamma function
1494 
1495 Cephes Math Library Release 2.8: June, 2000
1496 Copyright 1984, 1987, 2000 by Stephen L. Moshier
1497 *************************************************************************/
1498 double chisquarecdistribution(const double v, const double x);
1499 
1500 
1501 /*************************************************************************
1502 Inverse of complemented Chi-square distribution
1503 
1504 Finds the Chi-square argument x such that the integral
1505 from x to infinity of the Chi-square density is equal
1506 to the given cumulative probability y.
1507 
1508 This is accomplished using the inverse gamma integral
1509 function and the relation
1510 
1511  x/2 = igami( df/2, y );
1512 
1513 ACCURACY:
1514 
1515 See inverse incomplete gamma function
1516 
1517 
1518 Cephes Math Library Release 2.8: June, 2000
1519 Copyright 1984, 1987, 2000 by Stephen L. Moshier
1520 *************************************************************************/
1521 double invchisquaredistribution(const double v, const double y);
1522 
1523 /*************************************************************************
1524 Calculation of the value of the Legendre polynomial Pn.
1525 
1526 Parameters:
1527  n - degree, n>=0
1528  x - argument
1529 
1530 Result:
1531  the value of the Legendre polynomial Pn at x
1532 *************************************************************************/
1533 double legendrecalculate(const ae_int_t n, const double x);
1534 
1535 
1536 /*************************************************************************
1537 Summation of Legendre polynomials using Clenshaw's recurrence formula.
1538 
1539 This routine calculates
1540  c[0]*P0(x) + c[1]*P1(x) + ... + c[N]*PN(x)
1541 
1542 Parameters:
1543  n - degree, n>=0
1544  x - argument
1545 
1546 Result:
1547  the value of the Legendre polynomial at x
1548 *************************************************************************/
1549 double legendresum(const real_1d_array &c, const ae_int_t n, const double x);
1550 
1551 
1552 /*************************************************************************
1553 Representation of Pn as C[0] + C[1]*X + ... + C[N]*X^N
1554 
1555 Input parameters:
1556  N - polynomial degree, n>=0
1557 
1558 Output parameters:
1559  C - coefficients
1560 *************************************************************************/
1562 
1563 /*************************************************************************
1564 Beta function
1565 
1566 
1567  - -
1568  | (a) | (b)
1569 beta( a, b ) = -----------.
1570  -
1571  | (a+b)
1572 
1573 For large arguments the logarithm of the function is
1574 evaluated using lgam(), then exponentiated.
1575 
1576 ACCURACY:
1577 
1578  Relative error:
1579 arithmetic domain # trials peak rms
1580  IEEE 0,30 30000 8.1e-14 1.1e-14
1581 
1582 Cephes Math Library Release 2.0: April, 1987
1583 Copyright 1984, 1987 by Stephen L. Moshier
1584 *************************************************************************/
1585 double beta(const double a, const double b);
1586 
1587 /*************************************************************************
1588 Calculation of the value of the Chebyshev polynomials of the
1589 first and second kinds.
1590 
1591 Parameters:
1592  r - polynomial kind, either 1 or 2.
1593  n - degree, n>=0
1594  x - argument, -1 <= x <= 1
1595 
1596 Result:
1597  the value of the Chebyshev polynomial at x
1598 *************************************************************************/
1599 double chebyshevcalculate(const ae_int_t r, const ae_int_t n, const double x);
1600 
1601 
1602 /*************************************************************************
1603 Summation of Chebyshev polynomials using Clenshaw's recurrence formula.
1604 
1605 This routine calculates
1606  c[0]*T0(x) + c[1]*T1(x) + ... + c[N]*TN(x)
1607 or
1608  c[0]*U0(x) + c[1]*U1(x) + ... + c[N]*UN(x)
1609 depending on the R.
1610 
1611 Parameters:
1612  r - polynomial kind, either 1 or 2.
1613  n - degree, n>=0
1614  x - argument
1615 
1616 Result:
1617  the value of the Chebyshev polynomial at x
1618 *************************************************************************/
1619 double chebyshevsum(const real_1d_array &c, const ae_int_t r, const ae_int_t n, const double x);
1620 
1621 
1622 /*************************************************************************
1623 Representation of Tn as C[0] + C[1]*X + ... + C[N]*X^N
1624 
1625 Input parameters:
1626  N - polynomial degree, n>=0
1627 
1628 Output parameters:
1629  C - coefficients
1630 *************************************************************************/
1632 
1633 
1634 /*************************************************************************
1635 Conversion of a series of Chebyshev polynomials to a power series.
1636 
1637 Represents A[0]*T0(x) + A[1]*T1(x) + ... + A[N]*Tn(x) as
1638 B[0] + B[1]*X + ... + B[N]*X^N.
1639 
1640 Input parameters:
1641  A - Chebyshev series coefficients
1642  N - degree, N>=0
1643 
1644 Output parameters
1645  B - power series coefficients
1646 *************************************************************************/
1648 
1649 /*************************************************************************
1650 Student's t distribution
1651 
1652 Computes the integral from minus infinity to t of the Student
1653 t distribution with integer k > 0 degrees of freedom:
1654 
1655  t
1656  -
1657  | |
1658  - | 2 -(k+1)/2
1659  | ( (k+1)/2 ) | ( x )
1660  ---------------------- | ( 1 + --- ) dx
1661  - | ( k )
1662  sqrt( k pi ) | ( k/2 ) |
1663  | |
1664  -
1665  -inf.
1666 
1667 Relation to incomplete beta integral:
1668 
1669  1 - stdtr(k,t) = 0.5 * incbet( k/2, 1/2, z )
1670 where
1671  z = k/(k + t**2).
1672 
1673 For t < -2, this is the method of computation. For higher t,
1674 a direct method is derived from integration by parts.
1675 Since the function is symmetric about t=0, the area under the
1676 right tail of the density is found by calling the function
1677 with -t instead of t.
1678 
1679 ACCURACY:
1680 
1681 Tested at random 1 <= k <= 25. The "domain" refers to t.
1682  Relative error:
1683 arithmetic domain # trials peak rms
1684  IEEE -100,-2 50000 5.9e-15 1.4e-15
1685  IEEE -2,100 500000 2.7e-15 4.9e-17
1686 
1687 Cephes Math Library Release 2.8: June, 2000
1688 Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier
1689 *************************************************************************/
1690 double studenttdistribution(const ae_int_t k, const double t);
1691 
1692 
1693 /*************************************************************************
1694 Functional inverse of Student's t distribution
1695 
1696 Given probability p, finds the argument t such that stdtr(k,t)
1697 is equal to p.
1698 
1699 ACCURACY:
1700 
1701 Tested at random 1 <= k <= 100. The "domain" refers to p:
1702  Relative error:
1703 arithmetic domain # trials peak rms
1704  IEEE .001,.999 25000 5.7e-15 8.0e-16
1705  IEEE 10^-6,.001 25000 2.0e-12 2.9e-14
1706 
1707 Cephes Math Library Release 2.8: June, 2000
1708 Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier
1709 *************************************************************************/
1710 double invstudenttdistribution(const ae_int_t k, const double p);
1711 
1712 /*************************************************************************
1713 Binomial distribution
1714 
1715 Returns the sum of the terms 0 through k of the Binomial
1716 probability density:
1717 
1718  k
1719  -- ( n ) j n-j
1720  > ( ) p (1-p)
1721  -- ( j )
1722  j=0
1723 
1724 The terms are not summed directly; instead the incomplete
1725 beta integral is employed, according to the formula
1726 
1727 y = bdtr( k, n, p ) = incbet( n-k, k+1, 1-p ).
1728 
1729 The arguments must be positive, with p ranging from 0 to 1.
1730 
1731 ACCURACY:
1732 
1733 Tested at random points (a,b,p), with p between 0 and 1.
1734 
1735  a,b Relative error:
1736 arithmetic domain # trials peak rms
1737  For p between 0.001 and 1:
1738  IEEE 0,100 100000 4.3e-15 2.6e-16
1739 
1740 Cephes Math Library Release 2.8: June, 2000
1741 Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier
1742 *************************************************************************/
1743 double binomialdistribution(const ae_int_t k, const ae_int_t n, const double p);
1744 
1745 
1746 /*************************************************************************
1747 Complemented binomial distribution
1748 
1749 Returns the sum of the terms k+1 through n of the Binomial
1750 probability density:
1751 
1752  n
1753  -- ( n ) j n-j
1754  > ( ) p (1-p)
1755  -- ( j )
1756  j=k+1
1757 
1758 The terms are not summed directly; instead the incomplete
1759 beta integral is employed, according to the formula
1760 
1761 y = bdtrc( k, n, p ) = incbet( k+1, n-k, p ).
1762 
1763 The arguments must be positive, with p ranging from 0 to 1.
1764 
1765 ACCURACY:
1766 
1767 Tested at random points (a,b,p).
1768 
1769  a,b Relative error:
1770 arithmetic domain # trials peak rms
1771  For p between 0.001 and 1:
1772  IEEE 0,100 100000 6.7e-15 8.2e-16
1773  For p between 0 and .001:
1774  IEEE 0,100 100000 1.5e-13 2.7e-15
1775 
1776 Cephes Math Library Release 2.8: June, 2000
1777 Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier
1778 *************************************************************************/
1779 double binomialcdistribution(const ae_int_t k, const ae_int_t n, const double p);
1780 
1781 
1782 /*************************************************************************
1783 Inverse binomial distribution
1784 
1785 Finds the event probability p such that the sum of the
1786 terms 0 through k of the Binomial probability density
1787 is equal to the given cumulative probability y.
1788 
1789 This is accomplished using the inverse beta integral
1790 function and the relation
1791 
1792 1 - p = incbi( n-k, k+1, y ).
1793 
1794 ACCURACY:
1795 
1796 Tested at random points (a,b,p).
1797 
1798  a,b Relative error:
1799 arithmetic domain # trials peak rms
1800  For p between 0.001 and 1:
1801  IEEE 0,100 100000 2.3e-14 6.4e-16
1802  IEEE 0,10000 100000 6.6e-12 1.2e-13
1803  For p between 10^-6 and 0.001:
1804  IEEE 0,100 100000 2.0e-12 1.3e-14
1805  IEEE 0,10000 100000 1.5e-12 3.2e-14
1806 
1807 Cephes Math Library Release 2.8: June, 2000
1808 Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier
1809 *************************************************************************/
1810 double invbinomialdistribution(const ae_int_t k, const ae_int_t n, const double y);
1811 
1812 /*************************************************************************
1813 Airy function
1814 
1815 Solution of the differential equation
1816 
1817 y"(x) = xy.
1818 
1819 The function returns the two independent solutions Ai, Bi
1820 and their first derivatives Ai'(x), Bi'(x).
1821 
1822 Evaluation is by power series summation for small x,
1823 by rational minimax approximations for large x.
1824 
1825 
1826 
1827 ACCURACY:
1828 Error criterion is absolute when function <= 1, relative
1829 when function > 1, except * denotes relative error criterion.
1830 For large negative x, the absolute error increases as x^1.5.
1831 For large positive x, the relative error increases as x^1.5.
1832 
1833 Arithmetic domain function # trials peak rms
1834 IEEE -10, 0 Ai 10000 1.6e-15 2.7e-16
1835 IEEE 0, 10 Ai 10000 2.3e-14* 1.8e-15*
1836 IEEE -10, 0 Ai' 10000 4.6e-15 7.6e-16
1837 IEEE 0, 10 Ai' 10000 1.8e-14* 1.5e-15*
1838 IEEE -10, 10 Bi 30000 4.2e-15 5.3e-16
1839 IEEE -10, 10 Bi' 30000 4.9e-15 7.3e-16
1840 
1841 Cephes Math Library Release 2.8: June, 2000
1842 Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
1843 *************************************************************************/
1844 void airy(const double x, double &ai, double &aip, double &bi, double &bip);
1845 }
1846 
1848 //
1849 // THIS SECTION CONTAINS COMPUTATIONAL CORE DECLARATIONS (FUNCTIONS)
1850 //
1852 namespace alglib_impl
1853 {
1854 double gammafunction(double x, ae_state *_state);
1855 double lngamma(double x, double* sgngam, ae_state *_state);
1856 double errorfunction(double x, ae_state *_state);
1857 double errorfunctionc(double x, ae_state *_state);
1858 double normaldistribution(double x, ae_state *_state);
1859 double inverf(double e, ae_state *_state);
1860 double invnormaldistribution(double y0, ae_state *_state);
1861 double incompletegamma(double a, double x, ae_state *_state);
1862 double incompletegammac(double a, double x, ae_state *_state);
1863 double invincompletegammac(double a, double y0, ae_state *_state);
1864 double ellipticintegralk(double m, ae_state *_state);
1865 double ellipticintegralkhighprecision(double m1, ae_state *_state);
1866 double incompleteellipticintegralk(double phi, double m, ae_state *_state);
1867 double ellipticintegrale(double m, ae_state *_state);
1868 double incompleteellipticintegrale(double phi, double m, ae_state *_state);
1869 double hermitecalculate(ae_int_t n, double x, ae_state *_state);
1870 double hermitesum(/* Real */ ae_vector* c,
1871  ae_int_t n,
1872  double x,
1873  ae_state *_state);
1875  /* Real */ ae_vector* c,
1876  ae_state *_state);
1877 double dawsonintegral(double x, ae_state *_state);
1878 void sinecosineintegrals(double x,
1879  double* si,
1880  double* ci,
1881  ae_state *_state);
1883  double* shi,
1884  double* chi,
1885  ae_state *_state);
1886 double poissondistribution(ae_int_t k, double m, ae_state *_state);
1887 double poissoncdistribution(ae_int_t k, double m, ae_state *_state);
1888 double invpoissondistribution(ae_int_t k, double y, ae_state *_state);
1889 double besselj0(double x, ae_state *_state);
1890 double besselj1(double x, ae_state *_state);
1891 double besseljn(ae_int_t n, double x, ae_state *_state);
1892 double bessely0(double x, ae_state *_state);
1893 double bessely1(double x, ae_state *_state);
1894 double besselyn(ae_int_t n, double x, ae_state *_state);
1895 double besseli0(double x, ae_state *_state);
1896 double besseli1(double x, ae_state *_state);
1897 double besselk0(double x, ae_state *_state);
1898 double besselk1(double x, ae_state *_state);
1899 double besselkn(ae_int_t nn, double x, ae_state *_state);
1900 double incompletebeta(double a, double b, double x, ae_state *_state);
1901 double invincompletebeta(double a, double b, double y, ae_state *_state);
1902 double fdistribution(ae_int_t a, ae_int_t b, double x, ae_state *_state);
1903 double fcdistribution(ae_int_t a, ae_int_t b, double x, ae_state *_state);
1905  ae_int_t b,
1906  double y,
1907  ae_state *_state);
1908 void fresnelintegral(double x, double* c, double* s, ae_state *_state);
1910  double m,
1911  double* sn,
1912  double* cn,
1913  double* dn,
1914  double* ph,
1915  ae_state *_state);
1916 double psi(double x, ae_state *_state);
1917 double exponentialintegralei(double x, ae_state *_state);
1918 double exponentialintegralen(double x, ae_int_t n, ae_state *_state);
1919 double laguerrecalculate(ae_int_t n, double x, ae_state *_state);
1920 double laguerresum(/* Real */ ae_vector* c,
1921  ae_int_t n,
1922  double x,
1923  ae_state *_state);
1925  /* Real */ ae_vector* c,
1926  ae_state *_state);
1927 double chisquaredistribution(double v, double x, ae_state *_state);
1928 double chisquarecdistribution(double v, double x, ae_state *_state);
1929 double invchisquaredistribution(double v, double y, ae_state *_state);
1930 double legendrecalculate(ae_int_t n, double x, ae_state *_state);
1931 double legendresum(/* Real */ ae_vector* c,
1932  ae_int_t n,
1933  double x,
1934  ae_state *_state);
1936  /* Real */ ae_vector* c,
1937  ae_state *_state);
1938 double beta(double a, double b, ae_state *_state);
1940  ae_int_t n,
1941  double x,
1942  ae_state *_state);
1943 double chebyshevsum(/* Real */ ae_vector* c,
1944  ae_int_t r,
1945  ae_int_t n,
1946  double x,
1947  ae_state *_state);
1949  /* Real */ ae_vector* c,
1950  ae_state *_state);
1951 void fromchebyshev(/* Real */ ae_vector* a,
1952  ae_int_t n,
1953  /* Real */ ae_vector* b,
1954  ae_state *_state);
1955 double studenttdistribution(ae_int_t k, double t, ae_state *_state);
1956 double invstudenttdistribution(ae_int_t k, double p, ae_state *_state);
1958  ae_int_t n,
1959  double p,
1960  ae_state *_state);
1962  ae_int_t n,
1963  double p,
1964  ae_state *_state);
1966  ae_int_t n,
1967  double y,
1968  ae_state *_state);
1969 void airy(double x,
1970  double* ai,
1971  double* aip,
1972  double* bi,
1973  double* bip,
1974  ae_state *_state);
1975 
1976 }
1977 #endif
1978 
double exponentialintegralen(const double x, const ae_int_t n)
double invpoissondistribution(const ae_int_t k, const double y)
double invstudenttdistribution(const ae_int_t k, const double p)
double gammafunction(const double x)
double laguerrecalculate(const ae_int_t n, const double x)
double besselk1(double x, ae_state *_state)
void fresnelintegral(double x, double *c, double *s, ae_state *_state)
double poissoncdistribution(ae_int_t k, double m, ae_state *_state)
double besselyn(ae_int_t n, double x, ae_state *_state)
void jacobianellipticfunctions(double u, double m, double *sn, double *cn, double *dn, double *ph, ae_state *_state)
double laguerresum(const real_1d_array &c, const ae_int_t n, const double x)
double hermitecalculate(ae_int_t n, double x, ae_state *_state)
double exponentialintegralen(double x, ae_int_t n, ae_state *_state)
double binomialdistribution(const ae_int_t k, const ae_int_t n, const double p)
double hermitesum(const real_1d_array &c, const ae_int_t n, const double x)
double invincompletegammac(double a, double y0, ae_state *_state)
double psi(double x, ae_state *_state)
double invpoissondistribution(ae_int_t k, double y, ae_state *_state)
double chebyshevcalculate(const ae_int_t r, const ae_int_t n, const double x)
double invfdistribution(const ae_int_t a, const ae_int_t b, const double y)
double inverf(double e, ae_state *_state)
double besseli1(const double x)
double exponentialintegralei(const double x)
double chisquaredistribution(double v, double x, ae_state *_state)
double incompleteellipticintegralk(const double phi, const double m)
double errorfunctionc(const double x)
double incompleteellipticintegralk(double phi, double m, ae_state *_state)
double gammafunction(double x, ae_state *_state)
void fresnelintegral(const double x, double &c, double &s)
double besseli0(double x, ae_state *_state)
void chebyshevcoefficients(ae_int_t n, ae_vector *c, ae_state *_state)
double fcdistribution(ae_int_t a, ae_int_t b, double x, ae_state *_state)
void fromchebyshev(const real_1d_array &a, const ae_int_t n, real_1d_array &b)
double lngamma(const double x, double &sgngam)
double chisquarecdistribution(double v, double x, ae_state *_state)
double besseli0(const double x)
double invincompletebeta(const double a, const double b, const double y)
void laguerrecoefficients(const ae_int_t n, real_1d_array &c)
void airy(double x, double *ai, double *aip, double *bi, double *bip, ae_state *_state)
double ellipticintegrale(const double m)
double poissondistribution(const ae_int_t k, const double m)
double besselj0(double x, ae_state *_state)
void jacobianellipticfunctions(const double u, const double m, double &sn, double &cn, double &dn, double &ph)
double binomialcdistribution(const ae_int_t k, const ae_int_t n, const double p)
double exponentialintegralei(double x, ae_state *_state)
double inverf(const double e)
double psi(const double x)
double errorfunction(double x, ae_state *_state)
double incompleteellipticintegrale(double phi, double m, ae_state *_state)
double binomialdistribution(ae_int_t k, ae_int_t n, double p, ae_state *_state)
double chisquarecdistribution(const double v, const double x)
double fcdistribution(const ae_int_t a, const ae_int_t b, const double x)
double studenttdistribution(ae_int_t k, double t, ae_state *_state)
double fdistribution(ae_int_t a, ae_int_t b, double x, ae_state *_state)
void fromchebyshev(ae_vector *a, ae_int_t n, ae_vector *b, ae_state *_state)
double normaldistribution(const double x)
double dawsonintegral(double x, ae_state *_state)
double chebyshevsum(const real_1d_array &c, const ae_int_t r, const ae_int_t n, const double x)
double besselk1(const double x)
double incompletebeta(double a, double b, double x, ae_state *_state)
void sinecosineintegrals(const double x, double &si, double &ci)
double hermitecalculate(const ae_int_t n, const double x)
void airy(const double x, double &ai, double &aip, double &bi, double &bip)
double legendrecalculate(const ae_int_t n, const double x)
double invincompletegammac(const double a, const double y0)
double besselkn(const ae_int_t nn, const double x)
double laguerrecalculate(ae_int_t n, double x, ae_state *_state)
double normaldistribution(double x, ae_state *_state)
double hermitesum(ae_vector *c, ae_int_t n, double x, ae_state *_state)
double chebyshevcalculate(ae_int_t r, ae_int_t n, double x, ae_state *_state)
double besselj1(const double x)
double invincompletebeta(double a, double b, double y, ae_state *_state)
double ellipticintegralk(const double m)
double invchisquaredistribution(const double v, const double y)
double incompletegammac(const double a, const double x)
void legendrecoefficients(const ae_int_t n, real_1d_array &c)
alglib_impl::ae_int_t ae_int_t
Definition: ap.h:965
double chisquaredistribution(const double v, const double x)
double invnormaldistribution(const double y0)
void hyperbolicsinecosineintegrals(double x, double *shi, double *chi, ae_state *_state)
double bessely0(const double x)
double dawsonintegral(const double x)
double invnormaldistribution(double y0, ae_state *_state)
double incompletegammac(double a, double x, ae_state *_state)
double besseljn(const ae_int_t n, const double x)
double beta(const double a, const double b)
void hermitecoefficients(const ae_int_t n, real_1d_array &c)
double invstudenttdistribution(ae_int_t k, double p, ae_state *_state)
void chebyshevcoefficients(const ae_int_t n, real_1d_array &c)
double invbinomialdistribution(const ae_int_t k, const ae_int_t n, const double y)
double errorfunction(const double x)
void sinecosineintegrals(double x, double *si, double *ci, ae_state *_state)
double studenttdistribution(const ae_int_t k, const double t)
double bessely0(double x, ae_state *_state)
double poissondistribution(ae_int_t k, double m, ae_state *_state)
void laguerrecoefficients(ae_int_t n, ae_vector *c, ae_state *_state)
double beta(double a, double b, ae_state *_state)
double errorfunctionc(double x, ae_state *_state)
void legendrecoefficients(ae_int_t n, ae_vector *c, ae_state *_state)
double fdistribution(const ae_int_t a, const ae_int_t b, const double x)
double besselyn(const ae_int_t n, const double x)
double bessely1(const double x)
double besselj1(double x, ae_state *_state)
double binomialcdistribution(ae_int_t k, ae_int_t n, double p, ae_state *_state)
double incompletegamma(double a, double x, ae_state *_state)
double ellipticintegrale(double m, ae_state *_state)
double besselk0(double x, ae_state *_state)
double chebyshevsum(ae_vector *c, ae_int_t r, ae_int_t n, double x, ae_state *_state)
double legendrecalculate(ae_int_t n, double x, ae_state *_state)
double incompletegamma(const double a, const double x)
double invchisquaredistribution(double v, double y, ae_state *_state)
double besseljn(ae_int_t n, double x, ae_state *_state)
double incompleteellipticintegrale(const double phi, const double m)
double besselj0(const double x)
double besselk0(const double x)
double incompletebeta(const double a, const double b, const double x)
void hyperbolicsinecosineintegrals(const double x, double &shi, double &chi)
double ellipticintegralkhighprecision(const double m1)
double invbinomialdistribution(ae_int_t k, ae_int_t n, double y, ae_state *_state)
double ellipticintegralk(double m, ae_state *_state)
void hermitecoefficients(ae_int_t n, ae_vector *c, ae_state *_state)
ptrdiff_t ae_int_t
Definition: ap.h:185
double invfdistribution(ae_int_t a, ae_int_t b, double y, ae_state *_state)
double ellipticintegralkhighprecision(double m1, ae_state *_state)
double besselkn(ae_int_t nn, double x, ae_state *_state)
double poissoncdistribution(const ae_int_t k, const double m)
double lngamma(double x, double *sgngam, ae_state *_state)
double laguerresum(ae_vector *c, ae_int_t n, double x, ae_state *_state)
double bessely1(double x, ae_state *_state)
double legendresum(const real_1d_array &c, const ae_int_t n, const double x)
double besseli1(double x, ae_state *_state)
double legendresum(ae_vector *c, ae_int_t n, double x, ae_state *_state)
Page URL: http://wiki.math.ethz.ch/bin/view/Concepts/WebHome
21 August 2020
© 2020 Eidgenössische Technische Hochschule Zürich