[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Changeset] Add the hypot function
From: |
David Bateman |
Subject: |
[Changeset] Add the hypot function |
Date: |
Sun, 23 Mar 2008 22:31:11 +0100 |
User-agent: |
Thunderbird 2.0.0.12 (X11/20080306) |
Matlab maps the libm (or in their case the fdlib implementation of libm)
hypot function into matlab, that calculates sqrt(x.^2+y.^2) but in a way
that avoids overflows for large numbers. This patch adds the hypot
function to matlab and fixes a small issue with the atan2 function when
it is called with sparse matrices.
D.
# HG changeset patch
# User David Bateman <address@hidden>
# Date 1206307680 -3600
# Node ID a469b39eb915ffe9bcb1c6ed135c1a47963c59ff
# Parent a14add2f1a1c3000c4cfb83b1fb1bd98ed629654
Add the hypot function
diff --git a/src/ChangeLog b/src/ChangeLog
--- a/src/ChangeLog
+++ b/src/ChangeLog
@@ -1,3 +1,9 @@ 2008-03-21 John W. Eaton <address@hidden
+2008-03-23 David Bateman <address@hidden>
+
+ * data.cc (map_s_s): Fix for sparse/sparse mappers that resulted
+ in an empty sparse matrix being returned.
+ (Fhypot): New function based on the libm hypot function.
+
2008-03-21 John W. Eaton <address@hidden>
* ov-cell.h (octave_cell::is_constant): Return true.
diff --git a/src/data.cc b/src/data.cc
--- a/src/data.cc
+++ b/src/data.cc
@@ -353,7 +353,7 @@ map_s_s (d_dd_fcn f, const SparseMatrix&
}
else
{
- octave_idx_type nz = y.nnz ();
+ octave_idx_type nz = x.nnz () + y.nnz ();
retval = SparseMatrix (nr, nc, nz);
octave_idx_type ii = 0;
retval.cidx (ii) = 0;
@@ -403,6 +403,7 @@ map_s_s (d_dd_fcn f, const SparseMatrix&
retval.ridx (ii++) = r;
}
}
+ retval.cidx (j + 1) = ii;
}
retval.maybe_compress (false);
@@ -533,6 +534,156 @@ and @var{x}. The result is in range -pi
%! assert (size (atan2 (rand (2, 3, 4), 1), [2, 3, 4])
%! assert (size (atan2 (1, rand (2, 3, 4)), [2, 3, 4])
%! assert (size (atan2 (1, 2), [1, 1])
+*/
+
+DEFUN (hypot, args, ,
+ "-*- texinfo -*-\n\
address@hidden {Mapping Function} {} hypot (@var{x}, @var{y})\n\
+Compute square-root of the squares of @var{x} and @var{y}\n\
+element-by-element. This equivalent to @code{sqrt (@var{x}.^ 2 + @var{y}\n\
+.^ 2)}, but calculated in a manner that avoids overflows for large\n\
+values of @var{x} or @var{y}.\n\
address@hidden deftypefn")
+{
+ octave_value retval;
+
+ int nargin = args.length ();
+
+ if (nargin == 2 && args(0).is_defined () && args(1).is_defined ())
+ {
+ if (args(0).is_integer_type () || args(1).is_integer_type ())
+ error ("hypot: not defined for integer types");
+ else
+ {
+ octave_value arg_x = args(0);
+ octave_value arg_y = args(1);
+
+ dim_vector x_dims = arg_x.dims ();
+ dim_vector y_dims = arg_y.dims ();
+
+ bool x_is_scalar = x_dims.all_ones ();
+ bool y_is_scalar = y_dims.all_ones ();
+
+ if (y_is_scalar && x_is_scalar)
+ {
+ double x;
+ if (arg_x.is_complex_type ())
+ x = abs (arg_x.complex_value ());
+ else
+ x = arg_x.double_value ();
+
+ if (! error_state)
+ {
+ double y;
+ if (arg_y.is_complex_type ())
+ y = abs (arg_y.complex_value ());
+ else
+ y = arg_y.double_value ();
+
+ if (! error_state)
+ retval = hypot (x, y);
+ }
+ }
+ else if (y_is_scalar)
+ {
+ NDArray x;
+ if (arg_x.is_complex_type ())
+ x = arg_x.complex_array_value ().abs ();
+ else
+ x = arg_x.array_value ();
+
+ if (! error_state)
+ {
+ double y;
+ if (arg_y.is_complex_type ())
+ y = abs (arg_y.complex_value ());
+ else
+ y = arg_y.double_value ();
+
+ if (! error_state)
+ retval = map_m_d (hypot, x, y);
+ }
+ }
+ else if (x_is_scalar)
+ {
+ double x;
+ if (arg_x.is_complex_type ())
+ x = abs (arg_x.complex_value ());
+ else
+ x = arg_x.double_value ();
+
+ if (! error_state)
+ {
+ NDArray y;
+ if (arg_y.is_complex_type ())
+ y = arg_y.complex_array_value ().abs ();
+ else
+ y = arg_y.array_value ();
+
+ if (! error_state)
+ retval = map_d_m (hypot, x, y);
+ }
+ }
+ else if (y_dims == x_dims)
+ {
+ if (arg_x.is_sparse_type () && arg_y.is_sparse_type ())
+ {
+ SparseMatrix x;
+ if (arg_x.is_complex_type ())
+ x = arg_x.sparse_complex_matrix_value ().abs ();
+ else
+ x = arg_x.sparse_matrix_value ();
+
+ if (! error_state)
+ {
+ SparseMatrix y;
+ if (arg_y.is_complex_type ())
+ y = arg_y.sparse_complex_matrix_value ().abs ();
+ else
+ y = arg_y.sparse_matrix_value ();
+
+ if (! error_state)
+ retval = map_s_s (hypot, x, y);
+ }
+ }
+ else
+ {
+ NDArray x;
+ if (arg_x.is_complex_type ())
+ x = arg_x.complex_array_value ().abs ();
+ else
+ x = arg_x.array_value ();
+
+ if (! error_state)
+ {
+ NDArray y;
+ if (arg_y.is_complex_type ())
+ y = arg_y.complex_array_value ().abs ();
+ else
+ y = arg_y.array_value ();
+
+ if (! error_state)
+ retval = map_m_m (hypot, x, y);
+ }
+ }
+ }
+ else
+ error ("hypot: nonconformant matrices");
+ }
+ }
+ else
+ print_usage ();
+
+ return retval;
+}
+
+/*
+%! assert (size (hypot (zeros (0, 2), zeros (0, 2))), [0, 2])
+%! assert (size (hypot (rand (2, 3, 4), zeros (2, 3, 4))), [2, 3, 4])
+%! assert (size (hypot (rand (2, 3, 4), 1), [2, 3, 4])
+%! assert (size (hypot (1, rand (2, 3, 4)), [2, 3, 4])
+%! assert (size (hypot (1, 2), [1, 1])
+%! assert (hypot (1:10,1:10), sqrt(2) * [1:10]);
*/
DEFUN (fmod, args, ,
- [Changeset] Add the hypot function,
David Bateman <=