@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 \
@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 \
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
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_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
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)
@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 \
-# Doxyfile 0.49-991106
+# Doxyfile 0.49-991205
# This file describes the settings to be used by doxygen for a project
#
# 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
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
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
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>
-# Doxyfile 0.49-991106
+# Doxyfile 0.49-991205
# This file describes the settings to be used by doxygen for a project
#
# 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
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
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
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>
@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 \
@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 \
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)
@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 \
// 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) {
s->push_back(combine_ex_with_coeff_to_pair((*cit2).rest.subs(ls,lr),
(*cit2).coeff));
++cit2;
- }
+ }
return s;
}
++cit;
$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
#include "function.h"
#include "ex.h"
+#include "utils.h"
#include "debugmsg.h"
#ifndef NO_GINAC_NAMESPACE
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}
#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 {
* 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)) {
}
// 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());
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);
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);
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(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 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
}
-/** 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.
{
if (level == 1)
return this->hold();
-
+
// Construct a new series with evaluated coefficients
epvector new_seq;
new_seq.reserve(seq.size());
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
*/
{
ex e;
epvector::const_iterator it = seq.begin(), itend = seq.end();
-
+
while (it != itend) {
if (is_order_function(it->rest)) {
if (!no_order)
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));
if (!coeff.is_zero())
seq.push_back(expair(coeff, numeric(n)));
}
-
+
// Higher-order terms, if present
deriv = deriv.diff(s);
if (!deriv.is_zero())
nul.push_back(expair(Order(exONE()), exZERO()));
return series(var, point, nul);
}
-
+
// Series addition
epvector new_seq;
epvector::const_iterator a = seq.begin();
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) {
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
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();
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));
}
{
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))
// Series multiplication
epvector new_seq;
-
+
const symbol *s = static_cast<symbol *>(var.bp);
int a_max = degree(*s);
int b_max = other.degree(*s);
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_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)
}
-/** 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
{
}
*/
+/** 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();
int i;
const symbol *s = static_cast<symbol *>(var.bp);
int ldeg = ldegree(*s);
-
+
// Calculate coefficients of powered series
exvector co;
co.reserve(deg);
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;
// 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);
}
/** @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
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);
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
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
@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 \