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