Added integral class.
authorJens Vollinga <vollinga@thep.physik.uni-mainz.de>
Tue, 12 Oct 2004 13:32:11 +0000 (13:32 +0000)
committerJens Vollinga <vollinga@thep.physik.uni-mainz.de>
Tue, 12 Oct 2004 13:32:11 +0000 (13:32 +0000)
ginac/Makefile.am
ginac/basic.cpp
ginac/basic.h
ginac/ex.h
ginac/ginac.h
ginac/indexed.cpp
ginac/integral.cpp [new file with mode: 0644]
ginac/integral.h [new file with mode: 0644]
ginac/tinfos.h
ginac/wildcard.cpp
ginac/wildcard.h

index c7285c9..2db8a60 100644 (file)
@@ -5,10 +5,10 @@ libginac_la_SOURCES = add.cpp archive.cpp basic.cpp clifford.cpp color.cpp \
   constant.cpp ex.cpp expair.cpp expairseq.cpp exprseq.cpp fail.cpp \
   fderivative.cpp function.cpp idx.cpp indexed.cpp inifcns.cpp \
   inifcns_trans.cpp inifcns_gamma.cpp inifcns_nstdsums.cpp \
-  lst.cpp matrix.cpp mul.cpp ncmul.cpp normal.cpp numeric.cpp operators.cpp \
-  power.cpp registrar.cpp relational.cpp remember.cpp pseries.cpp print.cpp \
-  structure.cpp symbol.cpp symmetry.cpp tensor.cpp utils.cpp wildcard.cpp \
-  input_parser.yy input_lexer.ll \
+  integral.cpp lst.cpp matrix.cpp mul.cpp ncmul.cpp normal.cpp numeric.cpp \
+  operators.cpp power.cpp registrar.cpp relational.cpp remember.cpp \
+  pseries.cpp print.cpp structure.cpp symbol.cpp symmetry.cpp tensor.cpp \
+  utils.cpp wildcard.cpp input_parser.yy input_lexer.ll \
   input_lexer.h remember.h tostring.h utils.h
 libginac_la_LDFLAGS = -version-info $(LT_CURRENT):$(LT_REVISION):$(LT_AGE) \
   -release $(LT_RELEASE)
@@ -16,9 +16,9 @@ ginacincludedir = $(includedir)/ginac
 ginacinclude_HEADERS = ginac.h add.h archive.h assertion.h basic.h class_info.h \
   clifford.h color.h constant.h container.h ex.h expair.h expairseq.h exprseq.h \
   fail.h fderivative.h flags.h function.h hash_map.h idx.h indexed.h inifcns.h \
-  lst.h matrix.h mul.h ncmul.h normal.h numeric.h operators.h power.h print.h \
-  pseries.h ptr.h registrar.h relational.h structure.h symbol.h symmetry.h \
-  tensor.h tinfos.h version.h wildcard.h
+  integral.h lst.h matrix.h mul.h ncmul.h normal.h numeric.h operators.h \
+  power.h print.h pseries.h ptr.h registrar.h relational.h structure.h \
+  symbol.h symmetry.h tensor.h tinfos.h version.h wildcard.h
 AM_LFLAGS = -Pginac_yy -olex.yy.c
 AM_YFLAGS = -p ginac_yy -d
 EXTRA_DIST = function.pl input_parser.h version.h.in
index 98ba54d..c37ed3e 100644 (file)
@@ -482,6 +482,20 @@ ex basic::evalm() const
                return map(map_evalm);
 }
 
+/** Function object to be applied by basic::eval_integ(). */
+struct eval_integ_map_function : public map_function {
+       ex operator()(const ex & e) { return eval_integ(e); }
+} map_eval_integ;
+
+/** Evaluate integrals, if result is known. */
+ex basic::eval_integ() const
+{
+       if (nops() == 0)
+               return *this;
+       else
+               return map(map_eval_integ);
+}
+
 /** Perform automatic symbolic evaluations on indexed expression that
  *  contains this object as the base expression. */
 ex basic::eval_indexed(const basic & i) const
index a5b71a5..9f88d7b 100644 (file)
@@ -136,6 +136,7 @@ public: // only const functions please (may break reference counting)
        virtual ex eval(int level = 0) const;
        virtual ex evalf(int level = 0) const;
        virtual ex evalm() const;
+       virtual ex eval_integ() const;
 protected:
        virtual ex eval_ncmul(const exvector & v) const;
 public:
index b274b49..1c1e799 100644 (file)
@@ -117,6 +117,7 @@ public:
        ex evalf(int level = 0) const { return bp->evalf(level); }
        ex evalm() const { return bp->evalm(); }
        ex eval_ncmul(const exvector & v) const { return bp->eval_ncmul(v); }
+       ex eval_integ() const { return bp->eval_integ(); }
 
        // printing
        void print(const print_context & c, unsigned level = 0) const;
@@ -807,6 +808,9 @@ inline ex evalf(const ex & thisex, int level = 0)
 inline ex evalm(const ex & thisex)
 { return thisex.evalm(); }
 
+inline ex eval_integ(const ex & thisex)
+{ return thisex.eval_integ(); }
+
 inline ex diff(const ex & thisex, const symbol & s, unsigned nth = 1)
 { return thisex.diff(s, nth); }
 
index 2ba0ad7..54af6c4 100644 (file)
@@ -34,6 +34,7 @@
 
 #include "constant.h"
 #include "fail.h"
+#include "integral.h"
 #include "lst.h"
 #include "matrix.h"
 #include "numeric.h"
index d693373..29dd86a 100644 (file)
@@ -36,6 +36,7 @@
 #include "lst.h"
 #include "archive.h"
 #include "utils.h"
+#include "integral.h"
 
 namespace GiNaC {
 
@@ -522,6 +523,13 @@ exvector power::get_free_indices() const
                return basis_indices;
 }
 
+exvector integral::get_free_indices() const
+{
+       if (a.get_free_indices().size() || b.get_free_indices().size())
+               throw (std::runtime_error("integral::get_free_indices: boundary values should not have free indices"));
+       return f.get_free_indices();
+}
+
 /** Rename dummy indices in an expression.
  *
  *  @param e Expression to work on
diff --git a/ginac/integral.cpp b/ginac/integral.cpp
new file mode 100644 (file)
index 0000000..467cf29
--- /dev/null
@@ -0,0 +1,441 @@
+/** @file integral.cpp
+ *
+ *  Implementation of GiNaC's symbolic  integral. */
+
+/*
+ *  GiNaC Copyright (C) 1999-2004 Johannes Gutenberg University Mainz, Germany
+ *
+ *  This program is free software; you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation; either version 2 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program; if not, write to the Free Software
+ *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+ */
+
+#include "integral.h"
+#include "numeric.h"
+#include "symbol.h"
+#include "add.h"
+#include "mul.h"
+#include "power.h"
+#include "inifcns.h"
+#include "wildcard.h"
+#include "archive.h"
+#include "registrar.h"
+#include "utils.h"
+#include "operators.h"
+#include "relational.h"
+
+using namespace std;
+
+namespace GiNaC {
+
+GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(integral, basic,
+  print_func<print_dflt>(&integral::do_print).
+  print_func<print_latex>(&integral::do_print_latex))
+
+
+//////////
+// default constructor
+//////////
+
+integral::integral()
+               : inherited(TINFO_integral),
+               x((new symbol())->setflag(status_flags::dynallocated))
+{}
+
+//////////
+// other constructors
+//////////
+
+// public
+
+integral::integral(const ex & x_, const ex & a_, const ex & b_, const ex & f_)
+               : inherited(TINFO_integral), x(x_), a(a_), b(b_), f(f_)
+{
+       if (!is_a<symbol>(x)) {
+               throw(std::invalid_argument("first argument of integral must be of type symbol"));
+       }
+}
+
+//////////
+// archiving
+//////////
+
+integral::integral(const archive_node & n, lst & sym_lst) : inherited(n, sym_lst)
+{
+       n.find_ex("x", x, sym_lst);
+       n.find_ex("a", a, sym_lst);
+       n.find_ex("b", b, sym_lst);
+       n.find_ex("f", f, sym_lst);
+}
+
+void integral::archive(archive_node & n) const
+{
+       inherited::archive(n);
+       n.add_ex("x", x);
+       n.add_ex("a", a);
+       n.add_ex("b", b);
+       n.add_ex("f", f);
+}
+
+DEFAULT_UNARCHIVE(integral)
+
+//////////
+// functions overriding virtual functions from base classes
+//////////
+
+void integral::do_print(const print_context & c, unsigned level) const
+{
+       c.s << "integral(";
+       x.print(c);
+       c.s << ",";
+       a.print(c);
+       c.s << ",";
+       b.print(c);
+       c.s << ",";
+       f.print(c);
+       c.s << ")";
+}
+
+void integral::do_print_latex(const print_latex & c, unsigned level) const
+{
+       string varname = ex_to<symbol>(x).get_name();
+       if (level > precedence())
+               c.s << "\\left(";
+       c.s << "\\int_{";
+       a.print(c);
+       c.s << "}^{";
+       b.print(c);
+       c.s << "} d";
+       if (varname.size() > 1)
+               c.s << "\\," << varname << "\\:";
+       else
+               c.s << varname << "\\,";
+       f.print(c,precedence());
+       if (level > precedence())
+               c.s << "\\right)";
+}
+
+int integral::compare_same_type(const basic & other) const
+{
+       GINAC_ASSERT(is_exactly_a<integral>(other));
+       const integral &o = static_cast<const integral &>(other);
+
+       int cmpval = x.compare(o.x);
+       if (cmpval)
+               return cmpval;
+       cmpval = a.compare(o.a);
+       if (cmpval)
+               return cmpval;
+       cmpval = b.compare(o.b);
+       if (cmpval)
+               return cmpval;
+       return f.compare(o.f);
+}
+
+ex integral::eval(int level) const
+{
+       if ((level==1) && (flags & status_flags::evaluated))
+               return *this;
+       if (level == -max_recursion_level)
+               throw(std::runtime_error("max recursion level reached"));
+
+       ex eintvar = (level==1) ? x : x.eval(level-1);
+       ex ea      = (level==1) ? a : a.eval(level-1);
+       ex eb      = (level==1) ? b : b.eval(level-1);
+       ex ef      = (level==1) ? f : f.eval(level-1);
+
+       if (!ef.has(eintvar) && !haswild(ef))
+               return eb*ef-ea*ef;
+
+       if (ea==eb)
+               return _ex0;
+
+       if (are_ex_trivially_equal(eintvar,x) && are_ex_trivially_equal(ea,a)
+                       && are_ex_trivially_equal(eb,b) && are_ex_trivially_equal(ef,f))
+               return this->hold();
+       return (new integral(eintvar, ea, eb, ef))
+               ->setflag(status_flags::dynallocated | status_flags::evaluated);
+}
+
+ex integral::evalf(int level) const
+{
+       ex ea;
+       ex eb;
+       ex ef;
+
+       if (level==1) {
+               ea = a;
+               eb = b;
+               ef = f;
+       } else if (level == -max_recursion_level) {
+               throw(runtime_error("max recursion level reached"));
+       } else {
+               ea = a.evalf(level-1);
+               eb = b.evalf(level-1);
+               ef = f.evalf(level-1);
+       }
+
+       // 12.34 is just an arbitrary number used to check whether a number
+       // results after subsituting a number for the integration variable.
+       if (is_exactly_a<numeric>(ea) && is_exactly_a<numeric>(eb) 
+                       && is_exactly_a<numeric>(ef.subs(x==12.34).evalf())) {
+               try {
+                       return adaptivesimpson(x, ea, eb, ef);
+               } catch (runtime_error &rte) {}
+       }
+
+       if (are_ex_trivially_equal(a, ea) && are_ex_trivially_equal(b, eb)
+                               && are_ex_trivially_equal(f, ef))
+                       return *this;
+               else
+                       return (new integral(x, ea, eb, ef))
+                               ->setflag(status_flags::dynallocated);
+}
+
+int integral::max_integration_level = 15;
+ex integral::relative_integration_error = power(10,-8).evalf();
+
+ex subsvalue(const ex & var, const ex & value, const ex & fun)
+{
+       ex result = fun.subs(var==value).evalf();
+       if (is_a<numeric>(result))
+               return result;
+       throw logic_error("integrant does not evaluate to numeric");
+}
+
+/** Numeric integration routine based upon the "Adaptive Quadrature" one
+  * in "Numerical Analysis" by Burden and Faires. Parameters are integration
+  * variable, left boundary, right boundary, function to be integrated and
+  * the relative integration error. The function should evalf into a number
+  * after substituting the integration variable by a number. Another thing
+  * to note is that this implementation is no good at integrating functions
+  * with discontinuities. */
+ex adaptivesimpson(const ex & x, const ex & a, const ex & b, const ex & f, const ex & error)
+{
+       // use lookup table to be potentially much faster.
+       static exmap lookup;
+       static symbol ivar("ivar");
+       ex lookupex = integral(ivar,a,b,f.subs(x==ivar));
+       exmap::iterator emi = lookup.find(lookupex);
+       if (emi!=lookup.end())
+               return emi->second;
+
+       ex app = 0;
+       int i = 1;
+       exvector avec(integral::max_integration_level+1);
+       exvector hvec(integral::max_integration_level+1);
+       exvector favec(integral::max_integration_level+1);
+       exvector fbvec(integral::max_integration_level+1);
+       exvector fcvec(integral::max_integration_level+1);
+       exvector svec(integral::max_integration_level+1);
+       exvector errorvec(integral::max_integration_level+1);
+       vector<int> lvec(integral::max_integration_level+1);
+
+       avec[i] = a;
+       hvec[i] = (b-a)/2;
+       favec[i] = subsvalue(x, a, f);
+       fcvec[i] = subsvalue(x, a+hvec[i], f);
+       fbvec[i] = subsvalue(x, b, f);
+       svec[i] = hvec[i]*(favec[i]+4*fcvec[i]+fbvec[i])/3;
+       lvec[i] = 1;
+       errorvec[i] = integral::relative_integration_error*svec[i];
+
+       while (i>0) {
+               ex fd = subsvalue(x, avec[i]+hvec[i]/2, f);
+               ex fe = subsvalue(x, avec[i]+3*hvec[i]/2, f);
+               ex s1 = hvec[i]*(favec[i]+4*fd+fcvec[i])/6;
+               ex s2 = hvec[i]*(fcvec[i]+4*fe+fbvec[i])/6;
+               ex nu1 = avec[i];
+               ex nu2 = favec[i];
+               ex nu3 = fcvec[i];
+               ex nu4 = fbvec[i];
+               ex nu5 = hvec[i];
+               // hopefully prevents a crash if the function is zero sometimes.
+               ex nu6 = max(errorvec[i], (s1+s2)*integral::relative_integration_error);
+               ex nu7 = svec[i];
+               int nu8 = lvec[i];
+               --i;
+               if (abs(ex_to<numeric>(s1+s2-nu7)) <= nu6)
+                       app+=(s1+s2);
+               else {
+                       if (nu8>=integral::max_integration_level)
+                               throw runtime_error("max integration level reached");
+                       ++i;
+                       avec[i] = nu1+nu5;
+                       favec[i] = nu3;
+                       fcvec[i] = fe;
+                       fbvec[i] = nu4;
+                       hvec[i] = nu5/2;
+                       errorvec[i]=nu6/2;
+                       svec[i] = s2;
+                       lvec[i] = nu8+1;
+                       ++i;
+                       avec[i] = nu1;
+                       favec[i] = nu2;
+                       fcvec[i] = fd;
+                       fbvec[i] = nu3;
+                       hvec[i] = hvec[i-1];
+                       errorvec[i]=errorvec[i-1];
+                       svec[i] = s1;
+                       lvec[i] = lvec[i-1];
+               }
+       }
+
+       lookup[lookupex]=app;
+       return app;
+}
+
+int integral::degree(const ex & s) const
+{
+       return ((b-a)*f).degree(s);
+}
+
+int integral::ldegree(const ex & s) const
+{
+       return ((b-a)*f).ldegree(s);
+}
+
+ex integral::eval_ncmul(const exvector & v) const
+{
+       return f.eval_ncmul(v);
+}
+
+size_t integral::nops() const
+{
+       return 4;
+}
+
+ex integral::op(size_t i) const
+{
+       GINAC_ASSERT(i<4);
+
+       switch(i) {
+               case(0):
+                       return x;
+               case(1):
+                       return a;
+               case(2):
+                       return b;
+               case(3):
+                       return f;
+       }
+}
+
+ex & integral::let_op(size_t i)
+{
+       ensure_if_modifiable();
+       switch(i) {
+               case(0):
+                       return x;
+               case(1):
+                       return a;
+               case(2):
+                       return b;
+               case(3):
+                       return f;
+       }
+}
+
+ex integral::expand(unsigned options) const
+{
+       if (options==0 && (flags & status_flags::expanded))
+               return *this;
+
+       ex newa = a.expand(options);
+       ex newb = b.expand(options);
+       ex newf = f.expand(options);
+
+       if (is_a<add>(newf)) {
+               exvector v;
+               v.reserve(newf.nops());
+               for (size_t i=0; i<newf.nops(); ++i)
+                       v.push_back(integral(x, newa, newb, newf.op(i)).expand(options));
+               return ex(add(v)).expand(options);
+       }
+
+       if (is_a<mul>(newf)) {
+               ex prefactor = 1;
+               ex rest = 1;
+               for (size_t i=0; i<newf.nops(); ++i)
+                       if (newf.op(i).has(x))
+                               rest *= newf.op(i);
+                       else
+                               prefactor *= newf.op(i);
+               if (prefactor != 1)
+                       return (prefactor*integral(x, newa, newb, rest)).expand(options);
+       }
+
+       if (are_ex_trivially_equal(a, newa) && are_ex_trivially_equal(b, newb)
+                       && are_ex_trivially_equal(f, newf)) {
+               if (options==0)
+                       this->setflag(status_flags::expanded);
+               return *this;
+       }
+
+       const basic & newint = (new integral(x, newa, newb, newf))
+               ->setflag(status_flags::dynallocated);
+       if (options == 0)
+               newint.setflag(status_flags::expanded);
+       return newint;
+}
+
+ex integral::derivative(const symbol & s) const
+{      if (s==x)
+               throw(logic_error("differentiation with respect to dummy variable"));
+       return b.diff(s)*f.subs(x==b)-a.diff(s)*f.subs(x==a)+integral(x, a, b, f.diff(s));
+}
+
+unsigned integral::return_type() const
+{
+       return f.return_type();
+}
+
+unsigned integral::return_type_tinfo() const
+{
+       return f.return_type_tinfo();
+}
+
+ex integral::conjugate() const
+{
+       ex conja = a.conjugate();
+       ex conjb = b.conjugate();
+       ex conjf = f.conjugate().subs(x.conjugate()==x);
+
+       if (are_ex_trivially_equal(a, conja) && are_ex_trivially_equal(b, conjb)
+                       && are_ex_trivially_equal(f, conjf))
+               return *this;
+
+       return (new integral(x, conja, conjb, conjf))
+               ->setflag(status_flags::dynallocated);
+}
+
+ex integral::eval_integ() const
+{
+       if (!(flags & status_flags::expanded))
+               return this->expand().eval_integ();
+       
+       if (f==x)
+               return b*b/2-a*a/2;
+       if (is_a<power>(f) && f.op(0)==x) {
+               if (f.op(1)==-1)
+                       return log(b/a);
+               if (!f.op(1).has(x)) {
+                       ex primit = power(x,f.op(1)+1)/(f.op(1)+1);
+                       return primit.subs(x==b)-primit.subs(x==a);
+               }
+       }
+
+       return *this;
+}
+
+} // namespace GiNaC
diff --git a/ginac/integral.h b/ginac/integral.h
new file mode 100644 (file)
index 0000000..ef2e662
--- /dev/null
@@ -0,0 +1,96 @@
+/** @file integral.h
+ *
+ *  Interface to GiNaC's symbolic  integral. */
+
+/*
+ *  GiNaC Copyright (C) 1999-2004 Johannes Gutenberg University Mainz, Germany
+ *
+ *  This program is free software; you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation; either version 2 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program; if not, write to the Free Software
+ *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+ */
+
+#ifndef __GINAC_INTEGRAL_H__
+#define __GINAC_INTEGRAL_H__
+
+#include "basic.h"
+#include "ex.h"
+
+namespace GiNaC {
+
+/** Product of expressions. */
+class integral : public basic
+{
+       GINAC_DECLARE_REGISTERED_CLASS(integral, basic)
+       
+       // other constructors
+public:
+       integral(const ex & x_, const ex & a_, const ex & b_, const ex & f_);
+       
+       // functions overriding virtual functions from base classes
+public:
+       unsigned precedence() const {return 45;}
+       ex eval(int level=0) const;
+       ex evalf(int level=0) const;
+       int degree(const ex & s) const;
+       int ldegree(const ex & s) const;
+       ex eval_ncmul(const exvector & v) const;
+       size_t nops() const;
+       ex op(size_t i) const;
+       ex & let_op(size_t i);
+       ex expand(unsigned options = 0) const;
+       exvector get_free_indices() const;
+       unsigned return_type() const;
+       unsigned return_type_tinfo() const;
+       ex conjugate() const;
+       ex eval_integ() const;
+protected:
+       ex derivative(const symbol & s) const;
+
+       // new virtual functions which can be overridden by derived classes
+       // none
+       
+       // non-virtual functions in this class
+protected:
+       void do_print(const print_context & c, unsigned level) const;
+       void do_print_latex(const print_latex & c, unsigned level) const;
+public:
+       static int max_integration_level;
+       static ex relative_integration_error;
+private:
+       ex x;
+       ex a;
+       ex b;
+       ex f;
+};
+
+// utility functions
+
+/** Specialization of is_exactly_a<integral>(obj) for integral objects. */
+template<> inline bool is_exactly_a<integral>(const basic & obj)
+{
+       return obj.tinfo()==TINFO_integral;
+}
+
+
+GiNaC::ex adaptivesimpson(
+       const GiNaC::ex &x,
+       const GiNaC::ex &a,
+       const GiNaC::ex &b,
+       const GiNaC::ex &f,
+       const GiNaC::ex &error = integral::relative_integration_error
+);
+
+} // namespace GiNaC
+
+#endif // ndef __GINAC_INTEGRAL_H__
index 37d681b..f6b0551 100644 (file)
@@ -38,6 +38,7 @@ const unsigned TINFO_exprseq       = 0x00030001U;
 const unsigned TINFO_function      = 0x00031001U;
 const unsigned TINFO_fderivative   = 0x00032001U;
 const unsigned TINFO_ncmul         = 0x00031002U;
+const unsigned TINFO_integral      = 0x00033001U;
 
 const unsigned TINFO_lst           = 0x00040001U;
 
index 7ebc023..1bcf2ac 100644 (file)
@@ -119,4 +119,14 @@ bool wildcard::match(const ex & pattern, lst & repl_lst) const
        return is_equal(ex_to<basic>(pattern));
 }
 
+bool haswild(const ex & x)
+{
+       if (is_a<wildcard>(x))
+               return true;
+       for (int i=0; i<x.nops(); ++i)
+               if (haswild(x.op(i)))
+                       return true;
+       return false;
+}
+
 } // namespace GiNaC
index 1842ce7..e52c3fa 100644 (file)
@@ -75,6 +75,9 @@ inline ex wild(unsigned label = 0)
        return wildcard(label);
 }
 
+/** Check whether x has a wildcard anywhere as a subexpression. */
+bool haswild(const ex & x);
+
 } // namespace GiNaC
 
 #endif // ndef __GINAC_WILDCARD_H__