]> www.ginac.de Git - ginac.git/blob - ginac/basic.cpp
- in the output, list delimiters are now { } and matrix delimiters are [ ]
[ginac.git] / ginac / basic.cpp
1 /** @file basic.cpp
2  *
3  *  Implementation of GiNaC's ABC. */
4
5 /*
6  *  GiNaC Copyright (C) 1999-2001 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., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
21  */
22
23 #include <iostream>
24 #include <stdexcept>
25 #ifdef DO_GINAC_ASSERT
26 #  include <typeinfo>
27 #endif
28
29 #include "basic.h"
30 #include "ex.h"
31 #include "numeric.h"
32 #include "power.h"
33 #include "symbol.h"
34 #include "lst.h"
35 #include "ncmul.h"
36 #include "relational.h"
37 #include "print.h"
38 #include "archive.h"
39 #include "utils.h"
40 #include "debugmsg.h"
41
42 namespace GiNaC {
43
44 GINAC_IMPLEMENT_REGISTERED_CLASS_NO_CTORS(basic, void)
45
46 //////////
47 // default ctor, dtor, copy ctor assignment operator and helpers
48 //////////
49
50 // public
51
52 basic::basic(const basic & other) : tinfo_key(TINFO_basic), flags(0), refcount(0)
53 {
54         debugmsg("basic copy ctor", LOGLEVEL_CONSTRUCT);
55         copy(other);
56 }
57
58 const basic & basic::operator=(const basic & other)
59 {
60         debugmsg("basic operator=", LOGLEVEL_ASSIGNMENT);
61         if (this != &other) {
62                 destroy(true);
63                 copy(other);
64         }
65         return *this;
66 }
67
68 // protected
69
70 // none (all conditionally inlined)
71
72 //////////
73 // other ctors
74 //////////
75
76 // none (all conditionally inlined)
77
78 //////////
79 // archiving
80 //////////
81
82 /** Construct object from archive_node. */
83 basic::basic(const archive_node &n, const lst &sym_lst) : flags(0), refcount(0)
84 {
85         debugmsg("basic ctor from archive_node", LOGLEVEL_CONSTRUCT);
86
87         // Reconstruct tinfo_key from class name
88         std::string class_name;
89         if (n.find_string("class", class_name))
90                 tinfo_key = find_tinfo_key(class_name);
91         else
92                 throw (std::runtime_error("archive node contains no class name"));
93 }
94
95 /** Unarchive the object. */
96 DEFAULT_UNARCHIVE(basic)
97
98 /** Archive the object. */
99 void basic::archive(archive_node &n) const
100 {
101         n.add_string("class", class_name());
102 }
103
104 //////////
105 // functions overriding virtual functions from bases classes
106 //////////
107
108 // none
109
110 //////////
111 // new virtual functions which can be overridden by derived classes
112 //////////
113
114 // public
115
116 /** Output to stream.
117  *  @param c print context object that describes the output formatting
118  *  @param level value that is used to identify the precedence or indentation
119  *               level for placing parentheses and formatting */
120 void basic::print(const print_context & c, unsigned level) const
121 {
122         debugmsg("basic print", LOGLEVEL_PRINT);
123
124         if (is_of_type(c, print_tree)) {
125
126                 c.s << std::string(level, ' ') << class_name()
127                     << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec
128                     << ", nops=" << nops()
129                     << std::endl;
130                 for (unsigned i=0; i<nops(); ++i)
131                         op(i).print(c, level + static_cast<const print_tree &>(c).delta_indent);
132
133         } else
134                 c.s << "[" << class_name() << " object]";
135 }
136
137 /** Little wrapper arount print to be called within a debugger.
138  *  This is needed because you cannot call foo.print(cout) from within the
139  *  debugger because it might not know what cout is.  This method can be
140  *  invoked with no argument and it will simply print to stdout.
141  *
142  *  @see basic::print */
143 void basic::dbgprint(void) const
144 {
145         this->print(std::cerr);
146         std::cerr << std::endl;
147 }
148
149 /** Little wrapper arount printtree to be called within a debugger.
150  *
151  *  @see basic::dbgprint
152  *  @see basic::printtree */
153 void basic::dbgprinttree(void) const
154 {
155         this->print(print_tree(std::cerr));
156 }
157
158 /** Return relative operator precedence (for parenthizing output). */
159 unsigned basic::precedence(void) const
160 {
161         return 70;
162 }
163
164 /** Create a new copy of this on the heap.  One can think of this as simulating
165  *  a virtual copy constructor which is needed for instance by the refcounted
166  *  construction of an ex from a basic. */
167 basic * basic::duplicate() const
168 {
169         debugmsg("basic duplicate",LOGLEVEL_DUPLICATE);
170         return new basic(*this);
171 }
172
173 /** Information about the object.
174  *
175  *  @see class info_flags */
176 bool basic::info(unsigned inf) const
177 {
178         // all possible properties are false for basic objects
179         return false;
180 }
181
182 /** Number of operands/members. */
183 unsigned basic::nops() const
184 {
185         // iterating from 0 to nops() on atomic objects should be an empty loop,
186         // and accessing their elements is a range error.  Container objects should
187         // override this.
188         return 0;
189 }
190
191 /** Return operand/member at position i. */
192 ex basic::op(int i) const
193 {
194         return (const_cast<basic *>(this))->let_op(i);
195 }
196
197 /** Return modifyable operand/member at position i. */
198 ex & basic::let_op(int i)
199 {
200         throw(std::out_of_range("op() out of range"));
201 }
202
203 ex basic::operator[](const ex & index) const
204 {
205         if (is_exactly_of_type(*index.bp,numeric))
206                 return op(static_cast<const numeric &>(*index.bp).to_int());
207         
208         throw(std::invalid_argument("non-numeric indices not supported by this type"));
209 }
210
211 ex basic::operator[](int i) const
212 {
213         return op(i);
214 }
215
216 /** Search ocurrences.  An object 'has' an expression if it is the expression
217  *  itself or one of the children 'has' it.  As a consequence (according to
218  *  the definition of children) given e=x+y+z, e.has(x) is true but e.has(x+y)
219  *  is false.  The expression can also contain wildcards. */
220 bool basic::has(const ex & other) const
221 {
222         GINAC_ASSERT(other.bp!=0);
223         lst repl_lst;
224         if (match(*other.bp, repl_lst))
225                 return true;
226         for (unsigned i=0; i<nops(); i++)
227                 if (op(i).has(other))
228                         return true;
229         
230         return false;
231 }
232
233 /** Construct new expression by applying the specified function to all
234  *  sub-expressions. */
235 ex basic::map(map_func f) const
236 {
237         unsigned num = nops();
238         if (num == 0)
239                 return *this;
240
241         basic *copy = duplicate();
242         copy->setflag(status_flags::dynallocated);
243         copy->clearflag(status_flags::hash_calculated);
244         ex e(*copy);
245         for (unsigned i=0; i<num; i++)
246                 e.let_op(i) = f(e.op(i));
247         return e.eval();
248 }
249
250 /** Return degree of highest power in object s. */
251 int basic::degree(const ex & s) const
252 {
253         return 0;
254 }
255
256 /** Return degree of lowest power in object s. */
257 int basic::ldegree(const ex & s) const
258 {
259         return 0;
260 }
261
262 /** Return coefficient of degree n in object s. */
263 ex basic::coeff(const ex & s, int n) const
264 {
265         return n==0 ? *this : _ex0();
266 }
267
268 /** Sort expression in terms of powers of some object(s).
269  *  @param s object(s) to sort in
270  *  @param distributed recursive or distributed form (only used when s is a list) */
271 ex basic::collect(const ex & s, bool distributed) const
272 {
273         ex x;
274         if (is_ex_of_type(s, lst)) {
275
276                 // List of objects specified
277                 if (s.nops() == 1)
278                         return collect(s.op(0));
279
280                 else if (distributed) {
281
282                         // Get lower/upper degree of all symbols in list
283                         int num = s.nops();
284                         struct sym_info {
285                                 ex sym;
286                                 int ldeg, deg;
287                                 int cnt;  // current degree, 'counter'
288                                 ex coeff; // coefficient for degree 'cnt'
289                         };
290                         sym_info *si = new sym_info[num];
291                         ex c = *this;
292                         for (int i=0; i<num; i++) {
293                                 si[i].sym = s.op(i);
294                                 si[i].ldeg = si[i].cnt = this->ldegree(si[i].sym);
295                                 si[i].deg = this->degree(si[i].sym);
296                                 c = si[i].coeff = c.coeff(si[i].sym, si[i].cnt);
297                         }
298
299                         while (true) {
300
301                                 // Calculate coeff*x1^c1*...*xn^cn
302                                 ex y = _ex1();
303                                 for (int i=0; i<num; i++) {
304                                         int cnt = si[i].cnt;
305                                         y *= power(si[i].sym, cnt);
306                                 }
307                                 x += y * si[num - 1].coeff;
308
309                                 // Increment counters
310                                 int n = num - 1;
311                                 while (true) {
312                                         si[n].cnt++;
313                                         if (si[n].cnt <= si[n].deg) {
314                                                 // Update coefficients
315                                                 ex c;
316                                                 if (n == 0)
317                                                         c = *this;
318                                                 else
319                                                         c = si[n - 1].coeff;
320                                                 for (int i=n; i<num; i++)
321                                                         c = si[i].coeff = c.coeff(si[i].sym, si[i].cnt);
322                                                 break;
323                                         }
324                                         if (n == 0)
325                                                 goto done;
326                                         si[n].cnt = si[n].ldeg;
327                                         n--;
328                                 }
329                         }
330
331 done:           delete[] si;
332
333                 } else {
334
335                         // Recursive form
336                         x = *this;
337                         for (int n=s.nops()-1; n>=0; n--)
338                                 x = x.collect(s[n]);
339                 }
340
341         } else {
342
343                 // Only one object specified
344                 for (int n=this->ldegree(s); n<=this->degree(s); ++n)
345                         x += this->coeff(s,n)*power(s,n);
346         }
347         
348         // correct for lost fractional arguments and return
349         return x + (*this - x).expand();
350 }
351
352 /** Perform automatic non-interruptive symbolic evaluation on expression. */
353 ex basic::eval(int level) const
354 {
355         // There is nothing to do for basic objects:
356         return this->hold();
357 }
358
359 /** Evaluate object numerically. */
360 ex basic::evalf(int level) const
361 {
362         // There is nothing to do for basic objects:
363         return *this;
364 }
365
366 /** Evaluate sums and products of matrices. */
367 ex basic::evalm(void) const
368 {
369         unsigned num = nops();
370         if (num == 0)
371                 return *this;
372         else
373                 return map(GiNaC::evalm);
374 }
375
376 /** Perform automatic symbolic evaluations on indexed expression that
377  *  contains this object as the base expression. */
378 ex basic::eval_indexed(const basic & i) const
379  // this function can't take a "const ex & i" because that would result
380  // in an infinite eval() loop
381 {
382         // There is nothing to do for basic objects
383         return i.hold();
384 }
385
386 /** Add two indexed expressions. They are guaranteed to be of class indexed
387  *  (or a subclass) and their indices are compatible. This function is used
388  *  internally by simplify_indexed().
389  *
390  *  @param self First indexed expression; it's base object is *this
391  *  @param other Second indexed expression
392  *  @return sum of self and other 
393  *  @see ex::simplify_indexed() */
394 ex basic::add_indexed(const ex & self, const ex & other) const
395 {
396         return self + other;
397 }
398
399 /** Multiply an indexed expression with a scalar. This function is used
400  *  internally by simplify_indexed().
401  *
402  *  @param self Indexed expression; it's base object is *this
403  *  @param other Numeric value
404  *  @return product of self and other
405  *  @see ex::simplify_indexed() */
406 ex basic::scalar_mul_indexed(const ex & self, const numeric & other) const
407 {
408         return self * other;
409 }
410
411 /** Try to contract two indexed expressions that appear in the same product. 
412  *  If a contraction exists, the function overwrites one or both of the
413  *  expressions and returns true. Otherwise it returns false. It is
414  *  guaranteed that both expressions are of class indexed (or a subclass)
415  *  and that at least one dummy index has been found. This functions is
416  *  used internally by simplify_indexed().
417  *
418  *  @param self Pointer to first indexed expression; it's base object is *this
419  *  @param other Pointer to second indexed expression
420  *  @param v The complete vector of factors
421  *  @return true if the contraction was successful, false otherwise
422  *  @see ex::simplify_indexed() */
423 bool basic::contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const
424 {
425         // Do nothing
426         return false;
427 }
428
429 /** Check whether the expression matches a given pattern. For every wildcard
430  *  object in the pattern, an expression of the form "wildcard == matching_expression"
431  *  is added to repl_lst. */
432 bool basic::match(const ex & pattern, lst & repl_lst) const
433 {
434 /*
435         Sweet sweet shapes, sweet sweet shapes,
436         Thats the key thing, right right.
437         Feed feed face, feed feed shapes,
438         But who is the king tonight?
439         Who is the king tonight?
440         Pattern is the thing, the key thing-a-ling,
441         But who is the king of pattern?
442         But who is the king, the king thing-a-ling,
443         Who is the king of Pattern?
444         Bog is the king, the king thing-a-ling,
445         Bog is the king of Pattern.
446         Ba bu-bu-bu-bu bu-bu-bu-bu-bu-bu bu-bu
447         Bog is the king of Pattern.
448 */
449
450         if (is_ex_exactly_of_type(pattern, wildcard)) {
451
452                 // Wildcard matches anything, but check whether we already have found
453                 // a match for that wildcard first (if so, it the earlier match must
454                 // be the same expression)
455                 for (unsigned i=0; i<repl_lst.nops(); i++) {
456                         if (repl_lst.op(i).op(0).is_equal(pattern))
457                                 return is_equal(*repl_lst.op(i).op(1).bp);
458                 }
459                 repl_lst.append(pattern == *this);
460                 return true;
461
462         } else {
463
464                 // Expression must be of the same type as the pattern
465                 if (tinfo() != pattern.bp->tinfo())
466                         return false;
467
468                 // Number of subexpressions must match
469                 if (nops() != pattern.nops())
470                         return false;
471
472                 // No subexpressions? Then just compare the objects (there can't be
473                 // wildcards in the pattern)
474                 if (nops() == 0)
475                         return is_equal(*pattern.bp);
476
477                 // Otherwise the subexpressions must match one-to-one
478                 for (unsigned i=0; i<nops(); i++)
479                         if (!op(i).match(pattern.op(i), repl_lst))
480                                 return false;
481
482                 // Looks similar enough, match found
483                 return true;
484         }
485 }
486
487 /** Substitute a set of objects by arbitrary expressions. The ex returned
488  *  will already be evaluated. */
489 ex basic::subs(const lst & ls, const lst & lr, bool no_pattern) const
490 {
491         GINAC_ASSERT(ls.nops() == lr.nops());
492
493         if (no_pattern) {
494                 for (unsigned i=0; i<ls.nops(); i++) {
495                         if (is_equal(*ls.op(i).bp))
496                                 return lr.op(i);
497                 }
498         } else {
499                 for (unsigned i=0; i<ls.nops(); i++) {
500                         lst repl_lst;
501                         if (match(*ls.op(i).bp, repl_lst))
502                                 return lr.op(i).bp->subs(repl_lst, true); // avoid infinite recursion when re-substituting the wildcards
503                 }
504         }
505
506         return *this;
507 }
508
509 /** Default interface of nth derivative ex::diff(s, n).  It should be called
510  *  instead of ::derivative(s) for first derivatives and for nth derivatives it
511  *  just recurses down.
512  *
513  *  @param s symbol to differentiate in
514  *  @param nth order of differentiation
515  *  @see ex::diff */
516 ex basic::diff(const symbol & s, unsigned nth) const
517 {
518         // trivial: zeroth derivative
519         if (nth==0)
520                 return ex(*this);
521         
522         // evaluate unevaluated *this before differentiating
523         if (!(flags & status_flags::evaluated))
524                 return ex(*this).diff(s, nth);
525         
526         ex ndiff = this->derivative(s);
527         while (!ndiff.is_zero() &&    // stop differentiating zeros
528                nth>1) {
529                 ndiff = ndiff.diff(s);
530                 --nth;
531         }
532         return ndiff;
533 }
534
535 /** Return a vector containing the free indices of an expression. */
536 exvector basic::get_free_indices(void) const
537 {
538         return exvector(); // return an empty exvector
539 }
540
541 ex basic::simplify_ncmul(const exvector & v) const
542 {
543         return simplified_ncmul(v);
544 }
545
546 // protected
547
548 /** Default implementation of ex::diff(). It simply throws an error message.
549  *
550  *  @exception logic_error (differentiation not supported by this type)
551  *  @see ex::diff */
552 ex basic::derivative(const symbol & s) const
553 {
554         throw(std::logic_error("differentiation not supported by this type"));
555 }
556
557 /** Returns order relation between two objects of same type.  This needs to be
558  *  implemented by each class. It may never return anything else than 0,
559  *  signalling equality, or +1 and -1 signalling inequality and determining
560  *  the canonical ordering.  (Perl hackers will wonder why C++ doesn't feature
561  *  the spaceship operator <=> for denoting just this.) */
562 int basic::compare_same_type(const basic & other) const
563 {
564         return compare_pointers(this, &other);
565 }
566
567 /** Returns true if two objects of same type are equal.  Normally needs
568  *  not be reimplemented as long as it wasn't overwritten by some parent
569  *  class, since it just calls compare_same_type().  The reason why this
570  *  function exists is that sometimes it is easier to determine equality
571  *  than an order relation and then it can be overridden. */
572 bool basic::is_equal_same_type(const basic & other) const
573 {
574         return this->compare_same_type(other)==0;
575 }
576
577 unsigned basic::return_type(void) const
578 {
579         return return_types::commutative;
580 }
581
582 unsigned basic::return_type_tinfo(void) const
583 {
584         return tinfo();
585 }
586
587 /** Compute the hash value of an object and if it makes sense to store it in
588  *  the objects status_flags, do so.  The method inherited from class basic
589  *  computes a hash value based on the type and hash values of possible
590  *  members.  For this reason it is well suited for container classes but
591  *  atomic classes should override this implementation because otherwise they
592  *  would all end up with the same hashvalue. */
593 unsigned basic::calchash(void) const
594 {
595         unsigned v = golden_ratio_hash(tinfo());
596         for (unsigned i=0; i<nops(); i++) {
597                 v = rotate_left_31(v);
598                 v ^= (const_cast<basic *>(this))->op(i).gethash();
599         }
600         
601         // mask out numeric hashes:
602         v &= 0x7FFFFFFFU;
603         
604         // store calculated hash value only if object is already evaluated
605         if (flags & status_flags::evaluated) {
606                 setflag(status_flags::hash_calculated);
607                 hashvalue = v;
608         }
609
610         return v;
611 }
612
613 /** Expand expression, i.e. multiply it out and return the result as a new
614  *  expression. */
615 ex basic::expand(unsigned options) const
616 {
617         return this->setflag(status_flags::expanded);
618 }
619
620
621 //////////
622 // non-virtual functions in this class
623 //////////
624
625 // public
626
627 /** Substitute objects in an expression (syntactic substitution) and return
628  *  the result as a new expression.  There are two valid types of
629  *  replacement arguments: 1) a relational like object==ex and 2) a list of
630  *  relationals lst(object1==ex1,object2==ex2,...), which is converted to
631  *  subs(lst(object1,object2,...),lst(ex1,ex2,...)). */
632 ex basic::subs(const ex & e, bool no_pattern) const
633 {
634         if (e.info(info_flags::relation_equal)) {
635                 return subs(lst(e), no_pattern);
636         }
637         if (!e.info(info_flags::list)) {
638                 throw(std::invalid_argument("basic::subs(ex): argument must be a list"));
639         }
640         lst ls;
641         lst lr;
642         for (unsigned i=0; i<e.nops(); i++) {
643                 ex r = e.op(i);
644                 if (!r.info(info_flags::relation_equal)) {
645                         throw(std::invalid_argument("basic::subs(ex): argument must be a list of equations"));
646                 }
647                 ls.append(r.op(0));
648                 lr.append(r.op(1));
649         }
650         return subs(ls, lr, no_pattern);
651 }
652
653 /** Compare objects to establish canonical ordering.
654  *  All compare functions return: -1 for *this less than other, 0 equal,
655  *  1 greater. */
656 int basic::compare(const basic & other) const
657 {
658         unsigned hash_this = gethash();
659         unsigned hash_other = other.gethash();
660         
661         if (hash_this<hash_other) return -1;
662         if (hash_this>hash_other) return 1;
663         
664         unsigned typeid_this = tinfo();
665         unsigned typeid_other = other.tinfo();
666         
667         if (typeid_this<typeid_other) {
668 //              std::cout << "hash collision, different types: " 
669 //                        << *this << " and " << other << std::endl;
670 //              this->print(print_tree(std::cout));
671 //              std::cout << " and ";
672 //              other.print(print_tree(std::cout));
673 //              std::cout << std::endl;
674                 return -1;
675         }
676         if (typeid_this>typeid_other) {
677 //              std::cout << "hash collision, different types: " 
678 //                        << *this << " and " << other << std::endl;
679 //              this->print(print_tree(std::cout));
680 //              std::cout << " and ";
681 //              other.print(print_tree(std::cout));
682 //              std::cout << std::endl;
683                 return 1;
684         }
685         
686         GINAC_ASSERT(typeid(*this)==typeid(other));
687         
688 //      int cmpval = compare_same_type(other);
689 //      if ((cmpval!=0) && (hash_this<0x80000000U)) {
690 //              std::cout << "hash collision, same type: " 
691 //                        << *this << " and " << other << std::endl;
692 //              this->print(print_tree(std::cout));
693 //              std::cout << " and ";
694 //              other.print(print_tree(std::cout));
695 //              std::cout << std::endl;
696 //      }
697 //      return cmpval;
698         
699         return compare_same_type(other);
700 }
701
702 /** Test for equality.
703  *  This is only a quick test, meaning objects should be in the same domain.
704  *  You might have to .expand(), .normal() objects first, depending on the
705  *  domain of your computation, to get a more reliable answer.
706  *
707  *  @see is_equal_same_type */
708 bool basic::is_equal(const basic & other) const
709 {
710         if (this->gethash()!=other.gethash())
711                 return false;
712         if (this->tinfo()!=other.tinfo())
713                 return false;
714         
715         GINAC_ASSERT(typeid(*this)==typeid(other));
716         
717         return this->is_equal_same_type(other);
718 }
719
720 // protected
721
722 /** Stop further evaluation.
723  *
724  *  @see basic::eval */
725 const basic & basic::hold(void) const
726 {
727         return this->setflag(status_flags::evaluated);
728 }
729
730 /** Ensure the object may be modified without hurting others, throws if this
731  *  is not the case. */
732 void basic::ensure_if_modifiable(void) const
733 {
734         if (this->refcount>1)
735                 throw(std::runtime_error("cannot modify multiply referenced object"));
736 }
737
738 //////////
739 // global variables
740 //////////
741
742 int max_recursion_level = 1024;
743
744 } // namespace GiNaC