# HG changeset patch
# User Jaroslav Hajek <highegg@gmail.com>
# Date 1253280429 -7200
# Node ID e087d7c77ff949dd4d16cedfe6ce676fad804373
# Parent  ecdb275bd41b9015aa7f067119292beed5adfe17
improve linspace in liboctave

diff --git a/liboctave/Array.h b/liboctave/Array.h
--- a/liboctave/Array.h
+++ b/liboctave/Array.h
@@ -679,6 +679,39 @@
   static void instantiation_guard ();
 };
 
+// This is a simple wrapper template that will subclass an Array<T> type or any
+// later type derived from it and override the default non-const operator() to
+// not check for the array's uniqueness. It is, however, the user's
+// responsibility to ensure the array is actually unaliased whenever elements
+// are accessed.
+
+template<class ArrayClass>
+class NoAlias : public ArrayClass
+{
+  typedef typename ArrayClass::element_type T;
+public:
+  NoAlias () : ArrayClass () { }
+
+  // FIXME: this would be simpler once C++0x is available
+  template <class X>
+    explicit NoAlias (X x) : ArrayClass (x) { }
+
+  template <class X, class Y>
+    explicit NoAlias (X x, Y y) : ArrayClass (x, y) { }
+
+  template <class X, class Y, class Z>
+    explicit NoAlias (X x, Y y, Z z) : ArrayClass (x, y, z) { }
+
+  T& operator () (octave_idx_type n)
+    { return ArrayClass::xelem (n); }
+  T& operator () (octave_idx_type i, octave_idx_type j)
+    { return ArrayClass::xelem (i, j); }
+  T& operator () (octave_idx_type i, octave_idx_type j, octave_idx_type k)
+    { return ArrayClass::xelem (i, j, k); }
+  T& operator () (const Array<octave_idx_type>& ra_idx)
+    { return ArrayClass::xelem (ra_idx); }
+};
+
 #endif
 
 /*
diff --git a/liboctave/CMatrix.cc b/liboctave/CMatrix.cc
--- a/liboctave/CMatrix.cc
+++ b/liboctave/CMatrix.cc
@@ -4066,6 +4066,39 @@
   return result;
 }
 
+ComplexMatrix linspace (const ComplexColumnVector& x1, 
+                        const ComplexColumnVector& x2, 
+                        octave_idx_type n)
+
+{
+  if (n < 1) n = 1;
+
+  octave_idx_type m = x1.length ();
+
+  if (x2.length () != m)
+    (*current_liboctave_error_handler) ("linspace: vectors must be of equal length");
+
+  NoAlias<ComplexMatrix> retval;
+
+  retval.clear (m, n);
+  for (octave_idx_type i = 0; i < m; i++)
+    retval(i, 0) = x1(i);
+
+  // The last column is not needed while using delta.
+  Complex *delta = &retval(0, 1); 
+  for (octave_idx_type i = 0; i < m; i++)
+    delta[i] = (x2(i) - x1(i)) / (n - 1.0);
+
+  for (octave_idx_type j = 1; j < n-1; j++)
+    for (octave_idx_type i = 0; i < m; i++)
+      retval(i, j) = retval(i, j-1) + delta[i];
+
+  for (octave_idx_type i = 0; i < m; i++)
+    retval(i, n-1) = x2(i);
+
+  return retval;
+}
+
 MS_CMP_OPS (ComplexMatrix, Complex)
 MS_BOOL_OPS (ComplexMatrix, Complex)
 
diff --git a/liboctave/CMatrix.h b/liboctave/CMatrix.h
--- a/liboctave/CMatrix.h
+++ b/liboctave/CMatrix.h
@@ -408,6 +408,11 @@
 extern OCTAVE_API ComplexMatrix max (const ComplexMatrix& m, const Complex& c);
 extern OCTAVE_API ComplexMatrix max (const ComplexMatrix& a, const ComplexMatrix& b);
 
+extern OCTAVE_API ComplexMatrix linspace (const ComplexColumnVector& x1, 
+                                          const ComplexColumnVector& x2, 
+                                          octave_idx_type n);
+
+
 MS_CMP_OP_DECLS (ComplexMatrix, Complex, OCTAVE_API)
 MS_BOOL_OP_DECLS (ComplexMatrix, Complex, OCTAVE_API)
 
diff --git a/liboctave/CRowVector.cc b/liboctave/CRowVector.cc
--- a/liboctave/CRowVector.cc
+++ b/liboctave/CRowVector.cc
@@ -482,22 +482,15 @@
 ComplexRowVector
 linspace (const Complex& x1, const Complex& x2, octave_idx_type n)
 {
-  ComplexRowVector retval;
+  if (n < 1) n = 1;
+
+  NoAlias<ComplexRowVector> retval (n);
 
-  if (n > 0)
-    {
-      retval.resize (n);
-      Complex delta = (x2 - x1) / (n - 1.0);
-      retval.elem (0) = x1;
-      for (octave_idx_type i = 1; i < n-1; i++)
-	retval.elem (i) = x1 + 1.0 * i * delta;
-      retval.elem (n-1) = x2;
-    }
-  else
-    {
-      retval.resize (1);
-      retval.elem (0) = x2;
-    }
+  Complex delta = (x2 - x1) / (n - 1.0);
+  Complex y = retval(0) = x1;
+  for (octave_idx_type i = 1; i < n-1; i++)
+    retval(i) = y += delta;
+  retval(n-1) = x2;
 
   return retval;
 }
diff --git a/liboctave/ChangeLog b/liboctave/ChangeLog
--- a/liboctave/ChangeLog
+++ b/liboctave/ChangeLog
@@ -1,3 +1,19 @@
+2009-09-18  Jaroslav Hajek  <highegg@gmail.com>
+
+	* Array.h (NoAlias): New template class.
+	* dRowVector.cc (linspace): Rewrite.
+	* fRowVector.cc (linspace): Rewrite.
+	* CRowVector.cc (linspace): Rewrite.
+	* fCRowVector.cc (linspace): Rewrite.
+	* dMatrix.cc (linspace): New method.
+	* dMatrix.h (linspace): Declare it.
+	* fMatrix.cc (linspace): New method.
+	* fMatrix.h (linspace): Declare it.
+	* CMatrix.cc (linspace): New method.
+	* CMatrix.h (linspace): Declare it.
+	* fCMatrix.cc (linspace): New method.
+	* fCMatrix.h (linspace): Declare it.
+
 2009-09-17  John W. Eaton  <jwe@octave.org>
 
 	* oct-types.h.in: Delete.
diff --git a/liboctave/dMatrix.cc b/liboctave/dMatrix.cc
--- a/liboctave/dMatrix.cc
+++ b/liboctave/dMatrix.cc
@@ -3395,6 +3395,39 @@
   return result;
 }
 
+Matrix linspace (const ColumnVector& x1, 
+                 const ColumnVector& x2, 
+                 octave_idx_type n)
+
+{
+  if (n < 1) n = 1;
+
+  octave_idx_type m = x1.length ();
+
+  if (x2.length () != m)
+    (*current_liboctave_error_handler) ("linspace: vectors must be of equal length");
+
+  NoAlias<Matrix> retval;
+
+  retval.clear (m, n);
+  for (octave_idx_type i = 0; i < m; i++)
+    retval(i, 0) = x1(i);
+
+  // The last column is not needed while using delta.
+  double *delta = &retval(0, 1); 
+  for (octave_idx_type i = 0; i < m; i++)
+    delta[i] = (x2(i) - x1(i)) / (n - 1);
+
+  for (octave_idx_type j = 1; j < n-1; j++)
+    for (octave_idx_type i = 0; i < m; i++)
+      retval(i, j) = retval(i, j-1) + delta[i];
+
+  for (octave_idx_type i = 0; i < m; i++)
+    retval(i, n-1) = x2(i);
+
+  return retval;
+}
+
 MS_CMP_OPS (Matrix, double)
 MS_BOOL_OPS (Matrix, double)
 
diff --git a/liboctave/dMatrix.h b/liboctave/dMatrix.h
--- a/liboctave/dMatrix.h
+++ b/liboctave/dMatrix.h
@@ -357,6 +357,10 @@
 extern OCTAVE_API Matrix max (const Matrix& m, double d);
 extern OCTAVE_API Matrix max (const Matrix& a, const Matrix& b);
 
+extern OCTAVE_API Matrix linspace (const ColumnVector& x1, 
+                                   const ColumnVector& x2, 
+                                   octave_idx_type n);
+
 MS_CMP_OP_DECLS (Matrix, double, OCTAVE_API)
 MS_BOOL_OP_DECLS (Matrix, double, OCTAVE_API)
 
diff --git a/liboctave/dRowVector.cc b/liboctave/dRowVector.cc
--- a/liboctave/dRowVector.cc
+++ b/liboctave/dRowVector.cc
@@ -311,22 +311,15 @@
 RowVector
 linspace (double x1, double x2, octave_idx_type n)
 {
-  RowVector retval;
+  if (n < 1) n = 1;
+
+  NoAlias<RowVector> retval (n);
 
-  if (n > 1)
-    {
-      retval.resize (n);
-      double delta = (x2 - x1) / (n - 1);
-      retval.elem (0) = x1;
-      for (octave_idx_type i = 1; i < n-1; i++)
-	retval.elem (i) = x1 + i * delta;
-      retval.elem (n-1) = x2;
-    }
-  else
-    {
-      retval.resize (1);
-      retval.elem (0) = x2;
-    }
+  double delta = (x2 - x1) / (n - 1);
+  double y = retval(0) = x1;
+  for (octave_idx_type i = 1; i < n-1; i++)
+    retval(i) = y += delta;
+  retval(n-1) = x2;
 
   return retval;
 }
diff --git a/liboctave/fCMatrix.cc b/liboctave/fCMatrix.cc
--- a/liboctave/fCMatrix.cc
+++ b/liboctave/fCMatrix.cc
@@ -4059,6 +4059,39 @@
   return result;
 }
 
+FloatComplexMatrix linspace (const FloatComplexColumnVector& x1, 
+                             const FloatComplexColumnVector& x2, 
+                             octave_idx_type n)
+
+{
+  if (n < 1) n = 1;
+
+  octave_idx_type m = x1.length ();
+
+  if (x2.length () != m)
+    (*current_liboctave_error_handler) ("linspace: vectors must be of equal length");
+
+  NoAlias<FloatComplexMatrix> retval;
+
+  retval.clear (m, n);
+  for (octave_idx_type i = 0; i < m; i++)
+    retval(i, 0) = x1(i);
+
+  // The last column is not needed while using delta.
+  FloatComplex *delta = &retval(0, 1); 
+  for (octave_idx_type i = 0; i < m; i++)
+    delta[i] = (x2(i) - x1(i)) / (n - 1.0f);
+
+  for (octave_idx_type j = 1; j < n-1; j++)
+    for (octave_idx_type i = 0; i < m; i++)
+      retval(i, j) = retval(i, j-1) + delta[i];
+
+  for (octave_idx_type i = 0; i < m; i++)
+    retval(i, n-1) = x2(i);
+
+  return retval;
+}
+
 MS_CMP_OPS (FloatComplexMatrix, FloatComplex)
 MS_BOOL_OPS (FloatComplexMatrix, FloatComplex)
 
diff --git a/liboctave/fCMatrix.h b/liboctave/fCMatrix.h
--- a/liboctave/fCMatrix.h
+++ b/liboctave/fCMatrix.h
@@ -408,6 +408,10 @@
 extern OCTAVE_API FloatComplexMatrix max (const FloatComplexMatrix& m, const FloatComplex& c);
 extern OCTAVE_API FloatComplexMatrix max (const FloatComplexMatrix& a, const FloatComplexMatrix& b);
 
+extern OCTAVE_API FloatComplexMatrix linspace (const FloatComplexColumnVector& x1, 
+                                               const FloatComplexColumnVector& x2, 
+                                               octave_idx_type n);
+
 MS_CMP_OP_DECLS (FloatComplexMatrix, FloatComplex, OCTAVE_API)
 MS_BOOL_OP_DECLS (FloatComplexMatrix, FloatComplex, OCTAVE_API)
 
diff --git a/liboctave/fCRowVector.cc b/liboctave/fCRowVector.cc
--- a/liboctave/fCRowVector.cc
+++ b/liboctave/fCRowVector.cc
@@ -482,22 +482,15 @@
 FloatComplexRowVector
 linspace (const FloatComplex& x1, const FloatComplex& x2, octave_idx_type n)
 {
-  FloatComplexRowVector retval;
+  if (n < 1) n = 1;
+
+  NoAlias<FloatComplexRowVector> retval (n);
 
-  if (n > 0)
-    {
-      retval.resize (n);
-      FloatComplex delta = (x2 - x1) / static_cast<float> (n - 1.0);
-      retval.elem (0) = x1;
-      for (octave_idx_type i = 1; i < n-1; i++)
-	retval.elem (i) = x1 +  static_cast<float> (1.0) * i * delta;
-      retval.elem (n-1) = x2;
-    }
-  else
-    {
-      retval.resize (1);
-      retval.elem (0) = x2;
-    }
+  FloatComplex delta = (x2 - x1) / (n - 1.0f);
+  FloatComplex y = retval(0) = x1;
+  for (octave_idx_type i = 1; i < n-1; i++)
+    retval(i) = y += delta;
+  retval(n-1) = x2;
 
   return retval;
 }
diff --git a/liboctave/fMatrix.cc b/liboctave/fMatrix.cc
--- a/liboctave/fMatrix.cc
+++ b/liboctave/fMatrix.cc
@@ -3394,6 +3394,39 @@
   return result;
 }
 
+FloatMatrix linspace (const FloatColumnVector& x1, 
+                      const FloatColumnVector& x2, 
+                      octave_idx_type n)
+
+{
+  if (n < 1) n = 1;
+
+  octave_idx_type m = x1.length ();
+
+  if (x2.length () != m)
+    (*current_liboctave_error_handler) ("linspace: vectors must be of equal length");
+
+  NoAlias<FloatMatrix> retval;
+
+  retval.clear (m, n);
+  for (octave_idx_type i = 0; i < m; i++)
+    retval(i, 0) = x1(i);
+
+  // The last column is not needed while using delta.
+  float *delta = &retval(0, 1); 
+  for (octave_idx_type i = 0; i < m; i++)
+    delta[i] = (x2(i) - x1(i)) / (n - 1);
+
+  for (octave_idx_type j = 1; j < n-1; j++)
+    for (octave_idx_type i = 0; i < m; i++)
+      retval(i, j) = retval(i, j-1) + delta[i];
+
+  for (octave_idx_type i = 0; i < m; i++)
+    retval(i, n-1) = x2(i);
+
+  return retval;
+}
+
 MS_CMP_OPS (FloatMatrix, float)
 MS_BOOL_OPS (FloatMatrix, float)
 
diff --git a/liboctave/fMatrix.h b/liboctave/fMatrix.h
--- a/liboctave/fMatrix.h
+++ b/liboctave/fMatrix.h
@@ -358,6 +358,11 @@
 extern OCTAVE_API FloatMatrix max (const FloatMatrix& m, float d);
 extern OCTAVE_API FloatMatrix max (const FloatMatrix& a, const FloatMatrix& b);
 
+extern OCTAVE_API FloatMatrix linspace (const FloatColumnVector& x1, 
+                                        const FloatColumnVector& x2, 
+                                        octave_idx_type n);
+
+
 MS_CMP_OP_DECLS (FloatMatrix, float, OCTAVE_API)
 MS_BOOL_OP_DECLS (FloatMatrix, float, OCTAVE_API)
 
diff --git a/liboctave/fRowVector.cc b/liboctave/fRowVector.cc
--- a/liboctave/fRowVector.cc
+++ b/liboctave/fRowVector.cc
@@ -311,22 +311,15 @@
 FloatRowVector
 linspace (float x1, float x2, octave_idx_type n)
 {
-  FloatRowVector retval;
+  if (n < 1) n = 1;
+
+  NoAlias<FloatRowVector> retval (n);
 
-  if (n > 1)
-    {
-      retval.resize (n);
-      float delta = (x2 - x1) / (n - 1);
-      retval.elem (0) = x1;
-      for (octave_idx_type i = 1; i < n-1; i++)
-	retval.elem (i) = x1 + i * delta;
-      retval.elem (n-1) = x2;
-    }
-  else
-    {
-      retval.resize (1);
-      retval.elem (0) = x2;
-    }
+  float delta = (x2 - x1) / (n - 1);
+  float y = retval(0) = x1;
+  for (octave_idx_type i = 1; i < n-1; i++)
+    retval(i) = y += delta;
+  retval(n-1) = x2;
 
   return retval;
 }