From: Richard Kreckel <Richard.Kreckel@uni-mainz.de>
Date: Sat, 22 Jul 2000 04:47:04 +0000 (+0000)
Subject: - Fixed the bug that broke xloop's po_redux.
X-Git-Tag: release_0-6-3~3
X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=commitdiff_plain;h=fe9dbfb9947b24149b3ce7dd9285f27ab286cbd7

- Fixed the bug that broke xloop's po_redux.
- Interface Change: No argument 'bool branchcut' for series methods,
  instead 'unsigned options'.  Sorry.
: ----------------------------------------------------------------------
---

diff --git a/ginac/add.h b/ginac/add.h
index 7db48695..78464d3e 100644
--- a/ginac/add.h
+++ b/ginac/add.h
@@ -70,7 +70,7 @@ public:
     int ldegree(const symbol & s) const;
     ex coeff(const symbol & s, int n=1) const;
     ex eval(int level=0) const;
-    ex series(const relational & r, int order, bool branchcut = true) const;
+    ex series(const relational & r, int order, unsigned options = 0) const;
     ex normal(lst &sym_lst, lst &repl_lst, int level=0) const;
     numeric integer_content(void) const;
     ex smod(const numeric &xi) const;
diff --git a/ginac/basic.h b/ginac/basic.h
index 2263081c..43fc4958 100644
--- a/ginac/basic.h
+++ b/ginac/basic.h
@@ -138,7 +138,7 @@ public: // only const functions please (may break reference counting)
     virtual ex collect(const symbol & s) const;
     virtual ex eval(int level = 0) const;
     virtual ex evalf(int level = 0) const;
-    virtual ex series(const relational & r, int order, bool branchcut = true) const;
+    virtual ex series(const relational & r, int order, unsigned options = 0) const;
     virtual ex subs(const lst & ls, const lst & lr) const;
     virtual ex normal(lst &sym_lst, lst &repl_lst, int level = 0) const;
     virtual ex to_rational(lst &repl_lst) const;
diff --git a/ginac/ex.h b/ginac/ex.h
index 05f31573..216a4646 100644
--- a/ginac/ex.h
+++ b/ginac/ex.h
@@ -252,7 +252,7 @@ public:
     ex eval(int level = 0) const;
     ex evalf(int level = 0) const;
     ex diff(const symbol & s, unsigned nth = 1) const;
-    ex series(const ex & r, int order, bool branchcut = true) const;
+    ex series(const ex & r, int order, unsigned options = 0) const;
     ex subs(const lst & ls, const lst & lr) const;
     ex subs(const ex & e) const;
     exvector get_indices(void) const;
@@ -395,8 +395,8 @@ inline ex evalf(const ex & thisex, int level = 0)
 inline ex diff(const ex & thisex, const symbol & s, unsigned nth = 1)
 { return thisex.diff(s, nth); }
 
-inline ex series(const ex & thisex, const ex & r, int order, bool branchcut = true)
-{ return thisex.series(r, order, branchcut); }
+inline ex series(const ex & thisex, const ex & r, int order, unsigned options = 0)
+{ return thisex.series(r, order, options); }
 
 inline ex subs(const ex & thisex, const ex & e)
 { return thisex.subs(e); }
diff --git a/ginac/flags.h b/ginac/flags.h
index 8b9caed4..77e3ec01 100644
--- a/ginac/flags.h
+++ b/ginac/flags.h
@@ -33,6 +33,19 @@ public:
          };
 };
 
+class series_options {
+public:
+    enum { suppress_branchcut        = 0x0001
+         };
+};
+
+class determinant_options {
+public:
+    enum { laplace,
+           bareiss
+         };
+};
+
 class status_flags {
 public:
     enum { dynallocated              = 0x0001,
@@ -124,7 +137,7 @@ public:
     enum { delete_never, // let table grow undefinitely, not recommmended, but currently default
            delete_lru,   // least recently used
            delete_lfu,   // least frequently used
-           delete_cyclic // first one in list (oldest)
+           delete_cyclic // first (oldest) one in list
     };
 };
 
diff --git a/ginac/function.pl b/ginac/function.pl
index a69f58fe..f67c4b96 100755
--- a/ginac/function.pl
+++ b/ginac/function.pl
@@ -95,7 +95,7 @@ $typedef_derivative_funcp=generate(
 'const ex &','');
 
 $typedef_series_funcp=generate(
-'typedef ex (* series_funcp_${N})(${SEQ1}, const relational &, int, bool);'."\n",
+'typedef ex (* series_funcp_${N})(${SEQ1}, const relational &, int, unsigned);'."\n",
 'const ex &','');
 
 $eval_func_interface=generate('    function_options & eval_func(eval_funcp_${N} e);'."\n",'','');
@@ -145,9 +145,9 @@ $series_switch_statement=generate(
     <<'END_OF_SERIES_SWITCH_STATEMENT','seq[${N}-1]','');
     case ${N}:
         try {
-            res = ((series_funcp_${N})(registered_functions()[serial].series_f))(${SEQ1},r,order,branchcut);
+            res = ((series_funcp_${N})(registered_functions()[serial].series_f))(${SEQ1},r,order,options);
         } catch (do_taylor) {
-            res = basic::series(r, order, branchcut);
+            res = basic::series(r, order, options);
         }
         return res;
         break;
@@ -419,7 +419,7 @@ public:
     ex expand(unsigned options=0) const;
     ex eval(int level=0) const;
     ex evalf(int level=0) const;
-    ex series(const relational & r, int order, bool branchcut = true) const;
+    ex series(const relational & r, int order, unsigned options = 0) const;
     ex thisexprseq(const exvector & v) const;
     ex thisexprseq(exvector * vp) const;
 protected:
@@ -889,7 +889,7 @@ ex function::thisexprseq(exvector * vp) const
 
 /** Implementation of ex::series for functions.
  *  \@see ex::series */
-ex function::series(const relational & r, int order, bool branchcut = true) const
+ex function::series(const relational & r, int order, unsigned options = 0) const
 {
     GINAC_ASSERT(serial<registered_functions().size());
 
diff --git a/ginac/inifcns.cpp b/ginac/inifcns.cpp
index 14e53b5e..e0d1e91e 100644
--- a/ginac/inifcns.cpp
+++ b/ginac/inifcns.cpp
@@ -110,7 +110,7 @@ static ex csgn_eval(const ex & x)
 static ex csgn_series(const ex & arg,
                       const relational & rel,
                       int order,
-                      bool branchcut)
+                      unsigned options)
 {
     const ex arg_pt = arg.subs(rel);
     if (arg_pt.info(info_flags::numeric)) {
@@ -179,7 +179,7 @@ static ex Li2_deriv(const ex & x, unsigned deriv_param)
     return -log(1-x)/x;
 }
 
-static ex Li2_series(const ex &x, const relational &rel, int order, bool branchcut)
+static ex Li2_series(const ex &x, const relational &rel, int order, unsigned options)
 {
     const ex x_pt = x.subs(rel);
     if (x_pt.info(info_flags::numeric)) {
@@ -237,7 +237,8 @@ static ex Li2_series(const ex &x, const relational &rel, int order, bool branchc
             return ser.series(rel, order);
         }
         // third special case: x real, >=1 (branch cut)
-        if (ex_to_numeric(x_pt).is_real() && ex_to_numeric(x_pt)>1) {
+        if (!(options & series_options::suppress_branchcut) &&
+            ex_to_numeric(x_pt).is_real() && ex_to_numeric(x_pt)>1) {
             // method:
             // This is the branch cut: assemble the primitive series manually
             // and then add the corresponding complex step function.
@@ -250,7 +251,7 @@ static ex Li2_series(const ex &x, const relational &rel, int order, bool branchc
             // compute the intermediate terms:
             ex replarg = series(Li2(x), *s==foo, order);
             for (unsigned i=1; i<replarg.nops()-1; ++i)
-                seq.push_back(expair((replarg.op(i)/power(*s-foo,i)).series(foo==point,1,branchcut).op(0).subs(foo==*s),i));
+                seq.push_back(expair((replarg.op(i)/power(*s-foo,i)).series(foo==point,1,options).op(0).subs(foo==*s),i));
             // append an order term:
             seq.push_back(expair(Order(_ex1()), replarg.nops()-1));
             return pseries(rel, seq);
@@ -341,7 +342,7 @@ static ex Order_eval(const ex & x)
 	return Order(x).hold();
 }
 
-static ex Order_series(const ex & x, const relational & r, int order, bool branchcut)
+static ex Order_series(const ex & x, const relational & r, int order, unsigned options)
 {
 	// Just wrap the function into a pseries object
 	epvector new_seq;
diff --git a/ginac/inifcns_gamma.cpp b/ginac/inifcns_gamma.cpp
index 26b5bd78..bf60c1e5 100644
--- a/ginac/inifcns_gamma.cpp
+++ b/ginac/inifcns_gamma.cpp
@@ -88,7 +88,7 @@ static ex lgamma_deriv(const ex & x, unsigned deriv_param)
 static ex lgamma_series(const ex & arg,
                         const relational & rel,
                         int order,
-                        bool branchcut)
+                        unsigned options)
 {
     // method:
     // Taylor series where there is no pole falls back to psi function
@@ -184,7 +184,7 @@ static ex tgamma_deriv(const ex & x, unsigned deriv_param)
 static ex tgamma_series(const ex & arg,
                         const relational & rel,
                         int order,
-                        bool branchcut)
+                        unsigned options)
 {
     // method:
     // Taylor series where there is no pole falls back to psi function
@@ -283,7 +283,7 @@ static ex beta_series(const ex & arg1,
                       const ex & arg2,
                       const relational & rel,
                       int order,
-                      bool branchcut)
+                      unsigned options)
 {
     // method:
     // Taylor series where there is no pole of one of the tgamma functions
@@ -392,7 +392,7 @@ static ex psi1_deriv(const ex & x, unsigned deriv_param)
 static ex psi1_series(const ex & arg,
                       const relational & rel,
                       int order,
-                      bool branchcut)
+                      unsigned options)
 {
     // method:
     // Taylor series where there is no pole falls back to polygamma function
@@ -516,7 +516,7 @@ static ex psi2_series(const ex & n,
                       const ex & arg,
                       const relational & rel,
                       int order,
-                      bool branchcut)
+                      unsigned options)
 {
     // method:
     // Taylor series where there is no pole falls back to polygamma function
diff --git a/ginac/inifcns_trans.cpp b/ginac/inifcns_trans.cpp
index 75c92b76..8341651c 100644
--- a/ginac/inifcns_trans.cpp
+++ b/ginac/inifcns_trans.cpp
@@ -147,7 +147,7 @@ static ex log_deriv(const ex & x, unsigned deriv_param)
 static ex log_series(const ex &arg,
                      const relational &rel,
                      int order,
-                     bool branchcut)
+                     unsigned options)
 {
     GINAC_ASSERT(is_ex_exactly_of_type(rel.lhs(),symbol));
     ex arg_pt;
@@ -171,7 +171,7 @@ static ex log_series(const ex &arg,
         // Return a plain n*log(x) for the x^n part and series expand the
         // other part.  Add them together and reexpand again in order to have
         // one unnested pseries object.  All this also works for negative n.
-        const pseries argser = ex_to_pseries(arg.series(rel, order, branchcut));
+        const pseries argser = ex_to_pseries(arg.series(rel, order, options));
         const symbol *s = static_cast<symbol *>(rel.lhs().bp);
         const ex point = rel.rhs();
         const int n = argser.ldegree(*s);
@@ -179,12 +179,13 @@ static ex log_series(const ex &arg,
         seq.push_back(expair(n*log(*s-point), _ex0()));
         if (!argser.is_terminating() || argser.nops()!=1) {
             // in this case n more terms are needed
-            ex newarg = ex_to_pseries(arg.series(rel, order+n, branchcut)).shift_exponents(-n).convert_to_poly(true);
-            return pseries(rel, seq).add_series(ex_to_pseries(log(newarg).series(rel, order, branchcut)));
+            ex newarg = ex_to_pseries(arg.series(rel, order+n, options)).shift_exponents(-n).convert_to_poly(true);
+            return pseries(rel, seq).add_series(ex_to_pseries(log(newarg).series(rel, order, options)));
         } else  // it was a monomial
             return pseries(rel, seq);
     }
-    if (branchcut && arg_pt.info(info_flags::negative)) {
+    if (!(options & series_options::suppress_branchcut) &&
+         arg_pt.info(info_flags::negative)) {
         // method:
         // This is the branch cut: assemble the primitive series manually and
         // then add the corresponding complex step function.
@@ -444,7 +445,7 @@ static ex tan_deriv(const ex & x, unsigned deriv_param)
 static ex tan_series(const ex &x,
                      const relational &rel,
                      int order,
-                     bool branchcut)
+                     unsigned options)
 {
     // method:
     // Taylor series where there is no pole falls back to tan_deriv.
@@ -804,7 +805,7 @@ static ex tanh_deriv(const ex & x, unsigned deriv_param)
 static ex tanh_series(const ex &x,
                       const relational &rel,
                       int order,
-                      bool branchcut)
+                      unsigned options)
 {
     // method:
     // Taylor series where there is no pole falls back to tanh_deriv.
diff --git a/ginac/matrix.cpp b/ginac/matrix.cpp
index fc9d7e26..feae32fe 100644
--- a/ginac/matrix.cpp
+++ b/ginac/matrix.cpp
@@ -382,8 +382,8 @@ matrix matrix::mul(const matrix & other) const
     
     exvector prod(row*other.col);
     
-    for (unsigned r1=0; r1<row; ++r1) {
-        for (unsigned c=0; c<col; ++c) {
+    for (unsigned r1=0; r1<rows(); ++r1) {
+        for (unsigned c=0; c<cols(); ++c) {
             if (m[r1*col+c].is_zero())
                 continue;
             for (unsigned r2=0; r2<other.col; ++r2)
@@ -397,13 +397,13 @@ matrix matrix::mul(const matrix & other) const
 /** operator() to access elements.
  *
  *  @param ro row of element
- *  @param co column of element 
+ *  @param co column of element
  *  @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)
         throw (std::range_error("matrix::operator(): index out of range"));
-    
+
     return m[ro*col+co];
 }
 
@@ -626,9 +626,8 @@ matrix matrix::inverse(void) const
             throw (std::runtime_error("matrix::inverse(): singular matrix"));
         }
         if (indx != 0) {  // swap rows r and indx of matrix tmp
-            for (unsigned i=0; i<col; ++i) {
+            for (unsigned i=0; i<col; ++i)
                 tmp.m[r1*col+i].swap(tmp.m[indx*col+i]);
-            }
         }
         ex a1 = cpy.m[r1*col+r1];
         for (unsigned c=0; c<col; ++c) {
@@ -648,6 +647,16 @@ matrix matrix::inverse(void) const
     return tmp;
 }
 
+// superfluous helper function, to be removed:
+void matrix::swap(unsigned r1, unsigned c1, unsigned r2 ,unsigned c2)
+{
+    ensure_if_modifiable();
+    
+    ex tmp = (*this)(r1,c1);
+    set(r1,c1,(*this)(r2,c2));
+    set(r2,c2,tmp);
+}
+
 
 /** Solve a set of equations for an m x n matrix by fraction-free Gaussian
  *  elimination.  Based on algorithm 9.1 from 'Algorithms for Computer Algebra'
@@ -678,30 +687,30 @@ matrix matrix::fraction_free_elim(const matrix & vars,
     for (unsigned k=0; (k<n)&&(r<m); ++k) {
         // find a nonzero pivot
         unsigned p;
-        for (p=r; (p<m)&&(a.m[p*a.cols()+k].is_zero()); ++p) {}
+        for (p=r; (p<m)&&(a(p,k).is_zero()); ++p) {}
         // pivot is in row p
         if (p<m) {
             if (p!=r) {
                 // swap rows p and r
                 for (unsigned j=k; j<n; ++j)
-                    a.m[p*a.cols()+j].swap(a.m[r*a.cols()+j]);
-                a.m[p*a.cols()].swap(a.m[r*a.cols()]);
+                    a.swap(p,j,r,j);
+                b.swap(p,0,r,0);
                 // keep track of sign changes due to row exchange
-                sign = -sign;
+                sign *= -1;
             }
             for (unsigned i=r+1; i<m; ++i) {
                 for (unsigned j=k+1; j<n; ++j) {
-                    a.set(i,j,(a.m[r*a.cols()+k]*a.m[i*a.cols()+j]
-                              -a.m[r*a.cols()+j]*a.m[i*a.cols()+k])/divisor);
-                    a.set(i,j,a.m[i*a.cols()+j].normal());
+                    a.set(i,j,(a(r,k)*a(i,j)
+                              -a(r,j)*a(i,k))/divisor);
+                    a.set(i,j,a(i,j).normal());
                 }
-                b.set(i,0,(a.m[r*a.cols()+k]*b.m[i*b.cols()]
-                          -b.m[r*b.cols()]*a.m[i*a.cols()+k])/divisor);
-                b.set(i,0,b.m[i*b.cols()].normal());
+                b.set(i,0,(a(r,k)*b(i,0)
+                          -b(r,0)*a(i,k))/divisor);
+                b.set(i,0,b(i,0).normal());
                 a.set(i,k,0);
             }
-            divisor = a.m[r*a.cols()+k];
-            r++;
+            divisor = a(r,k);
+            ++r;
         }
     }
     
@@ -711,7 +720,7 @@ matrix matrix::fraction_free_elim(const matrix & vars,
     for (unsigned r=0; r<m; ++r) {
         int zero_in_this_row=0;
         for (unsigned c=0; c<n; ++c) {
-            if (a.m[r*a.cols()+c].is_zero())
+            if (a(r,c).is_zero())
                zero_in_this_row++;
             else
                 break;
@@ -726,65 +735,65 @@ matrix matrix::fraction_free_elim(const matrix & vars,
     unsigned last_assigned_sol = n+1;
     for (int r=m-1; r>=0; --r) {
         unsigned first_non_zero = 1;
-        while ((first_non_zero<=n)&&(a.m[r*a.cols()+(first_non_zero-1)].is_zero()))
+        while ((first_non_zero<=n)&&(a(r,first_non_zero-1).is_zero()))
             first_non_zero++;
         if (first_non_zero>n) {
             // row consists only of zeroes, corresponding rhs must be 0 as well
-            if (!b.m[r*b.cols()].is_zero()) {
+            if (!b(r,0).is_zero()) {
                 throw (std::runtime_error("matrix::fraction_free_elim(): singular matrix"));
             }
         } else {
             // assign solutions for vars between first_non_zero+1 and
             // last_assigned_sol-1: free parameters
             for (unsigned c=first_non_zero; c<last_assigned_sol-1; ++c)
-                sol.set(c,0,vars.m[c*vars.cols()]);
-            ex e = b.m[r*b.cols()];
+                sol.set(c,0,vars(c,0));
+            ex e = b(r,0);
             for (unsigned c=first_non_zero; c<n; ++c)
-                e -= a.m[r*a.cols()+c]*sol.m[c*sol.cols()];
+                e -= a(r,c)*sol(c,0);
             sol.set(first_non_zero-1,0,
-                    (e/a.m[r*a.cols()+(first_non_zero-1)]).normal());
+                    (e/a(r,first_non_zero-1)).normal());
             last_assigned_sol = first_non_zero;
         }
     }
     // assign solutions for vars between 1 and
     // last_assigned_sol-1: free parameters
     for (unsigned c=0; c<last_assigned_sol-1; ++c)
-        sol.set(c,0,vars.m[c*vars.cols()]);
+        sol.set(c,0,vars(c,0));
     
 #ifdef DO_GINAC_ASSERT
     // test solution with echelon matrix
     for (unsigned r=0; r<m; ++r) {
         ex e = 0;
         for (unsigned c=0; c<n; ++c)
-            e += a.m[r*a.cols()+c]*sol.m[c*sol.cols()];
-        if (!(e-b.m[r*b.cols()]).normal().is_zero()) {
+            e += a(r,c)*sol(c,0);
+        if (!(e-b(r,0)).normal().is_zero()) {
             cout << "e=" << e;
-            cout << "b(" << r <<",0)=" << b.m[r*b.cols()] << endl;
-            cout << "diff=" << (e-b.m[r*b.cols()]).normal() << endl;
+            cout << "b(" << r <<",0)=" << b(r,0) << endl;
+            cout << "diff=" << (e-b(r,0)).normal() << endl;
         }
-        GINAC_ASSERT((e-b.m[r*b.cols()]).normal().is_zero());
+        GINAC_ASSERT((e-b(r,0)).normal().is_zero());
     }
     
     // test solution with original matrix
     for (unsigned r=0; r<m; ++r) {
         ex e = 0;
         for (unsigned c=0; c<n; ++c)
-            e += this->m[r*cols()+c]*sol.m[c*sol.cols()];
+            e += this->m[r*cols()+c]*sol(c,0);
         try {
-            if (!(e-rhs.m[r*rhs.cols()]).normal().is_zero()) {
+            if (!(e-rhs(r,0)).normal().is_zero()) {
                 cout << "e==" << e << endl;
                 e.printtree(cout);
                 ex en = e.normal();
                 cout << "e.normal()=" << en << endl;
                 en.printtree(cout);
-                cout << "rhs(" << r <<",0)=" << rhs.m[r*rhs.cols()] << endl;
-                cout << "diff=" << (e-rhs.m[r*rhs.cols()]).normal() << endl;
+                cout << "rhs(" << r <<",0)=" << rhs(r,0) << endl;
+                cout << "diff=" << (e-rhs(r,0)).normal() << endl;
             }
         } catch (...) {
-            ex xxx = e - rhs.m[r*rhs.cols()];
+            ex xxx = e - rhs(r,0);
             cerr << "xxx=" << xxx << endl << endl;
         }
-        GINAC_ASSERT((e-rhs.m[r*rhs.cols()]).normal().is_zero());
+        GINAC_ASSERT((e-rhs(r,0)).normal().is_zero());
     }
 #endif // def DO_GINAC_ASSERT
     
diff --git a/ginac/matrix.h b/ginac/matrix.h
index f4665278..cef3dc36 100644
--- a/ginac/matrix.h
+++ b/ginac/matrix.h
@@ -27,6 +27,11 @@
 #include "basic.h"
 #include "ex.h"
 
+namespace std {
+    // forward declaration, so <stdexcept> need not be included:
+    class range_error;
+}
+
 #ifndef NO_NAMESPACE_GINAC
 namespace GiNaC {
 #endif // ndef NO_NAMESPACE_GINAC
@@ -100,6 +105,7 @@ protected:
     int division_free_elimination(void);
     int fraction_free_elimination(bool det = false);
     int pivot(unsigned ro, bool symbolic=true);
+    void swap(unsigned r1, unsigned c1, unsigned r2 ,unsigned c2);
     
 // member variables
 protected:
@@ -109,10 +115,12 @@ protected:
     static unsigned precedence;
 };
 
+
 // global constants
 extern const matrix some_matrix;
 extern const type_info & typeid_matrix;
 
+
 // wrapper functions around member functions
 
 inline unsigned nops(const matrix & m)
diff --git a/ginac/mul.h b/ginac/mul.h
index b8409f4d..bfb92366 100644
--- a/ginac/mul.h
+++ b/ginac/mul.h
@@ -71,7 +71,7 @@ public:
     ex coeff(const symbol & s, int n = 1) const;
     ex eval(int level=0) const;
     ex evalf(int level=0) const;
-    ex series(const relational & s, int order, bool branchcut = true) const;
+    ex series(const relational & s, int order, unsigned options = 0) const;
     ex normal(lst &sym_lst, lst &repl_lst, int level = 0) const;
     numeric integer_content(void) const;
     ex smod(const numeric &xi) const;
diff --git a/ginac/power.h b/ginac/power.h
index 3d8be419..afc59f33 100644
--- a/ginac/power.h
+++ b/ginac/power.h
@@ -73,7 +73,7 @@ public:
     ex coeff(const symbol & s, int n = 1) const;
     ex eval(int level=0) const;
     ex evalf(int level=0) const;
-    ex series(const relational & s, int order, bool branchcut = true) const;
+    ex series(const relational & s, int order, unsigned options = 0) const;
     ex subs(const lst & ls, const lst & lr) const;
     ex normal(lst &sym_lst, lst &repl_lst, int level = 0) const;
     ex to_rational(lst &repl_lst) const;
diff --git a/ginac/pseries.cpp b/ginac/pseries.cpp
index 38499480..f55c0fed 100644
--- a/ginac/pseries.cpp
+++ b/ginac/pseries.cpp
@@ -484,7 +484,7 @@ bool pseries::is_terminating(void) const
 
 /** Default implementation of ex::series(). This performs Taylor expansion.
  *  @see ex::series */
-ex basic::series(const relational & r, int order, bool branchcut) const
+ex basic::series(const relational & r, int order, unsigned options) const
 {
     epvector seq;
     numeric fac(1);
@@ -518,7 +518,7 @@ ex basic::series(const relational & r, int order, bool branchcut) const
 
 /** Implementation of ex::series() for symbols.
  *  @see ex::series */
-ex symbol::series(const relational & r, int order, bool branchcut) const
+ex symbol::series(const relational & r, int order, unsigned options) const
 {
     epvector seq;
     const ex point = r.rhs();
@@ -615,12 +615,12 @@ ex pseries::add_series(const pseries &other) const
 /** Implementation of ex::series() for sums. This performs series addition when
  *  adding pseries objects.
  *  @see ex::series */
-ex add::series(const relational & r, int order, bool branchcut) const
+ex add::series(const relational & r, int order, unsigned options) const
 {
     ex acc; // Series accumulator
     
     // Get first term from overall_coeff
-    acc = overall_coeff.series(r, order, branchcut);
+    acc = overall_coeff.series(r, order, options);
     
     // Add remaining terms
     epvector::const_iterator it = seq.begin();
@@ -630,7 +630,7 @@ ex add::series(const relational & r, int order, bool branchcut) const
         if (is_ex_exactly_of_type(it->rest, pseries))
             op = it->rest;
         else
-            op = it->rest.series(r, order, branchcut);
+            op = it->rest.series(r, order, options);
         if (!it->coeff.is_equal(_ex1()))
             op = ex_to_pseries(op).mul_const(ex_to_numeric(it->coeff));
         
@@ -720,12 +720,12 @@ ex pseries::mul_series(const pseries &other) const
 /** Implementation of ex::series() for product. This performs series
  *  multiplication when multiplying series.
  *  @see ex::series */
-ex mul::series(const relational & r, int order, bool branchcut) const
+ex mul::series(const relational & r, int order, unsigned options) const
 {
     ex acc; // Series accumulator
     
     // Get first term from overall_coeff
-    acc = overall_coeff.series(r, order, branchcut);
+    acc = overall_coeff.series(r, order, options);
     
     // Multiply with remaining terms
     epvector::const_iterator it = seq.begin();
@@ -738,7 +738,7 @@ ex mul::series(const relational & r, int order, bool branchcut) const
             acc = ex_to_pseries(acc).mul_const(ex_to_numeric(f));
             continue;
         } else if (!is_ex_exactly_of_type(op, pseries))
-            op = op.series(r, order, branchcut);
+            op = op.series(r, order, options);
         if (!it->coeff.is_equal(_ex1()))
             op = ex_to_pseries(op).power_const(ex_to_numeric(it->coeff), order);
 
@@ -810,20 +810,20 @@ pseries pseries::shift_exponents(int deg) const
 /** Implementation of ex::series() for powers. This performs Laurent expansion
  *  of reciprocals of series at singularities.
  *  @see ex::series */
-ex power::series(const relational & r, int order, bool branchcut) const
+ex power::series(const relational & r, int order, unsigned options) const
 {
     ex e;
     if (!is_ex_exactly_of_type(basis, pseries)) {
         // Basis is not a series, may there be a singulary?
         if (!exponent.info(info_flags::negint))
-            return basic::series(r, order, branchcut);
+            return basic::series(r, order, options);
         
         // Expression is of type something^(-int), check for singularity
         if (!basis.subs(r).is_zero())
-            return basic::series(r, order, branchcut);
+            return basic::series(r, order, options);
         
         // Singularity encountered, expand basis into series
-        e = basis.series(r, order, branchcut);
+        e = basis.series(r, order, options);
     } else {
         // Basis is a series
         e = basis;
@@ -835,7 +835,7 @@ ex power::series(const relational & r, int order, bool branchcut) const
 
 
 /** Re-expansion of a pseries object. */
-ex pseries::series(const relational & r, int order, bool branchcut) const
+ex pseries::series(const relational & r, int order, unsigned options) const
 {
     const ex p = r.rhs();
     GINAC_ASSERT(is_ex_exactly_of_type(r.lhs(),symbol));
@@ -859,7 +859,7 @@ ex pseries::series(const relational & r, int order, bool branchcut) const
             return pseries(r, new_seq);
         }
     } else
-        return convert_to_poly().series(r, order, branchcut);
+        return convert_to_poly().series(r, order, options);
 }
 
 
@@ -870,9 +870,9 @@ ex pseries::series(const relational & r, int order, bool branchcut) const
  *
  *  @param r  expansion relation, lhs holds variable and rhs holds point
  *  @param order  truncation order of series calculations
- *  @param branchcut  when set to false, branch cuts are not honored
+ *  @param options  of class series_options
  *  @return an expression holding a pseries object */
-ex ex::series(const ex & r, int order, bool branchcut) const
+ex ex::series(const ex & r, int order, unsigned options) const
 {
     GINAC_ASSERT(bp!=0);
     ex e;
@@ -886,7 +886,7 @@ ex ex::series(const ex & r, int order, bool branchcut) const
         throw (std::logic_error("ex::series(): expansion point has unknown type"));
     
     try {
-        e = bp->series(rel_, order, branchcut);
+        e = bp->series(rel_, order, options);
     } catch (std::exception &x) {
         throw (std::logic_error(std::string("unable to compute series (") + x.what() + ")"));
     }
diff --git a/ginac/pseries.h b/ginac/pseries.h
index f7c32aa5..48975724 100644
--- a/ginac/pseries.h
+++ b/ginac/pseries.h
@@ -67,7 +67,7 @@ public:
     ex collect(const symbol &s) const;
     ex eval(int level=0) const;
     ex evalf(int level=0) const;
-    ex series(const relational & r, int order, bool branchcut = true) const;
+    ex series(const relational & r, int order, unsigned options = 0) const;
     ex subs(const lst & ls, const lst & lr) const;
     ex normal(lst &sym_lst, lst &repl_lst, int level = 0) const;
     ex expand(unsigned options = 0) const;
diff --git a/ginac/symbol.h b/ginac/symbol.h
index 1a51ffb4..bd37892f 100644
--- a/ginac/symbol.h
+++ b/ginac/symbol.h
@@ -81,7 +81,7 @@ public:
     int ldegree(const symbol & s) const;
     ex coeff(const symbol & s, int n = 1) const;
     ex eval(int level = 0) const;
-    ex series(const relational & s, int order, bool branchcut = true) const;
+    ex series(const relational & s, int order, unsigned options = 0) const;
     ex normal(lst &sym_lst, lst &repl_lst, int level = 0) const;
     ex to_rational(lst &repl_lst) const;
     ex subs(const lst & ls, const lst & lr) const;