# HG changeset patch
# User Kai T. Ohlhus <k.ohlhus@gmail.com>
# Date 1643708945 -32400
#      Tue Feb 01 18:49:05 2022 +0900
# Node ID 6f02a336900969ae926be3226eaa2206702bfb29
# Parent  a711560715a81caead4f303311bf1bae0ff4d016
pow2: make C++ function

diff -r a711560715a8 -r 6f02a3369009 libinterp/corefcn/data.cc
--- a/libinterp/corefcn/data.cc	Mon Jan 31 17:28:28 2022 +0100
+++ b/libinterp/corefcn/data.cc	Tue Feb 01 18:49:05 2022 +0900
@@ -569,6 +569,139 @@ and @var{e} is an integer.  If
 %!assert <*42583> (all (log2 (pow2 (-1074:1023)) == -1074:1023))
 */
 
+template <typename T>
+void
+map_2_xpow2 (Array<T>& y, const Array<T>& f, const Array<T>& e)
+{
+  if (f.numel () == e.numel ())
+    {
+      y = Array<T>(f.dims ());
+      for (octave_idx_type i = 0; i < y.numel (); i++)
+        y.xelem (i) = ldexp (f.xelem (i), static_cast<int> (e.xelem (i)));
+    }
+  else if (f.numel () == 1)
+    {
+      y = Array<T>(e.dims ());
+      for (octave_idx_type i = 0; i < y.numel (); i++)
+        y.xelem (i) = ldexp (f.xelem (0), static_cast<int> (e.xelem (i)));
+    }
+  else if (e.numel () == 1)
+    {
+      y = Array<T>(f.dims ());
+      for (octave_idx_type i = 0; i < y.numel (); i++)
+        y.xelem (i) = ldexp (f.xelem (i), static_cast<int> (e.xelem (0)));
+    }
+  else
+    error ("pow2: Input dimensions must agree.");
+}
+
+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.length () == 1)
+    return ovl (args(0).exp2 ());
+
+  // 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.
+
+  // Check input.
+  if (! args(0).is_single_type () && ! args(0).is_double_type ())
+    err_wrong_type_arg ("pow2", args(0));
+
+  if (! args(1).is_single_type () && ! args(1).is_double_type ())
+    err_wrong_type_arg ("pow2", args(1));
+
+  if (args(0).iscomplex () || args(1).iscomplex ())
+    warning ("pow2: Imaginary part is ignored.");
+
+  // For Matlab compatibility: single type dominates.
+  if (args(0).is_single_type () || args(1).is_single_type ())
+    {
+      FloatNDArray y;
+      FloatNDArray f = args(0).float_array_value ();
+      FloatNDArray e = args(0).float_array_value ();
+      map_2_xpow2 (y, f, e);
+      return ovl (y);
+    }
+  else
+    {
+      NDArray y;
+      NDArray f = args(0).array_value ();
+      NDArray e = args(1).array_value ();
+      map_2_xpow2 (y, f, e);
+      return ovl (y);
+    }
+}
+
+/*
+%!test
+%! fcns = {@double, @single, @complex, @sparse};
+%! 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
+
+% 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];
+%! assert (pow2 (f, e), z);
+
+%!test
+%! f = 1;
+%! e = [1 + 1i, 1];
+%! z = [2, 2];
+%! assert (pow2 (f, e), z);
+
+%!test
+%! f = str2num (char ({'1/2','pi/4','-3/4','1/2','1-eps/2','1/2'}'));
+%! e = [1 2 2 -51 1024 -1021]';
+%! z = str2num (char ({'1','pi','-3','eps','realmax','realmin'}'));
+%! assert (pow2 (f, e), z);
+
+%!error <Invalid call> pow2 ()
+*/
+
 DEFUN (rem, args, ,
        doc: /* -*- texinfo -*-
 @deftypefn {} {} rem (@var{x}, @var{y})
diff -r a711560715a8 -r 6f02a3369009 libinterp/octave-value/ov-base.cc
--- a/libinterp/octave-value/ov-base.cc	Mon Jan 31 17:28:28 2022 +0100
+++ b/libinterp/octave-value/ov-base.cc	Tue Feb 01 18:49:05 2022 +0900
@@ -1094,6 +1094,7 @@ octave_base_value::get_umap_name (unary_
     "erfi",
     "dawson",
     "exp",
+    "exp2",
     "expm1",
     "isfinite",
     "fix",
diff -r a711560715a8 -r 6f02a3369009 libinterp/octave-value/ov-base.h
--- a/libinterp/octave-value/ov-base.h	Mon Jan 31 17:28:28 2022 +0100
+++ b/libinterp/octave-value/ov-base.h	Tue Feb 01 18:49:05 2022 +0900
@@ -800,6 +800,7 @@ public:
     umap_erfi,
     umap_dawson,
     umap_exp,
+    umap_exp2,
     umap_expm1,
     umap_isfinite,
     umap_fix,
diff -r a711560715a8 -r 6f02a3369009 libinterp/octave-value/ov-complex.cc
--- a/libinterp/octave-value/ov-complex.cc	Mon Jan 31 17:28:28 2022 +0100
+++ b/libinterp/octave-value/ov-complex.cc	Tue Feb 01 18:49:05 2022 +0900
@@ -525,6 +525,7 @@ octave_complex::map (unary_mapper_t umap
     SCALAR_MAPPER (cos, std::cos);
     SCALAR_MAPPER (cosh, std::cosh);
     SCALAR_MAPPER (exp, std::exp);
+    SCALAR_MAPPER (exp2, octave::math::exp2);
     SCALAR_MAPPER (expm1, octave::math::expm1);
     SCALAR_MAPPER (fix, octave::math::fix);
     SCALAR_MAPPER (floor, octave::math::floor);
diff -r a711560715a8 -r 6f02a3369009 libinterp/octave-value/ov-cx-mat.cc
--- a/libinterp/octave-value/ov-cx-mat.cc	Mon Jan 31 17:28:28 2022 +0100
+++ b/libinterp/octave-value/ov-cx-mat.cc	Tue Feb 01 18:49:05 2022 +0900
@@ -818,6 +818,7 @@ octave_complex_matrix::map (unary_mapper
     ARRAY_MAPPER (cos, Complex, std::cos);
     ARRAY_MAPPER (cosh, Complex, std::cosh);
     ARRAY_MAPPER (exp, Complex, std::exp);
+    ARRAY_MAPPER (exp2, Complex, octave::math::exp2);
     ARRAY_MAPPER (expm1, Complex, octave::math::expm1);
     ARRAY_MAPPER (fix, Complex, octave::math::fix);
     ARRAY_MAPPER (floor, Complex, octave::math::floor);
diff -r a711560715a8 -r 6f02a3369009 libinterp/octave-value/ov-cx-sparse.cc
--- a/libinterp/octave-value/ov-cx-sparse.cc	Mon Jan 31 17:28:28 2022 +0100
+++ b/libinterp/octave-value/ov-cx-sparse.cc	Tue Feb 01 18:49:05 2022 +0900
@@ -970,6 +970,7 @@ octave_sparse_complex_matrix::map (unary
     ARRAY_MAPPER (cos, Complex, std::cos);
     ARRAY_MAPPER (cosh, Complex, std::cosh);
     ARRAY_MAPPER (exp, Complex, std::exp);
+    ARRAY_MAPPER (exp2, Complex, octave::math::exp2);
     ARRAY_MAPPER (expm1, Complex, octave::math::expm1);
     ARRAY_MAPPER (fix, Complex, octave::math::fix);
     ARRAY_MAPPER (floor, Complex, octave::math::floor);
diff -r a711560715a8 -r 6f02a3369009 libinterp/octave-value/ov-float.cc
--- a/libinterp/octave-value/ov-float.cc	Mon Jan 31 17:28:28 2022 +0100
+++ b/libinterp/octave-value/ov-float.cc	Tue Feb 01 18:49:05 2022 +0900
@@ -386,6 +386,7 @@ octave_float_scalar::map (unary_mapper_t
     SCALAR_MAPPER (cos, ::cosf);
     SCALAR_MAPPER (cosh, ::coshf);
     SCALAR_MAPPER (exp, ::expf);
+    SCALAR_MAPPER (exp2, ::exp2f);
     SCALAR_MAPPER (expm1, octave::math::expm1);
     SCALAR_MAPPER (fix, octave::math::fix);
     SCALAR_MAPPER (floor, std::floor);
diff -r a711560715a8 -r 6f02a3369009 libinterp/octave-value/ov-flt-complex.cc
--- a/libinterp/octave-value/ov-flt-complex.cc	Mon Jan 31 17:28:28 2022 +0100
+++ b/libinterp/octave-value/ov-flt-complex.cc	Tue Feb 01 18:49:05 2022 +0900
@@ -467,6 +467,7 @@ octave_float_complex::map (unary_mapper_
     SCALAR_MAPPER (cos, std::cos);
     SCALAR_MAPPER (cosh, std::cosh);
     SCALAR_MAPPER (exp, std::exp);
+    SCALAR_MAPPER (exp2, octave::math::exp2);
     SCALAR_MAPPER (expm1, octave::math::expm1);
     SCALAR_MAPPER (fix, octave::math::fix);
     SCALAR_MAPPER (floor, octave::math::floor);
diff -r a711560715a8 -r 6f02a3369009 libinterp/octave-value/ov-flt-cx-mat.cc
--- a/libinterp/octave-value/ov-flt-cx-mat.cc	Mon Jan 31 17:28:28 2022 +0100
+++ b/libinterp/octave-value/ov-flt-cx-mat.cc	Tue Feb 01 18:49:05 2022 +0900
@@ -771,6 +771,7 @@ octave_float_complex_matrix::map (unary_
     ARRAY_MAPPER (cos, FloatComplex, std::cos);
     ARRAY_MAPPER (cosh, FloatComplex, std::cosh);
     ARRAY_MAPPER (exp, FloatComplex, std::exp);
+    ARRAY_MAPPER (exp2, FloatComplex, octave::math::exp2);
     ARRAY_MAPPER (expm1, FloatComplex, octave::math::expm1);
     ARRAY_MAPPER (fix, FloatComplex, octave::math::fix);
     ARRAY_MAPPER (floor, FloatComplex, octave::math::floor);
diff -r a711560715a8 -r 6f02a3369009 libinterp/octave-value/ov-flt-re-mat.cc
--- a/libinterp/octave-value/ov-flt-re-mat.cc	Mon Jan 31 17:28:28 2022 +0100
+++ b/libinterp/octave-value/ov-flt-re-mat.cc	Tue Feb 01 18:49:05 2022 +0900
@@ -835,6 +835,7 @@ octave_float_matrix::map (unary_mapper_t
     ARRAY_MAPPER (cos, float, ::cosf);
     ARRAY_MAPPER (cosh, float, ::coshf);
     ARRAY_MAPPER (exp, float, ::expf);
+    ARRAY_MAPPER (exp2, float, ::exp2f);
     ARRAY_MAPPER (expm1, float, octave::math::expm1);
     ARRAY_MAPPER (fix, float, octave::math::fix);
     ARRAY_MAPPER (floor, float, ::floorf);
diff -r a711560715a8 -r 6f02a3369009 libinterp/octave-value/ov-re-mat.cc
--- a/libinterp/octave-value/ov-re-mat.cc	Mon Jan 31 17:28:28 2022 +0100
+++ b/libinterp/octave-value/ov-re-mat.cc	Tue Feb 01 18:49:05 2022 +0900
@@ -962,6 +962,7 @@ octave_matrix::map (unary_mapper_t umap)
     ARRAY_MAPPER (cos, double, ::cos);
     ARRAY_MAPPER (cosh, double, ::cosh);
     ARRAY_MAPPER (exp, double, ::exp);
+    ARRAY_MAPPER (exp2, double, ::exp2);
     ARRAY_MAPPER (expm1, double, octave::math::expm1);
     ARRAY_MAPPER (fix, double, octave::math::fix);
     ARRAY_MAPPER (floor, double, ::floor);
diff -r a711560715a8 -r 6f02a3369009 libinterp/octave-value/ov-re-sparse.cc
--- a/libinterp/octave-value/ov-re-sparse.cc	Mon Jan 31 17:28:28 2022 +0100
+++ b/libinterp/octave-value/ov-re-sparse.cc	Tue Feb 01 18:49:05 2022 +0900
@@ -944,6 +944,7 @@ octave_sparse_matrix::map (unary_mapper_
     ARRAY_MAPPER (cos, double, ::cos);
     ARRAY_MAPPER (cosh, double, ::cosh);
     ARRAY_MAPPER (exp, double, ::exp);
+    ARRAY_MAPPER (exp2, double, ::exp2);
     ARRAY_MAPPER (expm1, double, octave::math::expm1);
     ARRAY_MAPPER (fix, double, octave::math::fix);
     ARRAY_MAPPER (floor, double, ::floor);
diff -r a711560715a8 -r 6f02a3369009 libinterp/octave-value/ov-scalar.cc
--- a/libinterp/octave-value/ov-scalar.cc	Mon Jan 31 17:28:28 2022 +0100
+++ b/libinterp/octave-value/ov-scalar.cc	Tue Feb 01 18:49:05 2022 +0900
@@ -407,6 +407,7 @@ octave_scalar::map (unary_mapper_t umap)
     SCALAR_MAPPER (cos, ::cos);
     SCALAR_MAPPER (cosh, ::cosh);
     SCALAR_MAPPER (exp, ::exp);
+    SCALAR_MAPPER (exp2, ::exp2);
     SCALAR_MAPPER (expm1, octave::math::expm1);
     SCALAR_MAPPER (fix, octave::math::fix);
     SCALAR_MAPPER (floor, std::floor);
diff -r a711560715a8 -r 6f02a3369009 libinterp/octave-value/ov.h
--- a/libinterp/octave-value/ov.h	Mon Jan 31 17:28:28 2022 +0100
+++ b/libinterp/octave-value/ov.h	Tue Feb 01 18:49:05 2022 +0900
@@ -1559,6 +1559,7 @@ public:
   MAPPER_FORWARD (erfi)
   MAPPER_FORWARD (dawson)
   MAPPER_FORWARD (exp)
+  MAPPER_FORWARD (exp2)
   MAPPER_FORWARD (expm1)
   MAPPER_FORWARD (isfinite)
   MAPPER_FORWARD (fix)
diff -r a711560715a8 -r 6f02a3369009 liboctave/numeric/lo-mappers.h
--- a/liboctave/numeric/lo-mappers.h	Mon Jan 31 17:28:28 2022 +0100
+++ b/liboctave/numeric/lo-mappers.h	Tue Feb 01 18:49:05 2022 +0900
@@ -98,6 +98,14 @@ namespace octave
     inline double exp2 (double x) { return std::exp2 (x); }
     inline float exp2 (float x) { return std::exp2f (x); }
 
+    // The C++ standard does not define std::exp2 for complex values.
+    template <typename T>
+    std::complex<T>
+    exp2 (const std::complex<T>& x)
+    {
+      return std::pow (std::complex<T> (2, 0), x);
+    }
+
     template <typename T>
     std::complex<T>
     ceil (const std::complex<T>& x)
diff -r a711560715a8 -r 6f02a3369009 scripts/specfun/module.mk
--- a/scripts/specfun/module.mk	Mon Jan 31 17:28:28 2022 +0100
+++ b/scripts/specfun/module.mk	Tue Feb 01 18:49:05 2022 +0900
@@ -19,7 +19,6 @@ FCN_FILE_DIRS += %reldir%
   %reldir%/nchoosek.m \
   %reldir%/nthroot.m \
   %reldir%/perms.m \
-  %reldir%/pow2.m \
   %reldir%/primes.m \
   %reldir%/reallog.m \
   %reldir%/realpow.m \
diff -r a711560715a8 -r 6f02a3369009 scripts/specfun/pow2.m
--- a/scripts/specfun/pow2.m	Mon Jan 31 17:28:28 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 ()
