27 #if defined (HAVE_CONFIG_H)
63 #if defined (HAVE_ACOSH)
75 #if defined (HAVE_ACOSHF)
87 #if defined (HAVE_COMPLEX_STD_ACOSH)
90 return log (x + sqrt (x + 1.0) * sqrt (x - 1.0));
97 #if defined (HAVE_COMPLEX_STD_ACOSH)
100 return log (x + sqrt (x + 1.0
f) * sqrt (x - 1.0
f));
107 #if defined (HAVE_ASINH)
119 #if defined (HAVE_ASINHF)
131 #if defined (HAVE_COMPLEX_STD_ASINH)
134 return log (x + sqrt (x*x + 1.0));
141 #if defined (HAVE_COMPLEX_STD_ASINH)
144 return log (x + sqrt (x*x + 1.0
f));
151 #if defined (HAVE_ATANH)
163 #if defined (HAVE_ATANHF)
175 #if defined (HAVE_COMPLEX_STD_ATANH)
178 return log ((1.0 + x) / (1.0 - x)) / 2.0;
185 #if defined (HAVE_COMPLEX_STD_ATANH)
188 return log ((1.0
f + x) / (1.0
f - x)) / 2.0f;
195 #if defined (HAVE_ERF)
207 #if defined (HAVE_ERFF)
226 Complex xd (x.real (), x.imag ());
234 #if defined (HAVE_ERFC)
246 #if defined (HAVE_ERFCF)
265 Complex xd (x.real (), x.imag ());
283 Complex xd (x.real (), x.imag ());
301 Complex xd (x.real (), x.imag ());
319 Complex xd (x.real (), x.imag ());
342 #if defined (HAVE_TGAMMA)
355 #if defined (HAVE_LGAMMA)
377 #if defined (HAVE_LGAMMA_R)
379 result = lgamma_r (x, &sgngam);
393 return result +
Complex (0., M_PI);
416 #if defined (HAVE_TGAMMA)
417 result = tgammaf (x);
429 #if defined (HAVE_LGAMMAF)
451 #if defined (HAVE_LGAMMAF_R)
453 result = lgammaf_r (x, &sgngam);
475 #if defined (HAVE_EXPM1)
480 double ax = fabs (x);
489 for (
int i = 2;
i < 7;
i++)
495 for (
int i = 0;
i < 4;
i++)
501 retval = (x > 0) ? s : -s / (1+s);
504 retval = exp (x) - 1;
517 double im = x.
imag ();
518 double u =
expm1 (x.real ());
519 double v = sin (im/2);
521 retval =
Complex (u*v + u + v, (u+1) * sin (im));
524 retval = std::exp (x) -
Complex (1);
532 #if defined (HAVE_EXPM1F)
546 for (
int i = 2;
i < 7;
i++)
552 for (
int i = 0;
i < 4;
i++)
558 retval = (x > 0) ? s : -s / (1+s);
561 retval = exp (x) - 1;
574 float im = x.
imag ();
575 float u =
expm1 (x.real ());
576 float v = sin (im/2);
589 #if defined (HAVE_LOG1P)
594 double ax = fabs (x);
599 double u = x / (2 +
x),
t = 1,
s = 0;
600 for (
int i = 2;
i < 12;
i += 2)
601 s += (
t *= u*u) / (
i+1);
603 retval = 2 * (
s + 1) * u;
619 if (fabs (r) < 0.5 && fabs (
i) < 0.5)
621 double u = 2*r + r*r +
i*
i;
631 template <
typename T>
635 static const T one_third = 0.3333333333333333333f;
641 return (x / (y*y) + y + y) / 3;
650 #if defined (HAVE_CBRT)
660 #if defined (HAVE_LOG1PF)
670 float u = x / (2 +
x),
t = 1.0
f,
s = 0;
671 for (
int i = 2;
i < 12;
i += 2)
672 s += (
t *= u*u) / (
i+1);
674 retval = 2 * (
s + 1.0f) * u;
690 if (fabs (r) < 0.5 && fabs (
i) < 0.5)
692 float u = 2*r + r*r +
i*
i;
705 #if defined (HAVE_CBRTF)
765 return x ==
static_cast<double> (
static_cast<long> (
x));
780 double zr = z.real ();
781 double zi = z.imag ();
792 if (zi == 0.0 && zr >= 0.0)
802 if ((static_cast<long> (alpha)) & 1)
812 if (ierr == 0 || ierr == 3)
814 tmp -= sin (M_PI * alpha) *
zbesy (z, alpha, kode, ierr);
840 double zr = z.real ();
841 double zi = z.imag ();
845 if (zr == 0.0 && zi == 0.0)
862 if (zi == 0.0 && zr >= 0.0)
873 if ((static_cast<long> (alpha - 0.5)) & 1)
883 if (ierr == 0 || ierr == 3)
885 tmp += sin (M_PI * alpha) *
zbesj (z, alpha, kode, ierr);
909 double zr = z.real ();
910 double zi = z.imag ();
921 if (zi == 0.0 && zr >= 0.0)
939 if (ierr == 0 || ierr == 3)
941 Complex tmp2 = (2.0 / M_PI) * sin (M_PI * alpha)
942 *
zbesk (z, alpha, kode, ierr);
947 tmp2 *= exp (-z -
std::abs (z.real ()));
974 double zr = z.real ();
975 double zi = z.imag ();
979 if (zr == 0.0 && zi == 0.0)
992 double rexpz = expz.real ();
993 double iexpz = expz.imag ();
995 double tmp = yr*rexpz - yi*iexpz;
997 yi = yr*iexpz + yi*rexpz;
1001 if (zi == 0.0 && zr >= 0.0)
1029 double zr = z.real ();
1030 double zi = z.imag ();
1032 F77_FUNC (
zbesh, ZBESH) (zr, zi, alpha, 2, 1, 1, &yr, &yi, nz,
ierr);
1038 double rexpz = expz.real ();
1039 double iexpz = expz.imag ();
1041 double tmp = yr*rexpz - yi*iexpz;
1043 yi = yr*iexpz + yi*rexpz;
1075 double zr = z.real ();
1076 double zi = z.imag ();
1078 F77_FUNC (
zbesh, ZBESH) (zr, zi, alpha, 2, 2, 1, &yr, &yi, nz,
ierr);
1084 double rexpz = expz.real ();
1085 double iexpz = expz.imag ();
1087 double tmp = yr*rexpz - yi*iexpz;
1089 yi = yr*iexpz + yi*rexpz;
1111 static inline Complex
1117 retval =
f (x, alpha, (scaled ? 2 : 1), ierr);
1170 if (x_nr != alpha_nr || x_nc != alpha_nc)
1171 (*current_liboctave_error_handler)
1172 (
"%s: the sizes of alpha and x must conform", fn);
1227 if (dv != alpha.
dims ())
1229 (
"%s: the sizes of alpha and x must conform", fn);
1261 #define SS_BESSEL(name, fcn) \
1263 name (double alpha, const Complex& x, bool scaled, octave_idx_type& ierr) \
1265 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
1268 #define SM_BESSEL(name, fcn) \
1270 name (double alpha, const ComplexMatrix& x, bool scaled, \
1271 Array<octave_idx_type>& ierr) \
1273 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
1276 #define MS_BESSEL(name, fcn) \
1278 name (const Matrix& alpha, const Complex& x, bool scaled, \
1279 Array<octave_idx_type>& ierr) \
1281 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
1284 #define MM_BESSEL(name, fcn) \
1286 name (const Matrix& alpha, const ComplexMatrix& x, bool scaled, \
1287 Array<octave_idx_type>& ierr) \
1289 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
1292 #define SN_BESSEL(name, fcn) \
1294 name (double alpha, const ComplexNDArray& x, bool scaled, \
1295 Array<octave_idx_type>& ierr) \
1297 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
1300 #define NS_BESSEL(name, fcn) \
1302 name (const NDArray& alpha, const Complex& x, bool scaled, \
1303 Array<octave_idx_type>& ierr) \
1305 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
1308 #define NN_BESSEL(name, fcn) \
1310 name (const NDArray& alpha, const ComplexNDArray& x, bool scaled, \
1311 Array<octave_idx_type>& ierr) \
1313 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
1316 #define RC_BESSEL(name, fcn) \
1318 name (const RowVector& alpha, const ComplexColumnVector& x, bool scaled, \
1319 Array<octave_idx_type>& ierr) \
1321 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
1324 #define ALL_BESSEL(name, fcn) \
1325 SS_BESSEL (name, fcn) \
1326 SM_BESSEL (name, fcn) \
1327 MS_BESSEL (name, fcn) \
1328 MM_BESSEL (name, fcn) \
1329 SN_BESSEL (name, fcn) \
1330 NS_BESSEL (name, fcn) \
1331 NN_BESSEL (name, fcn) \
1332 RC_BESSEL (name, fcn)
1404 return x ==
static_cast<float> (
static_cast<long> (
x));
1423 float expz = exp (
std::abs (z.imag ()));
1427 if (z.imag () == 0.0 && z.real () >= 0.0)
1437 if ((static_cast<long> (alpha)) & 1)
1446 *
cbesj (z, alpha, kode, ierr);
1448 if (ierr == 0 || ierr == 3)
1450 tmp -= sinf (static_cast<float> (M_PI) * alpha)
1451 *
cbesy (z, alpha, kode, ierr);
1478 if (z.real () == 0.0 && z.imag () == 0.0)
1489 float expz = exp (
std::abs (z.imag ()));
1493 if (z.imag () == 0.0 && z.real () >= 0.0)
1504 if ((static_cast<long> (alpha - 0.5)) & 1)
1513 *
cbesy (z, alpha, kode, ierr);
1515 if (ierr == 0 || ierr == 3)
1517 tmp += sinf (static_cast<float> (M_PI) * alpha)
1518 *
cbesj (z, alpha, kode, ierr);
1546 float expz = exp (
std::abs (z.real ()));
1550 if (z.imag () == 0.0 && z.real () >= 0.0)
1561 if (ierr == 0 || ierr == 3)
1564 * sinf (static_cast<float> (M_PI) * alpha)
1565 *
cbesk (z, alpha, kode, ierr);
1570 tmp2 *= exp (-z -
std::abs (z.real ()));
1598 if (z.real () == 0.0 && z.imag () == 0.0)
1611 float rexpz = expz.real ();
1612 float iexpz = expz.imag ();
1614 float tmp_r = y.real () * rexpz - y.imag () * iexpz;
1615 float tmp_i = y.real () * iexpz + y.imag () * rexpz;
1620 if (z.imag () == 0.0 && z.real () >= 0.0)
1654 float rexpz = expz.real ();
1655 float iexpz = expz.imag ();
1657 float tmp_r = y.real () * rexpz - y.imag () * iexpz;
1658 float tmp_i = y.real () * iexpz + y.imag () * rexpz;
1672 *
cbesh1 (z, alpha, kode, ierr);
1698 float rexpz = expz.real ();
1699 float iexpz = expz.imag ();
1701 float tmp_r = y.real () * rexpz - y.imag () * iexpz;
1702 float tmp_i = y.real () * iexpz + y.imag () * rexpz;
1716 *
cbesh2 (z, alpha, kode, ierr);
1727 static inline FloatComplex
1733 retval =
f (x, alpha, (scaled ? 2 : 1), ierr);
1758 const FloatComplex&
x,
1788 if (x_nr != alpha_nr || x_nc != alpha_nc)
1789 (*current_liboctave_error_handler)
1790 (
"%s: the sizes of alpha and x must conform", fn);
1846 if (dv != alpha.
dims ())
1848 (
"%s: the sizes of alpha and x must conform", fn);
1880 #define SS_BESSEL(name, fcn) \
1882 name (float alpha, const FloatComplex& x, bool scaled, octave_idx_type& ierr) \
1884 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
1887 #define SM_BESSEL(name, fcn) \
1888 FloatComplexMatrix \
1889 name (float alpha, const FloatComplexMatrix& x, bool scaled, \
1890 Array<octave_idx_type>& ierr) \
1892 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
1895 #define MS_BESSEL(name, fcn) \
1896 FloatComplexMatrix \
1897 name (const FloatMatrix& alpha, const FloatComplex& x, bool scaled, \
1898 Array<octave_idx_type>& ierr) \
1900 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
1903 #define MM_BESSEL(name, fcn) \
1904 FloatComplexMatrix \
1905 name (const FloatMatrix& alpha, const FloatComplexMatrix& x, bool scaled, \
1906 Array<octave_idx_type>& ierr) \
1908 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
1911 #define SN_BESSEL(name, fcn) \
1912 FloatComplexNDArray \
1913 name (float alpha, const FloatComplexNDArray& x, bool scaled, \
1914 Array<octave_idx_type>& ierr) \
1916 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
1919 #define NS_BESSEL(name, fcn) \
1920 FloatComplexNDArray \
1921 name (const FloatNDArray& alpha, const FloatComplex& x, bool scaled, \
1922 Array<octave_idx_type>& ierr) \
1924 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
1927 #define NN_BESSEL(name, fcn) \
1928 FloatComplexNDArray \
1929 name (const FloatNDArray& alpha, const FloatComplexNDArray& x, bool scaled, \
1930 Array<octave_idx_type>& ierr) \
1932 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
1935 #define RC_BESSEL(name, fcn) \
1936 FloatComplexMatrix \
1937 name (const FloatRowVector& alpha, const FloatComplexColumnVector& x, bool scaled, \
1938 Array<octave_idx_type>& ierr) \
1940 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
1943 #define ALL_BESSEL(name, fcn) \
1944 SS_BESSEL (name, fcn) \
1945 SM_BESSEL (name, fcn) \
1946 MS_BESSEL (name, fcn) \
1947 MM_BESSEL (name, fcn) \
1948 SN_BESSEL (name, fcn) \
1949 NS_BESSEL (name, fcn) \
1950 NN_BESSEL (name, fcn) \
1951 RC_BESSEL (name, fcn)
1978 double zr = z.real ();
1979 double zi = z.imag ();
1987 Complex expz = exp (- 2.0 / 3.0 * z * sqrt (z));
1989 double rexpz = expz.real ();
1990 double iexpz = expz.imag ();
1992 double tmp = ar*rexpz - ai*iexpz;
1994 ai = ar*iexpz + ai*rexpz;
1998 if (zi == 0.0 && (! scaled || zr >= 0.0))
2010 double zr = z.real ();
2011 double zi = z.imag ();
2021 double rexpz = expz.real ();
2022 double iexpz = expz.imag ();
2024 double tmp = ar*rexpz - ai*iexpz;
2026 ai = ar*iexpz + ai*rexpz;
2030 if (zi == 0.0 && (! scaled || zr >= 0.0))
2116 float ar = a.real ();
2117 float ai = a.imag ();
2121 FloatComplex expz = exp (- 2.0
f / 3.0
f * z * sqrt (z));
2123 float rexpz = expz.real ();
2124 float iexpz = expz.imag ();
2126 float tmp = ar*rexpz - ai*iexpz;
2128 ai = ar*iexpz + ai*rexpz;
2132 if (z.imag () == 0.0 && (! scaled || z.real () >= 0.0))
2148 float ar = a.real ();
2149 float ai = a.imag ();
2155 float rexpz = expz.real ();
2156 float iexpz = expz.imag ();
2158 float tmp = ar*rexpz - ai*iexpz;
2160 ai = ar*iexpz + ai*rexpz;
2164 if (z.imag () == 0.0 && (! scaled || z.real () >= 0.0))
2238 OCTAVE_NORETURN
static void
2246 (*current_liboctave_error_handler)
2247 (
"betainc: nonconformant arguments (x is %s, a is %s, b is %s)",
2248 d1_str.c_str (), d2_str.c_str (), d3_str.c_str ());
2251 OCTAVE_NORETURN
static void
2259 (*current_liboctave_error_handler)
2260 (
"betaincinv: nonconformant arguments (x is %s, a is %s, b is %s)",
2261 d1_str.c_str (), d2_str.c_str (), d3_str.c_str ());
2310 if (dv != b.
dims ())
2347 if (dv != b.
dims ())
2368 if (dv != a.
dims ())
2389 if (dv != a.
dims () || dv != b.
dims ())
2450 if (dv != b.
dims ())
2487 if (dv != b.
dims ())
2508 if (dv != a.
dims ())
2529 if (dv != a.
dims () || dv != b.
dims ())
2549 if (a < 0.0 || x < 0.0)
2550 (*current_liboctave_error_handler)
2551 (
"gammainc: A and X must be non-negative");
2615 if (nr != a_nr || nc != a_nc)
2616 (*current_liboctave_error_handler)
2617 (
"gammainc: nonconformant arguments (arg 1 is %dx%d, arg 2 is %dx%d)",
2684 if (dv != a.
dims ())
2689 (*current_liboctave_error_handler)
2690 (
"gammainc: nonconformant arguments (arg 1 is %s, arg 2 is %s)",
2691 x_str.c_str (), a_str.c_str ());
2712 if (a < 0.0 || x < 0.0)
2713 (*current_liboctave_error_handler)
2714 (
"gammainc: A and X must be non-negative");
2778 if (nr != a_nr || nc != a_nc)
2779 (*current_liboctave_error_handler)
2780 (
"gammainc: nonconformant arguments (arg 1 is %dx%d, arg 2 is %dx%d)",
2847 if (dv != a.
dims ())
2852 (*current_liboctave_error_handler)
2853 (
"gammainc: nonconformant arguments (arg 1 is %s, arg 2 is %s)",
2854 x_str.c_str (), a_str.c_str ());
2874 const double pi = 3.14159265358979323846;
2876 ? Complex (
std::log (-(1.0 + x)), pi)
2877 : Complex (
log1p (x)));
2882 const float pi = 3.14159265358979323846f;
2884 ? FloatComplex (
std::log (-(1.0
f + x)), pi)
2885 : FloatComplex (
log1p (x)));
2900 static const double a[] =
2902 -2.806989788730439e+01, 1.562324844726888e+02,
2903 -1.951109208597547e+02, 9.783370457507161e+01,
2904 -2.168328665628878e+01, 1.772453852905383e+00
2906 static const double b[] =
2908 -5.447609879822406e+01, 1.615858368580409e+02,
2909 -1.556989798598866e+02, 6.680131188771972e+01,
2910 -1.328068155288572e+01
2912 static const double c[] =
2914 -5.504751339936943e-03, -2.279687217114118e-01,
2915 -1.697592457770869e+00, -1.802933168781950e+00,
2916 3.093354679843505e+00, 2.077595676404383e+00
2918 static const double d[] =
2920 7.784695709041462e-03, 3.224671290700398e-01,
2921 2.445134137142996e+00, 3.754408661907416e+00
2924 static const double spi2 = 8.862269254527579e-01;
2925 static const double pbreak = 0.95150;
2926 double ax = fabs (x),
y;
2932 const double q = 0.5 *
x, r = q*q;
2933 const double yn = (((((a[0]*r + a[1])*r + a[2])*r + a[3])*r + a[4])*r + a[5])*q;
2934 const double yd = ((((b[0]*r + b[1])*r + b[2])*r + b[3])*r + b[4])*r + 1.0;
2940 const double q = sqrt (-2*
std::log (0.5*(1-ax)));
2941 const double yn = ((((c[0]*q + c[1])*q + c[2])*q + c[3])*q + c[4])*q + c[5];
2942 const double yd = (((d[0]*q + d[1])*q + d[2])*q + d[3])*q + 1.0;
2953 double u = (
erf (
y) -
x) * spi2 * exp (
y*
y);
2978 static const double a[] =
2980 -2.806989788730439e+01, 1.562324844726888e+02,
2981 -1.951109208597547e+02, 9.783370457507161e+01,
2982 -2.168328665628878e+01, 1.772453852905383e+00
2984 static const double b[] =
2986 -5.447609879822406e+01, 1.615858368580409e+02,
2987 -1.556989798598866e+02, 6.680131188771972e+01,
2988 -1.328068155288572e+01
2990 static const double c[] =
2992 -5.504751339936943e-03, -2.279687217114118e-01,
2993 -1.697592457770869e+00, -1.802933168781950e+00,
2994 3.093354679843505e+00, 2.077595676404383e+00
2996 static const double d[] =
2998 7.784695709041462e-03, 3.224671290700398e-01,
2999 2.445134137142996e+00, 3.754408661907416e+00
3002 static const double spi2 = 8.862269254527579e-01;
3003 static const double pbreak_lo = 0.04850;
3004 static const double pbreak_hi = 1.95150;
3008 if (x >= pbreak_lo && x <= pbreak_hi)
3011 const double q = 0.5*(1-
x), r = q*q;
3012 const double yn = (((((a[0]*r + a[1])*r + a[2])*r + a[3])*r + a[4])*r + a[5])*q;
3013 const double yd = ((((b[0]*r + b[1])*r + b[2])*r + b[3])*r + b[4])*r + 1.0;
3016 else if (x > 0.0 && x < 2.0)
3019 const double q = (x < 1
3023 const double yn = ((((c[0]*q + c[1])*q + c[2])*q + c[3])*q + c[4])*q + c[5];
3025 const double yd = (((d[0]*q + d[1])*q + d[2])*q + d[3])*q + 1.0;
3042 double u = (
erf (y) - (1-
x)) * spi2 * exp (y*y);
3078 double acu = 0.1E-14, ai, cx;
3081 double pp, psq, qq, rx, temp, term,
value, xx;
3088 if ((p <= 0.0 || q <= 0.0) || (x < 0.0 || 1.0 < x))
3096 if (x == 0.0 || x == 1.0)
3125 ns =
static_cast<int> (qq + cx * psq);
3138 term *= temp * rx / (pp + ai);
3142 if (temp <= acu && temp <= acu * value)
3145 + (qq - 1.0) *
std::log (cx) - beta) / pp;
3149 value = 1.0 -
value;
3195 double a, acu, adj, fpu, g,
h;
3198 double pp, prev, qq, r,
s, sae = -37.0, sq,
t, tx,
value,
w, xin, ycur, yprev;
3202 fpu =
pow (10.0, sae);
3207 if (p <= 0.0 || q <= 0.0)
3208 (*current_liboctave_error_handler) (
"betaincinv: wrong parameters");
3209 if (y < 0.0 || 1.0 < y)
3210 (*current_liboctave_error_handler) (
"betaincinv: wrong parameter Y");
3212 if (y == 0.0 || y == 1.0)
3236 ycur = r - (2.30753 + 0.27061 * r) / (1.0 + (0.99229 + 0.04481 * r) * r);
3238 if (1.0 < pp && 1.0 < qq)
3240 r = (ycur * ycur - 3.0) / 6.0;
3241 s = 1.0 / (pp + pp - 1.0);
3242 t = 1.0 / (qq + qq - 1.0);
3244 w = ycur * sqrt (h + r) / h - (t -
s) * (r + 5.0 / 6.0 - 2.0 / (3.0 * h));
3245 value = pp / (pp + qq * exp (w + w));
3250 t = 1.0 / (9.0 * qq);
3251 t = r *
pow (1.0 - t + ycur * sqrt (t), 3);
3255 value = 1.0 - exp ((
std::log ((1.0 - a) * qq) + beta) / qq);
3259 t = (4.0 * pp + r - 2.0) /
t;
3263 value = exp ((
std::log (a * pp) + beta) / pp);
3267 value = 1.0 - 2.0 / (t + 1.0);
3291 iex =
std::max (- 5.0 / pp / pp - 1.0 /
pow (a, 0.2) - 13.0, sae);
3293 acu =
pow (10.0, iex);
3297 ycur =
betain (value, pp, qq, beta, err);
3305 ycur = (ycur -
a) * exp (beta + r *
std::log (xin)
3308 if (ycur * yprev <= 0.0)
3326 if (0.0 <= tx && tx <= 1.0)
3338 value = 1.0 -
value;
3343 if (ycur * ycur <= acu)
3347 value = 1.0 -
value;
3352 if (tx != 0.0 && tx != 1.0)
3371 value = 1.0 -
value;
3415 if (dv != b.
dims ())
3452 if (dv != b.
dims ())
3473 if (dv != a.
dims ())
3495 if (dv != a.
dims () && dv != b.
dims ())
3511 ellipj (
double u,
double m,
double& sn,
double& cn,
double& dn,
double&
err)
3513 static const int Nmax = 16;
3514 double m1,
t=0, si_u, co_u, se_u, ta_u,
b,
c[Nmax],
a[Nmax], phi;
3519 (*current_liboctave_warning_with_id_handler)
3520 (
"Octave:ellipj-invalid-m",
3521 "ellipj: invalid M value, required value 0 <= M <= 1");
3528 double sqrt_eps = sqrt (std::numeric_limits<double>::epsilon ());
3534 t = 0.25*m*(u - si_u*co_u);
3535 sn = si_u - t * co_u;
3536 cn = co_u + t * si_u;
3537 dn = 1 - 0.5*m*si_u*si_u;
3539 else if ((1 - m) < sqrt_eps)
3547 sn = ta_u + 0.25*m1*(si_u*co_u -
u)*se_u*se_u;
3548 cn = se_u - 0.25*m1*(si_u*co_u -
u)*ta_u*se_u;
3549 dn = se_u + 0.25*m1*(si_u*co_u +
u)*ta_u*se_u;
3558 for (n = 1; n < Nmax; ++n)
3560 a[n] = (a[n - 1] +
b)/2;
3561 c[n] = (a[n - 1] -
b)/2;
3562 b = sqrt (a[n - 1]*b);
3563 if (c[n]/a[n] < std::numeric_limits<double>::epsilon ())
break;
3571 for (ii = 1; n > 0; ii = ii*2, --n) ;
3573 for (n = Nn; n > 0; --n)
3575 phi = (
asin ((c[n]/a[n])* sin (phi)) + phi)/2;
3579 dn = sqrt (1 - m*sn*sn);
3584 ellipj (
const Complex&
u,
double m, Complex& sn, Complex& cn, Complex& dn,
3587 double m1 = 1 -
m, ss1, cc1, dd1;
3589 ellipj (u.imag (), m1, ss1, cc1, dd1,
err);
3600 double ss, cc, dd, ddd;
3603 ddd = cc1*cc1 + m*ss*ss*ss1*ss1;
3604 sn =
Complex (ss*dd1/ddd, cc*dd*ss1*cc1/ddd);
3605 cn =
Complex (cc*cc1/ddd, -ss*dd*ss1*dd1/ddd);
3606 dn =
Complex (dd*cc1*dd1/ddd, -m*ss*cc*ss1/ddd);
3610 static const double pi = 3.14159265358979323846;
3612 template <
typename T>
3633 template <
typename T>
3641 static const T dg_coeff[10] =
3643 -0.83333333333333333e-1, 0.83333333333333333e-2,
3644 -0.39682539682539683e-2, 0.41666666666666667e-2,
3645 -0.75757575757575758e-2, 0.21092796092796093e-1,
3646 -0.83333333333333333e-1, 0.4432598039215686,
3647 -0.3053954330270122e+1, 0.125318899521531e+2
3650 T overz2 = T (1.0) / (zc * zc);
3655 p += dg_coeff[
k] * overz2k;
3656 p +=
xlog (zc) - T (0.5) / zc;
3660 template <
typename T>
3664 static const double euler_mascheroni =
3665 0.577215664901532860606512090082402431042;
3677 p =
psi (1 - z) - (pi / tan (pi * z));
3682 p = - euler_mascheroni;
3690 p += 1.0 / (2 *
k - 1);
3692 p = - euler_mascheroni - 2 *
std::log (2) + 2 * (
p);
3702 const signed char n = 10 - z;
3703 for (
signed char k = n - 1;
k >= 0;
k--)
3717 template <
typename T>
3723 typedef typename std::complex<T>::value_type
P;
3728 std::complex<T> dgam (0.0, 0.0);
3730 dgam = std::complex<T> (
psi (z_r), 0.0);
3732 dgam =
psi (P (1.0) - z)- (
P (pi) / tan (P (pi) * z));
3736 std::complex<T> z_m = z;
3739 unsigned char n = 8 - z_ra;
3740 z_m = z + std::complex<T> (n, 0.0);
3745 std::complex<T> z_p = z +
P (n - 1);
3746 for (
unsigned char k = n;
k > 0;
k--, z_p -= 1.0)
3747 dgam -=
P (1.0) / z_p;
3763 Complex
psi (
const Complex& z) {
return xpsi (z); }
3764 FloatComplex
psi (
const FloatComplex& z) {
return xpsi (z); }
3766 template <
typename T>
3789 template <
typename T>
3795 fortran_psifn<T> (z, n, &ans, &
ierr);
3804 ans = ans / (
pow (-1.0, n + 1) /
gamma (
double (n+1)));
3821 #if defined (OCTAVE_USE_DEPRECATED_FUNCTIONS)
void resize(octave_idx_type nr, octave_idx_type nc, const FloatComplex &rfv=FloatComplex(0))
F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_REAL const F77_REAL F77_REAL &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_DBLE F77_DBLE &F77_RET_T const F77_REAL F77_REAL &F77_RET_T F77_REAL F77_REAL &F77_RET_T F77_DBLE F77_DBLE &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE const F77_INT F77_INT * ierr
Octave interface to the compression and uncompression libraries.
subroutine dpsifn(X, N, KODE, M, ANS, NZ, IERR)
subroutine xerf(x, result)
std::string str(char sep= 'x') const
static FloatComplex cbesh1(const FloatComplex &z, float alpha, int kode, octave_idx_type &ierr)
subroutine dlgams(X, DLGAM, SGNGAM)
subroutine xdgamma(x, result)
static OCTAVE_NORETURN void err_betaincinv_nonconformant(const dim_vector &d1, const dim_vector &d2, const dim_vector &d3)
subroutine cairy(Z, ID, KODE, AI, NZ, IERR)
subroutine zairy(ZR, ZI, ID, KODE, AIR, AII, NZ, IERR)
std::complex< double > erfi(std::complex< double > z, double relerr=0)
octave_idx_type numel(void) const
Number of elements in the array.
identity matrix If supplied two scalar respectively For allows like xample val
F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_REAL const F77_REAL F77_REAL &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_DBLE F77_DBLE &F77_RET_T const F77_REAL F77_REAL &F77_RET_T F77_REAL F77_REAL &F77_RET_T F77_DBLE F77_DBLE &F77_RET_T const F77_DBLE const F77_DBLE * f
FloatComplexMatrix besselh1(const FloatRowVector &alpha, const FloatComplexColumnVector &x, bool scaled, Array< octave_idx_type > &ierr)
static Complex zbesk(const Complex &z, double alpha, int kode, octave_idx_type &ierr)
FloatComplex log1p(const FloatComplex &x)
OCTAVE_NORETURN liboctave_error_handler current_liboctave_error_handler
Complex besselk(double alpha, const Complex &x, bool scaled, octave_idx_type &ierr)
double betaincinv(double y, double p, double q)
Complex rc_lgamma(double x)
OCTAVE_EXPORT octave_value_list or N dimensional array whose elements are all equal to the base of natural logarithms The constant ex $e satisfies the equation log(e)
Complex besselh1(double alpha, const Complex &x, bool scaled, octave_idx_type &ierr)
subroutine cbiry(Z, ID, KODE, BI, IERR)
static T xlog(const T &x)
std::complex< double > Dawson(std::complex< double > z, double relerr=0)
octave_value real(void) const
double lo_ieee_nan_value(void)
FloatComplexMatrix besselh2(const FloatRowVector &alpha, const FloatComplexColumnVector &x, bool scaled, Array< octave_idx_type > &ierr)
OCTAVE_EXPORT octave_value_list return the number of command line arguments passed to Octave If called with the optional argument the function t
subroutine xderfc(x, result)
Complex asin(const Complex &x)
#define F77_XFCN(f, F, args)
subroutine xerfc(x, result)
subroutine xdacosh(x, result)
std::complex< double > erf(std::complex< double > z, double relerr=0)
FloatComplexMatrix besseli(const FloatRowVector &alpha, const FloatComplexColumnVector &x, bool scaled, Array< octave_idx_type > &ierr)
octave_idx_type rows(void) const
octave_idx_type numel(int n=0) const
Number of elements that a matrix with this dimensions would have.
octave_value imag(void) const
F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_REAL const F77_REAL F77_REAL &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_DBLE F77_DBLE &F77_RET_T const F77_REAL F77_REAL &F77_RET_T F77_REAL F77_REAL &F77_RET_T F77_DBLE F77_DBLE &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE * d
static Complex zbesj(const Complex &z, double alpha, int kode, octave_idx_type &ierr)
#define ALL_BESSEL(name, fcn)
FloatComplexNDArray biry(const FloatComplexNDArray &z, bool deriv, bool scaled, Array< octave_idx_type > &ierr)
subroutine xgamma(x, result)
calling an anonymous function involves an overhead quite comparable to the overhead of an m file function Passing a handle to a built in function is because the interpreter is not involved in the internal loop For a
const dim_vector & dims(void) const
Return a const-reference so that dims ()(i) works efficiently.
subroutine xdatanh(x, result)
FloatNDArray gammainc(const FloatNDArray &x, const FloatNDArray &a)
subroutine algams(X, ALGAM, SGNGAM)
static void fortran_psifn(const T z, const octave_idx_type n, T *ans, octave_idx_type *ierr)
Complex besseli(double alpha, const Complex &x, bool scaled, octave_idx_type &ierr)
FloatComplex expm1(const FloatComplex &x)
FloatComplexMatrix besselk(const FloatRowVector &alpha, const FloatComplexColumnVector &x, bool scaled, Array< octave_idx_type > &ierr)
nd deftypefn *octave_map m
static double do_erfinv(double x, bool refine)
static FloatComplex cbesk(const FloatComplex &z, float alpha, int kode, octave_idx_type &ierr)
std::complex< double > w(std::complex< double > z, double relerr=0)
static Complex zbesi(const Complex &z, double alpha, int kode, octave_idx_type &ierr)
static double do_erfcinv(double x, bool refine)
Array< float > betainc(const Array< float > &x, const Array< float > &a, const Array< float > &b)
subroutine zbiry(ZR, ZI, ID, KODE, BIR, BII, IERR)
static Complex zbesh2(const Complex &z, double alpha, int kode, octave_idx_type &ierr)
subroutine xacosh(x, result)
subroutine psifn(X, N, KODE, M, ANS, NZ, IERR)
void resize(const dim_vector &dv, const T &rfv)
static Complex bessel_return_value(const Complex &val, octave_idx_type ierr)
subroutine xgammainc(a, x, result)
static T lanczos_approximation_psi(const T zc)
octave_int< T > pow(const octave_int< T > &a, const octave_int< T > &b)
F77_RET_T F77_FUNC(xstopx, XSTOPX) const
std::complex< double > erfcx(std::complex< double > z, double relerr=0)
void fortran_psifn< float >(const float z, const octave_idx_type n, float *ans, octave_idx_type *ierr)
subroutine zbesh(ZR, ZI, FNU, KODE, M, N, CYR, CYI, NZ, IERR)
the sparsity preserving column transformation such that that defines the pivoting threshold can be given in which case it defines the c
subroutine xdasinh(x, result)
Array< double > betaincinv(const Array< double > &x, const Array< double > &a, const Array< double > &b)
With real return the complex result
subroutine cbesh(Z, FNU, KODE, M, N, CY, NZ, IERR)
subroutine xdbetai(x, a, b, result)
static FloatComplex cbesy(const FloatComplex &z, float alpha, int kode, octave_idx_type &ierr)
charNDArray max(char d, const charNDArray &m)
SparseMatrix atan2(const double &x, const SparseMatrix &y)
static Complex zbesy(const Complex &z, double alpha, int kode, octave_idx_type &ierr)
FloatComplex erfc(const FloatComplex &x)
static FloatComplex cbesh2(const FloatComplex &z, float alpha, int kode, octave_idx_type &ierr)
double betainc(double x, double a, double b)
void resize(octave_idx_type nr, octave_idx_type nc, const Complex &rfv=Complex(0))
Complex(* dptr)(const Complex &, double, int, octave_idx_type &)
static FloatComplex cbesi(const FloatComplex &z, float alpha, int kode, octave_idx_type &ierr)
=val(i)}if ode{val(i)}occurs in table i
OCTAVE_EXPORT octave_value_list return the value of the option it must match the dimension of the state and the relative tolerance must also be a vector of the same length tem it must match the dimension of the state and the absolute tolerance must also be a vector of the same length The local error test applied at each integration step is xample roup abs(local error in x(i))<
bool negative_sign(double x)
OCTAVE_EXPORT octave_value_list or N dimensional array whose elements are all equal to the IEEE symbol zero divided by nd tex zero divided by nd ifnottex and any operation involving another NaN value(5+NaN).Note that NaN always compares not equal to NaN(NaN!
subroutine xbetai(x, a, b, result)
issues an error eealso double
static FloatComplex cbesj(const FloatComplex &z, float alpha, int kode, octave_idx_type &ierr)
static Complex zbesh1(const Complex &z, double alpha, int kode, octave_idx_type &ierr)
Complex rc_log1p(double x)
the element is set to zero In other the statement xample y
OCTAVE_EXPORT octave_value_list or N dimensional array whose elements are all equal to the IEEE symbol NaN(Not a Number).NaN is the result of operations which do not produce a well defined 0 result.Common operations which produce a NaN are arithmetic with infinity ex($\infty-\infty $)
subroutine xsgammainc(a, x, result)
Complex besselh2(double alpha, const Complex &x, bool scaled, octave_idx_type &ierr)
void ellipj(double u, double m, double &sn, double &cn, double &dn, double &err)
#define F77_CONST_CMPLX_ARG(x)
void fortran_psifn< double >(const double z, const octave_idx_type n, double *ans, octave_idx_type *ierr)
double gammainc(double x, double a, bool &err)
std::complex< float > FloatComplex
OCTAVE_EXPORT octave_value_list error nd deftypefn *const octave_scalar_map err
FloatComplex(* fptr)(const FloatComplex &, float, int, octave_idx_type &)
Complex airy(const Complex &z, bool deriv, bool scaled, octave_idx_type &ierr)
FloatComplexMatrix besselj(const FloatRowVector &alpha, const FloatComplexColumnVector &x, bool scaled, Array< octave_idx_type > &ierr)
FloatComplexMatrix bessely(const FloatRowVector &alpha, const FloatComplexColumnVector &x, bool scaled, Array< octave_idx_type > &ierr)
std::complex< double > Complex
const T * fortran_vec(void) const
subroutine xatanh(x, result)
FloatComplexNDArray airy(const FloatComplexNDArray &z, bool deriv, bool scaled, Array< octave_idx_type > &ierr)
static bool is_integer_value(double x)
ColumnVector real(const ComplexColumnVector &a)
octave_idx_type cols(void) const
double betain(double x, double p, double q, double beta, bool &err)
Vector representing the dimensions (size) of an Array.
If this string is the system will ring the terminal sometimes it is useful to be able to print the original representation of the string
static OCTAVE_NORETURN void err_betainc_nonconformant(const dim_vector &d1, const dim_vector &d2, const dim_vector &d3)
F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_REAL const F77_REAL F77_REAL &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_DBLE F77_DBLE &F77_RET_T const F77_REAL F77_REAL &F77_RET_T F77_REAL F77_REAL &F77_RET_T F77_DBLE F77_DBLE &F77_RET_T const F77_DBLE * x
Complex besselj(double alpha, const Complex &x, bool scaled, octave_idx_type &ierr)
Complex bessely(double alpha, const Complex &x, bool scaled, octave_idx_type &ierr)
static Complex do_bessel(dptr f, const char *, double alpha, const Complex &x, bool scaled, octave_idx_type &ierr)
subroutine xasinh(x, result)
subroutine xderf(x, result)
Complex biry(const Complex &z, bool deriv, bool scaled, octave_idx_type &ierr)
std::complex< double > erfc(std::complex< double > z, double relerr=0)