Mercurial > octave

changeset 30716:a7c4fb821d64

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
pow2: convert to C++ and make Matlab compatible (bug #61968). * libinterp/corefcn/pow2.cc: new C++ implementation. * libinterp/corefcn/module.mk: add new function to build system. * scripts/specfun/pow2.m: remove m-file implementation. * scripts/specfun/module.mk: remove "pow2.m" from build system.
author Kai T. Ohlhus <k.ohlhus@gmail.com>
date Sun, 06 Feb 2022 23:36:46 +0900
parents 00db74e19df4
children 555676697eef
files libinterp/corefcn/module.mk libinterp/corefcn/pow2.cc scripts/specfun/module.mk scripts/specfun/pow2.m
diffstat 4 files changed, 313 insertions(+), 75 deletions(-) [+]
line wrap: on
line diff
--- a/libinterp/corefcn/module.mk	Sun Feb 06 15:06:54 2022 +0100
+++ b/libinterp/corefcn/module.mk	Sun Feb 06 23:36:46 2022 +0900
@@ -229,6 +229,7 @@
   %reldir%/ordschur.cc \
   %reldir%/pager.cc \
   %reldir%/pinv.cc \
+  %reldir%/pow2.cc \
   %reldir%/pr-flt-fmt.cc \
   %reldir%/pr-output.cc \
   %reldir%/procstream.cc \
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libinterp/corefcn/pow2.cc	Sun Feb 06 23:36:46 2022 +0900
@@ -0,0 +1,312 @@
+////////////////////////////////////////////////////////////////////////
+//
+// Copyright (C) 2022 The Octave Project Developers
+//
+// See the file COPYRIGHT.md in the top-level directory of this
+// distribution or <https://octave.org/copyright/>.
+//
+// This file is part of Octave.
+//
+// Octave is free software: you can redistribute it and/or modify it
+// under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// Octave is distributed in the hope that it will be useful, but
+// WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with Octave; see the file COPYING.  If not, see
+// <https://www.gnu.org/licenses/>.
+//
+////////////////////////////////////////////////////////////////////////
+
+#if defined (HAVE_CONFIG_H)
+#  include "config.h"
+#endif
+
+#include <cmath>
+
+#include "defun.h"
+#include "error.h"
+#include "errwarn.h"
+
+template <typename T>
+void
+map_2_xldexp (Array<T>& y, const Array<T>& f, const Array<T>& e)
+{
+  if (f.numel () == e.numel () || e.numel () == 1)
+    y = Array<T>(f.dims ());
+  else if (f.numel () == 1)
+    y = Array<T>(e.dims ());
+  else
+    error ("pow2: Input dimensions must agree.");
+
+  octave_idx_type f_inc = (f.numel () == 1) ? 0 : 1;
+  octave_idx_type e_inc = (e.numel () == 1) ? 0 : 1;
+
+  for (octave_idx_type i = 0; i < y.numel (); i++)
+    y.xelem (i) = std::ldexp (f.xelem (i * f_inc),
+                              static_cast<int> (e.xelem (i * e_inc)));
+}
+
+void
+map_2_xldexp_sparse (SparseMatrix& y, const SparseMatrix& f,
+                                      const SparseMatrix& e)
+{
+  if (e.numel () == 1)
+    {
+      int ee = static_cast<int> (e.data (0));
+      for (octave_idx_type i = 0; i < y.nnz (); i++)
+        y.data (i) = std::ldexp (f.data (i), ee);
+    }
+  else if (f.numel () == e.numel ())
+    {
+      octave_idx_type col = 1;
+      for (octave_idx_type i = 0; i < y.nnz (); i++)
+        {
+          // Determine current column.
+          while (i >= f.cidx (col))
+            col++;
+          int ee = static_cast<int> (e.elem (f.ridx (i), col - 1));
+          y.data (i) = std::ldexp (f.data (i), ee);
+        }
+    }
+  else
+    error ("pow2: Input dimensions must agree.");
+}
+
+OCTAVE_NAMESPACE_BEGIN
+
+DEFUN (pow2, args, ,
+       doc: /* -*- texinfo -*-
+@deftypefn  {} {@var{y} =} pow2 (@var{x})
+@deftypefnx {} {@var{y} =} pow2 (@var{f}, @var{e})
+With one input argument, compute
+@tex
+$y = 2^x$
+@end tex
+@ifnottex
+y = 2 .^ x
+@end ifnottex
+for each element of @var{x}.
+
+With two input arguments, return
+@tex
+$y = f \cdot 2^e$,
+@end tex
+@ifnottex
+y = f .* (2 .^ e).
+@end ifnottex
+where for complex inputs only the real part of both inputs is regarded
+and from @var{e} only the real integer part.  This calling form corresponds
+to C/C++ standard function ldexp().
+@seealso{log2, nextpow2, power}
+@end deftypefn */)
+{
+  if (args.length () < 1 || args.length () > 2)
+    print_usage ();
+
+  if (! args(0).isfloat ())
+    err_wrong_type_arg ("pow2", args(0));
+
+  // Call exp2(f) where possible for numerical more accurate results.
+  if (args.length () == 1)
+    {
+      if (args(0).iscomplex ())
+        {
+          // The C++ standard does not define exp2 for complex arguments.
+          // Therefore call `2.^x`.
+          octave_value retval = octave::binary_op (octave_value::op_el_pow,
+                                                   2, args(0));
+
+          // Preserve sparse datatype, but even for sparse input fill-up
+          // is unavoidable `2^0 == 1` thus cast only.
+          if (args(0).issparse ())
+            retval = octave_value (retval.sparse_complex_matrix_value ());
+
+          return ovl (retval);
+        }
+      else if (args(0).is_single_type ())
+        {
+          FloatNDArray x = args(0).float_array_value ();
+          FloatNDArray y (x.dims ());
+          for (octave_idx_type i = 0; i < y.numel (); i++)
+            y.xelem (i) = exp2f (x.xelem (i));
+          return ovl (y);
+        }
+      else
+        {
+          NDArray x = args(0).array_value ();
+          NDArray y (x.dims ());
+          for (octave_idx_type i = 0; i < y.numel (); i++)
+            y.xelem (i) = exp2 (x.xelem (i));
+
+          // Preserve sparse datatype, but even for sparse input fill-up
+          // is unavoidable `2^0 == 1` thus cast only.
+          if (args(0).issparse ())
+            return ovl (SparseMatrix (y));
+          else
+            return ovl (y);
+        }
+    }
+
+  // For Matlab compatibility, the two argument call `y = pow2 (f, e)`
+  // corresponds to std::ldexp() (see bug #61968).  The result y is
+  // computed quickly by adding the integer part of e to the floating-point
+  // exponent of f.
+
+  if (! args(1).isfloat ())
+    err_wrong_type_arg ("pow2", args(1));
+
+  if (args(0).iscomplex () || args(1).iscomplex ())
+    warning_with_id ("Octave:pow2:only-real",
+                     "pow2: Imaginary part is ignored.");
+
+  // Note: Matlab R2021a errors on `pow2 (sparse (f), single (e))`,
+  //       but sparsity in f determines output and can significantly
+  //       reduce computation, e.g. `N=1e5; pow2(speye(N),sparse(N,N))`.
+  if (args(0).issparse ())
+    {
+      SparseMatrix f = args(0).sparse_matrix_value ();
+
+      // Special case: return a sparse zero matrix in size of e.
+      if ((f.numel () == 1) && (f.nnz () == 0))
+        return ovl (SparseMatrix (args(1).rows (), args(1).columns ()));
+
+      // Only do sparse computation, if it pays off.  For scalar f fill-up
+      // is unavoidable even for sparse e because `f * 2^0 == f`.  Use dense
+      // code below in this case.
+      if (f.numel () > 1)
+        {
+          SparseMatrix e = args(1).sparse_matrix_value ();
+          SparseMatrix y = SparseMatrix (f);
+          map_2_xldexp_sparse (y, f, e);
+          return ovl (y);
+        }
+    }
+
+  if (args(0).is_single_type () || args(1).is_single_type ())
+    {
+      FloatNDArray f = args(0).float_array_value ();
+      FloatNDArray e = args(1).float_array_value ();
+      FloatNDArray y;
+      map_2_xldexp (y, f, e);
+      return ovl (y);
+    }
+  else
+    {
+      NDArray f = args(0).array_value ();
+      NDArray e = args(1).array_value ();
+      NDArray y;
+      map_2_xldexp (y, f, e);
+
+      // Preserve sparse datatype.
+      // Cases for efficient use of sparsity were treated above already.
+      if (args(0).issparse ())
+        return ovl (SparseMatrix (y));
+      else
+        return ovl (y);
+    }
+}
+
+/*
+## Call `y = pow2 (x)`
+
+%!test
+%! fcns = {@double, @single, @complex};
+%! x = [3, 0, -3];
+%! v = [8, 1, .125];
+%! for i = 1:numel (fcns)
+%!   fcn = fcns{i};
+%!   assert (pow2 (fcn (x)), fcn (v), sqrt (eps));
+%! endfor
+
+%!test
+%! fcns = {@double, @single, @complex, @sparse};
+%! x = [3, 1, -3];
+%! v = [8, 2, .125];
+%! for i = 1:numel (fcns)
+%!   fcn = fcns{i};
+%!   assert (pow2 (fcn (x)), fcn (v), sqrt (eps));
+%! endfor
+
+%!test
+%! fcns = {@double, @single, @complex, @sparse};
+%! x = [1, 1+1i, 1i];
+%! for i = 1:numel (fcns)
+%!   fcn = fcns{i};
+%!   assert (pow2 (fcn (x)), fcn (2) .^ fcn (x), sqrt (eps));
+%! endfor
+
+## Call `y = pow2 (f, e)`
+
+%!test
+%! fcns = {@double, @single, @complex, @sparse};
+%! f = [2 2];
+%! e = [2 2];
+%! z = [8 8];
+%! warning ("off", "Octave:pow2:only-real", "local");
+%! for i = 1:numel (fcns)
+%!   fcn = fcns{i};
+%!   assert (pow2 (fcn (f), fcn (e)), real (fcn (z)));
+%! endfor
+
+## Only integer part is taken into account.
+%!test
+%! f = 2;
+%! e = [2, 2.1, 2.2, 2.4, 2.5, 2.8];
+%! z = 8 .* ones (1, length (e));
+%! assert (pow2 (f, e), z);
+
+## Only real part is taken into account.
+%!test
+%! f = [1+1i, 1];
+%! e = 2;
+%! z = [4, 4];
+%! warning ("off", "Octave:pow2:only-real", "local");
+%! assert (pow2 (f, e), z);
+
+%!test
+%! f = 1;
+%! e = [1+1i, 1];
+%! z = [2, 2];
+%! warning ("off", "Octave:pow2:only-real", "local");
+%! assert (pow2 (f, e), z);
+
+%!test
+%! f = f = [1/2, pi/4, -3/4, 1/2, 1-eps()/2, 1/2];
+%! e = [1, 2, 2, -51, 1024, -1021];
+%! z = [1, pi, -3, eps(), realmax(), realmin()];
+%! assert (pow2 (f, e), z);
+
+## Tests for sparsity.
+%!assert (pow2 (sparse (0), ones  (3)), sparse (3, 3));
+%!assert (pow2 (sparse (1), ones  (3)), 2 .* sparse (ones (3)));
+%!assert (pow2 (sparse (1), speye (3)), sparse (ones (3) + eye (3)));
+%!assert (pow2 (sparse (3, 3), ones (3)), sparse (3, 3));
+%!assert (pow2 (speye (3), ones (3)), 2 .* speye (3));
+%!assert (pow2 (speye (3), 1),        2 .* speye (3));
+
+%!test
+%! f = speye (3);
+%! e = sparse (3, 3);
+%! e(1,1) = 1;
+%! e(1,3) = 1;
+%! z = f;
+%! z(1,1) = 2;
+%! assert (pow2 (f, e), z);
+
+## Large sparse matrix (only few real elements).
+%!test
+%! ## FIXME: `N = 1e5` would be a better test, but `assert` fills-up somehow.
+%! N = 1e3;
+%! assert (pow2 (speye  (N), sparse (N,N)), speye (N));
+%! assert (pow2 (sparse (0), speye  (N)),   sparse(N,N));
+
+%!error <Invalid call> pow2 ()
+*/
+
+OCTAVE_NAMESPACE_END
--- a/scripts/specfun/module.mk	Sun Feb 06 15:06:54 2022 +0100
+++ b/scripts/specfun/module.mk	Sun Feb 06 23:36:46 2022 +0900
@@ -19,7 +19,6 @@
   %reldir%/nchoosek.m \
   %reldir%/nthroot.m \
   %reldir%/perms.m \
-  %reldir%/pow2.m \
   %reldir%/primes.m \
   %reldir%/reallog.m \
   %reldir%/realpow.m \
--- a/scripts/specfun/pow2.m	Sun Feb 06 15:06:54 2022 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,74 +0,0 @@
-########################################################################
-##
-## Copyright (C) 1995-2022 The Octave Project Developers
-##
-## See the file COPYRIGHT.md in the top-level directory of this
-## distribution or <https://octave.org/copyright/>.
-##
-## This file is part of Octave.
-##
-## Octave is free software: you can redistribute it and/or modify it
-## under the terms of the GNU General Public License as published by
-## the Free Software Foundation, either version 3 of the License, or
-## (at your option) any later version.
-##
-## Octave is distributed in the hope that it will be useful, but
-## WITHOUT ANY WARRANTY; without even the implied warranty of
-## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-## GNU General Public License for more details.
-##
-## You should have received a copy of the GNU General Public License
-## along with Octave; see the file COPYING.  If not, see
-## <https://www.gnu.org/licenses/>.
-##
-########################################################################
-
-## -*- texinfo -*-
-## @deftypefn  {} {} pow2 (@var{x})
-## @deftypefnx {} {} pow2 (@var{f}, @var{e})
-## With one input argument, compute
-## @tex
-## $2^x$
-## @end tex
-## @ifnottex
-## 2 .^ x
-## @end ifnottex
-## for each element of @var{x}.
-##
-## With two input arguments, return
-## @tex
-## $f \cdot 2^e$.
-## @end tex
-## @ifnottex
-## f .* (2 .^ e).
-## @end ifnottex
-## @seealso{log2, nextpow2, power}
-## @end deftypefn
-
-function y = pow2 (f, e)
-
-  if (nargin < 1)
-    print_usage ();
-  endif
-
-  if (nargin == 1)
-    y = 2 .^ f;
-  else
-    y = f .* (2 .^ e);
-  endif
-
-endfunction
-
-
-%!test
-%! x = [3, 0, -3];
-%! v = [8, 1, .125];
-%! assert (pow2 (x), v, sqrt (eps));
-
-%!test
-%! x = [3, 0, -3, 4, 0, -4, 5, 0, -5];
-%! y = [-2, -2, -2, 1, 1, 1, 3, 3, 3];
-%! z = x .* (2 .^ y);
-%! assert (pow2 (x,y), z, sqrt (eps));
-
-%!error <Invalid call> pow2 ()