7becb3e172da764458a6a28251c1526e8beec511
[ginac.git] / ginac / ex.cpp
1 /** @file ex.cpp
2  *
3  *  Implementation of GiNaC's light-weight expression handles. */
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 <iostream>
24 #include <stdexcept>
25
26 #include "ex.h"
27 #include "add.h"
28 #include "mul.h"
29 #include "ncmul.h"
30 #include "numeric.h"
31 #include "matrix.h"
32 #include "power.h"
33 #include "lst.h"
34 #include "relational.h"
35 #include "input_lexer.h"
36 #include "utils.h"
37
38 namespace GiNaC {
39
40 //////////
41 // other constructors
42 //////////
43
44 // none (all inlined)
45
46 //////////
47 // non-virtual functions in this class
48 //////////
49
50 // public
51         
52 /** Print expression to stream. The formatting of the output is determined
53  *  by the kind of print_context object that is passed. Possible formattings
54  *  include ginsh-parsable output (the default), tree-like output for
55  *  debugging, and C++ source.
56  *  @see print_context */
57 void ex::print(const print_context & c, unsigned level) const
58 {
59         bp->print(c, level);
60 }
61
62 /** Little wrapper arount print to be called within a debugger. */
63 void ex::dbgprint() const
64 {
65         bp->dbgprint();
66 }
67
68 /** Little wrapper arount printtree to be called within a debugger. */
69 void ex::dbgprinttree() const
70 {
71         bp->dbgprinttree();
72 }
73
74 ex ex::expand(unsigned options) const
75 {
76         if (options == 0 && (bp->flags & status_flags::expanded)) // The "expanded" flag only covers the standard options; someone might want to re-expand with different options
77                 return *this;
78         else
79                 return bp->expand(options);
80 }
81
82 /** Compute partial derivative of an expression.
83  *
84  *  @param s  symbol by which the expression is derived
85  *  @param nth  order of derivative (default 1)
86  *  @return partial derivative as a new expression */
87 ex ex::diff(const symbol & s, unsigned nth) const
88 {
89         if (!nth)
90                 return *this;
91         else
92                 return bp->diff(s, nth);
93 }
94
95 /** Check whether expression matches a specified pattern. */
96 bool ex::match(const ex & pattern) const
97 {
98         exmap repl_lst;
99         return bp->match(pattern, repl_lst);
100 }
101
102 /** Find all occurrences of a pattern. The found matches are appended to
103  *  the "found" list. If the expression itself matches the pattern, the
104  *  children are not further examined. This function returns true when any
105  *  matches were found. */
106 bool ex::find(const ex & pattern, exset& found) const
107 {
108         if (match(pattern)) {
109                 found.insert(*this);
110                 return true;
111         }
112         bool any_found = false;
113         for (size_t i=0; i<nops(); i++)
114                 if (op(i).find(pattern, found))
115                         any_found = true;
116         return any_found;
117 }
118
119 /** Substitute objects in an expression (syntactic substitution) and return
120  *  the result as a new expression. */
121 ex ex::subs(const lst & ls, const lst & lr, unsigned options) const
122 {
123         GINAC_ASSERT(ls.nops() == lr.nops());
124
125         // Convert the lists to a map
126         exmap m;
127         for (lst::const_iterator its = ls.begin(), itr = lr.begin(); its != ls.end(); ++its, ++itr) {
128                 m.insert(std::make_pair(*its, *itr));
129
130                 // Search for products and powers in the expressions to be substituted
131                 // (for an optimization in expairseq::subs())
132                 if (is_exactly_a<mul>(*its) || is_exactly_a<power>(*its))
133                         options |= subs_options::pattern_is_product;
134         }
135         if (!(options & subs_options::pattern_is_product))
136                 options |= subs_options::pattern_is_not_product;
137
138         return bp->subs(m, options);
139 }
140
141 /** Substitute objects in an expression (syntactic substitution) and return
142  *  the result as a new expression.  There are two valid types of
143  *  replacement arguments: 1) a relational like object==ex and 2) a list of
144  *  relationals lst(object1==ex1,object2==ex2,...). */
145 ex ex::subs(const ex & e, unsigned options) const
146 {
147         if (e.info(info_flags::relation_equal)) {
148
149                 // Argument is a relation: convert it to a map
150                 exmap m;
151                 const ex & s = e.op(0);
152                 m.insert(std::make_pair(s, e.op(1)));
153
154                 if (is_exactly_a<mul>(s) || is_exactly_a<power>(s))
155                         options |= subs_options::pattern_is_product;
156                 else
157                         options |= subs_options::pattern_is_not_product;
158
159                 return bp->subs(m, options);
160
161         } else if (e.info(info_flags::list)) {
162
163                 // Argument is a list: convert it to a map
164                 exmap m;
165                 GINAC_ASSERT(is_a<lst>(e));
166                 for (lst::const_iterator it = ex_to<lst>(e).begin(); it != ex_to<lst>(e).end(); ++it) {
167                         ex r = *it;
168                         if (!r.info(info_flags::relation_equal))
169                                 throw(std::invalid_argument("basic::subs(ex): argument must be a list of equations"));
170                         const ex & s = r.op(0);
171                         m.insert(std::make_pair(s, r.op(1)));
172
173                         // Search for products and powers in the expressions to be substituted
174                         // (for an optimization in expairseq::subs())
175                         if (is_exactly_a<mul>(s) || is_exactly_a<power>(s))
176                                 options |= subs_options::pattern_is_product;
177                 }
178                 if (!(options & subs_options::pattern_is_product))
179                         options |= subs_options::pattern_is_not_product;
180
181                 return bp->subs(m, options);
182
183         } else
184                 throw(std::invalid_argument("ex::subs(ex): argument must be a relation_equal or a list"));
185 }
186
187 /** Traverse expression tree with given visitor, preorder traversal. */
188 void ex::traverse_preorder(visitor & v) const
189 {
190         accept(v);
191
192         size_t n = nops();
193         for (size_t i = 0; i < n; ++i)
194                 op(i).traverse_preorder(v);
195 }
196
197 /** Traverse expression tree with given visitor, postorder traversal. */
198 void ex::traverse_postorder(visitor & v) const
199 {
200         size_t n = nops();
201         for (size_t i = 0; i < n; ++i)
202                 op(i).traverse_postorder(v);
203
204         accept(v);
205 }
206
207 /** Return modifyable operand/member at position i. */
208 ex & ex::let_op(size_t i)
209 {
210         makewriteable();
211         return bp->let_op(i);
212 }
213
214 ex & ex::operator[](const ex & index)
215 {
216         makewriteable();
217         return (*bp)[index];
218 }
219
220 ex & ex::operator[](size_t i)
221 {
222         makewriteable();
223         return (*bp)[i];
224 }
225
226 /** Left hand side of relational expression. */
227 ex ex::lhs() const
228 {
229         if (!is_a<relational>(*this))
230                 throw std::runtime_error("ex::lhs(): not a relation");
231         return bp->op(0);
232 }
233
234 /** Right hand side of relational expression. */
235 ex ex::rhs() const
236 {
237         if (!is_a<relational>(*this))
238                 throw std::runtime_error("ex::rhs(): not a relation");
239         return bp->op(1);
240 }
241
242 /** Check whether expression is a polynomial. */
243 bool ex::is_polynomial(const ex & vars) const
244 {
245         if (is_a<lst>(vars)) {
246                 const lst & varlst = ex_to<lst>(vars);
247                 for (lst::const_iterator i=varlst.begin(); i!=varlst.end(); ++i)
248                         if (!bp->is_polynomial(*i))
249                                 return false;
250                 return true;
251         }
252         else
253                 return bp->is_polynomial(vars);
254 }
255
256 /** Check whether expression is zero or zero matrix. */
257 bool ex::is_zero_matrix() const
258 {
259         if (is_zero())
260                 return  true;
261         else {
262                 ex e = evalm();
263                 return is_a<matrix>(e) && ex_to<matrix>(e).is_zero_matrix();
264         }
265 }
266
267 // private
268
269 /** Make this ex writable (if more than one ex handle the same basic) by 
270  *  unlinking the object and creating an unshared copy of it. */
271 void ex::makewriteable()
272 {
273         GINAC_ASSERT(bp->flags & status_flags::dynallocated);
274         bp.makewritable();
275         GINAC_ASSERT(bp->get_refcount() == 1);
276 }
277
278 /** Share equal objects between expressions.
279  *  @see ex::compare(const ex &) */
280 void ex::share(const ex & other) const
281 {
282         if ((bp->flags | other.bp->flags) & status_flags::not_shareable)
283                 return;
284
285         if (bp->get_refcount() <= other.bp->get_refcount())
286                 bp = other.bp;
287         else
288                 other.bp = bp;
289 }
290
291 /** Helper function for the ex-from-basic constructor. This is where GiNaC's
292  *  automatic evaluator and memory management are implemented.
293  *  @see ex::ex(const basic &) */
294 ptr<basic> ex::construct_from_basic(const basic & other)
295 {
296         if (!(other.flags & status_flags::evaluated)) {
297
298                 // The object is not yet evaluated, so call eval() to evaluate
299                 // the top level. This will return either
300                 //  a) the original object with status_flags::evaluated set (when the
301                 //     eval() implementation calls hold())
302                 // or
303                 //  b) a different expression.
304                 //
305                 // eval() returns an ex, not a basic&, so this will go through
306                 // construct_from_basic() a second time. In case a) we end up in
307                 // the "else" branch below. In case b) we end up here again and
308                 // apply eval() once more. The recursion stops when eval() calls
309                 // hold() or returns an object that already has its "evaluated"
310                 // flag set, such as a symbol or a numeric.
311                 const ex & tmpex = other.eval(1);
312
313                 // Eventually, the eval() recursion goes through the "else" branch
314                 // below, which assures that the object pointed to by tmpex.bp is
315                 // allocated on the heap (either it was already on the heap or it
316                 // is a heap-allocated duplicate of another object).
317                 GINAC_ASSERT(tmpex.bp->flags & status_flags::dynallocated); 
318
319                 // If the original object is not referenced but heap-allocated,
320                 // it means that eval() hit case b) above. The original object is
321                 // no longer needed (it evaluated into something different), so we
322                 // delete it (because nobody else will).
323                 if ((other.get_refcount() == 0) && (other.flags & status_flags::dynallocated))
324                         delete &other; // yes, you can apply delete to a const pointer
325
326                 // We can't return a basic& here because the tmpex is destroyed as
327                 // soon as we leave the function, which would deallocate the
328                 // evaluated object.
329                 return tmpex.bp;
330
331         } else {
332
333                 // The easy case: making an "ex" out of an evaluated object.
334                 if (other.flags & status_flags::dynallocated) {
335
336                         // The object is already heap-allocated, so we can just make
337                         // another reference to it.
338                         return ptr<basic>(const_cast<basic &>(other));
339
340                 } else {
341
342                         // The object is not heap-allocated, so we create a duplicate
343                         // on the heap.
344                         basic *bp = other.duplicate();
345                         bp->setflag(status_flags::dynallocated);
346                         GINAC_ASSERT(bp->get_refcount() == 0);
347                         return bp;
348                 }
349         }
350 }
351
352 basic & ex::construct_from_int(int i)
353 {
354         switch (i) {  // prefer flyweights over new objects
355         case -12:
356                 return *const_cast<numeric *>(_num_12_p);
357         case -11:
358                 return *const_cast<numeric *>(_num_11_p);
359         case -10:
360                 return *const_cast<numeric *>(_num_10_p);
361         case -9:
362                 return *const_cast<numeric *>(_num_9_p);
363         case -8:
364                 return *const_cast<numeric *>(_num_8_p);
365         case -7:
366                 return *const_cast<numeric *>(_num_7_p);
367         case -6:
368                 return *const_cast<numeric *>(_num_6_p);
369         case -5:
370                 return *const_cast<numeric *>(_num_5_p);
371         case -4:
372                 return *const_cast<numeric *>(_num_4_p);
373         case -3:
374                 return *const_cast<numeric *>(_num_3_p);
375         case -2:
376                 return *const_cast<numeric *>(_num_2_p);
377         case -1:
378                 return *const_cast<numeric *>(_num_1_p);
379         case 0:
380                 return *const_cast<numeric *>(_num0_p);
381         case 1:
382                 return *const_cast<numeric *>(_num1_p);
383         case 2:
384                 return *const_cast<numeric *>(_num2_p);
385         case 3:
386                 return *const_cast<numeric *>(_num3_p);
387         case 4:
388                 return *const_cast<numeric *>(_num4_p);
389         case 5:
390                 return *const_cast<numeric *>(_num5_p);
391         case 6:
392                 return *const_cast<numeric *>(_num6_p);
393         case 7:
394                 return *const_cast<numeric *>(_num7_p);
395         case 8:
396                 return *const_cast<numeric *>(_num8_p);
397         case 9:
398                 return *const_cast<numeric *>(_num9_p);
399         case 10:
400                 return *const_cast<numeric *>(_num10_p);
401         case 11:
402                 return *const_cast<numeric *>(_num11_p);
403         case 12:
404                 return *const_cast<numeric *>(_num12_p);
405         default:
406                 basic *bp = new numeric(i);
407                 bp->setflag(status_flags::dynallocated);
408                 GINAC_ASSERT(bp->get_refcount() == 0);
409                 return *bp;
410         }
411 }
412         
413 basic & ex::construct_from_uint(unsigned int i)
414 {
415         switch (i) {  // prefer flyweights over new objects
416         case 0:
417                 return *const_cast<numeric *>(_num0_p);
418         case 1:
419                 return *const_cast<numeric *>(_num1_p);
420         case 2:
421                 return *const_cast<numeric *>(_num2_p);
422         case 3:
423                 return *const_cast<numeric *>(_num3_p);
424         case 4:
425                 return *const_cast<numeric *>(_num4_p);
426         case 5:
427                 return *const_cast<numeric *>(_num5_p);
428         case 6:
429                 return *const_cast<numeric *>(_num6_p);
430         case 7:
431                 return *const_cast<numeric *>(_num7_p);
432         case 8:
433                 return *const_cast<numeric *>(_num8_p);
434         case 9:
435                 return *const_cast<numeric *>(_num9_p);
436         case 10:
437                 return *const_cast<numeric *>(_num10_p);
438         case 11:
439                 return *const_cast<numeric *>(_num11_p);
440         case 12:
441                 return *const_cast<numeric *>(_num12_p);
442         default:
443                 basic *bp = new numeric(i);
444                 bp->setflag(status_flags::dynallocated);
445                 GINAC_ASSERT(bp->get_refcount() == 0);
446                 return *bp;
447         }
448 }
449         
450 basic & ex::construct_from_long(long i)
451 {
452         switch (i) {  // prefer flyweights over new objects
453         case -12:
454                 return *const_cast<numeric *>(_num_12_p);
455         case -11:
456                 return *const_cast<numeric *>(_num_11_p);
457         case -10:
458                 return *const_cast<numeric *>(_num_10_p);
459         case -9:
460                 return *const_cast<numeric *>(_num_9_p);
461         case -8:
462                 return *const_cast<numeric *>(_num_8_p);
463         case -7:
464                 return *const_cast<numeric *>(_num_7_p);
465         case -6:
466                 return *const_cast<numeric *>(_num_6_p);
467         case -5:
468                 return *const_cast<numeric *>(_num_5_p);
469         case -4:
470                 return *const_cast<numeric *>(_num_4_p);
471         case -3:
472                 return *const_cast<numeric *>(_num_3_p);
473         case -2:
474                 return *const_cast<numeric *>(_num_2_p);
475         case -1:
476                 return *const_cast<numeric *>(_num_1_p);
477         case 0:
478                 return *const_cast<numeric *>(_num0_p);
479         case 1:
480                 return *const_cast<numeric *>(_num1_p);
481         case 2:
482                 return *const_cast<numeric *>(_num2_p);
483         case 3:
484                 return *const_cast<numeric *>(_num3_p);
485         case 4:
486                 return *const_cast<numeric *>(_num4_p);
487         case 5:
488                 return *const_cast<numeric *>(_num5_p);
489         case 6:
490                 return *const_cast<numeric *>(_num6_p);
491         case 7:
492                 return *const_cast<numeric *>(_num7_p);
493         case 8:
494                 return *const_cast<numeric *>(_num8_p);
495         case 9:
496                 return *const_cast<numeric *>(_num9_p);
497         case 10:
498                 return *const_cast<numeric *>(_num10_p);
499         case 11:
500                 return *const_cast<numeric *>(_num11_p);
501         case 12:
502                 return *const_cast<numeric *>(_num12_p);
503         default:
504                 basic *bp = new numeric(i);
505                 bp->setflag(status_flags::dynallocated);
506                 GINAC_ASSERT(bp->get_refcount() == 0);
507                 return *bp;
508         }
509 }
510         
511 basic & ex::construct_from_ulong(unsigned long i)
512 {
513         switch (i) {  // prefer flyweights over new objects
514         case 0:
515                 return *const_cast<numeric *>(_num0_p);
516         case 1:
517                 return *const_cast<numeric *>(_num1_p);
518         case 2:
519                 return *const_cast<numeric *>(_num2_p);
520         case 3:
521                 return *const_cast<numeric *>(_num3_p);
522         case 4:
523                 return *const_cast<numeric *>(_num4_p);
524         case 5:
525                 return *const_cast<numeric *>(_num5_p);
526         case 6:
527                 return *const_cast<numeric *>(_num6_p);
528         case 7:
529                 return *const_cast<numeric *>(_num7_p);
530         case 8:
531                 return *const_cast<numeric *>(_num8_p);
532         case 9:
533                 return *const_cast<numeric *>(_num9_p);
534         case 10:
535                 return *const_cast<numeric *>(_num10_p);
536         case 11:
537                 return *const_cast<numeric *>(_num11_p);
538         case 12:
539                 return *const_cast<numeric *>(_num12_p);
540         default:
541                 basic *bp = new numeric(i);
542                 bp->setflag(status_flags::dynallocated);
543                 GINAC_ASSERT(bp->get_refcount() == 0);
544                 return *bp;
545         }
546 }
547         
548 basic & ex::construct_from_double(double d)
549 {
550         basic *bp = new numeric(d);
551         bp->setflag(status_flags::dynallocated);
552         GINAC_ASSERT(bp->get_refcount() == 0);
553         return *bp;
554 }
555
556 //////////
557 // static member variables
558 //////////
559
560 // none
561
562 //////////
563 // functions which are not member functions
564 //////////
565
566 // none
567
568 //////////
569 // global functions
570 //////////
571
572 // none
573
574
575 } // namespace GiNaC