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