Mercurial > octave-nkf
changeset 14979:8b4606e3be32
maint: periodic merge of default to jit
| author | Max Brister <max@2bass.com> |
|---|---|
| date | Tue, 03 Jul 2012 11:05:16 -0500 |
| parents | f649b66ef1af (current diff) 619fedc6ea61 (diff) |
| children | bb1f3a9bb122 |
| files | |
| diffstat | 18 files changed, 1000 insertions(+), 692 deletions(-) [+] |
line wrap: on
line diff
--- a/NEWS Wed Jun 27 23:43:06 2012 -0500 +++ b/NEWS Tue Jul 03 11:05:16 2012 -0500 @@ -74,9 +74,10 @@ ** Other new functions added in 3.8.0: - colorcube lines splinefit - erfcinv rgbplot tetramesh - findfigs shrinkfaces + betaincinv lines tetramesh + colorcube rgbplot + erfcinv shrinkfaces + findfigs splinefit ** Deprecated functions.
--- a/doc/interpreter/arith.txi Wed Jun 27 23:43:06 2012 -0500 +++ b/doc/interpreter/arith.txi Tue Jul 03 11:05:16 2012 -0500 @@ -255,6 +255,8 @@ @DOCSTRING(betainc) +@DOCSTRING(betaincinv) + @DOCSTRING(betaln) @DOCSTRING(bincoeff)
--- a/doc/interpreter/oop.txi Wed Jun 27 23:43:06 2012 -0500 +++ b/doc/interpreter/oop.txi Tue Jul 03 11:05:16 2012 -0500 @@ -586,7 +586,7 @@ @item @tab a >= b @tab ge (a, b) @tab Greater than or equal to operator @tab @item @tab a == b @tab eq (a, b) @tab Equal to operator @tab @item @tab a != b @tab ne (a, b) @tab Not equal to operator @tab -@item @tab a \& b @tab and (a, b) @tab Logical and operator @tab +@item @tab a & b @tab and (a, b) @tab Logical and operator @tab @item @tab a | b @tab or (a, b) @tab Logical or operator @tab @item @tab ! b @tab not (a) @tab Logical not operator @tab @item @tab a' @tab ctranspose (a) @tab Complex conjugate transpose operator @tab
--- a/liboctave/lo-ieee.h Wed Jun 27 23:43:06 2012 -0500 +++ b/liboctave/lo-ieee.h Tue Jul 03 11:05:16 2012 -0500 @@ -118,4 +118,26 @@ #define lo_ieee_signbit(x) (sizeof (x) == sizeof (float) ? \ __lo_ieee_float_signbit (x) : __lo_ieee_signbit (x)) +#ifdef __cplusplus + +template <typename T> +struct octave_numeric_limits +{ + static T NA (void) { return static_cast<T> (0); } +}; + +template <> +struct octave_numeric_limits<double> +{ + static double NA (void) { return octave_NA; } +}; + +template <> +struct octave_numeric_limits<float> +{ + static float NA (void) { return octave_Float_NA; } +}; + #endif + +#endif
--- a/liboctave/lo-specfun.cc Wed Jun 27 23:43:06 2012 -0500 +++ b/liboctave/lo-specfun.cc Tue Jul 03 11:05:16 2012 -0500 @@ -615,7 +615,7 @@ { FloatComplex retval; - float r = x.real (), i = x.imag(); + float r = x.real (), i = x.imag (); if (fabs (r) < 0.5 && fabs (i) < 0.5) { @@ -2158,6 +2158,28 @@ d1_str.c_str (), d2_str.c_str (), d3_str.c_str ()); } +static void +gripe_betaincinv_nonconformant (octave_idx_type r1, octave_idx_type c1, octave_idx_type r2, octave_idx_type c2, octave_idx_type r3, + octave_idx_type c3) +{ + (*current_liboctave_error_handler) + ("betaincinv: nonconformant arguments (x is %dx%d, a is %dx%d, b is %dx%d)", + r1, c1, r2, c2, r3, c3); +} + +static void +gripe_betaincinv_nonconformant (const dim_vector& d1, const dim_vector& d2, + const dim_vector& d3) +{ + std::string d1_str = d1.str (); + std::string d2_str = d2.str (); + std::string d3_str = d3.str (); + + (*current_liboctave_error_handler) + ("betaincinv: nonconformant arguments (x is %s, a is %s, b is %s)", + d1_str.c_str (), d2_str.c_str (), d3_str.c_str ()); +} + double betainc (double x, double a, double b) { @@ -2166,93 +2188,42 @@ return retval; } -Matrix -betainc (double x, double a, const Matrix& b) +Array<double> +betainc (double x, double a, const Array<double>& b) { - octave_idx_type nr = b.rows (); - octave_idx_type nc = b.cols (); - - Matrix retval (nr, nc); - - for (octave_idx_type j = 0; j < nc; j++) - for (octave_idx_type i = 0; i < nr; i++) - retval(i,j) = betainc (x, a, b(i,j)); - - return retval; -} - -Matrix -betainc (double x, const Matrix& a, double b) -{ - octave_idx_type nr = a.rows (); - octave_idx_type nc = a.cols (); - - Matrix retval (nr, nc); - - for (octave_idx_type j = 0; j < nc; j++) - for (octave_idx_type i = 0; i < nr; i++) - retval(i,j) = betainc (x, a(i,j), b); + dim_vector dv = b.dims (); + octave_idx_type nel = dv.numel (); + + Array<double> retval (dv); + + double *pretval = retval.fortran_vec (); + + for (octave_idx_type i = 0; i < nel; i++) + *pretval++ = betainc (x, a, b(i)); return retval; } -Matrix -betainc (double x, const Matrix& a, const Matrix& b) +Array<double> +betainc (double x, const Array<double>& a, double b) { - Matrix retval; - - octave_idx_type a_nr = a.rows (); - octave_idx_type a_nc = a.cols (); - - octave_idx_type b_nr = b.rows (); - octave_idx_type b_nc = b.cols (); - - if (a_nr == b_nr && a_nc == b_nc) - { - retval.resize (a_nr, a_nc); - - for (octave_idx_type j = 0; j < a_nc; j++) - for (octave_idx_type i = 0; i < a_nr; i++) - retval(i,j) = betainc (x, a(i,j), b(i,j)); - } - else - gripe_betainc_nonconformant (1, 1, a_nr, a_nc, b_nr, b_nc); + dim_vector dv = a.dims (); + octave_idx_type nel = dv.numel (); + + Array<double> retval (dv); + + double *pretval = retval.fortran_vec (); + + for (octave_idx_type i = 0; i < nel; i++) + *pretval++ = betainc (x, a(i), b); return retval; } -NDArray -betainc (double x, double a, const NDArray& b) +Array<double> +betainc (double x, const Array<double>& a, const Array<double>& b) { - dim_vector dv = b.dims (); - octave_idx_type nel = dv.numel (); - - NDArray retval (dv); - - for (octave_idx_type i = 0; i < nel; i++) - retval (i) = betainc (x, a, b(i)); - - return retval; -} - -NDArray -betainc (double x, const NDArray& a, double b) -{ - dim_vector dv = a.dims (); - octave_idx_type nel = dv.numel (); - - NDArray retval (dv); - - for (octave_idx_type i = 0; i < nel; i++) - retval (i) = betainc (x, a(i), b); - - return retval; -} - -NDArray -betainc (double x, const NDArray& a, const NDArray& b) -{ - NDArray retval; + Array<double> retval; dim_vector dv = a.dims (); if (dv == b.dims ()) @@ -2261,8 +2232,10 @@ retval.resize (dv); + double *pretval = retval.fortran_vec (); + for (octave_idx_type i = 0; i < nel; i++) - retval (i) = betainc (x, a(i), b(i)); + *pretval++ = betainc (x, a(i), b(i)); } else gripe_betainc_nonconformant (dim_vector (0, 0), dv, b.dims ()); @@ -2270,118 +2243,26 @@ return retval; } - -Matrix -betainc (const Matrix& x, double a, double b) -{ - octave_idx_type nr = x.rows (); - octave_idx_type nc = x.cols (); - - Matrix retval (nr, nc); - - for (octave_idx_type j = 0; j < nc; j++) - for (octave_idx_type i = 0; i < nr; i++) - retval(i,j) = betainc (x(i,j), a, b); - - return retval; -} - -Matrix -betainc (const Matrix& x, double a, const Matrix& b) +Array<double> +betainc (const Array<double>& x, double a, double b) { - Matrix retval; - - octave_idx_type nr = x.rows (); - octave_idx_type nc = x.cols (); - - octave_idx_type b_nr = b.rows (); - octave_idx_type b_nc = b.cols (); - - if (nr == b_nr && nc == b_nc) - { - retval.resize (nr, nc); - - for (octave_idx_type j = 0; j < nc; j++) - for (octave_idx_type i = 0; i < nr; i++) - retval(i,j) = betainc (x(i,j), a, b(i,j)); - } - else - gripe_betainc_nonconformant (nr, nc, 1, 1, b_nr, b_nc); + dim_vector dv = x.dims (); + octave_idx_type nel = dv.numel (); + + Array<double> retval (dv); + + double *pretval = retval.fortran_vec (); + + for (octave_idx_type i = 0; i < nel; i++) + *pretval++ = betainc (x(i), a, b); return retval; } -Matrix -betainc (const Matrix& x, const Matrix& a, double b) -{ - Matrix retval; - - octave_idx_type nr = x.rows (); - octave_idx_type nc = x.cols (); - - octave_idx_type a_nr = a.rows (); - octave_idx_type a_nc = a.cols (); - - if (nr == a_nr && nc == a_nc) - { - retval.resize (nr, nc); - - for (octave_idx_type j = 0; j < nc; j++) - for (octave_idx_type i = 0; i < nr; i++) - retval(i,j) = betainc (x(i,j), a(i,j), b); - } - else - gripe_betainc_nonconformant (nr, nc, a_nr, a_nc, 1, 1); - - return retval; -} - -Matrix -betainc (const Matrix& x, const Matrix& a, const Matrix& b) +Array<double> +betainc (const Array<double>& x, double a, const Array<double>& b) { - Matrix retval; - - octave_idx_type nr = x.rows (); - octave_idx_type nc = x.cols (); - - octave_idx_type a_nr = a.rows (); - octave_idx_type a_nc = a.cols (); - - octave_idx_type b_nr = b.rows (); - octave_idx_type b_nc = b.cols (); - - if (nr == a_nr && nr == b_nr && nc == a_nc && nc == b_nc) - { - retval.resize (nr, nc); - - for (octave_idx_type j = 0; j < nc; j++) - for (octave_idx_type i = 0; i < nr; i++) - retval(i,j) = betainc (x(i,j), a(i,j), b(i,j)); - } - else - gripe_betainc_nonconformant (nr, nc, a_nr, a_nc, b_nr, b_nc); - - return retval; -} - -NDArray -betainc (const NDArray& x, double a, double b) -{ - dim_vector dv = x.dims (); - octave_idx_type nel = dv.numel (); - - NDArray retval (dv); - - for (octave_idx_type i = 0; i < nel; i++) - retval (i) = betainc (x(i), a, b); - - return retval; -} - -NDArray -betainc (const NDArray& x, double a, const NDArray& b) -{ - NDArray retval; + Array<double> retval; dim_vector dv = x.dims (); if (dv == b.dims ()) @@ -2390,8 +2271,10 @@ retval.resize (dv); + double *pretval = retval.fortran_vec (); + for (octave_idx_type i = 0; i < nel; i++) - retval (i) = betainc (x(i), a, b(i)); + *pretval++ = betainc (x(i), a, b(i)); } else gripe_betainc_nonconformant (dv, dim_vector (0, 0), b.dims ()); @@ -2399,10 +2282,10 @@ return retval; } -NDArray -betainc (const NDArray& x, const NDArray& a, double b) +Array<double> +betainc (const Array<double>& x, const Array<double>& a, double b) { - NDArray retval; + Array<double> retval; dim_vector dv = x.dims (); if (dv == a.dims ()) @@ -2411,8 +2294,10 @@ retval.resize (dv); + double *pretval = retval.fortran_vec (); + for (octave_idx_type i = 0; i < nel; i++) - retval (i) = betainc (x(i), a(i), b); + *pretval++ = betainc (x(i), a(i), b); } else gripe_betainc_nonconformant (dv, a.dims (), dim_vector (0, 0)); @@ -2420,10 +2305,10 @@ return retval; } -NDArray -betainc (const NDArray& x, const NDArray& a, const NDArray& b) +Array<double> +betainc (const Array<double>& x, const Array<double>& a, const Array<double>& b) { - NDArray retval; + Array<double> retval; dim_vector dv = x.dims (); if (dv == a.dims () && dv == b.dims ()) @@ -2432,8 +2317,10 @@ retval.resize (dv); + double *pretval = retval.fortran_vec (); + for (octave_idx_type i = 0; i < nel; i++) - retval (i) = betainc (x(i), a(i), b(i)); + *pretval++ = betainc (x(i), a(i), b(i)); } else gripe_betainc_nonconformant (dv, a.dims (), b.dims ()); @@ -2449,93 +2336,42 @@ return retval; } -FloatMatrix -betainc (float x, float a, const FloatMatrix& b) +Array<float> +betainc (float x, float a, const Array<float>& b) { - octave_idx_type nr = b.rows (); - octave_idx_type nc = b.cols (); - - FloatMatrix retval (nr, nc); - - for (octave_idx_type j = 0; j < nc; j++) - for (octave_idx_type i = 0; i < nr; i++) - retval(i,j) = betainc (x, a, b(i,j)); - - return retval; -} - -FloatMatrix -betainc (float x, const FloatMatrix& a, float b) -{ - octave_idx_type nr = a.rows (); - octave_idx_type nc = a.cols (); - - FloatMatrix retval (nr, nc); - - for (octave_idx_type j = 0; j < nc; j++) - for (octave_idx_type i = 0; i < nr; i++) - retval(i,j) = betainc (x, a(i,j), b); + dim_vector dv = b.dims (); + octave_idx_type nel = dv.numel (); + + Array<float> retval (dv); + + float *pretval = retval.fortran_vec (); + + for (octave_idx_type i = 0; i < nel; i++) + *pretval++ = betainc (x, a, b(i)); return retval; } -FloatMatrix -betainc (float x, const FloatMatrix& a, const FloatMatrix& b) +Array<float> +betainc (float x, const Array<float>& a, float b) { - FloatMatrix retval; - - octave_idx_type a_nr = a.rows (); - octave_idx_type a_nc = a.cols (); - - octave_idx_type b_nr = b.rows (); - octave_idx_type b_nc = b.cols (); - - if (a_nr == b_nr && a_nc == b_nc) - { - retval.resize (a_nr, a_nc); - - for (octave_idx_type j = 0; j < a_nc; j++) - for (octave_idx_type i = 0; i < a_nr; i++) - retval(i,j) = betainc (x, a(i,j), b(i,j)); - } - else - gripe_betainc_nonconformant (1, 1, a_nr, a_nc, b_nr, b_nc); + dim_vector dv = a.dims (); + octave_idx_type nel = dv.numel (); + + Array<float> retval (dv); + + float *pretval = retval.fortran_vec (); + + for (octave_idx_type i = 0; i < nel; i++) + *pretval++ = betainc (x, a(i), b); return retval; } -FloatNDArray -betainc (float x, float a, const FloatNDArray& b) +Array<float> +betainc (float x, const Array<float>& a, const Array<float>& b) { - dim_vector dv = b.dims (); - octave_idx_type nel = dv.numel (); - - FloatNDArray retval (dv); - - for (octave_idx_type i = 0; i < nel; i++) - retval (i) = betainc (x, a, b(i)); - - return retval; -} - -FloatNDArray -betainc (float x, const FloatNDArray& a, float b) -{ - dim_vector dv = a.dims (); - octave_idx_type nel = dv.numel (); - - FloatNDArray retval (dv); - - for (octave_idx_type i = 0; i < nel; i++) - retval (i) = betainc (x, a(i), b); - - return retval; -} - -FloatNDArray -betainc (float x, const FloatNDArray& a, const FloatNDArray& b) -{ - FloatNDArray retval; + Array<float> retval; dim_vector dv = a.dims (); if (dv == b.dims ()) @@ -2544,8 +2380,10 @@ retval.resize (dv); + float *pretval = retval.fortran_vec (); + for (octave_idx_type i = 0; i < nel; i++) - retval (i) = betainc (x, a(i), b(i)); + *pretval++ = betainc (x, a(i), b(i)); } else gripe_betainc_nonconformant (dim_vector (0, 0), dv, b.dims ()); @@ -2553,118 +2391,26 @@ return retval; } - -FloatMatrix -betainc (const FloatMatrix& x, float a, float b) -{ - octave_idx_type nr = x.rows (); - octave_idx_type nc = x.cols (); - - FloatMatrix retval (nr, nc); - - for (octave_idx_type j = 0; j < nc; j++) - for (octave_idx_type i = 0; i < nr; i++) - retval(i,j) = betainc (x(i,j), a, b); - - return retval; -} - -FloatMatrix -betainc (const FloatMatrix& x, float a, const FloatMatrix& b) +Array<float> +betainc (const Array<float>& x, float a, float b) { - FloatMatrix retval; - - octave_idx_type nr = x.rows (); - octave_idx_type nc = x.cols (); - - octave_idx_type b_nr = b.rows (); - octave_idx_type b_nc = b.cols (); - - if (nr == b_nr && nc == b_nc) - { - retval.resize (nr, nc); - - for (octave_idx_type j = 0; j < nc; j++) - for (octave_idx_type i = 0; i < nr; i++) - retval(i,j) = betainc (x(i,j), a, b(i,j)); - } - else - gripe_betainc_nonconformant (nr, nc, 1, 1, b_nr, b_nc); + dim_vector dv = x.dims (); + octave_idx_type nel = dv.numel (); + + Array<float> retval (dv); + + float *pretval = retval.fortran_vec (); + + for (octave_idx_type i = 0; i < nel; i++) + *pretval++ = betainc (x(i), a, b); return retval; } -FloatMatrix -betainc (const FloatMatrix& x, const FloatMatrix& a, float b) -{ - FloatMatrix retval; - - octave_idx_type nr = x.rows (); - octave_idx_type nc = x.cols (); - - octave_idx_type a_nr = a.rows (); - octave_idx_type a_nc = a.cols (); - - if (nr == a_nr && nc == a_nc) - { - retval.resize (nr, nc); - - for (octave_idx_type j = 0; j < nc; j++) - for (octave_idx_type i = 0; i < nr; i++) - retval(i,j) = betainc (x(i,j), a(i,j), b); - } - else - gripe_betainc_nonconformant (nr, nc, a_nr, a_nc, 1, 1); - - return retval; -} - -FloatMatrix -betainc (const FloatMatrix& x, const FloatMatrix& a, const FloatMatrix& b) +Array<float> +betainc (const Array<float>& x, float a, const Array<float>& b) { - FloatMatrix retval; - - octave_idx_type nr = x.rows (); - octave_idx_type nc = x.cols (); - - octave_idx_type a_nr = a.rows (); - octave_idx_type a_nc = a.cols (); - - octave_idx_type b_nr = b.rows (); - octave_idx_type b_nc = b.cols (); - - if (nr == a_nr && nr == b_nr && nc == a_nc && nc == b_nc) - { - retval.resize (nr, nc); - - for (octave_idx_type j = 0; j < nc; j++) - for (octave_idx_type i = 0; i < nr; i++) - retval(i,j) = betainc (x(i,j), a(i,j), b(i,j)); - } - else - gripe_betainc_nonconformant (nr, nc, a_nr, a_nc, b_nr, b_nc); - - return retval; -} - -FloatNDArray -betainc (const FloatNDArray& x, float a, float b) -{ - dim_vector dv = x.dims (); - octave_idx_type nel = dv.numel (); - - FloatNDArray retval (dv); - - for (octave_idx_type i = 0; i < nel; i++) - retval (i) = betainc (x(i), a, b); - - return retval; -} - -FloatNDArray -betainc (const FloatNDArray& x, float a, const FloatNDArray& b) -{ - FloatNDArray retval; + Array<float> retval; dim_vector dv = x.dims (); if (dv == b.dims ()) @@ -2673,8 +2419,10 @@ retval.resize (dv); + float *pretval = retval.fortran_vec (); + for (octave_idx_type i = 0; i < nel; i++) - retval (i) = betainc (x(i), a, b(i)); + *pretval++ = betainc (x(i), a, b(i)); } else gripe_betainc_nonconformant (dv, dim_vector (0, 0), b.dims ()); @@ -2682,10 +2430,10 @@ return retval; } -FloatNDArray -betainc (const FloatNDArray& x, const FloatNDArray& a, float b) +Array<float> +betainc (const Array<float>& x, const Array<float>& a, float b) { - FloatNDArray retval; + Array<float> retval; dim_vector dv = x.dims (); if (dv == a.dims ()) @@ -2694,8 +2442,10 @@ retval.resize (dv); + float *pretval = retval.fortran_vec (); + for (octave_idx_type i = 0; i < nel; i++) - retval (i) = betainc (x(i), a(i), b); + *pretval++ = betainc (x(i), a(i), b); } else gripe_betainc_nonconformant (dv, a.dims (), dim_vector (0, 0)); @@ -2703,10 +2453,10 @@ return retval; } -FloatNDArray -betainc (const FloatNDArray& x, const FloatNDArray& a, const FloatNDArray& b) +Array<float> +betainc (const Array<float>& x, const Array<float>& a, const Array<float>& b) { - FloatNDArray retval; + Array<float> retval; dim_vector dv = x.dims (); if (dv == a.dims () && dv == b.dims ()) @@ -2715,8 +2465,10 @@ retval.resize (dv); + float *pretval = retval.fortran_vec (); + for (octave_idx_type i = 0; i < nel; i++) - retval (i) = betainc (x(i), a(i), b(i)); + *pretval++ = betainc (x(i), a(i), b(i)); } else gripe_betainc_nonconformant (dv, a.dims (), b.dims ()); @@ -3128,7 +2880,7 @@ (*current_liboctave_error_handler) ("gammainc: nonconformant arguments (arg 1 is %s, arg 2 is %s)", - x_str.c_str (), a_str. c_str ()); + x_str.c_str (), a_str.c_str ()); } done: @@ -3394,3 +3146,468 @@ { return erfcx_impl (x); } + +// +// Incomplete Beta function ratio +// +// Algorithm based on the one by John Burkardt. +// See http://people.sc.fsu.edu/~jburkardt/cpp_src/asa109/asa109.html +// +// The original code is distributed under the GNU LGPL v3 license. +// +// Reference: +// +// KL Majumder, GP Bhattacharjee, +// Algorithm AS 63: +// The incomplete Beta Integral, +// Applied Statistics, +// Volume 22, Number 3, 1973, pages 409-411. +// +double +betain (double x, double p, double q, double beta, bool& err) +{ + double acu = 0.1E-14, ai, cx; + bool indx; + int ns; + double pp, psq, qq, rx, temp, term, value, xx; + + value = x; + err = false; + + // Check the input arguments. + + if ((p <= 0.0 || q <= 0.0) || (x < 0.0 || 1.0 < x)) + { + err = true; + return value; + } + + // Special cases. + + if (x == 0.0 || x == 1.0) + { + return value; + } + + // Change tail if necessary and determine S. + + psq = p + q; + cx = 1.0 - x; + + if (p < psq * x) + { + xx = cx; + cx = x; + pp = q; + qq = p; + indx = true; + } + else + { + xx = x; + pp = p; + qq = q; + indx = false; + } + + term = 1.0; + ai = 1.0; + value = 1.0; + ns = (int) (qq + cx * psq); + + // Use the Soper reduction formula. + + rx = xx / cx; + temp = qq - ai; + if (ns == 0) + { + rx = xx; + } + + for ( ; ; ) + { + term = term * temp * rx / (pp + ai); + value = value + term; + temp = fabs (term); + + if (temp <= acu && temp <= acu * value) + { + value = value * exp (pp * log (xx) + + (qq - 1.0) * log (cx) - beta) / pp; + + if (indx) + { + value = 1.0 - value; + } + break; + } + + ai = ai + 1.0; + ns = ns - 1; + + if (0 <= ns) + { + temp = qq - ai; + if (ns == 0) + { + rx = xx; + } + } + else + { + temp = psq; + psq = psq + 1.0; + } + } + + return value; +} + +// +// Inverse of the incomplete Beta function +// +// Algorithm based on the one by John Burkardt. +// See http://people.sc.fsu.edu/~jburkardt/cpp_src/asa109/asa109.html +// +// The original code is distributed under the GNU LGPL v3 license. +// +// Reference: +// +// GW Cran, KJ Martin, GE Thomas, +// Remark AS R19 and Algorithm AS 109: +// A Remark on Algorithms AS 63: The Incomplete Beta Integral +// and AS 64: Inverse of the Incomplete Beta Integeral, +// Applied Statistics, +// Volume 26, Number 1, 1977, pages 111-114. +// +double +betaincinv (double y, double p, double q) { + double a, acu, adj, fpu, g, h; + int iex; + bool indx; + double pp, prev, qq, r, s, sae = -37.0, sq, t, tx, value, w, xin, ycur, yprev; + + double beta = lgamma (p) + lgamma (q) - lgamma (p + q); + bool err = false; + fpu = pow (10.0, sae); + value = y; + + // Test for admissibility of parameters. + + if (p <= 0.0 || q <= 0.0) + { + (*current_liboctave_error_handler) + ("betaincinv: wrong parameters"); + } + + if (y < 0.0 || 1.0 < y) + { + (*current_liboctave_error_handler) + ("betaincinv: wrong parameter Y"); + } + + if (y == 0.0 || y == 1.0) + { + return value; + } + + // Change tail if necessary. + + if (0.5 < y) + { + a = 1.0 - y; + pp = q; + qq = p; + indx = true; + } + else + { + a = y; + pp = p; + qq = q; + indx = false; + } + + // Calculate the initial approximation. + + r = sqrt (- log (a * a)); + + ycur = r - (2.30753 + 0.27061 * r) / (1.0 + (0.99229 + 0.04481 * r) * r); + + if (1.0 < pp && 1.0 < qq) + { + r = (ycur * ycur - 3.0) / 6.0; + s = 1.0 / (pp + pp - 1.0); + t = 1.0 / (qq + qq - 1.0); + h = 2.0 / (s + t); + w = ycur * sqrt (h + r) / h - (t - s) * (r + 5.0 / 6.0 - 2.0 / (3.0 * h)); + value = pp / (pp + qq * exp (w + w)); + } + else + { + r = qq + qq; + t = 1.0 / (9.0 * qq); + t = r * pow (1.0 - t + ycur * sqrt (t), 3); + + if (t <= 0.0) + { + value = 1.0 - exp ((log ((1.0 - a) * qq) + beta) / qq); + } + else + { + t = (4.0 * pp + r - 2.0) / t; + + if (t <= 1.0) + { + value = exp ((log (a * pp) + beta) / pp); + } + else + { + value = 1.0 - 2.0 / (t + 1.0); + } + } + } + + // Solve for X by a modified Newton-Raphson method, + // using the function BETAIN. + + r = 1.0 - pp; + t = 1.0 - qq; + yprev = 0.0; + sq = 1.0; + prev = 1.0; + + if (value < 0.0001) + { + value = 0.0001; + } + + if (0.9999 < value) + { + value = 0.9999; + } + + iex = std::max (- 5.0 / pp / pp - 1.0 / pow (a, 0.2) - 13.0, sae); + + acu = pow (10.0, iex); + + for ( ; ; ) + { + ycur = betain (value, pp, qq, beta, err); + + if (err) + { + return value; + } + + xin = value; + ycur = (ycur - a) * exp (beta + r * log (xin) + t * log (1.0 - xin)); + + if (ycur * yprev <= 0.0) + { + prev = std::max (sq, fpu); + } + + g = 1.0; + + for ( ; ; ) + { + for ( ; ; ) + { + adj = g * ycur; + sq = adj * adj; + + if (sq < prev) + { + tx = value - adj; + + if (0.0 <= tx && tx <= 1.0) + { + break; + } + } + g = g / 3.0; + } + + if (prev <= acu) + { + if (indx) + { + value = 1.0 - value; + } + return value; + } + + if (ycur * ycur <= acu) + { + if (indx) + { + value = 1.0 - value; + } + return value; + } + + if (tx != 0.0 && tx != 1.0) + { + break; + } + + g = g / 3.0; + } + + if (tx == value) + { + break; + } + + value = tx; + yprev = ycur; + } + + if (indx) + { + value = 1.0 - value; + } + + return value; +} + +Array<double> +betaincinv (double x, double a, const Array<double>& b) +{ + dim_vector dv = b.dims (); + octave_idx_type nel = dv.numel (); + + Array<double> retval (dv); + + double *pretval = retval.fortran_vec (); + + for (octave_idx_type i = 0; i < nel; i++) + *pretval++ = betaincinv (x, a, b(i)); + + return retval; +} + +Array<double> +betaincinv (double x, const Array<double>& a, double b) +{ + dim_vector dv = a.dims (); + octave_idx_type nel = dv.numel (); + + Array<double> retval (dv); + + double *pretval = retval.fortran_vec (); + + for (octave_idx_type i = 0; i < nel; i++) + *pretval++ = betaincinv (x, a(i), b); + + return retval; +} + +Array<double> +betaincinv (double x, const Array<double>& a, const Array<double>& b) +{ + Array<double> retval; + dim_vector dv = a.dims (); + + if (dv == b.dims ()) + { + octave_idx_type nel = dv.numel (); + + retval.resize (dv); + + double *pretval = retval.fortran_vec (); + + for (octave_idx_type i = 0; i < nel; i++) + *pretval++ = betaincinv (x, a(i), b(i)); + } + else + gripe_betaincinv_nonconformant (dim_vector (0, 0), dv, b.dims ()); + + return retval; +} + +Array<double> +betaincinv (const Array<double>& x, double a, double b) +{ + dim_vector dv = x.dims (); + octave_idx_type nel = dv.numel (); + + Array<double> retval (dv); + + double *pretval = retval.fortran_vec (); + + for (octave_idx_type i = 0; i < nel; i++) + *pretval++ = betaincinv (x(i), a, b); + + return retval; +} + +Array<double> +betaincinv (const Array<double>& x, double a, const Array<double>& b) +{ + Array<double> retval; + dim_vector dv = x.dims (); + + if (dv == b.dims ()) + { + octave_idx_type nel = dv.numel (); + + retval.resize (dv); + + double *pretval = retval.fortran_vec (); + + for (octave_idx_type i = 0; i < nel; i++) + *pretval++ = betaincinv (x(i), a, b(i)); + } + else + gripe_betaincinv_nonconformant (dv, dim_vector (0, 0), b.dims ()); + + return retval; +} + +Array<double> +betaincinv (const Array<double>& x, const Array<double>& a, double b) +{ + Array<double> retval; + dim_vector dv = x.dims (); + + if (dv == a.dims ()) + { + octave_idx_type nel = dv.numel (); + + retval.resize (dv); + + double *pretval = retval.fortran_vec (); + + for (octave_idx_type i = 0; i < nel; i++) + *pretval++ = betaincinv (x(i), a(i), b); + } + else + gripe_betaincinv_nonconformant (dv, a.dims (), dim_vector (0, 0)); + + return retval; +} + +Array<double> +betaincinv (const Array<double>& x, const Array<double>& a, const Array<double>& b) +{ + Array<double> retval; + dim_vector dv = x.dims (); + + if (dv == a.dims () && dv == b.dims ()) + { + octave_idx_type nel = dv.numel (); + + retval.resize (dv); + + double *pretval = retval.fortran_vec (); + + for (octave_idx_type i = 0; i < nel; i++) + *pretval++ = betaincinv (x(i), a(i), b(i)); + } + else + gripe_betaincinv_nonconformant (dv, a.dims (), b.dims ()); + + return retval; +}
--- a/liboctave/lo-specfun.h Wed Jun 27 23:43:06 2012 -0500 +++ b/liboctave/lo-specfun.h Tue Jul 03 11:05:16 2012 -0500 @@ -520,42 +520,24 @@ biry (const FloatComplexNDArray& z, bool deriv, bool scaled, Array<octave_idx_type>& ierr); extern OCTAVE_API double betainc (double x, double a, double b); -extern OCTAVE_API Matrix betainc (double x, double a, const Matrix& b); -extern OCTAVE_API Matrix betainc (double x, const Matrix& a, double b); -extern OCTAVE_API Matrix betainc (double x, const Matrix& a, const Matrix& b); - -extern OCTAVE_API NDArray betainc (double x, double a, const NDArray& b); -extern OCTAVE_API NDArray betainc (double x, const NDArray& a, double b); -extern OCTAVE_API NDArray betainc (double x, const NDArray& a, const NDArray& b); - -extern OCTAVE_API Matrix betainc (const Matrix& x, double a, double b); -extern OCTAVE_API Matrix betainc (const Matrix& x, double a, const Matrix& b); -extern OCTAVE_API Matrix betainc (const Matrix& x, const Matrix& a, double b); -extern OCTAVE_API Matrix betainc (const Matrix& x, const Matrix& a, const Matrix& b); - -extern OCTAVE_API NDArray betainc (const NDArray& x, double a, double b); -extern OCTAVE_API NDArray betainc (const NDArray& x, double a, const NDArray& b); -extern OCTAVE_API NDArray betainc (const NDArray& x, const NDArray& a, double b); -extern OCTAVE_API NDArray betainc (const NDArray& x, const NDArray& a, const NDArray& b); +extern OCTAVE_API Array<double> betainc (double x, double a, const Array<double>& b); +extern OCTAVE_API Array<double> betainc (double x, const Array<double>& a, double b); +extern OCTAVE_API Array<double> betainc (double x, const Array<double>& a, const Array<double>& b); +extern OCTAVE_API Array<double> betainc (const Array<double>& x, double a, double b); +extern OCTAVE_API Array<double> betainc (const Array<double>& x, double a, double b); +extern OCTAVE_API Array<double> betainc (const Array<double>& x, double a, const Array<double>& b); +extern OCTAVE_API Array<double> betainc (const Array<double>& x, const Array<double>& a, double b); +extern OCTAVE_API Array<double> betainc (const Array<double>& x, const Array<double>& a, const Array<double>& b); extern OCTAVE_API float betainc (float x, float a, float b); -extern OCTAVE_API FloatMatrix betainc (float x, float a, const FloatMatrix& b); -extern OCTAVE_API FloatMatrix betainc (float x, const FloatMatrix& a, float b); -extern OCTAVE_API FloatMatrix betainc (float x, const FloatMatrix& a, const FloatMatrix& b); - -extern OCTAVE_API FloatNDArray betainc (float x, float a, const FloatNDArray& b); -extern OCTAVE_API FloatNDArray betainc (float x, const FloatNDArray& a, float b); -extern OCTAVE_API FloatNDArray betainc (float x, const FloatNDArray& a, const FloatNDArray& b); - -extern OCTAVE_API FloatMatrix betainc (const FloatMatrix& x, float a, float b); -extern OCTAVE_API FloatMatrix betainc (const FloatMatrix& x, float a, const FloatMatrix& b); -extern OCTAVE_API FloatMatrix betainc (const FloatMatrix& x, const FloatMatrix& a, float b); -extern OCTAVE_API FloatMatrix betainc (const FloatMatrix& x, const FloatMatrix& a, const FloatMatrix& b); - -extern OCTAVE_API FloatNDArray betainc (const FloatNDArray& x, float a, float b); -extern OCTAVE_API FloatNDArray betainc (const FloatNDArray& x, float a, const FloatNDArray& b); -extern OCTAVE_API FloatNDArray betainc (const FloatNDArray& x, const FloatNDArray& a, float b); -extern OCTAVE_API FloatNDArray betainc (const FloatNDArray& x, const FloatNDArray& a, const FloatNDArray& b); +extern OCTAVE_API Array<float> betainc (float x, float a, const Array<float>& b); +extern OCTAVE_API Array<float> betainc (float x, const Array<float>& a, float b); +extern OCTAVE_API Array<float> betainc (float x, const Array<float>& a, const Array<float>& b); +extern OCTAVE_API Array<float> betainc (const Array<float>& x, float a, float b); +extern OCTAVE_API Array<float> betainc (const Array<float>& x, float a, float b); +extern OCTAVE_API Array<float> betainc (const Array<float>& x, float a, const Array<float>& b); +extern OCTAVE_API Array<float> betainc (const Array<float>& x, const Array<float>& a, float b); +extern OCTAVE_API Array<float> betainc (const Array<float>& x, const Array<float>& a, const Array<float>& b); extern OCTAVE_API double gammainc (double x, double a, bool& err); extern OCTAVE_API Matrix gammainc (double x, const Matrix& a); @@ -599,4 +581,14 @@ extern OCTAVE_API double erfcx (double x); extern OCTAVE_API float erfcx (float x); +extern OCTAVE_API double betaincinv (double x, double a, double b); +extern OCTAVE_API Array<double> betaincinv (double x, double a, const Array<double>& b); +extern OCTAVE_API Array<double> betaincinv (double x, const Array<double>& a, double b); +extern OCTAVE_API Array<double> betaincinv (double x, const Array<double>& a, const Array<double>& b); +extern OCTAVE_API Array<double> betaincinv (const Array<double>& x, double a, double b); +extern OCTAVE_API Array<double> betaincinv (const Array<double>& x, double a, double b); +extern OCTAVE_API Array<double> betaincinv (const Array<double>& x, double a, const Array<double>& b); +extern OCTAVE_API Array<double> betaincinv (const Array<double>& x, const Array<double>& a, double b); +extern OCTAVE_API Array<double> betaincinv (const Array<double>& x, const Array<double>& a, const Array<double>& b); + #endif
--- a/liboctave/lo-utils.cc Wed Jun 27 23:43:06 2012 -0500 +++ b/liboctave/lo-utils.cc Tue Jul 03 11:05:16 2012 -0500 @@ -195,10 +195,14 @@ return retval; } -static inline double -read_inf_nan_na (std::istream& is, char c0, char sign = '+') +// Note that the caller is responsible for repositioning the stream on +// failure. + +template <typename T> +T +read_inf_nan_na (std::istream& is, char c0) { - double d = 0.0; + T val = 0.0; switch (c0) { @@ -209,21 +213,12 @@ { char c2 = is.get (); if (c2 == 'f' || c2 == 'F') - d = sign == '-' ? -octave_Inf : octave_Inf; + val = std::numeric_limits<T>::infinity (); else - { - is.putback (c2); - is.putback (c1); - is.putback (c0); - is.setstate (std::ios::failbit); - } + is.setstate (std::ios::failbit); } else - { - is.putback (c1); - is.putback (c0); - is.setstate (std::ios::failbit); - } + is.setstate (std::ios::failbit); } break; @@ -234,19 +229,12 @@ { char c2 = is.get (); if (c2 == 'n' || c2 == 'N') - d = octave_NaN; + val = std::numeric_limits<T>::quiet_NaN (); else - { - is.putback (c2); - d = octave_NA; - } + val = octave_numeric_limits<T>::NA (); } else - { - is.putback (c1); - is.putback (c0); - is.setstate (std::ios::failbit); - } + is.setstate (std::ios::failbit); } break; @@ -254,74 +242,80 @@ abort (); } - return d; + return val; } // Read a double value. Discard any sign on NaN and NA. -template <> +template <typename T> double -octave_read_value (std::istream& is) +octave_read_fp_value (std::istream& is) { - double d = 0.0; + T val = 0.0; + + // FIXME -- resetting stream position is likely to fail unless we are + // reading from a file. + std::ios::streampos pos = is.tellg (); char c1 = ' '; while (isspace (c1)) c1 = is.get (); + bool neg = false; + switch (c1) { case '-': - { - char c2 = 0; - c2 = is.get (); - if (c2 == 'i' || c2 == 'I' || c2 == 'n' || c2 == 'N') - d = read_inf_nan_na (is, c2, c1); - else - { - is.putback (c2); - is.putback (c1); - is >> d; - } - } - break; + neg = true; + // fall through... case '+': { char c2 = 0; c2 = is.get (); if (c2 == 'i' || c2 == 'I' || c2 == 'n' || c2 == 'N') - d = read_inf_nan_na (is, c2, c1); + val = read_inf_nan_na<T> (is, c2); else { is.putback (c2); - is.putback (c1); - is >> d; + is >> val; } + + if (neg && ! is.fail ()) + val = -val; } break; case 'i': case 'I': case 'n': case 'N': - d = read_inf_nan_na (is, c1); + val = read_inf_nan_na<T> (is, c1); break; default: is.putback (c1); - is >> d; + is >> val; + break; } - return d; + std::ios::iostate status = is.rdstate (); + if (status & std::ios::failbit) + { + is.clear (); + is.seekg (pos); + is.setstate (status); + } + + return val; } -template <> -Complex -octave_read_value (std::istream& is) +template <typename T> +std::complex<T> +octave_read_cx_fp_value (std::istream& is) { - double re = 0.0, im = 0.0; + T re = 0.0, im = 0.0; - Complex cx = 0.0; + std::complex<T> cx = 0.0; char ch = ' '; @@ -330,16 +324,16 @@ if (ch == '(') { - re = octave_read_value<double> (is); + re = octave_read_value<T> (is); ch = is.get (); if (ch == ',') { - im = octave_read_value<double> (is); + im = octave_read_value<T> (is); ch = is.get (); if (ch == ')') - cx = Complex (re, im); + cx = std::complex<T> (re, im); else is.setstate (std::ios::failbit); } @@ -355,170 +349,26 @@ } return cx; - } -static inline float -read_float_inf_nan_na (std::istream& is, char c0, char sign = '+') +template <> OCTAVE_API double octave_read_value (std::istream& is) { - float d = 0.0; - - switch (c0) - { - case 'i': case 'I': - { - char c1 = is.get (); - if (c1 == 'n' || c1 == 'N') - { - char c2 = is.get (); - if (c2 == 'f' || c2 == 'F') - d = sign == '-' ? -octave_Float_Inf : octave_Float_Inf; - else - { - is.putback (c2); - is.putback (c1); - is.putback (c0); - is.setstate (std::ios::failbit); - } - } - else - { - is.putback (c1); - is.putback (c0); - is.setstate (std::ios::failbit); - } - } - break; - - case 'n': case 'N': - { - char c1 = is.get (); - if (c1 == 'a' || c1 == 'A') - { - char c2 = is.get (); - if (c2 == 'n' || c2 == 'N') - d = octave_Float_NaN; - else - { - is.putback (c2); - d = octave_Float_NA; - } - } - else - { - is.putback (c1); - is.putback (c0); - is.setstate (std::ios::failbit); - } - } - break; - - default: - abort (); - } - - return d; + return octave_read_fp_value<double> (is); } -// Read a float value. Discard any sign on NaN and NA. - -template <> -float -octave_read_value (std::istream& is) +template <> OCTAVE_API Complex octave_read_value (std::istream& is) { - float d = 0.0; - - char c1 = ' '; - - while (isspace (c1)) - c1 = is.get (); - - switch (c1) - { - case '-': - { - char c2 = 0; - c2 = is.get (); - if (c2 == 'i' || c2 == 'I' || c2 == 'n' || c2 == 'N') - d = read_float_inf_nan_na (is, c2, c1); - else - { - is.putback (c2); - is.putback (c1); - is >> d; - } - } - break; - - case '+': - { - char c2 = 0; - c2 = is.get (); - if (c2 == 'i' || c2 == 'I' || c2 == 'n' || c2 == 'N') - d = read_float_inf_nan_na (is, c2, c1); - else - { - is.putback (c2); - is.putback (c1); - is >> d; - } - } - break; - - case 'i': case 'I': - case 'n': case 'N': - d = read_float_inf_nan_na (is, c1); - break; - - default: - is.putback (c1); - is >> d; - } - - return d; + return octave_read_cx_fp_value<double> (is); } -template <> -FloatComplex -octave_read_value (std::istream& is) +template <> OCTAVE_API float octave_read_value (std::istream& is) { - float re = 0.0, im = 0.0; - - FloatComplex cx = 0.0; - - char ch = ' '; - - while (isspace (ch)) - ch = is.get (); - - if (ch == '(') - { - re = octave_read_value<float> (is); - ch = is.get (); + return octave_read_fp_value<float> (is); +} - if (ch == ',') - { - im = octave_read_value<float> (is); - ch = is.get (); - - if (ch == ')') - cx = FloatComplex (re, im); - else - is.setstate (std::ios::failbit); - } - else if (ch == ')') - cx = re; - else - is.setstate (std::ios::failbit); - } - else - { - is.putback (ch); - cx = octave_read_value<float> (is); - } - - return cx; - +template <> OCTAVE_API FloatComplex octave_read_value (std::istream& is) +{ + return octave_read_cx_fp_value<float> (is); } void
--- a/scripts/help/unimplemented.m Wed Jun 27 23:43:06 2012 -0500 +++ b/scripts/help/unimplemented.m Tue Jul 03 11:05:16 2012 -0500 @@ -98,7 +98,6 @@ "bar3", "bar3h", "bench", - "betaincinv", "bicgstabl", "brush", "builddocsearchdb",
--- a/scripts/io/strread.m Wed Jun 27 23:43:06 2012 -0500 +++ b/scripts/io/strread.m Tue Jul 03 11:05:16 2012 -0500 @@ -301,6 +301,13 @@ ## Format conversion specifiers following literals w/o space/delim ## in between are separate now. Separate those w trailing literals idy2 = find (! cellfun ("isempty", strfind (fmt_words, "%"))); + + ## Check for unsupported format specifiers + errpat = '(\[.*\]|[cq]|[nfdu]8|[nfdu]16|[nfdu]32|[nfdu]64)'; + if (! all (cellfun ("isempty", regexp (fmt_words(idy2), errpat)))) + error ("strread: %q, %c, %[] or bit width format specifiers are not supported yet."); + endif + a = strfind (fmt_words(idy2), "%"); b = regexp (fmt_words(idy2), '[nfdus]', 'end'); for jj = 1:numel (a) @@ -931,3 +938,19 @@ %! assert (isempty (b)); %! assert (isempty (c)); +%% Unsupported format specifiers +%!test +%!error <format specifiers are not supported> strread ('a', '%c') +%!error <format specifiers are not supported> strread ('a', '%*c %d') +%!error <format specifiers are not supported> strread ('a', '%q') +%!error <format specifiers are not supported> strread ('a', '%*q %d') +%!error <format specifiers are not supported> strread ('a', '%[a]') +%!error <format specifiers are not supported> strread ('a', '%*[a] %d') +%!error <format specifiers are not supported> strread ('a', '%[^a]') +%!error <format specifiers are not supported> strread ('a', '%*[â] %d') +%!error <format specifiers are not supported> strread ('a', '%d8') +%!error <format specifiers are not supported> strread ('a', '%*d8 %s') +%!error <format specifiers are not supported> strread ('a', '%f64') +%!error <format specifiers are not supported> strread ('a', '%*f64 %s') +%!error <format specifiers are not supported> strread ('a', '%u32') +%!error <format specifiers are not supported> strread ('a', '%*u32 %d')
--- a/scripts/miscellaneous/perl.m Wed Jun 27 23:43:06 2012 -0500 +++ b/scripts/miscellaneous/perl.m Tue Jul 03 11:05:16 2012 -0500 @@ -17,13 +17,14 @@ ## <http://www.gnu.org/licenses/>. ## -*- texinfo -*- -## @deftypefn {Function File} {[@var{output}, @var{status}] =} perl (@var{scriptfile}) -## @deftypefnx {Function File} {[@var{output}, @var{status}] =} perl (@var{scriptfile}, @var{argument1}, @var{argument2}, @dots{}) -## Invoke Perl script @var{scriptfile} with possibly a list of -## command line arguments. -## Returns output in @var{output} and status -## in @var{status}. -## @seealso{system} +## @deftypefn {Function File} {@var{output} =} perl (@var{scriptfile}) +## @deftypefnx {Function File} {@var{output} =} perl (@var{scriptfile}, @var{argument1}, @var{argument2}, @dots{}) +## @deftypefnx {Function File} {[@var{output}, @var{status}] =} perl (@dots{}) +## Invoke Perl script @var{scriptfile}, possibly with a list of command line +## arguments. Return output in @var{output} and optional status in +## @var{status}. If @var{scriptfile} is not an absolute file name it is +## is searched for in the current directory and then in the Octave loadpath. +## @seealso{system, python} ## @end deftypefn function [output, status] = perl (scriptfile = "-e ''", varargin) @@ -33,8 +34,13 @@ ## initial value in the argument list of the function definition. if (ischar (scriptfile) - && ((nargin != 1 && iscellstr (varargin)) - || (nargin == 1 && ! isempty (scriptfile)))) + && ( (nargin == 1 && ! isempty (scriptfile)) + || (nargin != 1 && iscellstr (varargin)))) + if (! strcmp (scriptfile(1:2), "-e")) + ## Attempt to find file in loadpath. No effect for absolute filenames. + scriptfile = file_in_loadpath (scriptfile); + endif + [status, output] = system (cstrcat ("perl ", scriptfile, sprintf (" %s", varargin{:}))); else
--- a/scripts/miscellaneous/python.m Wed Jun 27 23:43:06 2012 -0500 +++ b/scripts/miscellaneous/python.m Tue Jul 03 11:05:16 2012 -0500 @@ -18,26 +18,28 @@ ## <http://www.gnu.org/licenses/>. ## -*- texinfo -*- -## @deftypefn {Function File} {[@var{output}, @var{status}] =} python (@var{scriptfile}) -## @deftypefnx {Function File} {[@var{output}, @var{status}] =} python (@var{scriptfile}, @var{argument1}, @var{argument2}, @dots{}) -## Invoke python script @var{scriptfile} with possibly a list of -## command line arguments. -## Returns output in @var{output} and status -## in @var{status}. -## @seealso{system} +## @deftypefn {Function File} {@var{output} =} python (@var{scriptfile}) +## @deftypefnx {Function File} {@var{output} =} python (@var{scriptfile}, @var{argument1}, @var{argument2}, @dots{}) +## @deftypefnx {Function File} {[@var{output}, @var{status}] =} python (@dots{}) +## Invoke Python script @var{scriptfile}, possibly with a list of command line +## arguments. Return output in @var{output} and optional status in +## @var{status}. If @var{scriptfile} is not an absolute file name it is +## is searched for in the current directory and then in the Octave loadpath. +## @seealso{system, perl} ## @end deftypefn ## Author: Carnë Draug <carandraug+dev@gmail.com> function [output, status] = python (scriptfile = "-c ''", varargin) - ## VARARGIN is intialized to {}(1x0) if no additional arguments are - ## supplied, so there is no need to check for it, or provide an - ## initial value in the argument list of the function definition. + if (ischar (scriptfile) + && ( (nargin == 1 && ! isempty (scriptfile)) + || (nargin != 1 && iscellstr (varargin)))) + if (! strcmp (scriptfile(1:2), "-c")) + ## Attempt to find file in loadpath. No effect for absolute filenames. + scriptfile = file_in_loadpath (scriptfile); + endif - if (ischar (scriptfile) - && ((nargin != 1 && iscellstr (varargin)) - || (nargin == 1 && ! isempty (scriptfile)))) [status, output] = system (cstrcat ("python ", scriptfile, sprintf (" %s", varargin{:}))); else
--- a/scripts/polynomial/ppval.m Wed Jun 27 23:43:06 2012 -0500 +++ b/scripts/polynomial/ppval.m Tue Jul 03 11:05:16 2012 -0500 @@ -68,7 +68,7 @@ ndv = length (dimvec); ## Offsets. - dx = (xi - x(idx)); + dx = (xi - x(idx))(:)'; dx = repmat (dx, [prod(d), 1]); dx = reshape (dx, dimvec); dx = shiftdim (dx, ndv - 1); @@ -119,4 +119,11 @@ %!assert (ppval (pp2, xi), [1.1 1.3 1.9 1.1;1.1 1.3 1.9 1.1], abserr) %!assert (ppval (pp2, xi'), [1.1 1.3 1.9 1.1;1.1 1.3 1.9 1.1], abserr) %!assert (size (ppval (pp2, [xi;xi])), [2 2 4]) - +%!test +%! breaks = [0, 1, 2, 3]; +%! coefs = rand (6, 4); +%! pp = mkpp (breaks, coefs, 2); +%! ret = zeros (2, 4, 2); +%! ret(:,:,1) = ppval (pp, breaks'); +%! ret(:,:,2) = ppval (pp, breaks'); +%! assert (ppval (pp, [breaks',breaks']), ret)
--- a/scripts/sparse/gmres.m Wed Jun 27 23:43:06 2012 -0500 +++ b/scripts/sparse/gmres.m Tue Jul 03 11:05:16 2012 -0500 @@ -174,7 +174,7 @@ x = x_old + V(:, 1:restart_it) * Y(1:restart_it); - resvec(iter) = presn; + resvec(iter+1) = presn; if (norm (x - x_old, inf) <= eps) flag = 3; # Stagnation: no change between iterations break; @@ -191,8 +191,7 @@ flag = 0; # Converged to solution within tolerance endif - resvec = resvec(1:iter-1); - it = [ceil(iter / restart), rem(iter, restart)]; + it = [floor(iter/restart), restart_it-1]; endif endfunction
--- a/scripts/specfun/beta.m Wed Jun 27 23:43:06 2012 -0500 +++ b/scripts/specfun/beta.m Tue Jul 03 11:05:16 2012 -0500 @@ -31,6 +31,7 @@ ## @end example ## ## @end ifnottex +## @seealso{betaln, betainc} ## @end deftypefn ## Author: KH <Kurt.Hornik@wu-wien.ac.at>
--- a/scripts/strings/strcat.m Wed Jun 27 23:43:06 2012 -0500 +++ b/scripts/strings/strcat.m Tue Jul 03 11:05:16 2012 -0500 @@ -23,11 +23,35 @@ ## horizontally. If the arguments are cells strings, @code{strcat} ## returns a cell string with the individual cells concatenated. ## For numerical input, each element is converted to the -## corresponding ASCII character. Trailing white space is eliminated. +## corresponding ASCII character. Trailing white space for each of +## the inputs (@var{s1}, @var{S2}, @dots{}) is eliminated before they +## are concatenated. +## ## For example: ## ## @example ## @group +## strcat ("|", " leading space is preserved", "|") +## @result{} | leading space is perserved| +## @end group +## @end example +## +## @example +## @group +## strcat ("|", "trailing space is eliminated ", "|") +## @result{} |trailing space is eliminated| +## @end group +## @end example +## +## @example +## @group +## strcat ("homogeneous space |", " ", "| is also eliminated") +## @result{} homogeneous space || is also eliminated +## @end group +## @end example +## +## @example +## @group ## s = [ "ab"; "cde" ]; ## strcat (s, s, s) ## @result{}
--- a/src/DLD-FUNCTIONS/betainc.cc Wed Jun 27 23:43:06 2012 -0500 +++ b/src/DLD-FUNCTIONS/betainc.cc Tue Jul 03 11:05:16 2012 -0500 @@ -32,6 +32,9 @@ #include "oct-obj.h" #include "utils.h" +// FIXME: These functions do not need to be dynamically loaded. They should +// be placed elsewhere in the Octave code hierarchy. + DEFUN_DLD (betainc, args, , "-*- texinfo -*-\n\ @deftypefn {Mapping Function} {} betainc (@var{x}, @var{a}, @var{b})\n\ @@ -94,7 +97,7 @@ } else { - FloatNDArray b = b_arg.float_array_value (); + Array<float> b = b_arg.float_array_value (); if (! error_state) retval = betainc (x, a, b); @@ -103,7 +106,7 @@ } else { - FloatNDArray a = a_arg.float_array_value (); + Array<float> a = a_arg.float_array_value (); if (! error_state) { @@ -116,7 +119,7 @@ } else { - FloatNDArray b = b_arg.float_array_value (); + Array<float> b = b_arg.float_array_value (); if (! error_state) retval = betainc (x, a, b); @@ -126,7 +129,7 @@ } else { - FloatNDArray x = x_arg.float_array_value (); + Array<float> x = x_arg.float_array_value (); if (a_arg.is_scalar_type ()) { @@ -143,7 +146,7 @@ } else { - FloatNDArray b = b_arg.float_array_value (); + Array<float> b = b_arg.float_array_value (); if (! error_state) retval = betainc (x, a, b); @@ -152,7 +155,7 @@ } else { - FloatNDArray a = a_arg.float_array_value (); + Array<float> a = a_arg.float_array_value (); if (! error_state) { @@ -165,7 +168,7 @@ } else { - FloatNDArray b = b_arg.float_array_value (); + Array<float> b = b_arg.float_array_value (); if (! error_state) retval = betainc (x, a, b); @@ -195,7 +198,7 @@ } else { - NDArray b = b_arg.array_value (); + Array<double> b = b_arg.array_value (); if (! error_state) retval = betainc (x, a, b); @@ -204,7 +207,7 @@ } else { - NDArray a = a_arg.array_value (); + Array<double> a = a_arg.array_value (); if (! error_state) { @@ -217,7 +220,7 @@ } else { - NDArray b = b_arg.array_value (); + Array<double> b = b_arg.array_value (); if (! error_state) retval = betainc (x, a, b); @@ -227,7 +230,7 @@ } else { - NDArray x = x_arg.array_value (); + Array<double> x = x_arg.array_value (); if (a_arg.is_scalar_type ()) { @@ -244,7 +247,7 @@ } else { - NDArray b = b_arg.array_value (); + Array<double> b = b_arg.array_value (); if (! error_state) retval = betainc (x, a, b); @@ -253,7 +256,7 @@ } else { - NDArray a = a_arg.array_value (); + Array<double> a = a_arg.array_value (); if (! error_state) { @@ -266,7 +269,7 @@ } else { - NDArray b = b_arg.array_value (); + Array<double> b = b_arg.array_value (); if (! error_state) retval = betainc (x, a, b); @@ -283,7 +286,7 @@ } /* -%% test/octave.test/arith/betainc-1.m +## Double precision %!test %! a = [1, 1.5, 2, 3]; %! b = [4, 3, 2, 1]; @@ -295,7 +298,7 @@ %! assert (v1, v2, sqrt (eps)); %! assert (v3, v4, sqrt (eps)); -%% Single precision +## Single precision %!test %! a = single ([1, 1.5, 2, 3]); %! b = single ([4, 3, 2, 1]); @@ -307,7 +310,7 @@ %! assert (v1, v2, sqrt (eps ("single"))); %! assert (v3, v4, sqrt (eps ("single"))); -%% Mixed double/single precision +## Mixed double/single precision %!test %! a = single ([1, 1.5, 2, 3]); %! b = [4, 3, 2, 1]; @@ -316,15 +319,172 @@ %! x = [.2, .4, .6, .8]; %! v3 = betainc (x, a, b); %! v4 = 1-betainc (1.-x, b, a); -%! assert (v1, v2, sqrt (eps ('single'))); -%! assert (v3, v4, sqrt (eps ('single'))); +%! assert (v1, v2, sqrt (eps ("single"))); +%! assert (v3, v4, sqrt (eps ("single"))); + +%!error betainc () +%!error betainc (1) +%!error betainc (1,2) +%!error betainc (1,2,3,4) +*/ + +DEFUN_DLD (betaincinv, args, , + "-*- texinfo -*-\n\ +@deftypefn {Mapping Function} {} betaincinv (@var{y}, @var{a}, @var{b})\n\ +Compute the inverse of the incomplete Beta function, i.e., @var{x} such that\n\ +\n\ +@example\n\ +@var{y} == betainc (@var{x}, @var{a}, @var{b}) \n\ +@end example\n\ +@seealso{betainc, beta, betaln}\n\ +@end deftypefn") +{ + octave_value retval; + + int nargin = args.length (); + + if (nargin == 3) + { + octave_value x_arg = args(0); + octave_value a_arg = args(1); + octave_value b_arg = args(2); + + if (x_arg.is_scalar_type ()) + { + double x = x_arg.double_value (); + + if (a_arg.is_scalar_type ()) + { + double a = a_arg.double_value (); -%% test/octave.test/arith/betainc-2.m -%!error betainc () + if (! error_state) + { + if (b_arg.is_scalar_type ()) + { + double b = b_arg.double_value (); + + if (! error_state) + retval = betaincinv (x, a, b); + } + else + { + Array<double> b = b_arg.array_value (); + + if (! error_state) + retval = betaincinv (x, a, b); + } + } + } + else + { + Array<double> a = a_arg.array_value (); + + if (! error_state) + { + if (b_arg.is_scalar_type ()) + { + double b = b_arg.double_value (); + + if (! error_state) + retval = betaincinv (x, a, b); + } + else + { + Array<double> b = b_arg.array_value (); + + if (! error_state) + retval = betaincinv (x, a, b); + } + } + } + } + else + { + Array<double> x = x_arg.array_value (); -%% test/octave.test/arith/betainc-3.m -%!error betainc> betainc (1) + if (a_arg.is_scalar_type ()) + { + double a = a_arg.double_value (); + + if (! error_state) + { + if (b_arg.is_scalar_type ()) + { + double b = b_arg.double_value (); + + if (! error_state) + retval = betaincinv (x, a, b); + } + else + { + Array<double> b = b_arg.array_value (); + + if (! error_state) + retval = betaincinv (x, a, b); + } + } + } + else + { + Array<double> a = a_arg.array_value (); + + if (! error_state) + { + if (b_arg.is_scalar_type ()) + { + double b = b_arg.double_value (); + + if (! error_state) + retval = betaincinv (x, a, b); + } + else + { + Array<double> b = b_arg.array_value (); -%% test/octave.test/arith/betainc-4.m -%!error betainc> betainc (1,2) + if (! error_state) + retval = betaincinv (x, a, b); + } + } + } + } + + // FIXME: It would be better to have an algorithm for betaincinv which accepted + // float inputs and returned float outputs. As it is, we do extra work + // to calculate betaincinv to double precision and then throw that precision + // away. + if (x_arg.is_single_type () || a_arg.is_single_type () || + b_arg.is_single_type ()) + { + retval = Array<float> (retval.array_value ()); + } + } + else + print_usage (); + + return retval; +} + +/* +%!assert (betaincinv ([0.875 0.6875], [1 2], 3), [0.5 0.5], sqrt (eps)) +%!assert (betaincinv (0.5, 3, 3), 0.5, sqrt (eps)) +%!assert (betaincinv (0.34375, 4, 3), 0.5, sqrt (eps)) +%!assert (betaincinv (0.2265625, 5, 3), 0.5, sqrt (eps)) +%!assert (betaincinv (0.14453125, 6, 3), 0.5, sqrt (eps)) +%!assert (betaincinv (0.08984375, 7, 3), 0.5, sqrt (eps)) +%!assert (betaincinv (0.0546875, 8, 3), 0.5, sqrt (eps)) +%!assert (betaincinv (0.03271484375, 9, 3), 0.5, sqrt (eps)) +%!assert (betaincinv (0.019287109375, 10, 3), 0.5, sqrt (eps)) + +## Test class single as well +%!assert (betaincinv ([0.875 0.6875], [1 2], single (3)), [0.5 0.5], sqrt (eps ("single"))) +%!assert (betaincinv (0.5, 3, single (3)), 0.5, sqrt (eps ("single"))) +%!assert (betaincinv (0.34375, 4, single (3)), 0.5, sqrt (eps ("single"))) + +## Extreme values +%!assert (betaincinv (0, 42, 42), 0, sqrt (eps)) +%!assert (betaincinv (1, 42, 42), 1, sqrt (eps)) + +%!error betaincinv () +%!error betaincinv (1, 2) */ +
--- a/src/syscalls.cc Wed Jun 27 23:43:06 2012 -0500 +++ b/src/syscalls.cc Tue Jul 03 11:05:16 2012 -0500 @@ -1594,8 +1594,10 @@ DEFUNX ("canonicalize_file_name", Fcanonicalize_file_name, args, , "-*- texinfo -*-\n\ -@deftypefn {Built-in Function} {[@var{cname}, @var{status}, @var{msg}]} canonicalize_file_name (@var{name})\n\ -Return the canonical name of file @var{name}.\n\ +@deftypefn {Built-in Function} {[@var{cname}, @var{status}, @var{msg}] =} canonicalize_file_name (@var{fname})\n\ +Return the canonical name of file @var{fname}. If the file does not exist\n\ +the empty string (\"\") is returned.\n\ +@seealso{make_absolute_filename, is_absolute_filename, is_rooted_relative_filename}\n\ @end deftypefn") { octave_value_list retval;
--- a/src/utils.cc Wed Jun 27 23:43:06 2012 -0500 +++ b/src/utils.cc Tue Jul 03 11:05:16 2012 -0500 @@ -862,8 +862,9 @@ DEFUN (make_absolute_filename, args, , "-*- texinfo -*-\n\ @deftypefn {Built-in Function} {} make_absolute_filename (@var{file})\n\ -Return the full name of @var{file}, relative to the current directory.\n\ -@seealso{is_absolute_filename, is_rooted_relative_filename, isdir}\n\ +Return the full name of @var{file} beginning from the root of the file\n\ +system. No check is done for the existence of @var{file}.\n\ +@seealso{canonicalize_file_name, is_absolute_filename, is_rooted_relative_filename, isdir}\n\ @end deftypefn") { octave_value retval = std::string ();