- Introduced exception do_taylor to signal Taylor expansion is ok for series

- 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.

diff --git a/Makefile.in b/Makefile.in
index a6e7630faad8f9497a4ed774e3bbbfe2c64c376a..561c21bd607fb3e285a557aa7b9a7ef77c9a919a 100644 (file)
--- a/Makefile.in
+++ b/Makefile.in
@@ -352,7 +352,7 @@ distdir: $(DISTFILES)
        @for file in $(DISTFILES); do \
          d=$(srcdir); \
          if test -d $$d/$$file; then \
        @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 \
          else \
            test -f $(distdir)/$$file \
            || ln $$d/$$file $(distdir)/$$file 2> /dev/null \

diff --git a/check/Makefile.in b/check/Makefile.in
index 9189758af0aaa38175e03e81d797a33a75415126..f9e5c2cd51b697a7f73fc74cbf4bd3cdb1b6d7be 100644 (file)
--- a/check/Makefile.in
+++ b/check/Makefile.in
@@ -245,7 +245,7 @@ distdir: $(DISTFILES)
        @for file in $(DISTFILES); do \
          d=$(srcdir); \
          if test -d $$d/$$file; then \
        @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 \
          else \
            test -f $(distdir)/$$file \
            || ln $$d/$$file $(distdir)/$$file 2> /dev/null \

diff --git a/configure b/configure
index 4a6e3d2ccd7d4499d9621f7754c045cde6de73c5..273451401438724a46e41d7847e22c2bc697e847 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
 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
 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_MAJOR_VERSION=0
 GINACLIB_MINOR_VERSION=4

-GINACLIB_MICRO_VERSION=0
+GINACLIB_MICRO_VERSION=1

 GINACLIB_INTERFACE_AGE=0
 GINACLIB_INTERFACE_AGE=0

-GINACLIB_BINARY_AGE=0
+GINACLIB_BINARY_AGE=1

 GINACLIB_VERSION=$GINACLIB_MAJOR_VERSION.$GINACLIB_MINOR_VERSION.$GINACLIB_MICRO_VERSION
 
 
 GINACLIB_VERSION=$GINACLIB_MAJOR_VERSION.$GINACLIB_MINOR_VERSION.$GINACLIB_MICRO_VERSION
 
 

diff --git a/configure.in b/configure.in
index 96be9d3670070b41dfadcc5d6e01e83b704c7146..8e10c5a2f45c76ae786cd819067dff75c139e5d2 100644 (file)
--- a/configure.in
+++ b/configure.in
@@ -4,7 +4,7 @@ AC_INIT(ginac/basic.cpp)
 AC_PREREQ(2.13)
 
 dnl Configure options
 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
 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_MAJOR_VERSION=0
 GINACLIB_MINOR_VERSION=4

-GINACLIB_MICRO_VERSION=0
+GINACLIB_MICRO_VERSION=1

 GINACLIB_INTERFACE_AGE=0
 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)
 GINACLIB_VERSION=$GINACLIB_MAJOR_VERSION.$GINACLIB_MINOR_VERSION.$GINACLIB_MICRO_VERSION
 
 AC_SUBST(GINACLIB_MAJOR_VERSION)

diff --git a/doc/Makefile.in b/doc/Makefile.in
index 2498724b46169e857f3535e247a85fdb07d06b0f..4a11c86ad2afd723d04634719f786a65dfbb99a6 100644 (file)
--- a/doc/Makefile.in
+++ b/doc/Makefile.in
@@ -267,7 +267,7 @@ distdir: $(DISTFILES)
        @for file in $(DISTFILES); do \
          d=$(srcdir); \
          if test -d $$d/$$file; then \
        @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 \
          else \
            test -f $(distdir)/$$file \
            || ln $$d/$$file $(distdir)/$$file 2> /dev/null \

diff --git a/doc/reference/DoxyfileHTML b/doc/reference/DoxyfileHTML
index 65d561d162f6113300050477d493744056ee18a8..d9260e296cfd8752feabb430dcffcc543bdec0d1 100644 (file)
--- a/doc/reference/DoxyfileHTML
+++ b/doc/reference/DoxyfileHTML
@@ -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
 #
 
 # 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.
 
 # 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
 
 # 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
 
 
 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
 
 
 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
 # 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
 
 
 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
 # Doxygen uses this value to replace tabs by spaces in code fragments.
 
 TAB_SIZE             = 4


@@ -218,13 +222,6 @@ EXAMPLE_PATTERNS     =
 
 IMAGE_PATH           =
 
 
 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>
 # 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>

diff --git a/doc/reference/DoxyfileTEX b/doc/reference/DoxyfileTEX
index dd4c8804a0a4cd699d256cae11d628c4f12db506..ab7f05d4280a4516ce684331b0c5a9f918970263 100644 (file)
--- a/doc/reference/DoxyfileTEX
+++ b/doc/reference/DoxyfileTEX
@@ -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
 #
 
 # 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.
 
 # 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
 
 # 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
 
 
 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
 
 
 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
 # 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
 
 
 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
 # Doxygen uses this value to replace tabs by spaces in code fragments.
 
 TAB_SIZE             = 4


@@ -218,13 +222,6 @@ EXAMPLE_PATTERNS     =
 
 IMAGE_PATH           =
 
 
 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>
 # 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>

diff --git a/doc/reference/Makefile.in b/doc/reference/Makefile.in
index 152e777b7be0ea18850bd9a5f4e76a60931b803f..0931f3e08f2ea0b3c78f4f702289f09dc716797b 100644 (file)
--- a/doc/reference/Makefile.in
+++ b/doc/reference/Makefile.in
@@ -163,7 +163,7 @@ distdir: $(DISTFILES)
        @for file in $(DISTFILES); do \
          d=$(srcdir); \
          if test -d $$d/$$file; then \
        @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 \
          else \
            test -f $(distdir)/$$file \
            || ln $$d/$$file $(distdir)/$$file 2> /dev/null \

diff --git a/doc/tutorial/Makefile.in b/doc/tutorial/Makefile.in
index c1715667d27d96ba834bc84e0feb37b64db15189..92f9ec15058afa604de0de38ab7dd63e5935b1e8 100644 (file)
--- a/doc/tutorial/Makefile.in
+++ b/doc/tutorial/Makefile.in
@@ -315,7 +315,7 @@ distdir: $(DISTFILES)
        @for file in $(DISTFILES); do \
          d=$(srcdir); \
          if test -d $$d/$$file; then \
        @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 \
          else \
            test -f $(distdir)/$$file \
            || ln $$d/$$file $(distdir)/$$file 2> /dev/null \

diff --git a/doc/tutorial/ginac.texi b/doc/tutorial/ginac.texi
index 5e47963a71ab70dc75cb5687036f3b0bc3a7e796..028f2647cf53cc843b3a9c7051d34d97ac81c203 100644 (file)
--- a/doc/tutorial/ginac.texi
+++ b/doc/tutorial/ginac.texi
@@ -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
 
 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)
 
 @example
 DECLARE_FUNCTION_1P(cos)

diff --git a/ginac/Makefile.in b/ginac/Makefile.in
index 0a0c330651cbedf6c266594c2e2ffc53dcfcd4ce..90ece770fd93abef1098997b863c27e3c27ec900 100644 (file)
--- a/ginac/Makefile.in
+++ b/ginac/Makefile.in
@@ -288,7 +288,7 @@ distdir: $(DISTFILES)
        @for file in $(DISTFILES); do \
          d=$(srcdir); \
          if test -d $$d/$$file; then \
        @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 \
          else \
            test -f $(distdir)/$$file \
            || ln $$d/$$file $(distdir)/$$file 2> /dev/null \

diff --git a/ginac/expairseq.cpp b/ginac/expairseq.cpp
index 2e50b0e47d0ab582d1503d7682f25ec789d93c28..8a9c09525f021c176d71667d3abd4f949ed341ef 100644 (file)
--- a/ginac/expairseq.cpp
+++ b/ginac/expairseq.cpp
@@ -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)
     // 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)) {
     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());
             // 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) {
             // 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;
                 s->push_back(combine_ex_with_coeff_to_pair((*cit2).rest.subs(ls,lr),
                                                            (*cit2).coeff));
                 ++cit2;

-           }
+            }

             return s;
         }
         ++cit;
             return s;
         }
         ++cit;

diff --git a/ginac/function.pl b/ginac/function.pl
index 41018c93428c8cee29b830411ab43f92fc6a5e6f..b55845a3a16a0b107eba6225267f5fb0cf37ad03 100755 (executable)
--- a/ginac/function.pl
+++ b/ginac/function.pl
@@ -110,7 +110,12 @@ END_OF_DIFF_SWITCH_STATEMENT
 $series_switch_statement=generate(
     <<'END_OF_SERIES_SWITCH_STATEMENT','seq[${N}-1]','');
     case ${N}:
 $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
 
         break;
 END_OF_SERIES_SWITCH_STATEMENT
 


@@ -377,6 +382,7 @@ $implementation=<<END_OF_IMPLEMENTATION;
 
 #include "function.h"
 #include "ex.h"
 
 #include "function.h"
 #include "ex.h"

+#include "utils.h"

 #include "debugmsg.h"
 
 #ifndef NO_GINAC_NAMESPACE
 #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);
     }
     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}
     switch (registered_functions()[serial].nparams) {
         // the following lines have been generated for max. ${maxargs} parameters
 ${series_switch_statement}

diff --git a/ginac/inifcns_gamma.cpp b/ginac/inifcns_gamma.cpp
index 0d59eb306b67dd2612d55f3d648296d6f4f11700..7f3b9f5c924fcb7700e8dc4a8fcad4f4f2bd6c04 100644 (file)
--- a/ginac/inifcns_gamma.cpp
+++ b/ginac/inifcns_gamma.cpp
@@ -27,9 +27,12 @@
 #include "inifcns.h"
 #include "ex.h"
 #include "constant.h"
 #include "inifcns.h"
 #include "ex.h"
 #include "constant.h"

+#include "series.h"

 #include "numeric.h"
 #include "power.h"
 #include "numeric.h"
 #include "power.h"

+#include "relational.h"

 #include "symbol.h"
 #include "symbol.h"

+#include "utils.h"

 
 #ifndef NO_GINAC_NAMESPACE
 namespace GiNaC {
 
 #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...
  *
  *  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)) {
 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 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());
             // 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 {
                 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);
                 // 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)
 {
 
 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);
 }
 
 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;
 }
 
     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)
 
 //////////
 // 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);
     // 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();
 }    
     }
     return psi(n, x).hold();
 }    

diff --git a/ginac/inifcns_zeta.cpp b/ginac/inifcns_zeta.cpp
index c2a298a9aaf1661dedf2ce2245677e4982574d4e..e343b6fbe2692dc6c2eb92953309b65c0e2ec081 100644 (file)
--- a/ginac/inifcns_zeta.cpp
+++ b/ginac/inifcns_zeta.cpp
@@ -73,19 +73,7 @@ static ex zeta_evalf(ex const & x)
     return zeta(ex_to_numeric(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
 
 #ifndef NO_GINAC_NAMESPACE
 } // namespace GiNaC

diff --git a/ginac/normal.cpp b/ginac/normal.cpp
index 3846f408f0162617b3303dcd164dd4ad5a69f73b..3549dfee1f54df61d4c1889119651231f38be47b 100644 (file)
--- a/ginac/normal.cpp
+++ b/ginac/normal.cpp
@@ -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.
 class gcdheu_failed {};
 
 /** Compute GCD of multivariate polynomials using the heuristic GCD algorithm.

diff --git a/ginac/series.cpp b/ginac/series.cpp
index 7f8d39281bd5baf0a4fd95f70d6d5869450ff513..e21b2ba42be6bf97bbb0068b7eb1a74e9eb70556 100644 (file)
--- a/ginac/series.cpp
+++ b/ginac/series.cpp
@@ -182,7 +182,7 @@ ex series::eval(int level) const
 {
     if (level == 1)
         return this->hold();
 {
     if (level == 1)
         return this->hold();

-
+    

     // Construct a new series with evaluated coefficients
     epvector new_seq;
     new_seq.reserve(seq.size());
     // 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);
 }
 
     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);
 }
 
 ex series::evalf(int level) const
 {
     return convert_to_poly().evalf(level);
 }
 

-

 /*
  *  Construct expression (polynomial) out of series
  */
 /*
  *  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();
 {
     ex e;
     epvector::const_iterator it = seq.begin(), itend = seq.end();

-
+    

     while (it != itend) {
         if (is_order_function(it->rest)) {
             if (!no_order)
     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)));
     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));
     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)));
     }
         if (!coeff.is_zero())
             seq.push_back(expair(coeff, numeric(n)));
     }

-
+    

     // Higher-order terms, if present
     deriv = deriv.diff(s);
     if (!deriv.is_zero())
     // 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);
     }
         nul.push_back(expair(Order(exONE()), exZERO()));
         return series(var, point, nul);
     }

-
+    

     // Series addition
     epvector new_seq;
     epvector::const_iterator a = seq.begin();
     // 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();
             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) {
         // 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();
             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
         // 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
 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);
     // 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();
     // 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));
             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));
     }
         // 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 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))
     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;
 
     // Series multiplication
     epvector new_seq;

-
+    

     const symbol *s = static_cast<symbol *>(var.bp);
     int a_max = degree(*s);
     int b_max = other.degree(*s);
     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 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)))
     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;
     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)
     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
 {
 /*
 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
 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);
     // 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();
     // 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);
     int i;
     const symbol *s = static_cast<symbol *>(var.bp);
     int ldeg = ldegree(*s);

-
+    

     // Calculate coefficients of powered series
     exvector co;
     co.reserve(deg);
     // 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));
     }
             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;
     // 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);
         // 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);
         // 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;
     }
         // 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);
 }
     // Power e
     return ex_to_series(e).power_const(ex_to_numeric(exponent), order);
 }

diff --git a/ginac/utils.h b/ginac/utils.h
index 18a5e7bd3f86c85c17216ec9ef20793546244fcf..7be44f73281ed2c8829176f41938a509b1932ec1 100644 (file)
--- a/ginac/utils.h
+++ b/ginac/utils.h
@@ -1,6 +1,7 @@
 /** @file utils.h
  *
 /** @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
 
 /*
  *  GiNaC Copyright (C) 1999 Johannes Gutenberg University Mainz, Germany


@@ -40,6 +41,10 @@ string ToString(T const & t)
     return buf;
 }
 
     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);
 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) {
 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
 }
 
 // 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) {
                         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
 }
 
 #ifndef NO_GINAC_NAMESPACE

diff --git a/ginsh/Makefile.in b/ginsh/Makefile.in
index 8bb6e4f5793cba54cbbe692df65450516e7be19a..342ed4627716727b87f22f1cc5189715d2f8a0eb 100644 (file)
--- a/ginsh/Makefile.in
+++ b/ginsh/Makefile.in
@@ -309,7 +309,7 @@ distdir: $(DISTFILES)
        @for file in $(DISTFILES); do \
          d=$(srcdir); \
          if test -d $$d/$$file; then \
        @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 \
          else \
            test -f $(distdir)/$$file \
            || ln $$d/$$file $(distdir)/$$file 2> /dev/null \