# HG changeset patch
# User Jaroslav Hajek <highegg@gmail.com>
# Date 1258453356 -3600
# Node ID 9b62f2d8de6dfcdb689e3e7d32faa9cbb39a310c
# Parent  64270d3ad46985ad67c613772c53d0736334f756
improve r->c mapper strategy

diff -r 64270d3ad469 -r 9b62f2d8de6d src/ChangeLog
--- a/src/ChangeLog	Tue Nov 17 08:12:09 2009 -0500
+++ b/src/ChangeLog	Tue Nov 17 11:22:36 2009 +0100
@@ -1,3 +1,10 @@
+2009-11-17  Jaroslav Hajek  <highegg@gmail.com>
+
+	* ov-re-mat.cc (do_rc_map): New static function.
+	(octave_matrix::map): Use it here.
+	* ov-flt-re-mat.cc (do_rc_map): New static function.
+	(octave_float_matrix::map): Use it here.
+
 2009-11-17  Jaroslav Hajek  <highegg@gmail.com>
 
 	* DLD-FUNCTIONS/balance.cc: Fix docs.
diff -r 64270d3ad469 -r 9b62f2d8de6d src/ov-flt-re-mat.cc
--- a/src/ov-flt-re-mat.cc	Tue Nov 17 08:12:09 2009 -0500
+++ b/src/ov-flt-re-mat.cc	Tue Nov 17 11:22:36 2009 +0100
@@ -680,6 +680,45 @@
   return retval;
 }
 
+// This uses a smarter strategy for doing the complex->real mappers.  We
+// allocate an array for a real result and keep filling it until a complex
+// result is produced.
+static octave_value
+do_rc_map (const FloatNDArray& a, FloatComplex (&fcn) (float))
+{
+  octave_idx_type n = a.numel ();
+  NoAlias<FloatNDArray> rr (a.dims ());
+
+  for (octave_idx_type i = 0; i < n; i++)
+    {
+      OCTAVE_QUIT;
+
+      FloatComplex tmp = fcn (a(i));
+      if (tmp.imag () == 0.0)
+        rr(i) = tmp.real ();
+      else
+        {
+          NoAlias<FloatComplexNDArray> rc (a.dims ());
+
+          for (octave_idx_type j = 0; j < i; j++)
+            rc(j) = rr(j);
+
+          rc(i) = tmp;
+
+          for (octave_idx_type j = i+1; j < n; j++)
+            {
+              OCTAVE_QUIT;
+
+              rc(j) = fcn (a(i));
+            }
+
+          return new octave_float_complex_matrix (rc);
+        }
+    }
+
+  return rr;
+}
+
 octave_value
 octave_float_matrix::map (unary_mapper_t umap) const
 {
@@ -706,18 +745,22 @@
     case umap_ ## UMAP: \
       return octave_value (matrix.map<TYPE> (FCN))
 
-      ARRAY_MAPPER (acos, FloatComplex, rc_acos);
-      ARRAY_MAPPER (acosh, FloatComplex, rc_acosh);
+#define RC_ARRAY_MAPPER(UMAP, TYPE, FCN) \
+    case umap_ ## UMAP: \
+      return do_rc_map (matrix, FCN)
+
+      RC_ARRAY_MAPPER (acos, FloatComplex, rc_acos);
+      RC_ARRAY_MAPPER (acosh, FloatComplex, rc_acosh);
       ARRAY_MAPPER (angle, float, ::arg);
       ARRAY_MAPPER (arg, float, ::arg);
-      ARRAY_MAPPER (asin, FloatComplex, rc_asin);
+      RC_ARRAY_MAPPER (asin, FloatComplex, rc_asin);
       ARRAY_MAPPER (asinh, float, ::asinhf);
       ARRAY_MAPPER (atan, float, ::atanf);
-      ARRAY_MAPPER (atanh, FloatComplex, rc_atanh);
+      RC_ARRAY_MAPPER (atanh, FloatComplex, rc_atanh);
       ARRAY_MAPPER (erf, float, ::erff);
       ARRAY_MAPPER (erfc, float, ::erfcf);
       ARRAY_MAPPER (gamma, float, xgamma);
-      ARRAY_MAPPER (lgamma, FloatComplex, rc_lgamma);
+      RC_ARRAY_MAPPER (lgamma, FloatComplex, rc_lgamma);
       ARRAY_MAPPER (ceil, float, ::ceilf);
       ARRAY_MAPPER (cos, float, ::cosf);
       ARRAY_MAPPER (cosh, float, ::coshf);
@@ -725,16 +768,16 @@
       ARRAY_MAPPER (expm1, float, ::expm1f);
       ARRAY_MAPPER (fix, float, ::fix);
       ARRAY_MAPPER (floor, float, ::floorf);
-      ARRAY_MAPPER (log, FloatComplex, rc_log);
-      ARRAY_MAPPER (log2, FloatComplex, rc_log2);
-      ARRAY_MAPPER (log10, FloatComplex, rc_log10);
-      ARRAY_MAPPER (log1p, FloatComplex, rc_log1p);
+      RC_ARRAY_MAPPER (log, FloatComplex, rc_log);
+      RC_ARRAY_MAPPER (log2, FloatComplex, rc_log2);
+      RC_ARRAY_MAPPER (log10, FloatComplex, rc_log10);
+      RC_ARRAY_MAPPER (log1p, FloatComplex, rc_log1p);
       ARRAY_MAPPER (round, float, xround);
       ARRAY_MAPPER (roundb, float, xroundb);
       ARRAY_MAPPER (signum, float, ::signum);
       ARRAY_MAPPER (sin, float, ::sinf);
       ARRAY_MAPPER (sinh, float, ::sinhf);
-      ARRAY_MAPPER (sqrt, FloatComplex, rc_sqrt);
+      RC_ARRAY_MAPPER (sqrt, FloatComplex, rc_sqrt);
       ARRAY_MAPPER (tan, float, ::tanf);
       ARRAY_MAPPER (tanh, float, ::tanhf);
       ARRAY_MAPPER (isna, bool, octave_is_NA);
diff -r 64270d3ad469 -r 9b62f2d8de6d src/ov-re-mat.cc
--- a/src/ov-re-mat.cc	Tue Nov 17 08:12:09 2009 -0500
+++ b/src/ov-re-mat.cc	Tue Nov 17 11:22:36 2009 +0100
@@ -708,6 +708,45 @@
   return retval;
 }
 
+// This uses a smarter strategy for doing the complex->real mappers.  We
+// allocate an array for a real result and keep filling it until a complex
+// result is produced.
+static octave_value
+do_rc_map (const NDArray& a, Complex (&fcn) (double))
+{
+  octave_idx_type n = a.numel ();
+  NoAlias<NDArray> rr (a.dims ());
+
+  for (octave_idx_type i = 0; i < n; i++)
+    {
+      OCTAVE_QUIT;
+
+      Complex tmp = fcn (a(i));
+      if (tmp.imag () == 0.0)
+        rr(i) = tmp.real ();
+      else
+        {
+          NoAlias<ComplexNDArray> rc (a.dims ());
+
+          for (octave_idx_type j = 0; j < i; j++)
+            rc(j) = rr(j);
+
+          rc(i) = tmp;
+
+          for (octave_idx_type j = i+1; j < n; j++)
+            {
+              OCTAVE_QUIT;
+
+              rc(j) = fcn (a(i));
+            }
+
+          return new octave_complex_matrix (rc);
+        }
+    }
+
+  return rr;
+}
+
 octave_value
 octave_matrix::map (unary_mapper_t umap) const
 {
@@ -734,18 +773,22 @@
     case umap_ ## UMAP: \
       return octave_value (matrix.map<TYPE> (FCN))
 
-      ARRAY_MAPPER (acos, Complex, rc_acos);
-      ARRAY_MAPPER (acosh, Complex, rc_acosh);
+#define RC_ARRAY_MAPPER(UMAP, TYPE, FCN) \
+    case umap_ ## UMAP: \
+      return do_rc_map (matrix, FCN)
+
+      RC_ARRAY_MAPPER (acos, Complex, rc_acos);
+      RC_ARRAY_MAPPER (acosh, Complex, rc_acosh);
       ARRAY_MAPPER (angle, double, ::arg);
       ARRAY_MAPPER (arg, double, ::arg);
-      ARRAY_MAPPER (asin, Complex, rc_asin);
+      RC_ARRAY_MAPPER (asin, Complex, rc_asin);
       ARRAY_MAPPER (asinh, double, ::asinh);
       ARRAY_MAPPER (atan, double, ::atan);
-      ARRAY_MAPPER (atanh, Complex, rc_atanh);
+      RC_ARRAY_MAPPER (atanh, Complex, rc_atanh);
       ARRAY_MAPPER (erf, double, ::erf);
       ARRAY_MAPPER (erfc, double, ::erfc);
       ARRAY_MAPPER (gamma, double, xgamma);
-      ARRAY_MAPPER (lgamma, Complex, rc_lgamma);
+      RC_ARRAY_MAPPER (lgamma, Complex, rc_lgamma);
       ARRAY_MAPPER (ceil, double, ::ceil);
       ARRAY_MAPPER (cos, double, ::cos);
       ARRAY_MAPPER (cosh, double, ::cosh);
@@ -753,16 +796,16 @@
       ARRAY_MAPPER (expm1, double, ::expm1);
       ARRAY_MAPPER (fix, double, ::fix);
       ARRAY_MAPPER (floor, double, ::floor);
-      ARRAY_MAPPER (log, Complex, rc_log);
-      ARRAY_MAPPER (log2, Complex, rc_log2);
-      ARRAY_MAPPER (log10, Complex, rc_log10);
-      ARRAY_MAPPER (log1p, Complex, rc_log1p);
+      RC_ARRAY_MAPPER (log, Complex, rc_log);
+      RC_ARRAY_MAPPER (log2, Complex, rc_log2);
+      RC_ARRAY_MAPPER (log10, Complex, rc_log10);
+      RC_ARRAY_MAPPER (log1p, Complex, rc_log1p);
       ARRAY_MAPPER (round, double, xround);
       ARRAY_MAPPER (roundb, double, xroundb);
       ARRAY_MAPPER (signum, double, ::signum);
       ARRAY_MAPPER (sin, double, ::sin);
       ARRAY_MAPPER (sinh, double, ::sinh);
-      ARRAY_MAPPER (sqrt, Complex, rc_sqrt);
+      RC_ARRAY_MAPPER (sqrt, Complex, rc_sqrt);
       ARRAY_MAPPER (tan, double, ::tan);
       ARRAY_MAPPER (tanh, double, ::tanh);
       ARRAY_MAPPER (isna, bool, octave_is_NA);