# HG changeset patch
# User Rik <rik@octave.org>
# Date 1399244926 25200
# Node ID 8cc66f09158440c60e6d4758866338a3cef0a58d
# Parent  87fafe6d05d3c14116757a808920477f7bbb3045
Add "native" option to prod() (bug #40349).

* data.cc (Fprod): Rewrite docstring to mention new "native", "double" options.
Borrow code from Fsum replacing sum with prod.  Add %!tests for new behavior.

* fCNDArray.h (dprod): New function prototype.

* fCNDArray.cc (FloatComplexNDArray::dprod): New function.

* fNDArray.h (dprod): New function prototype.

* fNDArray.cc (FloatNDArray::dprod): New function.

* intNDArray.h (prod): New function prototype.

* intNDArray.cc (intNDArray<T>::prod): New function.

* mx-inlines.cc (op_dble_prod): New inline functions for implementing double
product.

diff -r 87fafe6d05d3 -r 8cc66f091584 libinterp/corefcn/data.cc
--- a/libinterp/corefcn/data.cc	Sat May 03 17:28:49 2014 -0700
+++ b/libinterp/corefcn/data.cc	Sun May 04 16:08:46 2014 -0700
@@ -1379,12 +1379,144 @@
        "-*- texinfo -*-\n\
 @deftypefn  {Built-in Function} {} prod (@var{x})\n\
 @deftypefnx {Built-in Function} {} prod (@var{x}, @var{dim})\n\
-Product of elements along dimension @var{dim}.  If @var{dim} is\n\
-omitted, it defaults to the first non-singleton dimension.\n\
+@deftypefnx {Built-in Function} {} prod (@dots{}, \"native\")\n\
+@deftypefnx {Built-in Function} {} prod (@dots{}, \"double\")\n\
+Product of elements along dimension @var{dim}.\n\
+\n\
+If @var{dim} is omitted, it defaults to the first non-singleton dimension.\n\
+\n\
+The optional @qcode{\"type\"} input determines the class of the variable used for\n\
+calculations.  If the argument @qcode{\"native\"} is given, then the operation is\n\
+performed in the same type as the original argument, rather than the default double\n\
+type.\n\
+\n\
+For example:\n\
+\n\
+@example\n\
+@group\n\
+prod ([true, true])\n\
+   @result{} 1\n\
+prod ([true, true], \"native\")\n\
+   @result{} true\n\
+@end group\n\
+@end example\n\
+\n\
+On the contrary, if @qcode{\"double\"} is given, the operation is performed in\n\
+double precision even for single precision inputs.\n\
 @seealso{cumprod, sum}\n\
 @end deftypefn")
 {
-  DATA_REDUCTION (prod);
+  octave_value retval;
+
+  int nargin = args.length ();
+
+  bool isnative = false;
+  bool isdouble = false;
+
+  if (nargin > 1 && args(nargin - 1).is_string ())
+    {
+      std::string str = args(nargin - 1).string_value ();
+
+      if (! error_state)
+        {
+          if (str == "native")
+            isnative = true;
+          else if (str == "double")
+            isdouble = true;
+          else
+            error ("prod: unrecognized type argument '%s'", str.c_str ());
+          nargin --;
+        }
+    }
+
+  if (error_state)
+    return retval;
+
+  if (nargin == 1 || nargin == 2)
+    {
+      octave_value arg = args(0);
+
+      int dim = -1;
+      if (nargin == 2)
+        {
+          dim = args(1).int_value () - 1;
+          if (dim < 0)
+            error ("prod: invalid dimension DIM = %d", dim + 1);
+        }
+
+      if (! error_state)
+        {
+          switch (arg.builtin_type ())
+            {
+            case btyp_double:
+              if (arg.is_sparse_type ())
+                retval = arg.sparse_matrix_value ().prod (dim);
+              else
+                retval = arg.array_value ().prod (dim);
+              break;
+            case btyp_complex:
+              if (arg.is_sparse_type ())
+                retval = arg.sparse_complex_matrix_value ().prod (dim);
+              else
+                retval = arg.complex_array_value ().prod (dim);
+              break;
+            case btyp_float:
+              if (isdouble)
+                retval = arg.float_array_value ().dprod (dim);
+              else
+                retval = arg.float_array_value ().prod (dim);
+              break;
+            case btyp_float_complex:
+              if (isdouble)
+                retval = arg.float_complex_array_value ().dprod (dim);
+              else
+                retval = arg.float_complex_array_value ().prod (dim);
+              break;
+
+#define MAKE_INT_BRANCH(X) \
+            case btyp_ ## X: \
+              if (isnative) \
+                retval = arg.X ## _array_value ().prod (dim); \
+              else \
+                retval = arg.array_value ().prod (dim); \
+              break;
+            MAKE_INT_BRANCH (int8);
+            MAKE_INT_BRANCH (int16);
+            MAKE_INT_BRANCH (int32);
+            MAKE_INT_BRANCH (int64);
+            MAKE_INT_BRANCH (uint8);
+            MAKE_INT_BRANCH (uint16);
+            MAKE_INT_BRANCH (uint32);
+            MAKE_INT_BRANCH (uint64);
+#undef MAKE_INT_BRANCH
+
+            // GAGME: Accursed Matlab compatibility...
+            case btyp_char:
+              retval = arg.array_value (true).prod (dim);
+              break;
+            case btyp_bool:
+              if (arg.is_sparse_type ())
+                {
+                  if (isnative)
+                    retval = arg.sparse_bool_matrix_value ().all (dim);
+                  else
+                    retval = arg.sparse_matrix_value ().prod (dim);
+                }
+              else if (isnative)
+                retval = arg.bool_array_value ().all (dim);
+              else
+                retval = NDArray (arg.bool_array_value ().all (dim));
+              break;
+
+            default:
+              gripe_wrong_type_arg ("prod", arg);
+            }
+        }
+    }
+  else
+    print_usage ();
+
+  return retval;
 }
 
 /*
@@ -1398,6 +1530,13 @@
 %!assert (prod (single ([i, 2+i, -3+2i, 4])), single (-4 - 32i))
 %!assert (prod (single ([1, 2, 3; i, 2i, 3i; 1+i, 2+2i, 3+3i])), single ([-1+i, -8+8i, -27+27i]))
 
+%% Test sparse
+%!assert (prod (sparse ([1, 2, 3])), sparse (6))
+%!assert (prod (sparse ([-1; -2; -3])), sparse (-6))
+## Commented out until bug #42290 is fixed
+#%!assert (prod (sparse ([i, 2+i, -3+2i, 4])), sparse (-4 - 32i))
+#%!assert (prod (sparse ([1, 2, 3; i, 2i, 3i; 1+i, 2+2i, 3+3i])), sparse ([-1+i, -8+8i, -27+27i]))
+
 %!assert (prod ([1, 2; 3, 4], 1), [3, 8])
 %!assert (prod ([1, 2; 3, 4], 2), [2; 12])
 %!assert (prod (zeros (1, 0)), 1)
@@ -1428,7 +1567,24 @@
 %!assert (prod (zeros (0, 2, "single"), 1), single ([1, 1]))
 %!assert (prod (zeros (0, 2, "single"), 2), zeros (0, 1, "single"))
 
+%% Test "double" type argument
+%!assert (prod (single ([1, 2, 3]), "double"), 6)
+%!assert (prod (single ([-1; -2; -3]), "double"), -6)
+%!assert (prod (single ([i, 2+i, -3+2i, 4]), "double"), -4 - 32i)
+%!assert (prod (single ([1, 2, 3; i, 2i, 3i; 1+i, 2+2i, 3+3i]), "double"), [-1+i, -8+8i, -27+27i])
+
+%% Test "native" type argument
+%!assert (prod (uint8 ([1, 2, 3]), "native"), uint8 (6))
+%!assert (prod (uint8 ([-1; -2; -3]), "native"), uint8 (0))
+%!assert (prod (int8 ([1, 2, 3]), "native"), int8 (6))
+%!assert (prod (int8 ([-1; -2; -3]), "native"), int8 (-6))
+%!assert (prod ([true false; true true], "native"), [true false])
+%!assert (prod ([true false; true true], 2, "native"), [false; true])
+
+%% Test input validation
 %!error prod ()
+%!error prod (1,2,3)
+%!error <unrecognized type argument 'foobar'> prod (1, "foobar")
 */
 
 static bool
@@ -2807,7 +2963,8 @@
             MAKE_INT_BRANCH (uint32);
             MAKE_INT_BRANCH (uint64);
 #undef MAKE_INT_BRANCH
-              // GAGME: Accursed Matlab compatibility...
+            
+            // GAGME: Accursed Matlab compatibility...
             case btyp_char:
               if (isextra)
                 retval = arg.array_value (true).xsum (dim);
diff -r 87fafe6d05d3 -r 8cc66f091584 liboctave/array/fCNDArray.cc
--- a/liboctave/array/fCNDArray.cc	Sat May 03 17:28:49 2014 -0700
+++ b/liboctave/array/fCNDArray.cc	Sun May 04 16:08:46 2014 -0700
@@ -616,6 +616,12 @@
   return do_mx_red_op<FloatComplex, FloatComplex> (*this, dim, mx_inline_prod);
 }
 
+ComplexNDArray
+FloatComplexNDArray::dprod (int dim) const
+{
+  return do_mx_red_op<Complex, FloatComplex> (*this, dim, mx_inline_dprod);
+}
+
 FloatComplexNDArray
 FloatComplexNDArray::sum (int dim) const
 {
diff -r 87fafe6d05d3 -r 8cc66f091584 liboctave/array/fCNDArray.h
--- a/liboctave/array/fCNDArray.h	Sat May 03 17:28:49 2014 -0700
+++ b/liboctave/array/fCNDArray.h	Sun May 04 16:08:46 2014 -0700
@@ -83,6 +83,7 @@
   FloatComplexNDArray cumprod (int dim = -1) const;
   FloatComplexNDArray cumsum (int dim = -1) const;
   FloatComplexNDArray prod (int dim = -1) const;
+  ComplexNDArray dprod (int dim = -1) const;
   FloatComplexNDArray sum (int dim = -1) const;
   ComplexNDArray dsum (int dim = -1) const;
   FloatComplexNDArray sumsq (int dim = -1) const;
diff -r 87fafe6d05d3 -r 8cc66f091584 liboctave/array/fNDArray.cc
--- a/liboctave/array/fNDArray.cc	Sat May 03 17:28:49 2014 -0700
+++ b/liboctave/array/fNDArray.cc	Sun May 04 16:08:46 2014 -0700
@@ -627,6 +627,12 @@
   return do_mx_red_op<float, float> (*this, dim, mx_inline_prod);
 }
 
+NDArray
+FloatNDArray::dprod (int dim) const
+{
+  return do_mx_red_op<double, float> (*this, dim, mx_inline_dprod);
+}
+
 FloatNDArray
 FloatNDArray::sum (int dim) const
 {
diff -r 87fafe6d05d3 -r 8cc66f091584 liboctave/array/fNDArray.h
--- a/liboctave/array/fNDArray.h	Sat May 03 17:28:49 2014 -0700
+++ b/liboctave/array/fNDArray.h	Sun May 04 16:08:46 2014 -0700
@@ -90,6 +90,7 @@
   FloatNDArray cumprod (int dim = -1) const;
   FloatNDArray cumsum (int dim = -1) const;
   FloatNDArray prod (int dim = -1) const;
+  NDArray dprod (int dim = -1) const;
   FloatNDArray sum (int dim = -1) const;
   NDArray dsum (int dim = -1) const;
   FloatNDArray sumsq (int dim = -1) const;
diff -r 87fafe6d05d3 -r 8cc66f091584 liboctave/array/intNDArray.cc
--- a/liboctave/array/intNDArray.cc	Sat May 03 17:28:49 2014 -0700
+++ b/liboctave/array/intNDArray.cc	Sun May 04 16:08:46 2014 -0700
@@ -213,6 +213,13 @@
 
 template <class T>
 intNDArray<T>
+intNDArray<T>::prod (int dim) const
+{
+  return do_mx_red_op<T, T> (*this, dim, mx_inline_prod);
+}
+
+template <class T>
+intNDArray<T>
 intNDArray<T>::sum (int dim) const
 {
   return do_mx_red_op<T, T> (*this, dim, mx_inline_sum);
diff -r 87fafe6d05d3 -r 8cc66f091584 liboctave/array/intNDArray.h
--- a/liboctave/array/intNDArray.h	Sat May 03 17:28:49 2014 -0700
+++ b/liboctave/array/intNDArray.h	Sun May 04 16:08:46 2014 -0700
@@ -89,6 +89,7 @@
   intNDArray cummin (int dim = -1) const;
   intNDArray cummin (Array<octave_idx_type>& index, int dim = -1) const;
 
+  intNDArray prod (int dim) const;
   intNDArray sum (int dim) const;
   NDArray dsum (int dim) const;
   intNDArray cumsum (int dim) const;
diff -r 87fafe6d05d3 -r 8cc66f091584 liboctave/operators/mx-inlines.cc
--- a/liboctave/operators/mx-inlines.cc	Sat May 03 17:28:49 2014 -0700
+++ b/liboctave/operators/mx-inlines.cc	Sun May 04 16:08:46 2014 -0700
@@ -485,6 +485,14 @@
 #define OP_RED_SUMSQ(ac, el) ac += el*el
 #define OP_RED_SUMSQC(ac, el) ac += cabsq (el)
 
+inline void op_dble_prod (double& ac, float el)
+{ ac *= el; }
+inline void op_dble_prod (Complex& ac, const FloatComplex& el)
+{ ac *= el; } // FIXME: guaranteed?
+template <class T>
+inline void op_dble_prod (double& ac, const octave_int<T>& el)
+{ ac *= el.double_value (); }
+
 inline void op_dble_sum (double& ac, float el)
 { ac += el; }
 inline void op_dble_sum (Complex& ac, const FloatComplex& el)
@@ -514,6 +522,7 @@
 OP_RED_FCN (mx_inline_dsum, T, PROMOTE_DOUBLE(T), op_dble_sum, 0.0)
 OP_RED_FCN (mx_inline_count, bool, T, OP_RED_SUM, 0)
 OP_RED_FCN (mx_inline_prod, T, T, OP_RED_PROD, 1)
+OP_RED_FCN (mx_inline_dprod, T, PROMOTE_DOUBLE(T), op_dble_prod, 1)
 OP_RED_FCN (mx_inline_sumsq, T, T, OP_RED_SUMSQ, 0)
 OP_RED_FCN (mx_inline_sumsq, std::complex<T>, T, OP_RED_SUMSQC, 0)
 OP_RED_FCN (mx_inline_any, T, bool, OP_RED_ANYC, false)
@@ -539,6 +548,7 @@
 OP_RED_FCN2 (mx_inline_dsum, T, PROMOTE_DOUBLE(T), op_dble_sum, 0.0)
 OP_RED_FCN2 (mx_inline_count, bool, T, OP_RED_SUM, 0)
 OP_RED_FCN2 (mx_inline_prod, T, T, OP_RED_PROD, 1)
+OP_RED_FCN2 (mx_inline_dprod, T, PROMOTE_DOUBLE(T), op_dble_prod, 0.0)
 OP_RED_FCN2 (mx_inline_sumsq, T, T, OP_RED_SUMSQ, 0)
 OP_RED_FCN2 (mx_inline_sumsq, std::complex<T>, T, OP_RED_SUMSQC, 0)
 
@@ -612,6 +622,7 @@
 OP_RED_FCNN (mx_inline_dsum, T, PROMOTE_DOUBLE(T))
 OP_RED_FCNN (mx_inline_count, bool, T)
 OP_RED_FCNN (mx_inline_prod, T, T)
+OP_RED_FCNN (mx_inline_dprod, T, PROMOTE_DOUBLE(T))
 OP_RED_FCNN (mx_inline_sumsq, T, T)
 OP_RED_FCNN (mx_inline_sumsq, std::complex<T>, T)
 OP_RED_FCNN (mx_inline_any, T, bool)