From aff6beb8e799e6827c40975ed2f22b51976b1cb8 Mon Sep 17 00:00:00 2001
From: Richard Kreckel <Richard.Kreckel@uni-mainz.de>
Date: Mon, 13 Mar 2000 17:15:05 +0000
Subject: [PATCH] - More drastic performance boost on matrix stuff. - Addition
 of some comments.

---
 ginac/add.cpp       |  28 +++----
 ginac/basic.cpp     |   7 +-
 ginac/expairseq.cpp |   4 +-
 ginac/matrix.cpp    | 198 +++++++++++++++++++++++++-------------------
 ginac/matrix.h      |  10 +--
 ginac/mul.cpp       |  28 +++----
 ginac/power.cpp     |  25 +++---
 7 files changed, 164 insertions(+), 136 deletions(-)

diff --git a/ginac/add.cpp b/ginac/add.cpp
index 5db2045c..08d275a2 100644
--- a/ginac/add.cpp
+++ b/ginac/add.cpp
@@ -91,7 +91,7 @@ add::add(const ex & lh, const ex & rh)
 {
     debugmsg("add constructor from ex,ex",LOGLEVEL_CONSTRUCT);
     tinfo_key = TINFO_add;
-    overall_coeff=_ex0();
+    overall_coeff = _ex0();
     construct_from_2_ex(lh,rh);
     GINAC_ASSERT(is_canonical());
 }
@@ -100,7 +100,7 @@ add::add(const exvector & v)
 {
     debugmsg("add constructor from exvector",LOGLEVEL_CONSTRUCT);
     tinfo_key = TINFO_add;
-    overall_coeff=_ex0();
+    overall_coeff = _ex0();
     construct_from_exvector(v);
     GINAC_ASSERT(is_canonical());
 }
@@ -126,7 +126,7 @@ add::add(const epvector & v)
 {
     debugmsg("add constructor from epvector",LOGLEVEL_CONSTRUCT);
     tinfo_key = TINFO_add;
-    overall_coeff=_ex0();
+    overall_coeff = _ex0();
     construct_from_epvector(v);
     GINAC_ASSERT(is_canonical());
 }
@@ -135,7 +135,7 @@ add::add(const epvector & v, const ex & oc)
 {
     debugmsg("add constructor from epvector,ex",LOGLEVEL_CONSTRUCT);
     tinfo_key = TINFO_add;
-    overall_coeff=oc;
+    overall_coeff = oc;
     construct_from_epvector(v);
     GINAC_ASSERT(is_canonical());
 }
@@ -145,7 +145,7 @@ add::add(epvector * vp, const ex & oc)
     debugmsg("add constructor from epvector *,ex",LOGLEVEL_CONSTRUCT);
     tinfo_key = TINFO_add;
     GINAC_ASSERT(vp!=0);
-    overall_coeff=oc;
+    overall_coeff = oc;
     construct_from_epvector(*vp);
     delete vp;
     GINAC_ASSERT(is_canonical());
@@ -309,9 +309,9 @@ bool add::info(unsigned inf) const
 
 int add::degree(const symbol & s) const
 {
-    int deg=INT_MIN;
+    int deg = INT_MIN;
     if (!overall_coeff.is_equal(_ex0())) {
-        deg=0;
+        deg = 0;
     }
     int cur_deg;
     for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
@@ -323,13 +323,13 @@ int add::degree(const symbol & s) const
 
 int add::ldegree(const symbol & s) const
 {
-    int deg=INT_MAX;
+    int deg = INT_MAX;
     if (!overall_coeff.is_equal(_ex0())) {
-        deg=0;
+        deg = 0;
     }
     int cur_deg;
     for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
-        cur_deg=(*cit).rest.ldegree(s);
+        cur_deg = (*cit).rest.ldegree(s);
         if (cur_deg<deg) deg=cur_deg;
     }
     return deg;
@@ -365,7 +365,7 @@ ex add::eval(int level) const
         return (new add(evaled_seqp,overall_coeff))->
                    setflag(status_flags::dynallocated);
     }
-
+    
 #ifdef DO_GINAC_ASSERT
     for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
         GINAC_ASSERT(!is_ex_exactly_of_type((*cit).rest,add));
@@ -375,13 +375,13 @@ ex add::eval(int level) const
         GINAC_ASSERT(!is_ex_exactly_of_type((*cit).rest,numeric));
     }
 #endif // def DO_GINAC_ASSERT
-
+    
     if (flags & status_flags::evaluated) {
         GINAC_ASSERT(seq.size()>0);
         GINAC_ASSERT((seq.size()>1)||!overall_coeff.is_equal(_ex0()));
         return *this;
     }
-
+    
     int seq_size=seq.size();
     if (seq_size==0) {
         // +(;c) -> c
@@ -528,7 +528,7 @@ ex add::recombine_pair_to_ex(const expair & p) const
 
 ex add::expand(unsigned options) const
 {
-    epvector * vp=expandchildren(options);
+    epvector * vp = expandchildren(options);
     if (vp==0) {
         return *this;
     }
diff --git a/ginac/basic.cpp b/ginac/basic.cpp
index cf69c42c..8eec72b1 100644
--- a/ginac/basic.cpp
+++ b/ginac/basic.cpp
@@ -229,9 +229,9 @@ ex & basic::let_op(int i)
 
 ex basic::operator[](const ex & index) const
 {
-    if (is_exactly_of_type(*index.bp,numeric)) {
+    if (is_exactly_of_type(*index.bp,numeric))
         return op(static_cast<const numeric &>(*index.bp).to_int());
-    }
+    
     throw(std::invalid_argument("non-numeric indices not supported by this type"));
 }
 
@@ -240,6 +240,8 @@ ex basic::operator[](int i) const
     return op(i);
 }
 
+/** Search ocurrences.  An object  'has' an expression if it is the expression
+ *  itself or one of the children 'has' it. */
 bool basic::has(const ex & other) const
 {
     GINAC_ASSERT(other.bp!=0);
@@ -283,6 +285,7 @@ ex basic::eval(int level) const
     return this->hold();
 }
 
+/** Evaluate object numerically. */
 ex basic::evalf(int level) const
 {
     return *this;
diff --git a/ginac/expairseq.cpp b/ginac/expairseq.cpp
index d791abf2..e1b4852c 100644
--- a/ginac/expairseq.cpp
+++ b/ginac/expairseq.cpp
@@ -1545,8 +1545,8 @@ bool expairseq::is_canonical() const
 
 epvector * expairseq::expandchildren(unsigned options) const
 {
-    epvector::const_iterator last=seq.end();
-    epvector::const_iterator cit=seq.begin();
+    epvector::const_iterator last = seq.end();
+    epvector::const_iterator cit = seq.begin();
     while (cit!=last) {
         const ex & expanded_ex=(*cit).rest.expand(options);
         if (!are_ex_trivially_equal((*cit).rest,expanded_ex)) {
diff --git a/ginac/matrix.cpp b/ginac/matrix.cpp
index 8fbe8628..0c627919 100644
--- a/ginac/matrix.cpp
+++ b/ginac/matrix.cpp
@@ -250,14 +250,12 @@ ex matrix::eval(int level) const
     debugmsg("matrix eval",LOGLEVEL_MEMBER_FUNCTION);
     
     // check if we have to do anything at all
-    if ((level==1)&&(flags & status_flags::evaluated)) {
+    if ((level==1)&&(flags & status_flags::evaluated))
         return *this;
-    }
     
     // emergency break
-    if (level == -max_recursion_level) {
+    if (level == -max_recursion_level)
         throw (std::runtime_error("matrix::eval(): recursion limit exceeded"));
-    }
     
     // eval() entry by entry
     exvector m2(row*col);
@@ -278,9 +276,8 @@ ex matrix::evalf(int level) const
     debugmsg("matrix evalf",LOGLEVEL_MEMBER_FUNCTION);
         
     // check if we have to do anything at all
-    if (level==1) {
+    if (level==1)
         return *this;
-    }
     
     // emergency break
     if (level == -max_recursion_level) {
@@ -336,9 +333,8 @@ int matrix::compare_same_type(const basic & other) const
  *  @exception logic_error (incompatible matrices) */
 matrix matrix::add(const matrix & other) const
 {
-    if (col != other.col || row != other.row) {
+    if (col != other.col || row != other.row)
         throw (std::logic_error("matrix::add(): incompatible matrices"));
-    }
     
     exvector sum(this->m);
     exvector::iterator i;
@@ -351,14 +347,14 @@ matrix matrix::add(const matrix & other) const
     return matrix(row,col,sum);
 }
 
+
 /** Difference of matrices.
  *
  *  @exception logic_error (incompatible matrices) */
 matrix matrix::sub(const matrix & other) const
 {
-    if (col != other.col || row != other.row) {
+    if (col != other.col || row != other.row)
         throw (std::logic_error("matrix::sub(): incompatible matrices"));
-    }
     
     exvector dif(this->m);
     exvector::iterator i;
@@ -371,14 +367,14 @@ matrix matrix::sub(const matrix & other) const
     return matrix(row,col,dif);
 }
 
+
 /** Product of matrices.
  *
  *  @exception logic_error (incompatible matrices) */
 matrix matrix::mul(const matrix & other) const
 {
-    if (col != other.row) {
+    if (col != other.row)
         throw (std::logic_error("matrix::mul(): incompatible matrices"));
-    }
     
     exvector prod(row*other.col);
     for (unsigned i=0; i<row; ++i) {
@@ -391,6 +387,7 @@ matrix matrix::mul(const matrix & other) const
     return matrix(row, other.col, prod);
 }
 
+
 /** operator() to access elements.
  *
  *  @param ro row of element
@@ -398,27 +395,27 @@ matrix matrix::mul(const matrix & other) const
  *  @exception range_error (index out of range) */
 const ex & matrix::operator() (unsigned ro, unsigned co) const
 {
-    if (ro<0 || ro>=row || co<0 || co>=col) {
+    if (ro<0 || ro>=row || co<0 || co>=col)
         throw (std::range_error("matrix::operator(): index out of range"));
-    }
     
     return m[ro*col+co];
 }
 
+
 /** Set individual elements manually.
  *
  *  @exception range_error (index out of range) */
 matrix & matrix::set(unsigned ro, unsigned co, ex value)
 {
-    if (ro<0 || ro>=row || co<0 || co>=col) {
+    if (ro<0 || ro>=row || co<0 || co>=col)
         throw (std::range_error("matrix::set(): index out of range"));
-    }
     
     ensure_if_modifiable();
     m[ro*col+co] = value;
     return *this;
 }
 
+
 /** Transposed of an m x n matrix, producing a new n x m matrix object that
  *  represents the transposed. */
 matrix matrix::transpose(void) const
@@ -432,80 +429,72 @@ matrix matrix::transpose(void) const
     return matrix(col,row,trans);
 }
 
-/*  Leverrier algorithm for large matrices having at least one symbolic entry.
- *  This routine is only called internally by matrix::determinant(). The
- *  algorithm is very bad for symbolic matrices since it returns expressions
- *  that are quite hard to expand. */
-/*ex matrix::determinant_symbolic_leverrier(const matrix & M)
- *{
- *    GINAC_ASSERT(M.rows()==M.cols());  // cannot happen, just in case...
- *    
- *    matrix B(M);
- *    matrix I(M.row, M.col);
- *    ex c=B.trace();
- *    for (unsigned i=1; i<M.row; ++i) {
- *        for (unsigned j=0; j<M.row; ++j)
- *            I.m[j*M.col+j] = c;
- *        B = M.mul(B.sub(I));
- *        c = B.trace()/ex(i+1);
- *    }
- *    if (M.row%2) {
- *        return c;
- *    } else {
- *        return -c;
- *    }
- *}*/
 
 /** Determinant of square matrix.  This routine doesn't actually calculate the
  *  determinant, it only implements some heuristics about which algorithm to
- *  call.  When the parameter for normalization is explicitly turned off this
- *  method does not normalize its result at the end, which might imply that
- *  the symbolic 2x2 matrix [[a/(a-b),1],[b/(a-b),1]] is not immediatly
- *  recognized to be unity.  (This is Mathematica's default behaviour, it
- *  should be used with care.)
+ *  call.  If all the elements of the matrix are elements of an integral domain
+ *  the determinant is also in that integral domain and the result is expanded
+ *  only.  If one or more elements are from a quotient field the determinant is
+ *  usually also in that quotient field and the result is normalized before it
+ *  is returned.  This implies that the determinant of the symbolic 2x2 matrix
+ *  [[a/(a-b),1],[b/(a-b),1]] is returned as unity.  (In this respect, it
+ *  behaves like MapleV and unlike Mathematica.)
  *
- *  @param     normalized may be set to false if no normalization of the
- *             result is desired (i.e. to force Mathematica behavior, Maple
- *             does normalize the result).
  *  @return    the determinant as a new expression
  *  @exception logic_error (matrix not square) */
-ex matrix::determinant(bool normalized) const
+ex matrix::determinant(void) const
 {
-    if (row != col) {
+    if (row != col)
         throw (std::logic_error("matrix::determinant(): matrix not square"));
-    }
-
-    // check, if there are non-numeric entries in the matrix:
+    GINAC_ASSERT(row*col==m.capacity());
+    
+    ex det;
+    bool numeric_flag = true;
+    bool normal_flag = false;
     for (exvector::const_iterator r=m.begin(); r!=m.end(); ++r) {
-        if (!(*r).info(info_flags::numeric)) {
-            if (normalized)
-                // return determinant_symbolic_minor().normal();
-                return determinant_symbolic_minor().normal();
-            else
-                return determinant_symbolic_perm();
-        }
+        if (!(*r).info(info_flags::numeric))
+            numeric_flag = false;
+        if ((*r).info(info_flags::rational_function) &&
+            !(*r).info(info_flags::crational_polynomial))
+            normal_flag = true;
     }
-    // if it turns out that all elements are numeric
-    return determinant_numeric();
+    
+    if (numeric_flag)
+        det = determinant_numeric();
+    else
+        if (normal_flag)
+            det = determinant_symbolic_minor().normal();
+        else
+            det = determinant_symbolic_minor();
+    
+    return det;
 }
 
-/** Trace of a matrix.
+
+/** Trace of a matrix.  The result is normalized if it is in some quotient
+ *  field and expanded only otherwise.  This implies that the trace of the
+ *  symbolic 2x2 matrix [[a/(a-b),x],[y,b/(b-a)]] is recognized to be unity.
  *
  *  @return    the sum of diagonal elements
  *  @exception logic_error (matrix not square) */
 ex matrix::trace(void) const
 {
-    if (row != col) {
+    if (row != col)
         throw (std::logic_error("matrix::trace(): matrix not square"));
-    }
+    GINAC_ASSERT(row*col==m.capacity());
     
     ex tr;
     for (unsigned r=0; r<col; ++r)
         tr += m[r*col+r];
-
-    return tr;
+    
+    if (tr.info(info_flags::rational_function) &&
+        !tr.info(info_flags::crational_polynomial))
+        return tr.normal();
+    else
+        return tr.expand();
 }
 
+
 /** Characteristic Polynomial.  The characteristic polynomial of a matrix M is
  *  defined as the determiant of (M - lambda * 1) where 1 stands for the unit
  *  matrix of the same dimension as M.  This method returns the characteristic
@@ -516,9 +505,8 @@ ex matrix::trace(void) const
  *  @see       matrix::determinant() */
 ex matrix::charpoly(const ex & lambda) const
 {
-    if (row != col) {
+    if (row != col)
         throw (std::logic_error("matrix::charpoly(): matrix not square"));
-    }
     
     matrix M(*this);
     for (unsigned r=0; r<col; ++r)
@@ -527,6 +515,7 @@ ex matrix::charpoly(const ex & lambda) const
     return (M.determinant());
 }
 
+
 /** Inverse of this matrix.
  *
  *  @return    the inverted matrix
@@ -534,9 +523,8 @@ ex matrix::charpoly(const ex & lambda) const
  *  @exception runtime_error (singular matrix) */
 matrix matrix::inverse(void) const
 {
-    if (row != col) {
+    if (row != col)
         throw (std::logic_error("matrix::inverse(): matrix not square"));
-    }
     
     matrix tmp(row,col);
     // set tmp to the unit matrix
@@ -573,6 +561,7 @@ matrix matrix::inverse(void) const
     return tmp;
 }
 
+
 // superfluous helper function
 void matrix::ffe_swap(unsigned r1, unsigned c1, unsigned r2 ,unsigned c2)
 {
@@ -801,6 +790,7 @@ matrix matrix::old_solve(const matrix & v) const
     return matrix(v.row, v.col, sol);
 }
 
+
 // protected
 
 /** Determinant of purely numeric matrix, using pivoting.
@@ -826,29 +816,56 @@ ex matrix::determinant_numeric(void) const
             }
         }
     }
+    
     return det;
 }
 
+
+/*  Leverrier algorithm for large matrices having at least one symbolic entry.
+ *  This routine is only called internally by matrix::determinant(). The
+ *  algorithm is very bad for symbolic matrices since it returns expressions
+ *  that are quite hard to expand. */
+/*ex matrix::determinant_symbolic_leverrier(const matrix & M)
+ *{
+ *    GINAC_ASSERT(M.rows()==M.cols());  // cannot happen, just in case...
+ *    
+ *    matrix B(M);
+ *    matrix I(M.row, M.col);
+ *    ex c=B.trace();
+ *    for (unsigned i=1; i<M.row; ++i) {
+ *        for (unsigned j=0; j<M.row; ++j)
+ *            I.m[j*M.col+j] = c;
+ *        B = M.mul(B.sub(I));
+ *        c = B.trace()/ex(i+1);
+ *    }
+ *    if (M.row%2) {
+ *        return c;
+ *    } else {
+ *        return -c;
+ *    }
+ *}*/
+
+
 /** Recursive determinant for small matrices having at least one symbolic
  *  entry.  The basic algorithm, known as Laplace-expansion, is enhanced by
  *  some bookkeeping to avoid calculation of the same submatrices ("minors")
  *  more than once.  According to W.M.Gentleman and S.C.Johnson this algorithm
- *  is better than elimination schemes for sparse multivariate polynomials and
- *  also for dense univariate polynomials once the dimesion becomes larger
- *  than 7.
+ *  is better than elimination schemes for matrices of sparse multivariate
+ *  polynomials and also for matrices of dense univariate polynomials if the
+ *  matrix' dimesion is larger than 7.
  *
  *  @see matrix::determinant() */
 ex matrix::determinant_symbolic_minor(void) const
 {
-    // for small matrices the algorithm does not make sense:
+    // for small matrices the algorithm does not make any sense:
     if (this->row==1)
         return m[0];
     if (this->row==2)
-        return (m[0]*m[3]-m[2]*m[1]);
+        return (m[0]*m[3]-m[2]*m[1]).expand();
     if (this->row==3)
-        return ((m[4]*m[8]-m[5]*m[7])*m[0]-
-                (m[1]*m[8]-m[2]*m[7])*m[3]+
-                (m[1]*m[5]-m[4]*m[2])*m[6]);
+        return (m[0]*m[4]*m[8]-m[0]*m[5]*m[7]-
+                m[1]*m[3]*m[8]+m[2]*m[3]*m[7]+
+                m[1]*m[5]*m[6]-m[2]*m[4]*m[6]).expand();
     
     // This algorithm can best be understood by looking at a naive
     // implementation of Laplace-expansion, like this one:
@@ -873,7 +890,7 @@ ex matrix::determinant_symbolic_minor(void) const
     //     else
     //         det += m[r1*col] * minorM.determinant_symbolic_minor();
     // }
-    // return det;
+    // return det.expand();
     // What happens is that while proceeding down many of the minors are
     // computed more than once.  In particular, there are binomial(n,k)
     // kxk minors and each one is computed factorial(n-k) times.  Therefore
@@ -882,15 +899,18 @@ ex matrix::determinant_symbolic_minor(void) const
     // calculated in step c-1.  We therefore only have to store at most 
     // 2*binomial(n,n/2) minors.
     
+    // Unique flipper counter for partitioning into minors
+    vector<unsigned> Pkey;
+    Pkey.reserve(this->col);
+    // key for minor determinant (a subpartition of Pkey)
+    vector<unsigned> Mkey;
+    Mkey.reserve(this->col-1);
     // we store our subminors in maps, keys being the rows they arise from
     typedef map<vector<unsigned>,class ex> Rmap;
     typedef map<vector<unsigned>,class ex>::value_type Rmap_value;
-    Rmap A, B;
+    Rmap A;
+    Rmap B;
     ex det;
-    vector<unsigned> Pkey;    // Unique flipper counter for partitioning into minors
-    Pkey.reserve(this->col);
-    vector<unsigned> Mkey;    // key for minor determinant (a subpartition of Pkey)
-    Mkey.reserve(this->col-1);
     // initialize A with last column:
     for (unsigned r=0; r<this->col; ++r) {
         Pkey.erase(Pkey.begin(),Pkey.end());
@@ -922,7 +942,9 @@ ex matrix::determinant_symbolic_minor(void) const
                 else
                     det += m[Pkey[r]*this->col+c]*A[Mkey];
             }
-            // Store the new determinant at its place in B:
+            // prevent build-up of deep nesting of expressions saves some time:
+            det = det.expand();
+            // store the new determinant at its place in B:
             B.insert(Rmap_value(Pkey,det));
             // increment our strange flipper counter
             for (fc=this->col-c; fc>0; --fc) {
@@ -934,7 +956,7 @@ ex matrix::determinant_symbolic_minor(void) const
                 for (unsigned j=fc; j<this->col-c; ++j)
                     Pkey[j] = Pkey[j-1]+1;
         } while(fc);
-        // change the role of A and B:
+        // next column, so change the role of A and B:
         A = B;
         B.clear();
     }
@@ -942,6 +964,7 @@ ex matrix::determinant_symbolic_minor(void) const
     return det;
 }
 
+
 /** Determinant built by application of the full permutation group.  This
  *  routine is only called internally by matrix::determinant(). */
 ex matrix::determinant_symbolic_perm(void) const
@@ -965,6 +988,7 @@ ex matrix::determinant_symbolic_perm(void) const
     return det;
 }
 
+
 /** Perform the steps of an ordinary Gaussian elimination to bring the matrix
  *  into an upper echelon form.
  *
@@ -987,9 +1011,11 @@ int matrix::gauss_elimination(void)
                 this->m[r2*col+c] = _ex0();
         }
     }
+    
     return sign;
 }
 
+
 /** Partial pivoting method.
  *  Usual pivoting (symbolic==false) returns the index to the element with the
  *  largest absolute value in column ro and swaps the current row with the one
diff --git a/ginac/matrix.h b/ginac/matrix.h
index 3a6f0e53..8a70ada8 100644
--- a/ginac/matrix.h
+++ b/ginac/matrix.h
@@ -77,9 +77,9 @@ protected:
     
     // non-virtual functions in this class
 public:
-    unsigned rows() const            //! get number of rows.
+    unsigned rows(void) const        //! get number of rows.
         { return row; }
-    unsigned cols() const            //! get number of columns.
+    unsigned cols(void) const        //! get number of columns.
         { return col; }
     matrix add(const matrix & other) const;
     matrix sub(const matrix & other) const;
@@ -87,7 +87,7 @@ public:
     const ex & operator() (unsigned ro, unsigned co) const;
     matrix & set(unsigned ro, unsigned co, ex value);
     matrix transpose(void) const;
-    ex determinant(bool normalized=true) const;
+    ex determinant(void) const;
     ex trace(void) const;
     ex charpoly(const ex & lambda) const;
     matrix inverse(void) const;
@@ -144,8 +144,8 @@ inline unsigned cols(const matrix & m)
 inline matrix transpose(const matrix & m)
 { return m.transpose(); }
 
-inline ex determinant(const matrix & m, bool normalized=true)
-{ return m.determinant(normalized); }
+inline ex determinant(const matrix & m)
+{ return m.determinant(); }
 
 inline ex trace(const matrix & m)
 { return m.trace(); }
diff --git a/ginac/mul.cpp b/ginac/mul.cpp
index 976db46e..efabd3be 100644
--- a/ginac/mul.cpp
+++ b/ginac/mul.cpp
@@ -656,25 +656,25 @@ ex mul::expand(unsigned options) const
     intvector positions_of_adds;
     intvector number_of_add_operands;
 
-    epvector * expanded_seqp=expandchildren(options);
+    epvector * expanded_seqp = expandchildren(options);
 
     const epvector & expanded_seq = expanded_seqp==0 ? seq : *expanded_seqp;
 
     positions_of_adds.resize(expanded_seq.size());
     number_of_add_operands.resize(expanded_seq.size());
 
-    int number_of_adds=0;
-    int number_of_expanded_terms=1;
+    int number_of_adds = 0;
+    int number_of_expanded_terms = 1;
 
-    unsigned current_position=0;
-    epvector::const_iterator last=expanded_seq.end();
+    unsigned current_position = 0;
+    epvector::const_iterator last = expanded_seq.end();
     for (epvector::const_iterator cit=expanded_seq.begin(); cit!=last; ++cit) {
         if (is_ex_exactly_of_type((*cit).rest,add)&&
             (ex_to_numeric((*cit).coeff).is_equal(_num1()))) {
-            positions_of_adds[number_of_adds]=current_position;
-            const add & expanded_addref=ex_to_add((*cit).rest);
-            unsigned addref_nops=expanded_addref.nops();
-            number_of_add_operands[number_of_adds]=addref_nops;
+            positions_of_adds[number_of_adds] = current_position;
+            const add & expanded_addref = ex_to_add((*cit).rest);
+            unsigned addref_nops = expanded_addref.nops();
+            number_of_add_operands[number_of_adds] = addref_nops;
             number_of_expanded_terms *= addref_nops;
             number_of_adds++;
         }
@@ -759,11 +759,11 @@ ex mul::expand(unsigned options) const
 
 epvector * mul::expandchildren(unsigned options) const
 {
-    epvector::const_iterator last=seq.end();
-    epvector::const_iterator cit=seq.begin();
+    epvector::const_iterator last = seq.end();
+    epvector::const_iterator cit = seq.begin();
     while (cit!=last) {
-        const ex & factor=recombine_pair_to_ex(*cit);
-        const ex & expanded_factor=factor.expand(options);
+        const ex & factor = recombine_pair_to_ex(*cit);
+        const ex & expanded_factor = factor.expand(options);
         if (!are_ex_trivially_equal(factor,expanded_factor)) {
 
             // something changed, copy seq, eval and return it
@@ -771,7 +771,7 @@ epvector * mul::expandchildren(unsigned options) const
             s->reserve(seq.size());
 
             // copy parts of seq which are known not to have changed
-            epvector::const_iterator cit2=seq.begin();
+            epvector::const_iterator cit2 = seq.begin();
             while (cit2!=cit) {
                 s->push_back(*cit2);
                 ++cit2;
diff --git a/ginac/power.cpp b/ginac/power.cpp
index 8cb62238..c8f65dc4 100644
--- a/ginac/power.cpp
+++ b/ginac/power.cpp
@@ -554,8 +554,8 @@ unsigned power::return_type_tinfo(void) const
 
 ex power::expand(unsigned options) const
 {
-    ex expanded_basis=basis.expand(options);
-
+    ex expanded_basis = basis.expand(options);
+    
     if (!is_ex_exactly_of_type(exponent,numeric)||
         !ex_to_numeric(exponent).is_integer()) {
         if (are_ex_trivially_equal(basis,expanded_basis)) {
@@ -565,19 +565,19 @@ ex power::expand(unsigned options) const
                     setflag(status_flags::dynallocated);
         }
     }
-
+    
     // integer numeric exponent
     const numeric & num_exponent=ex_to_numeric(exponent);
     int int_exponent = num_exponent.to_int();
-
+    
     if (int_exponent > 0 && is_ex_exactly_of_type(expanded_basis,add)) {
         return expand_add(ex_to_add(expanded_basis), int_exponent);
     }
-
+    
     if (is_ex_exactly_of_type(expanded_basis,mul)) {
         return expand_mul(ex_to_mul(expanded_basis), num_exponent);
     }
-
+    
     // cannot expand further
     if (are_ex_trivially_equal(basis,expanded_basis)) {
         return this->hold();
@@ -600,12 +600,11 @@ ex power::expand(unsigned options) const
 ex power::expand_add(const add & a, int n) const
 {
     // expand a^n where a is an add and n is an integer
-
-    if (n==2) {
+    
+    if (n==2)
         return expand_add_2(a);
-    }
     
-    int m=a.nops();
+    int m = a.nops();
     exvector sum;
     sum.reserve((n+1)*(m-1));
     intvector k(m-1);
@@ -618,7 +617,7 @@ ex power::expand_add(const add & a, int n) const
         k_cum[l]=0;
         upper_limit[l]=n;
     }
-
+    
     while (1) {
         exvector term;
         term.reserve(m+1);
@@ -634,7 +633,7 @@ ex power::expand_add(const add & a, int n) const
                 term.push_back(power(b,k[l]));
             }
         }
-
+        
         const ex & b=a.op(l);
         GINAC_ASSERT(!is_ex_exactly_of_type(b,add));
         GINAC_ASSERT(!is_ex_exactly_of_type(b,power)||
@@ -645,7 +644,7 @@ ex power::expand_add(const add & a, int n) const
         } else {
             term.push_back(power(b,n-k_cum[m-2]));
         }
-
+        
         numeric f=binomial(numeric(n),numeric(k[0]));
         for (l=1; l<m-1; l++) {
             f=f*binomial(numeric(n-k_cum[l-1]),numeric(k[l]));
-- 
2.44.0