Mercurial > octave-nkf
changeset 7667:fb3a6c53c2b2
Allow negative zero imaginary part to be treated as zero for erf, erfc, gamma and lgamma mapper function
| author | David Bateman <dbateman@free.fr> |
|---|---|
| date | Fri, 28 Mar 2008 18:26:29 +0100 |
| parents | 693ac94c2854 |
| children | 7ae60e7eb6f8 |
| files | src/ChangeLog src/mappers.cc src/ov-complex.cc src/ov-complex.h src/ov-cx-mat.cc src/ov-cx-mat.h src/ov-cx-sparse.cc src/ov-cx-sparse.h src/ov-re-sparse.cc |
| diffstat | 9 files changed, 278 insertions(+), 2 deletions(-) [+] |
line wrap: on
line diff
--- a/src/ChangeLog Fri Mar 28 13:09:35 2008 -0400 +++ b/src/ChangeLog Fri Mar 28 18:26:29 2008 +0100 @@ -1,3 +1,27 @@ +2008-03-28 David Bateman <dbateman@free.fr> + + * ov-complex.cc (SCALAR_MAPPER, CD_SCALAR_MAPPER): New macro for + complex values with zero imaginary part. + (erf, erfc, gamma, lgamma): Use the new mappers to define these + mapper functions. + * ov-complex.h (erf, erfc, gamma, lgamma): Declare them. + * ov-cx-mat.cc (any_element_less_than, any_element_greater_than): + New static functions + (DARRAY_MAPPER, CD_ARRAY_MAPPER): New macro for complex values + with zero imaginary part. + (erf, erfc, gamma, lgamma): Use the new mappers to define these + mapper functions. + * ov-cx-mat.h (erf, erfc, gamma, lgamma): Declare them. + * ov-cx-sparse.cc (any_element_less_than, any_element_greater_than): + New static functions + (DSPARSE_MAPPER, CD_SPARSE_MAPPER): New macro for complex values + with zero imaginary part. + (erf, erfc, gamma, lgamma): Use the new mappers to define these + mapper functions. + * ov-cx-sparse.h (erf, erfc, gamma, lgamma): Declare them. + * ov-re-sparse.cc (CD_SPARSE_MAPPER): Use correct mapper functors. + * mapper.cc: Add tests for the above cases. + 2008-03-27 John W. Eaton <jwe@octave.org> * DLD-FUNCTIONS/max.cc: Rename from minmax.cc.
--- a/src/mappers.cc Fri Mar 28 13:09:35 2008 -0400 +++ b/src/mappers.cc Fri Mar 28 18:26:29 2008 +0100 @@ -323,6 +323,14 @@ return retval; } +/* + +%!test +%! a = -1i*sqrt(-1/(6.4187*6.4187)); +%! assert (erf(a), erf(real(a))); + +*/ + DEFUN (erfc, args, , "-*- texinfo -*-\n\ @deftypefn {Mapping Function} {} erfc (@var{z})\n\ @@ -347,6 +355,14 @@ return retval; } +/* + +%!test +%! a = -1i*sqrt(-1/(6.4187*6.4187)); +%! assert (erfc(a), erfc(real(a))); + +*/ + DEFUN (exp, args, , "-*- texinfo -*-\n\ @deftypefn {Mapping Function} {} exp (@var{x})\n\ @@ -467,6 +483,14 @@ return retval; } +/* + +%!test +%! a = -1i*sqrt(-1/(6.4187*6.4187)); +%! assert (gamma(a), gamma(real(a))); + +*/ + DEFUN (imag, args, , "-*- texinfo -*-\n\ @deftypefn {Mapping Function} {} imag (@var{z})\n\ @@ -753,6 +777,14 @@ return retval; } +/* + +%!test +%! a = -1i*sqrt(-1/(6.4187*6.4187)); +%! assert (lgamma(a), lgamma(real(a))); + +*/ + DEFUN (log, args, , "-*- texinfo -*-\n\ @deftypefn {Mapping Function} {} log (@var{x})\n\
--- a/src/ov-complex.cc Fri Mar 28 13:09:35 2008 -0400 +++ b/src/ov-complex.cc Fri Mar 28 18:26:29 2008 +0100 @@ -358,6 +358,42 @@ return octave_value (FCN (scalar)); \ } +#define SCALAR_MAPPER(MAP, FCN) \ + octave_value \ + octave_complex::MAP (void) const \ + { \ + if (scalar.imag () == 0) \ + return octave_value (FCN (scalar.real ())); \ + else \ + { \ + error ("%s: not defined for complex arguments", #MAP); \ + return octave_value (); \ + } \ + } + +#define CD_SCALAR_MAPPER(MAP, RFCN, CFCN, L1, L2) \ + octave_value \ + octave_complex::MAP (void) const \ + { \ + if (scalar.imag () == 0) \ + { \ + double re = scalar.real (); \ + return (re < L1 || re > L2 \ + ? octave_value (CFCN (scalar)) \ + : octave_value (RFCN (re))); \ + } \ + else \ + { \ + error ("%s: not defined for complex arguments", #MAP); \ + return octave_value (); \ + } \ + } + +SCALAR_MAPPER (erf, ::erf) +SCALAR_MAPPER (erfc, ::erfc) +SCALAR_MAPPER (gamma, xgamma) +CD_SCALAR_MAPPER (lgamma, xlgamma, xlgamma, 0.0, octave_Inf) + COMPLEX_MAPPER (abs, xabs) COMPLEX_MAPPER (acos, ::acos) COMPLEX_MAPPER (acosh, ::acosh)
--- a/src/ov-complex.h Fri Mar 28 13:09:35 2008 -0400 +++ b/src/ov-complex.h Fri Mar 28 18:26:29 2008 +0100 @@ -150,6 +150,10 @@ mxArray *as_mxArray (void) const; + octave_value erf (void) const; + octave_value erfc (void) const; + octave_value gamma (void) const; + octave_value lgamma (void) const; octave_value abs (void) const; octave_value acos (void) const; octave_value acosh (void) const;
--- a/src/ov-cx-mat.cc Fri Mar 28 13:09:35 2008 -0400 +++ b/src/ov-cx-mat.cc Fri Mar 28 18:26:29 2008 +0100 @@ -670,6 +670,40 @@ return x.real (); } +static bool +any_element_less_than (const NDArray& a, double val) +{ + octave_idx_type len = a.length (); + const double *m = a.fortran_vec (); + + for (octave_idx_type i = 0; i < len; i++) + { + OCTAVE_QUIT; + + if (m[i] < val) + return true; + } + + return false; +} + +static bool +any_element_greater_than (const NDArray& a, double val) +{ + octave_idx_type len = a.length (); + const double *m = a.fortran_vec (); + + for (octave_idx_type i = 0; i < len; i++) + { + OCTAVE_QUIT; + + if (m[i] > val) + return true; + } + + return false; +} + #define ARRAY_MAPPER(MAP, AMAP, FCN) \ octave_value \ octave_complex_matrix::MAP (void) const \ @@ -678,6 +712,56 @@ return matrix.map (cmap); \ } +#define DARRAY_MAPPER(MAP, AMAP, FCN) \ + octave_value \ + octave_complex_matrix::MAP (void) const \ + { \ + static ComplexNDArray::dmapper dmap = ximag; \ + NDArray m = matrix.map (dmap); \ + if (m.all_elements_are_zero ()) \ + { \ + dmap = xreal; \ + m = matrix.map (dmap); \ + static AMAP cmap = FCN; \ + return m.map (cmap); \ + } \ + else \ + { \ + error ("%s: not defined for complex arguments", #MAP); \ + return octave_value (); \ + } \ + } + +#define CD_ARRAY_MAPPER(MAP, RFCN, CFCN, L1, L2) \ + octave_value \ + octave_complex_matrix::MAP (void) const \ + { \ + static ComplexNDArray::dmapper idmap = ximag; \ + NDArray m = matrix.map (idmap); \ + if (m.all_elements_are_zero ()) \ + { \ + static ComplexNDArray::dmapper rdmap = xreal; \ + m = matrix.map (rdmap); \ + static NDArray::dmapper dmap = RFCN; \ + static NDArray::cmapper cmap = CFCN; \ + return (any_element_less_than (m, L1) \ + ? octave_value (m.map (cmap)) \ + : (any_element_greater_than (m, L2) \ + ? octave_value (m.map (cmap)) \ + : octave_value (m.map (dmap)))); \ + } \ + else \ + { \ + /*error ("%s: not defined for complex arguments", #MAP); */ \ + return octave_value (m); \ + } \ + } + +DARRAY_MAPPER (erf, NDArray::dmapper, ::erf) +DARRAY_MAPPER (erfc, NDArray::dmapper, ::erfc) +DARRAY_MAPPER (gamma, NDArray::dmapper, xgamma) +CD_ARRAY_MAPPER (lgamma, xlgamma, xlgamma, 0.0, octave_Inf) + ARRAY_MAPPER (abs, ComplexNDArray::dmapper, xabs) ARRAY_MAPPER (acos, ComplexNDArray::cmapper, ::acos) ARRAY_MAPPER (acosh, ComplexNDArray::cmapper, ::acosh)
--- a/src/ov-cx-mat.h Fri Mar 28 13:09:35 2008 -0400 +++ b/src/ov-cx-mat.h Fri Mar 28 18:26:29 2008 +0100 @@ -153,6 +153,10 @@ mxArray *as_mxArray (void) const; + octave_value erf (void) const; + octave_value erfc (void) const; + octave_value gamma (void) const; + octave_value lgamma (void) const; octave_value abs (void) const; octave_value acos (void) const; octave_value acosh (void) const;
--- a/src/ov-cx-sparse.cc Fri Mar 28 13:09:35 2008 -0400 +++ b/src/ov-cx-sparse.cc Fri Mar 28 18:26:29 2008 +0100 @@ -839,6 +839,44 @@ return x.real (); } +static bool +any_element_less_than (const SparseMatrix& a, double val) +{ + octave_idx_type len = a.nnz (); + + if (val > 0. && len != a.numel ()) + return true; + + for (octave_idx_type i = 0; i < len; i++) + { + OCTAVE_QUIT; + + if (a.data(i) < val) + return true; + } + + return false; +} + +static bool +any_element_greater_than (const SparseMatrix& a, double val) +{ + octave_idx_type len = a.nnz (); + + if (val < 0. && len != a.numel ()) + return true; + + for (octave_idx_type i = 0; i < len; i++) + { + OCTAVE_QUIT; + + if (a.data(i) > val) + return true; + } + + return false; +} + #define SPARSE_MAPPER(MAP, AMAP, FCN) \ octave_value \ octave_sparse_complex_matrix::MAP (void) const \ @@ -847,6 +885,56 @@ return matrix.map (cmap); \ } +#define DSPARSE_MAPPER(MAP, AMAP, FCN) \ + octave_value \ + octave_sparse_complex_matrix::MAP (void) const \ + { \ + static SparseComplexMatrix::dmapper dmap = ximag; \ + SparseMatrix m = matrix.map (dmap); \ + if (m.all_elements_are_zero ()) \ + { \ + dmap = xreal; \ + m = matrix.map (dmap); \ + static AMAP cmap = FCN; \ + return m.map (cmap); \ + } \ + else \ + { \ + error ("%s: not defined for complex arguments", #MAP); \ + return octave_value (); \ + } \ + } + +#define CD_SPARSE_MAPPER(MAP, RFCN, CFCN, L1, L2) \ + octave_value \ + octave_sparse_complex_matrix::MAP (void) const \ + { \ + static SparseComplexMatrix::dmapper idmap = ximag; \ + SparseMatrix m = matrix.map (idmap); \ + if (m.all_elements_are_zero ()) \ + { \ + static SparseComplexMatrix::dmapper rdmap = xreal; \ + m = matrix.map (rdmap); \ + static SparseMatrix::dmapper dmap = RFCN; \ + static SparseMatrix::cmapper cmap = CFCN; \ + return (any_element_less_than (m, L1) \ + ? octave_value (m.map (cmap)) \ + : (any_element_greater_than (m, L2) \ + ? octave_value (m.map (cmap)) \ + : octave_value (m.map (dmap)))); \ + } \ + else \ + { \ + error ("%s: not defined for complex arguments", #MAP); \ + return octave_value (); \ + } \ + } + +DSPARSE_MAPPER (erf, SparseMatrix::dmapper, ::erf) +DSPARSE_MAPPER (erfc, SparseMatrix::dmapper, ::erfc) +DSPARSE_MAPPER (gamma, SparseMatrix::dmapper, xgamma) +CD_SPARSE_MAPPER (lgamma, xlgamma, xlgamma, 0.0, octave_Inf) + SPARSE_MAPPER (abs, SparseComplexMatrix::dmapper, xabs) SPARSE_MAPPER (acos, SparseComplexMatrix::cmapper, ::acos) SPARSE_MAPPER (acosh, SparseComplexMatrix::cmapper, ::acosh)
--- a/src/ov-cx-sparse.h Fri Mar 28 13:09:35 2008 -0400 +++ b/src/ov-cx-sparse.h Fri Mar 28 18:26:29 2008 +0100 @@ -152,6 +152,10 @@ mxArray *as_mxArray (void) const; + octave_value erf (void) const; + octave_value erfc (void) const; + octave_value gamma (void) const; + octave_value lgamma (void) const; octave_value abs (void) const; octave_value acos (void) const; octave_value acosh (void) const;
--- a/src/ov-re-sparse.cc Fri Mar 28 13:09:35 2008 -0400 +++ b/src/ov-re-sparse.cc Fri Mar 28 18:26:29 2008 +0100 @@ -886,8 +886,8 @@ octave_value \ octave_sparse_matrix::MAP (void) const \ { \ - static NDArray::dmapper dmap = RFCN; \ - static NDArray::cmapper cmap = CFCN; \ + static SparseMatrix::dmapper dmap = RFCN; \ + static SparseMatrix::cmapper cmap = CFCN; \ \ return (any_element_less_than (matrix, L1) \ ? octave_value (matrix.map (cmap)) \