Implemented the Blitz++ trick to allow the initialization of lists and matrices
authorChristian Bauer <Christian.Bauer@uni-mainz.de>
Tue, 30 Sep 2003 19:58:59 +0000 (19:58 +0000)
committerChristian Bauer <Christian.Bauer@uni-mainz.de>
Tue, 30 Sep 2003 19:58:59 +0000 (19:58 +0000)
from comma-separated sequences of expressions, like this:

  lst l;
  l = x, 2, y, x+y;

  matrix M(3, 3);
  M = x, y, 0,
     -y, x, 0,
      0, 0, 1;

This is both faster and produces much smaller code than the old constructors
lst(ex, ex, ...) and matrix(unsigned, unsigned, lst), especially in the case
of matrices, and is now the recommended way to create these objects.

check/exam_indexed.cpp
doc/tutorial/ginac.texi
ginac/container.h
ginac/matrix.h

index 69816a7..9deb487 100644 (file)
@@ -241,19 +241,17 @@ static unsigned edyn_check()
 
        // Lorentz transformation matrix (boost along x axis)
        matrix L(4, 4);
-       L(0, 0) = gamma;
-       L(0, 1) = -beta*gamma;
-       L(1, 0) = -beta*gamma;
-       L(1, 1) = gamma;
-       L(2, 2) = 1; L(3, 3) = 1;
+       L =       gamma, -beta*gamma, 0, 0,
+           -beta*gamma,       gamma, 0, 0,
+                     0,           0, 1, 0,
+                     0,           0, 0, 1;
 
        // Electromagnetic field tensor
-       matrix F(4, 4, lst(
-                0, -Ex, -Ey, -Ez,
+       matrix F(4, 4);
+       F =  0, -Ex, -Ey, -Ez,
                Ex,   0, -Bz,  By,
                Ey,  Bz,   0, -Bx,
-               Ez, -By,  Bx,   0
-       ));
+               Ez, -By,  Bx,   0;
 
        // Indices
        symbol s_mu("mu"), s_nu("nu"), s_rho("rho"), s_sigma("sigma");
index d07518e..e1f45b3 100644 (file)
@@ -1338,14 +1338,27 @@ packages, but are sometimes used to supply a variable number of arguments of
 the same type to GiNaC methods such as @code{subs()} and some @code{matrix}
 constructors, so you should have a basic understanding of them.
 
-Lists of up to 16 expressions can be directly constructed from single
+Lists can be constructed by assigning a comma-separated sequence of
 expressions:
 
 @example
 @{
     symbol x("x"), y("y");
-    lst l(x, 2, y, x+y);
-    // now, l is a list holding the expressions 'x', '2', 'y', and 'x+y'
+    lst l;
+    l = x, 2, y, x+y;
+    // now, l is a list holding the expressions 'x', '2', 'y', and 'x+y',
+    // in that order
+    ...
+@end example
+
+There are also constructors that allow direct creation of lists of up to
+16 expressions, which is often more convenient but slightly less efficient:
+
+@example
+    ...
+    // This produces the same list 'l' as above:
+    // lst l(x, 2, y, x+y);
+    // lst l = lst(x, 2, y, x+y);
     ...
 @end example
 
@@ -1400,10 +1413,10 @@ the C++ standard library:
 @example
     ...
     // print the elements of the list (requires #include <iterator>)
-    copy(l.begin(), l.end(), ostream_iterator<ex>(cout, "\n"));
+    std::copy(l.begin(), l.end(), ostream_iterator<ex>(cout, "\n"));
 
     // sum up the elements of the list (requires #include <numeric>)
-    ex sum = accumulate(l.begin(), l.end(), ex(0));
+    ex sum = std::accumulate(l.begin(), l.end(), ex(0));
     cout << sum << endl;  // prints '2+2*x+2*y'
     ...
 @end example
@@ -1450,8 +1463,9 @@ You can bring the elements of a list into a canonical order with @code{sort()}:
 
 @example
     ...
-    lst l1(x, 2, y, x+y);
-    lst l2(2, x+y, x, y);
+    lst l1, l2;
+    l1 = x, 2, y, x+y;
+    l2 = 2, x+y, x, y;
     l1.sort();
     l2.sort();
     // l1 and l2 are now equal
@@ -1463,7 +1477,8 @@ elements with @code{unique()}:
 
 @example
     ...
-    lst l3(x, 2, 2, 2, y, x+y, y+x);
+    lst l3;
+    l3 = x, 2, 2, 2, y, x+y, y+x;
     l3.unique();        // l3 is now @{x, 2, y, x+y@}
 @}
 @end example
@@ -1559,16 +1574,39 @@ matrix with @math{m} rows and @math{n} columns are accessed with two
 second one in the range 0@dots{}@math{n-1}.
 
 There are a couple of ways to construct matrices, with or without preset
-elements:
+elements. The constructor
+
+@example
+matrix::matrix(unsigned r, unsigned c);
+@end example
+
+creates a matrix with @samp{r} rows and @samp{c} columns with all elements
+set to zero.
+
+The fastest way to create a matrix with preinitialized elements is to assign
+a list of comma-separated expressions to an empty matrix (see below for an
+example). But you can also specify the elements as a (flat) list with
+
+@example
+matrix::matrix(unsigned r, unsigned c, const lst & l);
+@end example
+
+The function
 
 @cindex @code{lst_to_matrix()}
+@example
+ex lst_to_matrix(const lst & l);
+@end example
+
+constructs a matrix from a list of lists, each list representing a matrix row.
+
+There is also a set of functions for creating some special types of
+matrices:
+
 @cindex @code{diag_matrix()}
 @cindex @code{unit_matrix()}
 @cindex @code{symbolic_matrix()}
 @example
-matrix::matrix(unsigned r, unsigned c);
-matrix::matrix(unsigned r, unsigned c, const lst & l);
-ex lst_to_matrix(const lst & l);
 ex diag_matrix(const lst & l);
 ex unit_matrix(unsigned x);
 ex unit_matrix(unsigned r, unsigned c);
@@ -1576,16 +1614,11 @@ ex symbolic_matrix(unsigned r, unsigned c, const string & base_name);
 ex symbolic_matrix(unsigned r, unsigned c, const string & base_name, const string & tex_base_name);
 @end example
 
-The first two functions are @code{matrix} constructors which create a matrix
-with @samp{r} rows and @samp{c} columns. The matrix elements can be
-initialized from a (flat) list of expressions @samp{l}. Otherwise they are
-all set to zero. The @code{lst_to_matrix()} function constructs a matrix
-from a list of lists, each list representing a matrix row. @code{diag_matrix()}
-constructs a diagonal matrix given the list of diagonal elements.
-@code{unit_matrix()} creates an @samp{x} by @samp{x} (or @samp{r} by @samp{c})
-unit matrix. And finally, @code{symbolic_matrix} constructs a matrix filled
-with newly generated symbols made of the specified base name and the
-position of each element in the matrix.
+@code{diag_matrix()} constructs a diagonal matrix given the list of diagonal
+elements. @code{unit_matrix()} creates an @samp{x} by @samp{x} (or @samp{r}
+by @samp{c}) unit matrix. And finally, @code{symbolic_matrix} constructs a
+matrix filled with newly generated symbols made of the specified base name
+and the position of each element in the matrix.
 
 Matrix elements can be accessed and set using the parenthesis (function call)
 operator:
@@ -1599,18 +1632,24 @@ It is also possible to access the matrix elements in a linear fashion with
 the @code{op()} method. But C++-style subscripting with square brackets
 @samp{[]} is not available.
 
-Here are a couple of examples of constructing matrices:
+Here are a couple of examples for constructing matrices:
 
 @example
 @{
     symbol a("a"), b("b");
 
     matrix M(2, 2);
-    M(0, 0) = a;
-    M(1, 1) = b;
+    M = a, 0,
+        0, b;
     cout << M << endl;
      // -> [[a,0],[0,b]]
 
+    matrix M2(2, 2);
+    M2(0, 0) = a;
+    M2(1, 1) = b;
+    cout << M2 << endl;
+     // -> [[a,0],[0,b]]
+
     cout << matrix(2, 2, lst(a, 0, 0, b)) << endl;
      // -> [[a,0],[0,b]]
 
@@ -1647,9 +1686,13 @@ and @math{C}:
 
 @example
 @{
-    matrix A(2, 2, lst(1, 2, 3, 4));
-    matrix B(2, 2, lst(-1, 0, 2, 1));
-    matrix C(2, 2, lst(8, 4, 2, 1));
+    matrix A(2, 2), B(2, 2), C(2, 2);
+    A =  1, 2,
+         3, 4;
+    B = -1, 0,
+         2, 1;
+    C =  8, 4,
+         2, 1;
 
     matrix result = A.mul(B).sub(C.mul_scalar(2));
     cout << result << endl;
@@ -2448,7 +2491,10 @@ and scalar products):
     symbol x("x"), y("y");
 
     // A is a 2x2 matrix, X is a 2x1 vector
-    matrix A(2, 2, lst(1, 2, 3, 4)), X(2, 1, lst(x, y));
+    matrix A(2, 2), X(2, 1);
+    A = 1, 2,
+        3, 4;
+    X = x, y;
 
     cout << indexed(A, i, i) << endl;
      // -> 5
@@ -4792,16 +4838,16 @@ GiNaC contains the following predefined mathematical functions:
 @item @code{Order(x)}
 @tab order term function in truncated power series
 @cindex @code{Order()}
-@item @code{Li(n,x)}
+@item @code{Li(n, x)}
 @tab polylogarithm
 @cindex @code{Li()}
-@item @code{S(n,p,x)}
+@item @code{S(n, p, x)}
 @tab Nielsen's generalized polylogarithm
 @cindex @code{S()}
-@item @code{H(m_lst,x)}
+@item @code{H(m_lst, x)}
 @tab harmonic polylogarithm
 @cindex @code{H()}
-@item @code{Li(m_lst,x_lst)}
+@item @code{Li(m_lst, x_lst)}
 @tab multiple polylogarithm
 @cindex @code{Li()}
 @item @code{mZeta(m_lst)}
@@ -4851,12 +4897,11 @@ let us solve the two equations @code{a*x+b*y==3} and @code{x-y==b}:
 @example
 @{
     symbol a("a"), b("b"), x("x"), y("y");
-    lst eqns;
-    eqns.append(a*x+b*y==3).append(x-y==b);
-    lst vars;
-    vars.append(x).append(y);
+    lst eqns, vars;
+    eqns = a*x+b*y==3, x-y==b;
+    vars = x, y;
     cout << lsolve(eqns, vars) << endl;
-    // -> @{x==(3+b^2)/(b+a),y==(3-b*a)/(b+a)@}
+     // -> @{x==(3+b^2)/(b+a),y==(3-b*a)/(b+a)@}
 @end example
 
 When the linear equations @code{eqns} are underdetermined, the solution
@@ -5199,7 +5244,8 @@ And the stored expressions can be retrieved by their name:
 
 @example
     // ...
-    lst syms(x, y);
+    lst syms;
+    syms = x, y;
 
     ex ex1 = a2.unarchive_ex(syms, "foo");
     ex ex2 = a2.unarchive_ex(syms, "the second one");
index 4b6d055..fb929f2 100644 (file)
@@ -67,6 +67,64 @@ template <>
 inline void container_storage<std::vector>::reserve(std::vector<ex> & v, size_t n) { v.reserve(n); }
 
 
+/** Helper template to allow initialization of containers via an overloaded
+ *  comma operator (idea stolen from Blitz++). */
+template <typename T, typename STLT>
+class container_init {
+public:
+       container_init(STLT & s) : stlt(s) {}
+
+       container_init<T, STLT> operator,(const T & x)
+       {
+               stlt.push_back(x);
+               return container_init<T, STLT>(stlt);
+       }
+
+       // The following specializations produce much tighter code than the
+       // general case above
+
+       container_init<T, STLT> operator,(int x)
+       {
+               stlt.push_back(x);
+               return container_init<T, STLT>(stlt);
+       }
+
+       container_init<T, STLT> operator,(unsigned int x)
+       {
+               stlt.push_back(x);
+               return container_init<T, STLT>(stlt);
+       }
+
+       container_init<T, STLT> operator,(long x)
+       {
+               stlt.push_back(x);
+               return container_init<T, STLT>(stlt);
+       }
+
+       container_init<T, STLT> operator,(unsigned long x)
+       {
+               stlt.push_back(x);
+               return container_init<T, STLT>(stlt);
+       }
+
+       container_init<T, STLT> operator,(double x)
+       {
+               stlt.push_back(x);
+               return container_init<T, STLT>(stlt);
+       }
+
+       container_init<T, STLT> operator,(const symbol & x)
+       {
+               stlt.push_back(T(x));
+               return container_init<T, STLT>(stlt);
+       }
+
+private:
+       container_init();
+       STLT & stlt;
+};
+
+
 /** Wrapper template for making GiNaC classes out of STL containers. */
 template <template <class> class C>
 class container : public basic, public container_storage<C> {
@@ -290,6 +348,15 @@ public:
                this->seq.push_back(p16);
        }
 
+       // First step of initialization of container with a comma-separated
+       // sequence of expressions. Subsequent steps are handled by
+       // container_init<>::operator,().
+       container_init<ex, STLT> operator=(const ex & x)
+       {
+               this->seq.push_back(x);
+               return container_init<ex, STLT>(seq);
+       }
+
        // functions overriding virtual functions from base classes
 public:
        bool info(unsigned inf) const { return inherited::info(inf); }
index 36e7284..fc6a80b 100644 (file)
 
 namespace GiNaC {
 
+
+/** Helper template to allow initialization of matrices via an overloaded
+ *  comma operator (idea stolen from Blitz++). */
+template <typename T, typename It>
+class matrix_init {
+public:
+       matrix_init(It i) : iter(i) {}
+
+       matrix_init<T, It> operator,(const T & x)
+       {
+               *iter = x;
+               return matrix_init<T, It>(++iter);
+       }
+
+       // The following specializations produce much tighter code than the
+       // general case above
+
+       matrix_init<T, It> operator,(int x)
+       {
+               *iter = T(x);
+               return matrix_init<T, It>(++iter);
+       }
+
+       matrix_init<T, It> operator,(unsigned int x)
+       {
+               *iter = T(x);
+               return matrix_init<T, It>(++iter);
+       }
+
+       matrix_init<T, It> operator,(long x)
+       {
+               *iter = T(x);
+               return matrix_init<T, It>(++iter);
+       }
+
+       matrix_init<T, It> operator,(unsigned long x)
+       {
+               *iter = T(x);
+               return matrix_init<T, It>(++iter);
+       }
+
+       matrix_init<T, It> operator,(double x)
+       {
+               *iter = T(x);
+               return matrix_init<T, It>(++iter);
+       }
+
+       matrix_init<T, It> operator,(const symbol & x)
+       {
+               *iter = T(x);
+               return matrix_init<T, It>(++iter);
+       }
+
+private:
+       matrix_init();
+       It iter;
+};
+
+
 /** Symbolic matrices. */
 class matrix : public basic
 {
@@ -40,6 +99,14 @@ public:
        matrix(unsigned r, unsigned c);
        matrix(unsigned r, unsigned c, const exvector & m2);
        matrix(unsigned r, unsigned c, const lst & l);
+
+       // First step of initialization of matrix with a comma-separated seqeuence
+       // of expressions. Subsequent steps are handled by matrix_init<>::operator,().
+       matrix_init<ex, exvector::iterator> operator=(const ex & x)
+       {
+               m[0] = x;
+               return matrix_init<ex, exvector::iterator>(++m.begin());
+       }
        
        // functions overriding virtual functions from base classes
 public: