- Introduced exception do_taylor to signal Taylor expansion is ok for series
authorRichard Kreckel <Richard.Kreckel@uni-mainz.de>
Fri, 10 Dec 1999 14:59:18 +0000 (14:59 +0000)
committerRichard Kreckel <Richard.Kreckel@uni-mainz.de>
Fri, 10 Dec 1999 14:59:18 +0000 (14:59 +0000)
- Finished a clean implementation of gamma function's series expansion
- Bumped up version number to 1.4.1 to reflect the changes and the
  soon-to-come second prerelease.

19 files changed:
Makefile.in
check/Makefile.in
configure
configure.in
doc/Makefile.in
doc/reference/DoxyfileHTML
doc/reference/DoxyfileTEX
doc/reference/Makefile.in
doc/tutorial/Makefile.in
doc/tutorial/ginac.texi
ginac/Makefile.in
ginac/expairseq.cpp
ginac/function.pl
ginac/inifcns_gamma.cpp
ginac/inifcns_zeta.cpp
ginac/normal.cpp
ginac/series.cpp
ginac/utils.h
ginsh/Makefile.in

index a6e7630..561c21b 100644 (file)
@@ -352,7 +352,7 @@ distdir: $(DISTFILES)
        @for file in $(DISTFILES); do \
          d=$(srcdir); \
          if test -d $$d/$$file; then \
-           cp -pr $$/$$file $(distdir)/$$file; \
+           cp -pr $$d/$$file $(distdir)/$$file; \
          else \
            test -f $(distdir)/$$file \
            || ln $$d/$$file $(distdir)/$$file 2> /dev/null \
index 9189758..f9e5c2c 100644 (file)
@@ -245,7 +245,7 @@ distdir: $(DISTFILES)
        @for file in $(DISTFILES); do \
          d=$(srcdir); \
          if test -d $$d/$$file; then \
-           cp -pr $$/$$file $(distdir)/$$file; \
+           cp -pr $$d/$$file $(distdir)/$$file; \
          else \
            test -f $(distdir)/$$file \
            || ln $$d/$$file $(distdir)/$$file 2> /dev/null \
index 4a6e3d2..2734514 100755 (executable)
--- a/configure
+++ b/configure
@@ -12,7 +12,7 @@ ac_help=
 ac_default_prefix=/usr/local
 # Any additions from configure.in:
 ac_help="$ac_help
-  --enable-help-doc       build HTML documentation [default=yes]"
+  --enable-html-doc       build HTML documentation [default=yes]"
 ac_help="$ac_help
   --enable-ps-doc         build PostScript documentation [default=yes]"
 ac_help="$ac_help
@@ -560,9 +560,9 @@ fi
 
 GINACLIB_MAJOR_VERSION=0
 GINACLIB_MINOR_VERSION=4
-GINACLIB_MICRO_VERSION=0
+GINACLIB_MICRO_VERSION=1
 GINACLIB_INTERFACE_AGE=0
-GINACLIB_BINARY_AGE=0
+GINACLIB_BINARY_AGE=1
 GINACLIB_VERSION=$GINACLIB_MAJOR_VERSION.$GINACLIB_MINOR_VERSION.$GINACLIB_MICRO_VERSION
 
 
index 96be9d3..8e10c5a 100644 (file)
@@ -4,7 +4,7 @@ AC_INIT(ginac/basic.cpp)
 AC_PREREQ(2.13)
 
 dnl Configure options
-AC_ARG_ENABLE(html-doc, [  --enable-help-doc       build HTML documentation [default=yes]], , enable_html_doc=yes)
+AC_ARG_ENABLE(html-doc, [  --enable-html-doc       build HTML documentation [default=yes]], , enable_html_doc=yes)
 AC_ARG_ENABLE(ps-doc,   [  --enable-ps-doc         build PostScript documentation [default=yes]], , enable_ps_doc=yes)
 
 dnl GiNaC version information
@@ -23,9 +23,9 @@ dnl (don't we all *love* autoconf?)...
 
 GINACLIB_MAJOR_VERSION=0
 GINACLIB_MINOR_VERSION=4
-GINACLIB_MICRO_VERSION=0
+GINACLIB_MICRO_VERSION=1
 GINACLIB_INTERFACE_AGE=0
-GINACLIB_BINARY_AGE=0
+GINACLIB_BINARY_AGE=1
 GINACLIB_VERSION=$GINACLIB_MAJOR_VERSION.$GINACLIB_MINOR_VERSION.$GINACLIB_MICRO_VERSION
 
 AC_SUBST(GINACLIB_MAJOR_VERSION)
index 2498724..4a11c86 100644 (file)
@@ -267,7 +267,7 @@ distdir: $(DISTFILES)
        @for file in $(DISTFILES); do \
          d=$(srcdir); \
          if test -d $$d/$$file; then \
-           cp -pr $$/$$file $(distdir)/$$file; \
+           cp -pr $$d/$$file $(distdir)/$$file; \
          else \
            test -f $(distdir)/$$file \
            || ln $$d/$$file $(distdir)/$$file 2> /dev/null \
index 65d561d..d9260e2 100644 (file)
@@ -1,4 +1,4 @@
-# Doxyfile 0.49-991106
+# Doxyfile 0.49-991205
 
 # This file describes the settings to be used by doxygen for a project
 #
@@ -46,7 +46,7 @@ QUIET                = NO
 # generated by doxygen. Possible values are YES and NO. If left blank
 # NO is used.
 
-WARNINGS             = NO
+WARNINGS             = YES
 
 # The DISABLE_INDEX tag can be used to turn on/off the condensed index at
 # top of each HTML page. The value NO (the default) enables the index and
@@ -120,12 +120,16 @@ INTERNAL_DOCS        = NO
 
 CLASS_DIAGRAMS       = YES
 
-# If the SOURCE_BROWSER tag is set to YES than the body of a member or
-# function will be appended as a block of code to the documentation of.
-# that member or function.
+# If the SOURCE_BROWSER tag is set to YES then a list of source files will
+# be generated. Documented entities will be cross-referenced with these sources.
 
 SOURCE_BROWSER       = YES
 
+# Setting the INLINE_SOURCES tag to YES will include the body
+# of functions and classes directly in the documentation.
+
+INLINE_SOURCES       = YES
+
 # If the CASE_SENSE_NAMES tag is set to NO (the default) then Doxygen
 # will only generate file names in lower case letters. If set to
 # YES upper case letters are also allowed. This is useful if you have
@@ -158,7 +162,7 @@ INHERIT_DOCS         = YES
 
 INLINE_INFO          = YES
 
-# the TAB_SIZE tag can be used to set the number of spaces in a tab
+# the TAB_SIZE tag can be used to set the number of spaces in a tab.
 # Doxygen uses this value to replace tabs by spaces in code fragments.
 
 TAB_SIZE             = 4
@@ -218,13 +222,6 @@ EXAMPLE_PATTERNS     =
 
 IMAGE_PATH           =
 
-# If the value of the IMAGE_PATH tag contains directories, you can use the
-# IMAGE_PATTERNS tag to specify one or more wildcard pattern (like *.gif 
-# and *.eps) to filter out the image files in the directories. If left 
-# blank all files are included.
-
-IMAGE_PATTERNS       =
-
 # The INPUT_FILTER tag can be used to specify a program that doxygen should
 # invoke to filter for each input file. Doxygen will invoke the filter program 
 # by executing (via popen()) the command <filter> <input-file>, where <filter>
index dd4c880..ab7f05d 100644 (file)
@@ -1,4 +1,4 @@
-# Doxyfile 0.49-991106
+# Doxyfile 0.49-991205
 
 # This file describes the settings to be used by doxygen for a project
 #
@@ -46,7 +46,7 @@ QUIET                = NO
 # generated by doxygen. Possible values are YES and NO. If left blank
 # NO is used.
 
-WARNINGS             = NO
+WARNINGS             = YES
 
 # The DISABLE_INDEX tag can be used to turn on/off the condensed index at
 # top of each HTML page. The value NO (the default) enables the index and
@@ -120,12 +120,16 @@ INTERNAL_DOCS        = NO
 
 CLASS_DIAGRAMS       = YES
 
-# If the SOURCE_BROWSER tag is set to YES than the body of a member or
-# function will be appended as a block of code to the documentation of.
-# that member or function.
+# If the SOURCE_BROWSER tag is set to YES then a list of source files will
+# be generated. Documented entities will be cross-referenced with these sources.
 
 SOURCE_BROWSER       = NO
 
+# Setting the INLINE_SOURCES tag to YES will include the body
+# of functions and classes directly in the documentation.
+
+INLINE_SOURCES       = NO
+
 # If the CASE_SENSE_NAMES tag is set to NO (the default) then Doxygen
 # will only generate file names in lower case letters. If set to
 # YES upper case letters are also allowed. This is useful if you have
@@ -158,7 +162,7 @@ INHERIT_DOCS         = YES
 
 INLINE_INFO          = YES
 
-# the TAB_SIZE tag can be used to set the number of spaces in a tab
+# the TAB_SIZE tag can be used to set the number of spaces in a tab.
 # Doxygen uses this value to replace tabs by spaces in code fragments.
 
 TAB_SIZE             = 4
@@ -218,13 +222,6 @@ EXAMPLE_PATTERNS     =
 
 IMAGE_PATH           =
 
-# If the value of the IMAGE_PATH tag contains directories, you can use the
-# IMAGE_PATTERNS tag to specify one or more wildcard pattern (like *.gif 
-# and *.eps) to filter out the image files in the directories. If left 
-# blank all files are included.
-
-IMAGE_PATTERNS       =
-
 # The INPUT_FILTER tag can be used to specify a program that doxygen should
 # invoke to filter for each input file. Doxygen will invoke the filter program 
 # by executing (via popen()) the command <filter> <input-file>, where <filter>
index 152e777..0931f3e 100644 (file)
@@ -163,7 +163,7 @@ distdir: $(DISTFILES)
        @for file in $(DISTFILES); do \
          d=$(srcdir); \
          if test -d $$d/$$file; then \
-           cp -pr $$/$$file $(distdir)/$$file; \
+           cp -pr $$d/$$file $(distdir)/$$file; \
          else \
            test -f $(distdir)/$$file \
            || ln $$d/$$file $(distdir)/$$file 2> /dev/null \
index c171566..92f9ec1 100644 (file)
@@ -315,7 +315,7 @@ distdir: $(DISTFILES)
        @for file in $(DISTFILES); do \
          d=$(srcdir); \
          if test -d $$d/$$file; then \
-           cp -pr $$/$$file $(distdir)/$$file; \
+           cp -pr $$d/$$file $(distdir)/$$file; \
          else \
            test -f $(distdir)/$$file \
            || ln $$d/$$file $(distdir)/$$file 2> /dev/null \
index 5e47963..028f264 100644 (file)
@@ -1566,9 +1566,15 @@ REGISTER_FUNCTION(cos, cos_eval_method, cos_evalf_method, cos_diff, NULL);
 
 The first argument is the function's name, the second, third and fourth
 bind the corresponding methods to this objects and the fifth is a slot
-for inserting a method for series expansion.  Also, the new function
-needs to be declared somewhere.  This may also be done by a convenient
-preprocessor macro:
+for inserting a method for series expansion.  (If set to @code{NULL} it
+defaults to simple Taylor expansion, which is correct if there are no
+poles involved.  The way GiNaC handles poles in case there are any is
+best understood by studying one of the examples, like the Gamma function
+for instance.  In essence the function first checks if there is a pole
+at the evaluation point and falls back to Taylor expansion if there
+isn't.  Then, the pole is regularized by some suitable transformation.)
+Also, the new function needs to be declared somewhere.  This may also be
+done by a convenient preprocessor macro:
 
 @example
 DECLARE_FUNCTION_1P(cos)
index 0a0c330..90ece77 100644 (file)
@@ -288,7 +288,7 @@ distdir: $(DISTFILES)
        @for file in $(DISTFILES); do \
          d=$(srcdir); \
          if test -d $$d/$$file; then \
-           cp -pr $$/$$file $(distdir)/$$file; \
+           cp -pr $$d/$$file $(distdir)/$$file; \
          else \
            test -f $(distdir)/$$file \
            || ln $$d/$$file $(distdir)/$$file 2> /dev/null \
index 2e50b0e..8a9c095 100644 (file)
@@ -1504,17 +1504,18 @@ epvector * expairseq::subschildren(lst const & ls, lst const & lr) const
     // returns a NULL pointer if nothing had to be substituted
     // returns a pointer to a newly created epvector otherwise
     // (which has to be deleted somewhere else)
-
+    GINAC_ASSERT(ls.nops()==lr.nops());
+    
     epvector::const_iterator last=seq.end();
     epvector::const_iterator cit=seq.begin();
     while (cit!=last) {
         ex const & subsed_ex=(*cit).rest.subs(ls,lr);
         if (!are_ex_trivially_equal((*cit).rest,subsed_ex)) {
-
+            
             // something changed, copy seq, subs and return it
             epvector *s=new epvector;
             s->reserve(seq.size());
-
+            
             // copy parts of seq which are known not to have changed
             epvector::const_iterator cit2=seq.begin();
             while (cit2!=cit) {
@@ -1530,7 +1531,7 @@ epvector * expairseq::subschildren(lst const & ls, lst const & lr) const
                 s->push_back(combine_ex_with_coeff_to_pair((*cit2).rest.subs(ls,lr),
                                                            (*cit2).coeff));
                 ++cit2;
-           }
+            }
             return s;
         }
         ++cit;
index 41018c9..b55845a 100755 (executable)
@@ -110,7 +110,12 @@ END_OF_DIFF_SWITCH_STATEMENT
 $series_switch_statement=generate(
     <<'END_OF_SERIES_SWITCH_STATEMENT','seq[${N}-1]','');
     case ${N}:
-        return ((series_funcp_${N})(registered_functions()[serial].s))(${SEQ1},s,point,order);
+        try {
+            res = ((series_funcp_${N})(registered_functions()[serial].s))(${SEQ1},s,point,order);
+        } catch (do_taylor) {
+            res = basic::series(s, point, order);
+        }
+        return res;
         break;
 END_OF_SERIES_SWITCH_STATEMENT
 
@@ -377,6 +382,7 @@ $implementation=<<END_OF_IMPLEMENTATION;
 
 #include "function.h"
 #include "ex.h"
+#include "utils.h"
 #include "debugmsg.h"
 
 #ifndef NO_GINAC_NAMESPACE
@@ -602,6 +608,7 @@ ex function::series(symbol const & s, ex const & point, int order) const
     if (registered_functions()[serial].s==0) {
         return basic::series(s, point, order);
     }
+    ex res;
     switch (registered_functions()[serial].nparams) {
         // the following lines have been generated for max. ${maxargs} parameters
 ${series_switch_statement}
index 0d59eb3..7f3b9f5 100644 (file)
 #include "inifcns.h"
 #include "ex.h"
 #include "constant.h"
+#include "series.h"
 #include "numeric.h"
 #include "power.h"
+#include "relational.h"
 #include "symbol.h"
+#include "utils.h"
 
 #ifndef NO_GINAC_NAMESPACE
 namespace GiNaC {
@@ -43,7 +46,7 @@ namespace GiNaC {
  *  arguments and that's it. Somebody ought to provide some good numerical
  *  evaluation some day...
  *
- *  @exception fail_numeric("complex_infinity") or something similar... */
+ *  @exception std::domain_error("gamma_eval(): simple pole") */
 static ex gamma_eval(ex const & x)
 {
     if (x.info(info_flags::numeric)) {
@@ -58,7 +61,7 @@ static ex gamma_eval(ex const & x)
         }
         // trap half integer arguments:
         if ((x*2).info(info_flags::integer)) {
-            // trap positive x=(n+1/2)
+            // trap positive x==(n+1/2)
             // gamma(n+1/2) -> Pi^(1/2)*(1*3*..*(2*n-1))/(2^n)
             if ((x*2).info(info_flags::posint)) {
                 numeric n = ex_to_numeric(x).sub(numHALF());
@@ -66,7 +69,7 @@ static ex gamma_eval(ex const & x)
                 coefficient = coefficient.div(numTWO().power(n));
                 return coefficient * pow(Pi,numHALF());
             } else {
-                // trap negative x=(-n+1/2)
+                // trap negative x==(-n+1/2)
                 // gamma(-n+1/2) -> Pi^(1/2)*(-2)^n/(1*3*..*(2*n-1))
                 numeric n = abs(ex_to_numeric(x).sub(numHALF()));
                 numeric coefficient = numeric(-2).power(n);
@@ -96,12 +99,21 @@ static ex gamma_diff(ex const & x, unsigned diff_param)
 
 static ex gamma_series(ex const & x, symbol const & s, ex const & point, int order)
 {
-       // FIXME: Only handle one special case for now...
-       if (x.is_equal(s) && point.is_zero()) {
-               ex e = 1 / s - EulerGamma + s * (pow(Pi, 2) / 12 + pow(EulerGamma, 2) / 2) + Order(pow(s, 2));
-               return e.series(s, point, order);
-       } else
-               throw(std::logic_error("don't know the series expansion of this particular gamma function"));
+    // method:
+    // Taylor series where there is no pole falls back to psi functions.
+    // On a pole at -n use the identity
+    //   series(GAMMA(x),x=-n,order) ==
+    //   series(GAMMA(x+n+1)/(x*(x+1)...*(x+n)),x=-n,order+1);
+    ex xpoint = x.subs(s==point);
+    if (!xpoint.info(info_flags::integer) || xpoint.info(info_flags::positive))
+        throw do_taylor();
+    // if we got here we have to care for a simple pole at -n:
+    numeric n = -ex_to_numeric(xpoint);
+    ex ser_numer = gamma(x+n+exONE());
+    ex ser_denom = exONE();
+    for (numeric p; p<=n; ++p)
+        ser_denom *= x+p;
+    return (ser_numer/ser_denom).series(s, point, order+1);
 }
 
 REGISTER_FUNCTION(gamma, gamma_eval, gamma_evalf, gamma_diff, gamma_series);
@@ -166,16 +178,7 @@ static ex beta_diff(ex const & x, ex const & y, unsigned diff_param)
     return retval;
 }
 
-static ex beta_series(ex const & x, ex const & y, symbol const & s, ex const & point, int order)
-{
-       if (x.is_equal(s) && point.is_zero()) {
-               ex e = 1 / s - EulerGamma + s * (pow(Pi, 2) / 12 + pow(EulerGamma, 2) / 2) + Order(pow(s, 2));
-               return e.series(s, point, order);
-       } else
-               throw(std::logic_error("don't know the series expansion of this particular beta function"));
-}
-
-REGISTER_FUNCTION(beta, beta_eval, beta_evalf, beta_diff, beta_series);
+REGISTER_FUNCTION(beta, beta_eval, beta_evalf, beta_diff, NULL);
 
 //////////
 // Psi-function (aka polygamma-function)
@@ -240,8 +243,15 @@ static ex psi2_eval(ex const & n, ex const & x)
     // psi(0,x) -> psi(x)
     if (n.is_zero())
         return psi(x);
-    if (n.info(info_flags::numeric) && x.info(info_flags::numeric)) {
-        // do some stuff...
+    // psi(-1,x) -> log(gamma(x))
+    if (n.is_equal(exMINUSONE()))
+        return log(gamma(x));
+    if (n.info(info_flags::numeric) && n.info(info_flags::posint) &&
+        x.info(info_flags::numeric)) {
+        numeric nn = ex_to_numeric(n);
+        numeric nx = ex_to_numeric(x);
+        if (x.is_equal(exONE()))
+            return numMINUSONE().power(nn+numONE())*factorial(nn)*zeta(ex(nn+numONE()));
     }
     return psi(n, x).hold();
 }    
index c2a298a..e343b6f 100644 (file)
@@ -73,19 +73,7 @@ static ex zeta_evalf(ex const & x)
     return zeta(ex_to_numeric(x));
 }
 
-static ex zeta_diff(ex const & x, unsigned diff_param)
-{
-    GINAC_ASSERT(diff_param==0);
-    
-    return exZERO();  // should return zeta(numONE(),x);
-}
-
-static ex zeta_series(ex const & x, symbol const & s, ex const & point, int order)
-{
-    throw(std::logic_error("don't know the series expansion of the zeta function"));
-}
-
-REGISTER_FUNCTION(zeta, zeta_eval, zeta_evalf, zeta_diff, zeta_series);
+REGISTER_FUNCTION(zeta, zeta_eval, zeta_evalf, NULL, NULL);
 
 #ifndef NO_GINAC_NAMESPACE
 } // namespace GiNaC
index 3846f40..3549dfe 100644 (file)
@@ -966,7 +966,7 @@ ex mul::smod(const numeric &xi) const
 }
 
 
-/** Exception thrown by heur_gcd() to signal failure */
+/** Exception thrown by heur_gcd() to signal failure. */
 class gcdheu_failed {};
 
 /** Compute GCD of multivariate polynomials using the heuristic GCD algorithm.
index 7f8d392..e21b2ba 100644 (file)
@@ -182,7 +182,7 @@ ex series::eval(int level) const
 {
     if (level == 1)
         return this->hold();
-
+    
     // Construct a new series with evaluated coefficients
     epvector new_seq;
     new_seq.reserve(seq.size());
@@ -194,12 +194,12 @@ ex series::eval(int level) const
     return (new series(var, point, new_seq))->setflag(status_flags::dynallocated | status_flags::evaluated);
 }
 
+/** Evaluate numerically.  The order term is dropped. */
 ex series::evalf(int level) const
 {
     return convert_to_poly().evalf(level);
 }
 
-
 /*
  *  Construct expression (polynomial) out of series
  */
@@ -211,7 +211,7 @@ ex series::convert_to_poly(bool no_order) const
 {
     ex e;
     epvector::const_iterator it = seq.begin(), itend = seq.end();
-
+    
     while (it != itend) {
         if (is_order_function(it->rest)) {
             if (!no_order)
@@ -238,7 +238,7 @@ ex basic::series(symbol const & s, ex const & point, int order) const
     ex coeff = deriv.subs(s == point);
     if (!coeff.is_zero())
         seq.push_back(expair(coeff, numeric(0)));
-
+    
     int n;
     for (n=1; n<order; n++) {
         fac = fac.mul(numeric(n));
@@ -251,7 +251,7 @@ ex basic::series(symbol const & s, ex const & point, int order) const
         if (!coeff.is_zero())
             seq.push_back(expair(coeff, numeric(n)));
     }
-
+    
     // Higher-order terms, if present
     deriv = deriv.diff(s);
     if (!deriv.is_zero())
@@ -274,7 +274,7 @@ ex series::add_series(const series &other) const
         nul.push_back(expair(Order(exONE()), exZERO()));
         return series(var, point, nul);
     }
-
+    
     // Series addition
     epvector new_seq;
     epvector::const_iterator a = seq.begin();
@@ -292,7 +292,7 @@ ex series::add_series(const series &other) const
             break;
         } else
             pow_a = ex_to_numeric((*a).coeff).to_int();
-
+        
         // If b is empty, fill up with elements from a and stop
         if (b == b_end) {
             while (a != a_end) {
@@ -302,7 +302,7 @@ ex series::add_series(const series &other) const
             break;
         } else
             pow_b = ex_to_numeric((*b).coeff).to_int();
-
+        
         // a and b are non-empty, compare powers
         if (pow_a < pow_b) {
             // a has lesser power, get coefficient from a
@@ -374,10 +374,10 @@ ex add::series(symbol const & s, ex const & point, int order) const
 ex add::series(symbol const & s, ex const & point, int order) const
 {
     ex acc; // Series accumulator
-
+    
     // Get first term from overall_coeff
     acc = overall_coeff.series(s,point,order);
-
+    
     // Add remaining terms
     epvector::const_iterator it = seq.begin();
     epvector::const_iterator itend = seq.end();
@@ -389,7 +389,7 @@ ex add::series(symbol const & s, ex const & point, int order) const
             op = it->rest.series(s, point, order);
         if (!it->coeff.is_equal(exONE()))
             op = ex_to_series(op).mul_const(ex_to_numeric(it->coeff));
-
+        
         // Series addition
         acc = ex_to_series(acc).add_series(ex_to_series(op));
     }
@@ -406,7 +406,7 @@ ex series::mul_const(const numeric &other) const
 {
     epvector new_seq;
     new_seq.reserve(seq.size());
-
+    
     epvector::const_iterator it = seq.begin(), itend = seq.end();
     while (it != itend) {
         if (!is_order_function(it->rest))
@@ -436,7 +436,7 @@ ex series::mul_series(const series &other) const
 
     // Series multiplication
     epvector new_seq;
-
+    
     const symbol *s = static_cast<symbol *>(var.bp);
     int a_max = degree(*s);
     int b_max = other.degree(*s);
@@ -444,7 +444,7 @@ ex series::mul_series(const series &other) const
     int b_min = other.ldegree(*s);
     int cdeg_min = a_min + b_min;
     int cdeg_max = a_max + b_max;
-
+    
     int higher_order_a = INT_MAX;
     int higher_order_b = INT_MAX;
     if (is_order_function(coeff(*s, a_max)))
@@ -454,7 +454,7 @@ ex series::mul_series(const series &other) const
     int higher_order_c = min(higher_order_a, higher_order_b);
     if (cdeg_max >= higher_order_c)
         cdeg_max = higher_order_c - 1;
-
+    
     for (int cdeg=cdeg_min; cdeg<=cdeg_max; cdeg++) {
         ex co = exZERO();
         // c(i)=a(0)b(i)+...+a(i)b(0)
@@ -473,8 +473,6 @@ ex series::mul_series(const series &other) const
 }
 
 
-/** Implementation of ex::series() for product. This performs series multiplication when multiplying series.
- *  @see ex::series */
 /*
 ex mul::series(symbol const & s, ex const & point, int order) const
 {
@@ -513,13 +511,16 @@ ex mul::series(symbol const & s, ex const & point, int order) const
 }
 */
 
+/** Implementation of ex::series() for product. This performs series
+ *  multiplication when multiplying series.
+ *  @see ex::series */
 ex mul::series(symbol const & s, ex const & point, int order) const
 {
     ex acc; // Series accumulator
-
+    
     // Get first term from overall_coeff
     acc = overall_coeff.series(s, point, order);
-
+    
     // Multiply with remaining terms
     epvector::const_iterator it = seq.begin();
     epvector::const_iterator itend = seq.end();
@@ -551,7 +552,7 @@ ex series::power_const(const numeric &p, int deg) const
     int i;
     const symbol *s = static_cast<symbol *>(var.bp);
     int ldeg = ldegree(*s);
-
+    
     // Calculate coefficients of powered series
     exvector co;
     co.reserve(deg);
@@ -572,7 +573,7 @@ ex series::power_const(const numeric &p, int deg) const
             all_sums_zero = false;
         co.push_back(co0 * sum / numeric(i));
     }
-
+    
     // Construct new series (of non-zero coefficients)
     epvector new_seq;
     bool higher_order = false;
@@ -600,18 +601,18 @@ ex power::series(symbol const & s, ex const & point, int order) const
         // Basis is not a series, may there be a singulary?
         if (!exponent.info(info_flags::negint))
             return basic::series(s, point, order);
-
+        
         // Expression is of type something^(-int), check for singularity
         if (!basis.subs(s == point).is_zero())
             return basic::series(s, point, order);
-
+        
         // Singularity encountered, expand basis into series
         e = basis.series(s, point, order);
     } else {
         // Basis is a series
         e = basis;
     }
-
+    
     // Power e
     return ex_to_series(e).power_const(ex_to_numeric(exponent), order);
 }
index 18a5e7b..7be44f7 100644 (file)
@@ -1,6 +1,7 @@
 /** @file utils.h
  *
- *  Interface to several small and furry utilities. */
+ *  Interface to several small and furry utilities needed within GiNaC but not
+ *  of interest to the user of the library. */
 
 /*
  *  GiNaC Copyright (C) 1999 Johannes Gutenberg University Mainz, Germany
@@ -40,6 +41,10 @@ string ToString(T const & t)
     return buf;
 }
 
+/** Exception thrown by classes which provide their own series expansion to
+ *  signal that ordinary Taylor expansion is safe. */
+class do_taylor {};
+
 unsigned log2(unsigned n);
 
 int compare_pointers(void const * a, void const * b);
@@ -91,18 +96,18 @@ template <class InputIterator1, class InputIterator2, class OutputIterator,
 OutputIterator mymerge(InputIterator1 first1, InputIterator1 last1,
                        InputIterator2 first2, InputIterator2 last2,
                        OutputIterator result, Compare comp) {
-  while (first1 != last1 && first2 != last2) {
-    if (comp(*first1, *first2)) {
-      *result = *first1;
-      ++first1;
-    }
-    else {
-      *result = *first2;
-      ++first2;
+    while (first1 != last1 && first2 != last2) {
+        if (comp(*first1, *first2)) {
+            *result = *first1;
+            ++first1;
+        }
+        else {
+            *result = *first2;
+            ++first2;
+        }
+        ++result;
     }
-    ++result;
-  }
-  return copy(first2, last2, copy(first1, last1, result));
+    return copy(first2, last2, copy(first1, last1, result));
 }
 
 // like merge(), but three lists with *last2<*first3
@@ -112,25 +117,25 @@ OutputIterator mymerge3(InputIterator1 first1, InputIterator1 last1,
                         InputIterator2 first2, InputIterator2 last2,
                         InputIterator3 first3, InputIterator3 last3,
                         OutputIterator result, Compare comp) {
-  while (first1 != last1 && first2 != last2) {
-    if (comp(*first1, *first2)) {
-      *result = *first1;
-      ++first1;
+    while (first1 != last1 && first2 != last2) {
+        if (comp(*first1, *first2)) {
+            *result = *first1;
+            ++first1;
+        }
+        else {
+            *result = *first2;
+            ++first2;
+        }
+        ++result;
     }
-    else {
-      *result = *first2;
-      ++first2;
+    
+    if (first1==last1) {
+        // list1 empty, copy rest of list2, then list3
+        return copy(first3, last3, copy(first2, last2, result));
+    } else {
+        // list2 empty, merge rest of list1 with list3
+        return mymerge(first1,last1,first3,last3,result,comp);
     }
-    ++result;
-  }
-
-  if (first1==last1) {
-    // list1 empty, copy rest of list2, then list3
-    return copy(first3, last3, copy(first2, last2, result));
-  } else {
-    // list2 empty, merge rest of list1 with list3
-    return mymerge(first1,last1,first3,last3,result,comp);
-  }
 }
 
 #ifndef NO_GINAC_NAMESPACE
index 8bb6e4f..342ed46 100644 (file)
@@ -309,7 +309,7 @@ distdir: $(DISTFILES)
        @for file in $(DISTFILES); do \
          d=$(srcdir); \
          if test -d $$d/$$file; then \
-           cp -pr $$/$$file $(distdir)/$$file; \
+           cp -pr $$d/$$file $(distdir)/$$file; \
          else \
            test -f $(distdir)/$$file \
            || ln $$d/$$file $(distdir)/$$file 2> /dev/null \