0f27df6964fd2ed88cb2cb1a5101f2a8734dd190
[ginac.git] / ginac / integral.cpp
1 /** @file integral.cpp
2  *
3  *  Implementation of GiNaC's symbolic  integral. */
4
5 /*
6  *  GiNaC Copyright (C) 1999-2016 Johannes Gutenberg University Mainz, Germany
7  *
8  *  This program is free software; you can redistribute it and/or modify
9  *  it under the terms of the GNU General Public License as published by
10  *  the Free Software Foundation; either version 2 of the License, or
11  *  (at your option) any later version.
12  *
13  *  This program is distributed in the hope that it will be useful,
14  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
16  *  GNU General Public License for more details.
17  *
18  *  You should have received a copy of the GNU General Public License
19  *  along with this program; if not, write to the Free Software
20  *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
21  */
22
23 #include "integral.h"
24 #include "numeric.h"
25 #include "symbol.h"
26 #include "add.h"
27 #include "mul.h"
28 #include "power.h"
29 #include "inifcns.h"
30 #include "wildcard.h"
31 #include "archive.h"
32 #include "registrar.h"
33 #include "utils.h"
34 #include "operators.h"
35 #include "relational.h"
36
37 using namespace std;
38
39 namespace GiNaC {
40
41 GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(integral, basic,
42   print_func<print_dflt>(&integral::do_print).
43   print_func<print_latex>(&integral::do_print_latex))
44
45
46 //////////
47 // default constructor
48 //////////
49
50 integral::integral()
51                 : x(dynallocate<symbol>())
52 {}
53
54 //////////
55 // other constructors
56 //////////
57
58 // public
59
60 integral::integral(const ex & x_, const ex & a_, const ex & b_, const ex & f_)
61                 :  x(x_), a(a_), b(b_), f(f_)
62 {
63         if (!is_a<symbol>(x)) {
64                 throw(std::invalid_argument("first argument of integral must be of type symbol"));
65         }
66 }
67
68 //////////
69 // archiving
70 //////////
71
72 void integral::read_archive(const archive_node& n, lst& sym_lst)
73 {
74         inherited::read_archive(n, sym_lst);
75         n.find_ex("x", x, sym_lst);
76         n.find_ex("a", a, sym_lst);
77         n.find_ex("b", b, sym_lst);
78         n.find_ex("f", f, sym_lst);
79 }
80
81 void integral::archive(archive_node & n) const
82 {
83         inherited::archive(n);
84         n.add_ex("x", x);
85         n.add_ex("a", a);
86         n.add_ex("b", b);
87         n.add_ex("f", f);
88 }
89
90 //////////
91 // functions overriding virtual functions from base classes
92 //////////
93
94 void integral::do_print(const print_context & c, unsigned level) const
95 {
96         c.s << "integral(";
97         x.print(c);
98         c.s << ",";
99         a.print(c);
100         c.s << ",";
101         b.print(c);
102         c.s << ",";
103         f.print(c);
104         c.s << ")";
105 }
106
107 void integral::do_print_latex(const print_latex & c, unsigned level) const
108 {
109         string varname = ex_to<symbol>(x).get_name();
110         if (level > precedence())
111                 c.s << "\\left(";
112         c.s << "\\int_{";
113         a.print(c);
114         c.s << "}^{";
115         b.print(c);
116         c.s << "} d";
117         if (varname.size() > 1)
118                 c.s << "\\," << varname << "\\:";
119         else
120                 c.s << varname << "\\,";
121         f.print(c,precedence());
122         if (level > precedence())
123                 c.s << "\\right)";
124 }
125
126 int integral::compare_same_type(const basic & other) const
127 {
128         GINAC_ASSERT(is_exactly_a<integral>(other));
129         const integral &o = static_cast<const integral &>(other);
130
131         int cmpval = x.compare(o.x);
132         if (cmpval)
133                 return cmpval;
134         cmpval = a.compare(o.a);
135         if (cmpval)
136                 return cmpval;
137         cmpval = b.compare(o.b);
138         if (cmpval)
139                 return cmpval;
140         return f.compare(o.f);
141 }
142
143 ex integral::eval() const
144 {
145         if (flags & status_flags::evaluated)
146                 return *this;
147
148         if (!f.has(x) && !haswild(f))
149                 return b*f-a*f;
150
151         if (a==b)
152                 return _ex0;
153
154         return this->hold();
155 }
156
157 ex integral::evalf(int level) const
158 {
159         ex ea;
160         ex eb;
161         ex ef;
162
163         if (level==1) {
164                 ea = a;
165                 eb = b;
166                 ef = f;
167         } else if (level == -max_recursion_level) {
168                 throw(runtime_error("max recursion level reached"));
169         } else {
170                 ea = a.evalf(level-1);
171                 eb = b.evalf(level-1);
172                 ef = f.evalf(level-1);
173         }
174
175         // 12.34 is just an arbitrary number used to check whether a number
176         // results after substituting a number for the integration variable.
177         if (is_exactly_a<numeric>(ea) && is_exactly_a<numeric>(eb) &&
178             is_exactly_a<numeric>(ef.subs(x==12.34).evalf())) {
179                         return adaptivesimpson(x, ea, eb, ef);
180         }
181
182         if (are_ex_trivially_equal(a, ea) && are_ex_trivially_equal(b, eb) &&
183             are_ex_trivially_equal(f, ef))
184                         return *this;
185                 else
186                         return dynallocate<integral>(x, ea, eb, ef);
187 }
188
189 int integral::max_integration_level = 15;
190 ex integral::relative_integration_error = 1e-8;
191
192 ex subsvalue(const ex & var, const ex & value, const ex & fun)
193 {
194         ex result = fun.subs(var==value).evalf();
195         if (is_a<numeric>(result))
196                 return result;
197         throw logic_error("integrand does not evaluate to numeric");
198 }
199
200 struct error_and_integral
201 {
202         error_and_integral(const ex &err, const ex &integ)
203                 :error(err), integral(integ){}
204         ex error;
205         ex integral;
206 };
207
208 struct error_and_integral_is_less
209 {
210         bool operator()(const error_and_integral &e1,const error_and_integral &e2) const
211         {
212                 int c = e1.integral.compare(e2.integral);
213                 if(c < 0)
214                         return true;
215                 if(c > 0)
216                         return false;
217                 return ex_is_less()(e1.error, e2.error);
218         }
219 };
220
221 typedef map<error_and_integral, ex, error_and_integral_is_less> lookup_map;
222
223 /** Numeric integration routine based upon the "Adaptive Quadrature" one
224   * in "Numerical Analysis" by Burden and Faires. Parameters are integration
225   * variable, left boundary, right boundary, function to be integrated and
226   * the relative integration error. The function should evalf into a number
227   * after substituting the integration variable by a number. Another thing
228   * to note is that this implementation is no good at integrating functions
229   * with discontinuities. */
230 ex adaptivesimpson(const ex & x, const ex & a_in, const ex & b_in, const ex & f, const ex & error)
231 {
232         // Check whether boundaries and error are numbers.
233         ex a = is_exactly_a<numeric>(a_in) ? a_in : a_in.evalf();
234         ex b = is_exactly_a<numeric>(b_in) ? b_in : b_in.evalf();
235         if(!is_exactly_a<numeric>(a) || !is_exactly_a<numeric>(b))
236                 throw std::runtime_error("For numerical integration the boundaries of the integral should evalf into numbers.");
237         if(!is_exactly_a<numeric>(error))
238                 throw std::runtime_error("For numerical integration the error should be a number.");
239
240         // Use lookup table to be potentially much faster.
241         static lookup_map lookup;
242         static symbol ivar("ivar");
243         ex lookupex = integral(ivar,a,b,f.subs(x==ivar));
244         lookup_map::iterator emi = lookup.find(error_and_integral(error, lookupex));
245         if (emi!=lookup.end())
246                 return emi->second;
247
248         ex app = 0;
249         int i = 1;
250         exvector avec(integral::max_integration_level+1);
251         exvector hvec(integral::max_integration_level+1);
252         exvector favec(integral::max_integration_level+1);
253         exvector fbvec(integral::max_integration_level+1);
254         exvector fcvec(integral::max_integration_level+1);
255         exvector svec(integral::max_integration_level+1);
256         exvector errorvec(integral::max_integration_level+1);
257         vector<int> lvec(integral::max_integration_level+1);
258
259         avec[i] = a;
260         hvec[i] = (b-a)/2;
261         favec[i] = subsvalue(x, a, f);
262         fcvec[i] = subsvalue(x, a+hvec[i], f);
263         fbvec[i] = subsvalue(x, b, f);
264         svec[i] = hvec[i]*(favec[i]+4*fcvec[i]+fbvec[i])/3;
265         lvec[i] = 1;
266         errorvec[i] = error*abs(svec[i]);
267
268         while (i>0) {
269                 ex fd = subsvalue(x, avec[i]+hvec[i]/2, f);
270                 ex fe = subsvalue(x, avec[i]+3*hvec[i]/2, f);
271                 ex s1 = hvec[i]*(favec[i]+4*fd+fcvec[i])/6;
272                 ex s2 = hvec[i]*(fcvec[i]+4*fe+fbvec[i])/6;
273                 ex nu1 = avec[i];
274                 ex nu2 = favec[i];
275                 ex nu3 = fcvec[i];
276                 ex nu4 = fbvec[i];
277                 ex nu5 = hvec[i];
278                 // hopefully prevents a crash if the function is zero sometimes.
279                 ex nu6 = max(errorvec[i], abs(s1+s2)*error);
280                 ex nu7 = svec[i];
281                 int nu8 = lvec[i];
282                 --i;
283                 if (abs(ex_to<numeric>(s1+s2-nu7)) <= nu6)
284                         app+=(s1+s2);
285                 else {
286                         if (nu8>=integral::max_integration_level)
287                                 throw runtime_error("max integration level reached");
288                         ++i;
289                         avec[i] = nu1+nu5;
290                         favec[i] = nu3;
291                         fcvec[i] = fe;
292                         fbvec[i] = nu4;
293                         hvec[i] = nu5/2;
294                         errorvec[i]=nu6/2;
295                         svec[i] = s2;
296                         lvec[i] = nu8+1;
297                         ++i;
298                         avec[i] = nu1;
299                         favec[i] = nu2;
300                         fcvec[i] = fd;
301                         fbvec[i] = nu3;
302                         hvec[i] = hvec[i-1];
303                         errorvec[i]=errorvec[i-1];
304                         svec[i] = s1;
305                         lvec[i] = lvec[i-1];
306                 }
307         }
308
309         lookup[error_and_integral(error, lookupex)]=app;
310         return app;
311 }
312
313 int integral::degree(const ex & s) const
314 {
315         return ((b-a)*f).degree(s);
316 }
317
318 int integral::ldegree(const ex & s) const
319 {
320         return ((b-a)*f).ldegree(s);
321 }
322
323 ex integral::eval_ncmul(const exvector & v) const
324 {
325         return f.eval_ncmul(v);
326 }
327
328 size_t integral::nops() const
329 {
330         return 4;
331 }
332
333 ex integral::op(size_t i) const
334 {
335         GINAC_ASSERT(i<4);
336
337         switch (i) {
338                 case 0:
339                         return x;
340                 case 1:
341                         return a;
342                 case 2:
343                         return b;
344                 case 3:
345                         return f;
346                 default:
347                         throw (std::out_of_range("integral::op() out of range"));
348         }
349 }
350
351 ex & integral::let_op(size_t i)
352 {
353         ensure_if_modifiable();
354         switch (i) {
355                 case 0:
356                         return x;
357                 case 1:
358                         return a;
359                 case 2:
360                         return b;
361                 case 3:
362                         return f;
363                 default:
364                         throw (std::out_of_range("integral::let_op() out of range"));
365         }
366 }
367
368 ex integral::expand(unsigned options) const
369 {
370         if (options==0 && (flags & status_flags::expanded))
371                 return *this;
372
373         ex newa = a.expand(options);
374         ex newb = b.expand(options);
375         ex newf = f.expand(options);
376
377         if (is_a<add>(newf)) {
378                 exvector v;
379                 v.reserve(newf.nops());
380                 for (size_t i=0; i<newf.nops(); ++i)
381                         v.push_back(integral(x, newa, newb, newf.op(i)).expand(options));
382                 return ex(add(v)).expand(options);
383         }
384
385         if (is_a<mul>(newf)) {
386                 ex prefactor = 1;
387                 ex rest = 1;
388                 for (size_t i=0; i<newf.nops(); ++i)
389                         if (newf.op(i).has(x))
390                                 rest *= newf.op(i);
391                         else
392                                 prefactor *= newf.op(i);
393                 if (prefactor != 1)
394                         return (prefactor*integral(x, newa, newb, rest)).expand(options);
395         }
396
397         if (are_ex_trivially_equal(a, newa) && are_ex_trivially_equal(b, newb) &&
398             are_ex_trivially_equal(f, newf)) {
399                 if (options==0)
400                         this->setflag(status_flags::expanded);
401                 return *this;
402         }
403
404         const integral & newint = dynallocate<integral>(x, newa, newb, newf);
405         if (options == 0)
406                 newint.setflag(status_flags::expanded);
407         return newint;
408 }
409
410 ex integral::derivative(const symbol & s) const
411 {
412         if (s==x)
413                 throw(logic_error("differentiation with respect to dummy variable"));
414         return b.diff(s)*f.subs(x==b)-a.diff(s)*f.subs(x==a)+integral(x, a, b, f.diff(s));
415 }
416
417 unsigned integral::return_type() const
418 {
419         return f.return_type();
420 }
421
422 return_type_t integral::return_type_tinfo() const
423 {
424         return f.return_type_tinfo();
425 }
426
427 ex integral::conjugate() const
428 {
429         ex conja = a.conjugate();
430         ex conjb = b.conjugate();
431         ex conjf = f.conjugate().subs(x.conjugate()==x);
432
433         if (are_ex_trivially_equal(a, conja) && are_ex_trivially_equal(b, conjb) &&
434             are_ex_trivially_equal(f, conjf))
435                 return *this;
436
437         return dynallocate<integral>(x, conja, conjb, conjf);
438 }
439
440 ex integral::eval_integ() const
441 {
442         if (!(flags & status_flags::expanded))
443                 return this->expand().eval_integ();
444         
445         if (f==x)
446                 return b*b/2-a*a/2;
447         if (is_a<power>(f) && f.op(0)==x) {
448                 if (f.op(1)==-1)
449                         return log(b/a);
450                 if (!f.op(1).has(x)) {
451                         ex primit = power(x,f.op(1)+1)/(f.op(1)+1);
452                         return primit.subs(x==b)-primit.subs(x==a);
453                 }
454         }
455
456         return *this;
457 }
458
459 GINAC_BIND_UNARCHIVER(integral);
460 } // namespace GiNaC