54 F77_FUNC (
zbesj, ZBESJ) (
const double&,
const double&,
const double&,
56 double*,
double*, octave_idx_type&,
60 F77_FUNC (
zbesy, ZBESY) (
const double&,
const double&,
const double&,
61 const octave_idx_type&,
const octave_idx_type&,
62 double*,
double*, octave_idx_type&,
double*,
63 double*, octave_idx_type&);
66 F77_FUNC (
zbesi, ZBESI) (
const double&,
const double&,
const double&,
67 const octave_idx_type&,
const octave_idx_type&,
68 double*,
double*, octave_idx_type&,
72 F77_FUNC (
zbesk, ZBESK) (
const double&,
const double&,
const double&,
73 const octave_idx_type&,
const octave_idx_type&,
74 double*,
double*, octave_idx_type&,
78 F77_FUNC (
zbesh, ZBESH) (
const double&,
const double&,
const double&,
79 const octave_idx_type&,
const octave_idx_type&,
80 const octave_idx_type&,
double*,
double*,
81 octave_idx_type&, octave_idx_type&);
85 const octave_idx_type&,
const octave_idx_type&,
86 FloatComplex*, octave_idx_type&, octave_idx_type&);
90 const octave_idx_type&,
const octave_idx_type&,
91 FloatComplex*, octave_idx_type&,
92 FloatComplex*, octave_idx_type&);
96 const octave_idx_type&,
const octave_idx_type&,
97 FloatComplex*, octave_idx_type&, octave_idx_type&);
101 const octave_idx_type&,
const octave_idx_type&,
102 FloatComplex*, octave_idx_type&, octave_idx_type&);
106 const octave_idx_type&,
const octave_idx_type&,
107 const octave_idx_type&, FloatComplex*,
108 octave_idx_type&, octave_idx_type&);
112 const octave_idx_type&,
const octave_idx_type&,
113 double&,
double&, octave_idx_type&,
117 F77_FUNC (
cairy, CAIRY) (
const float&,
const float&,
const octave_idx_type&,
118 const octave_idx_type&,
float&,
float&,
119 octave_idx_type&, octave_idx_type&);
123 const octave_idx_type&,
const octave_idx_type&,
124 double&,
double&, octave_idx_type&);
127 F77_FUNC (
cbiry, CBIRY) (
const float&,
const float&,
const octave_idx_type&,
128 const octave_idx_type&,
float&,
float&,
163 const double&,
double&);
167 const float&,
float&);
188 #if !defined (HAVE_ACOSH)
198 #if !defined (HAVE_ACOSHF)
208 #if !defined (HAVE_ASINH)
218 #if !defined (HAVE_ASINHF)
228 #if !defined (HAVE_ATANH)
238 #if !defined (HAVE_ATANHF)
248 #if !defined (HAVE_ERF)
258 #if !defined (HAVE_ERFF)
268 #if !defined (HAVE_ERFC)
278 #if !defined (HAVE_ERFCF)
380 #if defined (HAVE_TGAMMA)
393 #if defined (HAVE_LGAMMA)
415 #if defined (HAVE_LGAMMA_R)
417 result = lgamma_r (x, &sgngam);
431 return result +
Complex (0., M_PI);
452 #if defined (HAVE_TGAMMA)
453 result = tgammaf (x);
465 #if defined (HAVE_LGAMMAF)
487 #if defined (HAVE_LGAMMAF_R)
489 result = lgammaf_r (x, &sgngam);
508 #if !defined (HAVE_EXPM1)
514 double ax = fabs (x);
523 for (
int i = 2; i < 7; i++)
529 for (
int i = 0; i < 4; i++)
535 retval = (x > 0) ? s : -s / (1+s);
538 retval = exp (x) - 1;
551 double im = x.imag ();
552 double u =
expm1 (x.real ());
553 double v = sin (im/2);
555 retval =
Complex (u*v + u + v, (u+1) * sin (im));
558 retval = std::exp (x) -
Complex (1);
563 #if !defined (HAVE_EXPM1F)
578 for (
int i = 2; i < 7; i++)
584 for (
int i = 0; i < 4; i++)
590 retval = (x > 0) ? s : -s / (1+s);
593 retval = exp (x) - 1;
606 float im = x.imag ();
607 float u =
expm1 (x.real ());
608 float v = sin (im/2);
618 #if !defined (HAVE_LOG1P)
624 double ax = fabs (x);
629 double u = x / (2 +
x), t = 1, s = 0;
630 for (
int i = 2; i < 12; i += 2)
631 s += (t *= u*u) / (i+1);
633 retval = 2 * (s + 1) * u;
636 retval = gnulib::log (1 + x);
647 double r = x.real (), i = x.imag ();
649 if (fabs (r) < 0.5 && fabs (i) < 0.5)
651 double u = 2*r + r*r + i*i;
656 retval = std::log (
Complex (1) + x);
661 #if !defined (HAVE_CBRT)
664 static const double one_third = 0.3333333333333333333;
670 return (x / (y*y) + y + y) / 3;
677 #if !defined (HAVE_LOG1PF)
688 float u = x / (2 +
x), t = 1.0
f, s = 0;
689 for (
int i = 2; i < 12; i += 2)
690 s += (t *= u*u) / (i+1);
692 retval = 2 * (s + 1.0f) * u;
695 retval = gnulib::logf (1.0
f + x);
706 float r = x.real (), i = x.imag ();
708 if (fabs (r) < 0.5 && fabs (i) < 0.5)
710 float u = 2*r + r*r + i*i;
720 #if !defined (HAVE_CBRTF)
723 static const float one_third = 0.3333333333333333333f;
729 return (x / (y*y) + y + y) / 3;
784 return x ==
static_cast<double> (
static_cast<long> (
x));
799 double zr = z.real ();
800 double zi = z.imag ();
811 if (zi == 0.0 && zr >= 0.0)
821 if ((static_cast<long> (alpha)) & 1)
829 Complex tmp = cos (M_PI * alpha) *
zbesj (z, alpha, kode, ierr);
831 if (ierr == 0 || ierr == 3)
833 tmp -= sin (M_PI * alpha) *
zbesy (z, alpha, kode, ierr);
858 double zr = z.real ();
859 double zi = z.imag ();
863 if (zr == 0.0 && zi == 0.0)
880 if (zi == 0.0 && zr >= 0.0)
891 if ((static_cast<long> (alpha - 0.5)) & 1)
899 Complex tmp = cos (M_PI * alpha) *
zbesy (z, alpha, kode, ierr);
901 if (ierr == 0 || ierr == 3)
903 tmp += sin (M_PI * alpha) *
zbesj (z, alpha, kode, ierr);
926 double zr = z.real ();
927 double zi = z.imag ();
938 if (zi == 0.0 && zr >= 0.0)
956 if (ierr == 0 || ierr == 3)
958 Complex tmp2 = (2.0 / M_PI) * sin (M_PI * alpha)
959 *
zbesk (z, alpha, kode, ierr);
964 tmp2 *= exp (-z -
std::abs (z.real ()));
990 double zr = z.real ();
991 double zi = z.imag ();
995 if (zr == 0.0 && zi == 0.0)
1008 double rexpz =
real (expz);
1009 double iexpz =
imag (expz);
1011 double tmp = yr*rexpz - yi*iexpz;
1013 yi = yr*iexpz + yi*rexpz;
1017 if (zi == 0.0 && zr >= 0.0)
1045 double zr = z.real ();
1046 double zi = z.imag ();
1048 F77_FUNC (
zbesh, ZBESH) (zr, zi, alpha, 2, 1, 1, &yr, &yi, nz,
ierr);
1054 double rexpz =
real (expz);
1055 double iexpz =
imag (expz);
1057 double tmp = yr*rexpz - yi*iexpz;
1059 yi = yr*iexpz + yi*rexpz;
1071 Complex tmp = exp (M_PI * alpha * eye) *
zbesh1 (z, alpha, kode, ierr);
1091 double zr = z.real ();
1092 double zi = z.imag ();
1094 F77_FUNC (
zbesh, ZBESH) (zr, zi, alpha, 2, 2, 1, &yr, &yi, nz,
ierr);
1100 double rexpz =
real (expz);
1101 double iexpz =
imag (expz);
1103 double tmp = yr*rexpz - yi*iexpz;
1105 yi = yr*iexpz + yi*rexpz;
1117 Complex tmp = exp (-M_PI * alpha * eye) *
zbesh2 (z, alpha, kode, ierr);
1133 retval =
f (x, alpha, (scaled ? 2 : 1), ierr);
1151 retval(i,j) =
f (
x(i,j), alpha, (scaled ? 2 : 1),
ierr(i,j));
1169 retval(i,j) =
f (x, alpha(i,j), (scaled ? 2 : 1),
ierr(i,j));
1186 if (x_nr == alpha_nr && x_nc == alpha_nc)
1197 retval(i,j) =
f (
x(i,j), alpha(i,j), (scaled ? 2 : 1),
ierr(i,j));
1201 (
"%s: the sizes of alpha and x must conform", fn);
1217 retval(i) =
f (
x(i), alpha, (scaled ? 2 : 1),
ierr(i));
1233 retval(i) =
f (x, alpha(i), (scaled ? 2 : 1),
ierr(i));
1245 if (dv == alpha.
dims ())
1253 retval(i) =
f (
x(i), alpha(i), (scaled ? 2 : 1),
ierr(i));
1257 (
"%s: the sizes of alpha and x must conform", fn);
1276 retval(i,j) =
f (
x(i), alpha(j), (scaled ? 2 : 1),
ierr(i,j));
1281 #define SS_BESSEL(name, fcn) \
1283 name (double alpha, const Complex& x, bool scaled, octave_idx_type& ierr) \
1285 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
1288 #define SM_BESSEL(name, fcn) \
1290 name (double alpha, const ComplexMatrix& x, bool scaled, \
1291 Array<octave_idx_type>& ierr) \
1293 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
1296 #define MS_BESSEL(name, fcn) \
1298 name (const Matrix& alpha, const Complex& x, bool scaled, \
1299 Array<octave_idx_type>& ierr) \
1301 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
1304 #define MM_BESSEL(name, fcn) \
1306 name (const Matrix& alpha, const ComplexMatrix& x, bool scaled, \
1307 Array<octave_idx_type>& ierr) \
1309 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
1312 #define SN_BESSEL(name, fcn) \
1314 name (double alpha, const ComplexNDArray& x, bool scaled, \
1315 Array<octave_idx_type>& ierr) \
1317 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
1320 #define NS_BESSEL(name, fcn) \
1322 name (const NDArray& alpha, const Complex& x, bool scaled, \
1323 Array<octave_idx_type>& ierr) \
1325 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
1328 #define NN_BESSEL(name, fcn) \
1330 name (const NDArray& alpha, const ComplexNDArray& x, bool scaled, \
1331 Array<octave_idx_type>& ierr) \
1333 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
1336 #define RC_BESSEL(name, fcn) \
1338 name (const RowVector& alpha, const ComplexColumnVector& x, bool scaled, \
1339 Array<octave_idx_type>& ierr) \
1341 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
1344 #define ALL_BESSEL(name, fcn) \
1345 SS_BESSEL (name, fcn) \
1346 SM_BESSEL (name, fcn) \
1347 MS_BESSEL (name, fcn) \
1348 MM_BESSEL (name, fcn) \
1349 SN_BESSEL (name, fcn) \
1350 NS_BESSEL (name, fcn) \
1351 NN_BESSEL (name, fcn) \
1352 RC_BESSEL (name, fcn)
1421 return x ==
static_cast<float> (
static_cast<long> (
x));
1443 if (
imag (z) == 0.0 &&
real (z) >= 0.0)
1453 if ((static_cast<long> (alpha)) & 1)
1461 FloatComplex tmp = cosf (static_cast<float> (M_PI) * alpha)
1462 *
cbesj (z, alpha, kode, ierr);
1464 if (ierr == 0 || ierr == 3)
1466 tmp -= sinf (static_cast<float> (M_PI) * alpha)
1467 *
cbesy (z, alpha, kode, ierr);
1493 if (
real (z) == 0.0 &&
imag (z) == 0.0)
1507 if (
imag (z) == 0.0 &&
real (z) >= 0.0)
1518 if ((static_cast<long> (alpha - 0.5)) & 1)
1526 FloatComplex tmp = cosf (static_cast<float> (M_PI) * alpha)
1527 *
cbesy (z, alpha, kode, ierr);
1529 if (ierr == 0 || ierr == 3)
1531 tmp += sinf (static_cast<float> (M_PI) * alpha)
1532 *
cbesj (z, alpha, kode, ierr);
1562 if (
imag (z) == 0.0 &&
real (z) >= 0.0)
1573 if (ierr == 0 || ierr == 3)
1576 * sinf (static_cast<float> (M_PI) * alpha)
1577 *
cbesk (z, alpha, kode, ierr);
1582 tmp2 *= exp (-z -
std::abs (z.real ()));
1609 if (
real (z) == 0.0 &&
imag (z) == 0.0)
1621 float rexpz =
real (expz);
1622 float iexpz =
imag (expz);
1624 float tmp_r =
real (y) * rexpz -
imag (y) * iexpz;
1625 float tmp_i =
real (y) * iexpz +
imag (y) * rexpz;
1630 if (
imag (z) == 0.0 &&
real (z) >= 0.0)
1663 float rexpz =
real (expz);
1664 float iexpz =
imag (expz);
1666 float tmp_r =
real (y) * rexpz -
imag (y) * iexpz;
1667 float tmp_i =
real (y) * iexpz +
imag (y) * rexpz;
1680 FloatComplex tmp = exp (static_cast<float> (M_PI) * alpha * eye)
1681 *
cbesh1 (z, alpha, kode, ierr);
1706 float rexpz =
real (expz);
1707 float iexpz =
imag (expz);
1709 float tmp_r =
real (y) * rexpz -
imag (y) * iexpz;
1710 float tmp_i =
real (y) * iexpz +
imag (y) * rexpz;
1723 FloatComplex tmp = exp (-static_cast<float> (M_PI) * alpha * eye)
1724 *
cbesh2 (z, alpha, kode, ierr);
1735 static inline FloatComplex
1739 FloatComplex retval;
1741 retval =
f (x, alpha, (scaled ? 2 : 1), ierr);
1759 retval(i,j) =
f (
x(i,j), alpha, (scaled ? 2 : 1),
ierr(i,j));
1766 const FloatComplex&
x,
1778 retval(i,j) =
f (x, alpha(i,j), (scaled ? 2 : 1),
ierr(i,j));
1796 if (x_nr == alpha_nr && x_nc == alpha_nc)
1807 retval(i,j) =
f (
x(i,j), alpha(i,j), (scaled ? 2 : 1),
ierr(i,j));
1811 (
"%s: the sizes of alpha and x must conform", fn);
1827 retval(i) =
f (
x(i), alpha, (scaled ? 2 : 1),
ierr(i));
1843 retval(i) =
f (x, alpha(i), (scaled ? 2 : 1),
ierr(i));
1856 if (dv == alpha.
dims ())
1864 retval(i) =
f (
x(i), alpha(i), (scaled ? 2 : 1),
ierr(i));
1868 (
"%s: the sizes of alpha and x must conform", fn);
1887 retval(i,j) =
f (
x(i), alpha(j), (scaled ? 2 : 1),
ierr(i,j));
1892 #define SS_BESSEL(name, fcn) \
1894 name (float alpha, const FloatComplex& x, bool scaled, octave_idx_type& ierr) \
1896 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
1899 #define SM_BESSEL(name, fcn) \
1900 FloatComplexMatrix \
1901 name (float alpha, const FloatComplexMatrix& x, bool scaled, \
1902 Array<octave_idx_type>& ierr) \
1904 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
1907 #define MS_BESSEL(name, fcn) \
1908 FloatComplexMatrix \
1909 name (const FloatMatrix& alpha, const FloatComplex& x, bool scaled, \
1910 Array<octave_idx_type>& ierr) \
1912 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
1915 #define MM_BESSEL(name, fcn) \
1916 FloatComplexMatrix \
1917 name (const FloatMatrix& alpha, const FloatComplexMatrix& x, bool scaled, \
1918 Array<octave_idx_type>& ierr) \
1920 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
1923 #define SN_BESSEL(name, fcn) \
1924 FloatComplexNDArray \
1925 name (float alpha, const FloatComplexNDArray& x, bool scaled, \
1926 Array<octave_idx_type>& ierr) \
1928 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
1931 #define NS_BESSEL(name, fcn) \
1932 FloatComplexNDArray \
1933 name (const FloatNDArray& alpha, const FloatComplex& x, bool scaled, \
1934 Array<octave_idx_type>& ierr) \
1936 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
1939 #define NN_BESSEL(name, fcn) \
1940 FloatComplexNDArray \
1941 name (const FloatNDArray& alpha, const FloatComplexNDArray& x, bool scaled, \
1942 Array<octave_idx_type>& ierr) \
1944 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
1947 #define RC_BESSEL(name, fcn) \
1948 FloatComplexMatrix \
1949 name (const FloatRowVector& alpha, const FloatComplexColumnVector& x, bool scaled, \
1950 Array<octave_idx_type>& ierr) \
1952 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
1955 #define ALL_BESSEL(name, fcn) \
1956 SS_BESSEL (name, fcn) \
1957 SM_BESSEL (name, fcn) \
1958 MS_BESSEL (name, fcn) \
1959 MM_BESSEL (name, fcn) \
1960 SN_BESSEL (name, fcn) \
1961 NS_BESSEL (name, fcn) \
1962 NN_BESSEL (name, fcn) \
1963 RC_BESSEL (name, fcn)
1990 double zr = z.real ();
1991 double zi = z.imag ();
1999 Complex expz = exp (- 2.0 / 3.0 * z * sqrt (z));
2001 double rexpz =
real (expz);
2002 double iexpz =
imag (expz);
2004 double tmp = ar*rexpz - ai*iexpz;
2006 ai = ar*iexpz + ai*rexpz;
2010 if (zi == 0.0 && (! scaled || zr >= 0.0))
2022 double zr = z.real ();
2023 double zi = z.imag ();
2033 double rexpz =
real (expz);
2034 double iexpz =
imag (expz);
2036 double tmp = ar*rexpz - ai*iexpz;
2038 ai = ar*iexpz + ai*rexpz;
2042 if (zi == 0.0 && (! scaled || zr >= 0.0))
2061 retval(i,j) =
airy (z(i,j), deriv, scaled,
ierr(i,j));
2079 retval(i,j) =
biry (z(i,j), deriv, scaled,
ierr(i,j));
2095 retval(i) =
airy (z(i), deriv, scaled,
ierr(i));
2111 retval(i) =
biry (z(i), deriv, scaled,
ierr(i));
2124 float zr = z.real ();
2125 float zi = z.imag ();
2133 FloatComplex expz = exp (- 2.0
f / 3.0
f * z * sqrt (z));
2135 float rexpz =
real (expz);
2136 float iexpz =
imag (expz);
2138 float tmp = ar*rexpz - ai*iexpz;
2140 ai = ar*iexpz + ai*rexpz;
2144 if (zi == 0.0 && (! scaled || zr >= 0.0))
2156 float zr = z.real ();
2157 float zi = z.imag ();
2165 FloatComplex expz = exp (
std::abs (
real (2.0
f / 3.0
f * z * sqrt (z))));
2167 float rexpz =
real (expz);
2168 float iexpz =
imag (expz);
2170 float tmp = ar*rexpz - ai*iexpz;
2172 ai = ar*iexpz + ai*rexpz;
2176 if (zi == 0.0 && (! scaled || zr >= 0.0))
2195 retval(i,j) =
airy (z(i,j), deriv, scaled,
ierr(i,j));
2213 retval(i,j) =
biry (z(i,j), deriv, scaled,
ierr(i,j));
2229 retval(i) =
airy (z(i), deriv, scaled,
ierr(i));
2245 retval(i) =
biry (z(i), deriv, scaled,
ierr(i));
2254 std::string d1_str = d1.
str ();
2255 std::string d2_str = d2.
str ();
2256 std::string d3_str = d3.
str ();
2258 (*current_liboctave_error_handler)
2259 (
"betainc: nonconformant arguments (x is %s, a is %s, b is %s)",
2260 d1_str.c_str (), d2_str.c_str (), d3_str.c_str ());
2267 std::string d1_str = d1.
str ();
2268 std::string d2_str = d2.
str ();
2269 std::string d3_str = d3.
str ();
2271 (*current_liboctave_error_handler)
2272 (
"betaincinv: nonconformant arguments (x is %s, a is %s, b is %s)",
2273 d1_str.c_str (), d2_str.c_str (), d3_str.c_str ());
2295 *pretval++ =
betainc (x, a, b(i));
2311 *pretval++ =
betainc (x, a(i), b);
2322 if (dv == b.
dims ())
2331 *pretval++ =
betainc (x, a(i), b(i));
2361 if (dv == b.
dims ())
2370 *pretval++ =
betainc (
x(i), a, b(i));
2384 if (dv == a.
dims ())
2393 *pretval++ =
betainc (
x(i), a(i), b);
2407 if (dv == a.
dims () && dv == b.
dims ())
2416 *pretval++ =
betainc (
x(i), a(i), b(i));
2443 *pretval++ =
betainc (x, a, b(i));
2459 *pretval++ =
betainc (x, a(i), b);
2470 if (dv == b.
dims ())
2479 *pretval++ =
betainc (x, a(i), b(i));
2509 if (dv == b.
dims ())
2518 *pretval++ =
betainc (
x(i), a, b(i));
2532 if (dv == a.
dims ())
2541 *pretval++ =
betainc (
x(i), a(i), b);
2555 if (dv == a.
dims () && dv == b.
dims ())
2564 *pretval++ =
betainc (
x(i), a(i), b(i));
2581 if (a < 0.0 || x < 0.0)
2582 (*current_liboctave_error_handler)
2583 (
"gammainc: A and X must be non-negative");
2604 result(i,j) =
gammainc (x, a(i,j), err);
2631 result(i,j) =
gammainc (
x(i,j), a, err);
2656 if (nr == a_nr && nc == a_nc)
2665 result(i,j) =
gammainc (
x(i,j), a(i,j), err);
2675 (
"gammainc: nonconformant arguments (arg 1 is %dx%d, arg 2 is %dx%d)",
2676 nr, nc, a_nr, a_nc);
2696 result(i) =
gammainc (x, a(i), err);
2744 if (dv == a.
dims ())
2762 std::string x_str = dv.
str ();
2763 std::string a_str = a.
dims ().
str ();
2765 (*current_liboctave_error_handler)
2766 (
"gammainc: nonconformant arguments (arg 1 is %s, arg 2 is %s)",
2767 x_str.c_str (), a_str. c_str ());
2782 if (a < 0.0 || x < 0.0)
2783 (*current_liboctave_error_handler)
2784 (
"gammainc: A and X must be non-negative");
2805 result(i,j) =
gammainc (x, a(i,j), err);
2832 result(i,j) =
gammainc (
x(i,j), a, err);
2857 if (nr == a_nr && nc == a_nc)
2866 result(i,j) =
gammainc (
x(i,j), a(i,j), err);
2876 (
"gammainc: nonconformant arguments (arg 1 is %dx%d, arg 2 is %dx%d)",
2877 nr, nc, a_nr, a_nc);
2897 result(i) =
gammainc (x, a(i), err);
2945 if (dv == a.
dims ())
2963 std::string x_str = dv.
str ();
2964 std::string a_str = a.
dims ().
str ();
2966 (*current_liboctave_error_handler)
2967 (
"gammainc: nonconformant arguments (arg 1 is %s, arg 2 is %s)",
2968 x_str.c_str (), a_str.c_str ());
2979 const double pi = 3.14159265358979323846;
2981 ?
Complex (gnulib::log (-(1.0 + x)), pi)
2987 const float pi = 3.14159265358979323846f;
2989 ? FloatComplex (gnulib::logf (-(1.0
f + x)), pi)
2990 : FloatComplex (
log1pf (x)));
3003 static const double a[] =
3005 -2.806989788730439e+01, 1.562324844726888e+02,
3006 -1.951109208597547e+02, 9.783370457507161e+01,
3007 -2.168328665628878e+01, 1.772453852905383e+00
3009 static const double b[] =
3011 -5.447609879822406e+01, 1.615858368580409e+02,
3012 -1.556989798598866e+02, 6.680131188771972e+01,
3013 -1.328068155288572e+01
3015 static const double c[] =
3017 -5.504751339936943e-03, -2.279687217114118e-01,
3018 -1.697592457770869e+00, -1.802933168781950e+00,
3019 3.093354679843505e+00, 2.077595676404383e+00
3021 static const double d[] =
3023 7.784695709041462e-03, 3.224671290700398e-01,
3024 2.445134137142996e+00, 3.754408661907416e+00
3027 static const double spi2 = 8.862269254527579e-01;
3028 static const double pbreak = 0.95150;
3029 double ax = fabs (x), y;
3035 const double q = 0.5 *
x, r = q*q;
3036 const double yn = (((((a[0]*r + a[1])*r + a[2])*r + a[3])*r + a[4])*r + a[5])*q;
3037 const double yd = ((((b[0]*r + b[1])*r + b[2])*r + b[3])*r + b[4])*r + 1.0;
3043 const double q = sqrt (-2*gnulib::log (0.5*(1-ax)));
3044 const double yn = ((((c[0]*q + c[1])*q + c[2])*q + c[3])*q + c[4])*q + c[5];
3045 const double yd = (((d[0]*q + d[1])*q + d[2])*q + d[3])*q + 1.0;
3046 y = yn / yd *
signum (-x);
3056 double u = (
erf (y) -
x) * spi2 * exp (y*y);
3080 static const double a[] =
3082 -2.806989788730439e+01, 1.562324844726888e+02,
3083 -1.951109208597547e+02, 9.783370457507161e+01,
3084 -2.168328665628878e+01, 1.772453852905383e+00
3086 static const double b[] =
3088 -5.447609879822406e+01, 1.615858368580409e+02,
3089 -1.556989798598866e+02, 6.680131188771972e+01,
3090 -1.328068155288572e+01
3092 static const double c[] =
3094 -5.504751339936943e-03, -2.279687217114118e-01,
3095 -1.697592457770869e+00, -1.802933168781950e+00,
3096 3.093354679843505e+00, 2.077595676404383e+00
3098 static const double d[] =
3100 7.784695709041462e-03, 3.224671290700398e-01,
3101 2.445134137142996e+00, 3.754408661907416e+00
3104 static const double spi2 = 8.862269254527579e-01;
3105 static const double pbreak_lo = 0.04850;
3106 static const double pbreak_hi = 1.95150;
3110 if (x >= pbreak_lo && x <= pbreak_hi)
3113 const double q = 0.5*(1-
x), r = q*q;
3114 const double yn = (((((a[0]*r + a[1])*r + a[2])*r + a[3])*r + a[4])*r + a[5])*q;
3115 const double yd = ((((b[0]*r + b[1])*r + b[2])*r + b[3])*r + b[4])*r + 1.0;
3118 else if (x > 0.0 && x < 2.0)
3121 const double q = (x < 1
3122 ? sqrt (-2*gnulib::log (0.5*x))
3123 : sqrt (-2*gnulib::log (0.5*(2-x))));
3125 const double yn = ((((c[0]*q + c[1])*q + c[2])*q + c[3])*q + c[4])*q + c[5];
3127 const double yd = (((d[0]*q + d[1])*q + d[2])*q + d[3])*q + 1.0;
3144 double u = (
erf (y) - (1-
x)) * spi2 * exp (y*y);
3178 betain (
double x,
double p,
double q,
double beta,
bool& err)
3180 double acu = 0.1E-14, ai, cx;
3183 double pp, psq, qq, rx, temp, term, value, xx;
3190 if ((p <= 0.0 || q <= 0.0) || (x < 0.0 || 1.0 < x))
3198 if (x == 0.0 || x == 1.0)
3227 ns =
static_cast<int> (qq + cx * psq);
3240 term = term * temp * rx / (pp + ai);
3241 value = value + term;
3244 if (temp <= acu && temp <= acu * value)
3246 value = value * exp (pp * gnulib::log (xx)
3247 + (qq - 1.0) * gnulib::log (cx) - beta) / pp;
3251 value = 1.0 - value;
3297 double a, acu, adj, fpu, g, h;
3300 double pp, prev, qq, r, s, sae = -37.0, sq, t, tx, value,
w, xin, ycur, yprev;
3304 fpu =
pow (10.0, sae);
3309 if (p <= 0.0 || q <= 0.0)
3311 (*current_liboctave_error_handler)
3312 (
"betaincinv: wrong parameters");
3315 if (y < 0.0 || 1.0 < y)
3317 (*current_liboctave_error_handler)
3318 (
"betaincinv: wrong parameter Y");
3321 if (y == 0.0 || y == 1.0)
3345 r = sqrt (- gnulib::log (a * a));
3347 ycur = r - (2.30753 + 0.27061 * r) / (1.0 + (0.99229 + 0.04481 * r) * r);
3349 if (1.0 < pp && 1.0 < qq)
3351 r = (ycur * ycur - 3.0) / 6.0;
3352 s = 1.0 / (pp + pp - 1.0);
3353 t = 1.0 / (qq + qq - 1.0);
3355 w = ycur * sqrt (h + r) / h - (t - s) * (r + 5.0 / 6.0 - 2.0 / (3.0 * h));
3356 value = pp / (pp + qq * exp (w + w));
3361 t = 1.0 / (9.0 * qq);
3362 t = r *
pow (1.0 - t + ycur * sqrt (t), 3);
3366 value = 1.0 - exp ((gnulib::log ((1.0 - a) * qq) + beta) / qq);
3370 t = (4.0 * pp + r - 2.0) / t;
3374 value = exp ((gnulib::log (a * pp) + beta) / pp);
3378 value = 1.0 - 2.0 / (t + 1.0);
3402 iex =
std::max (- 5.0 / pp / pp - 1.0 /
pow (a, 0.2) - 13.0, sae);
3404 acu =
pow (10.0, iex);
3408 ycur =
betain (value, pp, qq, beta, err);
3416 ycur = (ycur - a) * exp (beta + r * gnulib::log (xin)
3417 + t * gnulib::log (1.0 - xin));
3419 if (ycur * yprev <= 0.0)
3437 if (0.0 <= tx && tx <= 1.0)
3449 value = 1.0 - value;
3454 if (ycur * ycur <= acu)
3458 value = 1.0 - value;
3463 if (tx != 0.0 && tx != 1.0)
3482 value = 1.0 - value;
3526 if (dv == b.
dims ())
3565 if (dv == b.
dims ())
3588 if (dv == a.
dims ())
3612 if (dv == a.
dims () && dv == b.
dims ())
3630 ellipj (
double u,
double m,
double& sn,
double& cn,
double& dn,
double& err)
3632 static const int Nmax = 16;
3633 double m1, t=0, si_u, co_u, se_u, ta_u, b, c[Nmax], a[Nmax], phi;
3638 (*current_liboctave_warning_with_id_handler)
3639 (
"Octave:ellipj-invalid-m",
"ellipj: expecting 0 <= M <= 1");
3646 double sqrt_eps = sqrt (std::numeric_limits<double>::epsilon ());
3652 t = 0.25*m*(u - si_u*co_u);
3653 sn = si_u - t * co_u;
3654 cn = co_u + t * si_u;
3655 dn = 1 - 0.5*m*si_u*si_u;
3657 else if ((1 - m) < sqrt_eps)
3665 sn = ta_u + 0.25*m1*(si_u*co_u - u)*se_u*se_u;
3666 cn = se_u - 0.25*m1*(si_u*co_u - u)*ta_u*se_u;
3667 dn = se_u + 0.25*m1*(si_u*co_u + u)*ta_u*se_u;
3676 for (n = 1; n < Nmax; ++n)
3678 a[n] = (a[n - 1] + b)/2;
3679 c[n] = (a[n - 1] - b)/2;
3680 b = sqrt (a[n - 1]*b);
3681 if (c[n]/a[n] < std::numeric_limits<double>::epsilon ())
break;
3689 for (ii = 1; n > 0; ii = ii*2, --n) ;
3691 for (n = Nn; n > 0; --n)
3693 phi = (
asin ((c[n]/a[n])* sin (phi)) + phi)/2;
3697 dn = sqrt (1 - m*sn*sn);
3705 double m1 = 1 - m, ss1, cc1, dd1;
3718 double ss, cc, dd, ddd;
3721 ddd = cc1*cc1 + m*ss*ss*ss1*ss1;
3722 sn =
Complex (ss*dd1/ddd, cc*dd*ss1*cc1/ddd);
3723 cn =
Complex (cc*cc1/ddd, -ss*dd*ss1*dd1/ddd);
3724 dn =
Complex (dd*cc1*dd1/ddd, -m*ss*cc*ss1/ddd);
static bool is_integer_value(double x)
void resize(octave_idx_type nr, octave_idx_type nc, const FloatComplex &rfv=FloatComplex(0))
double betain(double x, double p, double q, double beta, bool &err)
static FloatComplex cbesj(const FloatComplex &z, float alpha, int kode, octave_idx_type &ierr)
Complex besselk(double alpha, const Complex &x, bool scaled, octave_idx_type &ierr)
void resize(octave_idx_type nr, octave_idx_type nc, float rfv=0)
subroutine xerf(x, result)
static Complex zbesj(const Complex &z, double alpha, int kode, octave_idx_type &ierr)
std::string str(char sep= 'x') const
static FloatComplex cbesi(const FloatComplex &z, float alpha, int kode, octave_idx_type &ierr)
subroutine dlgams(X, DLGAM, SGNGAM)
subroutine xdgamma(x, result)
subroutine cairy(Z, ID, KODE, AI, NZ, IERR)
void resize(octave_idx_type nr, octave_idx_type nc, double rfv=0)
Complex besseli(double alpha, const Complex &x, bool scaled, octave_idx_type &ierr)
subroutine zairy(ZR, ZI, ID, KODE, AIR, AII, NZ, IERR)
std::complex< double > erfi(std::complex< double > z, double relerr=0)
static double do_erfinv(double x, bool refine)
Complex(* dptr)(const Complex &, double, int, octave_idx_type &)
Complex airy(const Complex &z, bool deriv, bool scaled, octave_idx_type &ierr)
F77_RET_T const double const double double const octave_idx_type octave_idx_type * ierr
FloatComplex(* fptr)(const FloatComplex &, float, int, octave_idx_type &)
bool xnegative_sign(double x)
subroutine cbiry(Z, ID, KODE, BI, IERR)
std::complex< double > Dawson(std::complex< double > z, double relerr=0)
double lo_ieee_nan_value(void)
subroutine xderfc(x, result)
static void gripe_betaincinv_nonconformant(const dim_vector &d1, const dim_vector &d2, const dim_vector &d3)
double betaincinv(double y, double p, double q)
#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)
double gammainc(double x, double a, bool &err)
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.
F77_RET_T const double const double double * d
static Complex do_bessel(dptr f, const char *, double alpha, const Complex &x, bool scaled, octave_idx_type &ierr)
#define ALL_BESSEL(name, fcn)
const dim_vector & dims(void) const
Return a const-reference so that dims ()(i) works efficiently.
subroutine xdatanh(x, result)
Complex rc_log1p(double x)
liboctave_error_handler current_liboctave_error_handler
subroutine algams(X, ALGAM, SGNGAM)
Complex rc_lgamma(double x)
F77_RET_T const double const double * f
Complex besselh1(double alpha, const Complex &x, bool scaled, octave_idx_type &ierr)
std::complex< double > w(std::complex< double > z, double relerr=0)
subroutine zbiry(ZR, ZI, ID, KODE, BIR, BII, IERR)
subroutine xacosh(x, result)
void resize(const dim_vector &dv, const T &rfv)
subroutine xgammainc(a, x, result)
octave_int< T > pow(const octave_int< T > &a, const octave_int< T > &b)
static double do_erfcinv(double x, bool refine)
Complex biry(const Complex &z, bool deriv, bool scaled, octave_idx_type &ierr)
static void gripe_betainc_nonconformant(const dim_vector &d1, const dim_vector &d2, const dim_vector &d3)
std::complex< double > erfcx(std::complex< double > z, double relerr=0)
F77_RET_T const double const double const octave_idx_type const octave_idx_type double double octave_idx_type octave_idx_type &F77_RET_T const double const double const octave_idx_type const octave_idx_type double double octave_idx_type double double octave_idx_type &F77_RET_T const double const double const octave_idx_type const octave_idx_type double double octave_idx_type octave_idx_type &F77_RET_T const double const double const octave_idx_type const octave_idx_type double double octave_idx_type octave_idx_type &F77_RET_T const double const double const octave_idx_type const octave_idx_type const octave_idx_type double double octave_idx_type octave_idx_type &F77_RET_T const float const octave_idx_type const octave_idx_type FloatComplex octave_idx_type octave_idx_type &F77_RET_T const float const octave_idx_type const octave_idx_type FloatComplex octave_idx_type FloatComplex octave_idx_type &F77_RET_T const float const octave_idx_type const octave_idx_type FloatComplex octave_idx_type octave_idx_type &F77_RET_T const float const octave_idx_type const octave_idx_type FloatComplex octave_idx_type octave_idx_type &F77_RET_T const float const octave_idx_type const octave_idx_type const octave_idx_type FloatComplex octave_idx_type octave_idx_type &F77_RET_T const double const octave_idx_type const octave_idx_type double double octave_idx_type octave_idx_type &F77_RET_T const float const octave_idx_type const octave_idx_type float float octave_idx_type octave_idx_type &F77_RET_T const double const octave_idx_type const octave_idx_type double double octave_idx_type &F77_RET_T const float const octave_idx_type const octave_idx_type float float octave_idx_type &F77_RET_T double &F77_RET_T float &F77_RET_T double &F77_RET_T float &F77_RET_T double &F77_RET_T float &F77_RET_T double &F77_RET_T float &F77_RET_T double &F77_RET_T float &F77_RET_T const double const double double &F77_RET_T const float const float float &F77_RET_T double &F77_RET_T float &F77_RET_T const double double &F77_RET_T const float float &F77_RET_T double double &F77_RET_T float float &double acosh(double x)
subroutine zbesh(ZR, ZI, FNU, KODE, M, N, CYR, CYI, NZ, IERR)
OCTAVE_API double D_NINT(double x)
double betainc(double x, double a, double b)
static Complex zbesh2(const Complex &z, double alpha, int kode, octave_idx_type &ierr)
subroutine xdasinh(x, result)
static FloatComplex cbesh1(const FloatComplex &z, float alpha, int kode, octave_idx_type &ierr)
subroutine cbesh(Z, FNU, KODE, M, N, CY, NZ, IERR)
subroutine xdbetai(x, a, b, result)
charNDArray max(char d, const charNDArray &m)
octave_idx_type length(void) const
Number of elements in the array.
static FloatComplex cbesy(const FloatComplex &z, float alpha, int kode, octave_idx_type &ierr)
SparseMatrix atan2(const double &x, const SparseMatrix &y)
Complex besselj(double alpha, const Complex &x, bool scaled, octave_idx_type &ierr)
void resize(octave_idx_type nr, octave_idx_type nc, const Complex &rfv=Complex(0))
Complex bessely(double alpha, const Complex &x, bool scaled, octave_idx_type &ierr)
static Complex zbesi(const Complex &z, double alpha, int kode, octave_idx_type &ierr)
static FloatComplex cbesh2(const FloatComplex &z, float alpha, int kode, octave_idx_type &ierr)
subroutine xbetai(x, a, b, result)
void ellipj(double u, double m, double &sn, double &cn, double &dn, double &err)
subroutine xsgammainc(a, x, result)
Complex besselh2(double alpha, const Complex &x, bool scaled, octave_idx_type &ierr)
Complex asin(const Complex &x)
ColumnVector imag(const ComplexColumnVector &a)
std::complex< float > FloatComplex
static Complex bessel_return_value(const Complex &val, octave_idx_type ierr)
static Complex zbesh1(const Complex &z, double alpha, int kode, octave_idx_type &ierr)
std::complex< double > Complex
const T * fortran_vec(void) const
subroutine xatanh(x, result)
ColumnVector real(const ComplexColumnVector &a)
octave_idx_type cols(void) const
static Complex zbesk(const Complex &z, double alpha, int kode, octave_idx_type &ierr)
static Complex zbesy(const Complex &z, double alpha, int kode, octave_idx_type &ierr)
static FloatComplex cbesk(const FloatComplex &z, float alpha, int kode, octave_idx_type &ierr)
F77_RET_T const double * x
subroutine xasinh(x, result)
subroutine xderf(x, result)
F77_RET_T F77_FUNC(zbesj, ZBESJ)(const double &
std::complex< double > erfc(std::complex< double > z, double relerr=0)