ginac/add.cpp

   1 /** @file add.cpp
   2  *
   3  *  Implementation of GiNaC's sums of expressions. */
   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
  26 #include "add.h"
  27 #include "mul.h"
  28 #include "archive.h"
  29 #include "debugmsg.h"
  30 #include "utils.h"
  31
  32 #ifndef NO_NAMESPACE_GINAC
  33 namespace GiNaC {
  34 #endif // ndef NO_NAMESPACE_GINAC
  35
  36 GINAC_IMPLEMENT_REGISTERED_CLASS(add, expairseq)
  37
  38 //////////
  39 // default constructor, destructor, copy constructor assignment operator and helpers
  40 //////////
  41
  42 // public
  43
  44 add::add()
  45 {
  46         debugmsg("add default constructor",LOGLEVEL_CONSTRUCT);
  47         tinfo_key = TINFO_add;
  48 }
  49
  50 // protected
  51
  52 void add::copy(const add & other)
  53 {
  54         inherited::copy(other);
  55 }
  56
  57 void add::destroy(bool call_parent)
  58 {
  59         if (call_parent) inherited::destroy(call_parent);
  60 }
  61
  62 //////////
  63 // other constructors
  64 //////////
  65
  66 // public
  67
  68 add::add(const ex & lh, const ex & rh)
  69 {
  70         debugmsg("add constructor from ex,ex",LOGLEVEL_CONSTRUCT);
  71         tinfo_key = TINFO_add;
  72         overall_coeff = _ex0();
  73         construct_from_2_ex(lh,rh);
  74         GINAC_ASSERT(is_canonical());
  75 }
  76
  77 add::add(const exvector & v)
  78 {
  79         debugmsg("add constructor from exvector",LOGLEVEL_CONSTRUCT);
  80         tinfo_key = TINFO_add;
  81         overall_coeff = _ex0();
  82         construct_from_exvector(v);
  83         GINAC_ASSERT(is_canonical());
  84 }
  85
  86 /*
  87 add::add(const epvector & v, bool do_not_canonicalize)
  88 {
  89         debugmsg("add constructor from epvector,bool",LOGLEVEL_CONSTRUCT);
  90         tinfo_key = TINFO_add;
  91         if (do_not_canonicalize) {
  92                 seq=v;
  93 #ifdef EXPAIRSEQ_USE_HASHTAB
  94                 combine_same_terms(); // to build hashtab
  95 #endif // def EXPAIRSEQ_USE_HASHTAB
  96         } else {
  97                 construct_from_epvector(v);
  98         }
  99         GINAC_ASSERT(is_canonical());
 100 }
 101 */
 102
 103 add::add(const epvector & v)
 104 {
 105         debugmsg("add constructor from epvector",LOGLEVEL_CONSTRUCT);
 106         tinfo_key = TINFO_add;
 107         overall_coeff = _ex0();
 108         construct_from_epvector(v);
 109         GINAC_ASSERT(is_canonical());
 110 }
 111
 112 add::add(const epvector & v, const ex & oc)
 113 {
 114         debugmsg("add constructor from epvector,ex",LOGLEVEL_CONSTRUCT);
 115         tinfo_key = TINFO_add;
 116         overall_coeff = oc;
 117         construct_from_epvector(v);
 118         GINAC_ASSERT(is_canonical());
 119 }
 120
 121 add::add(epvector * vp, const ex & oc)
 122 {
 123         debugmsg("add constructor from epvector *,ex",LOGLEVEL_CONSTRUCT);
 124         tinfo_key = TINFO_add;
 125         GINAC_ASSERT(vp!=0);
 126         overall_coeff = oc;
 127         construct_from_epvector(*vp);
 128         delete vp;
 129         GINAC_ASSERT(is_canonical());
 130 }
 131
 132 //////////
 133 // archiving
 134 //////////
 135
 136 /** Construct object from archive_node. */
 137 add::add(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst)
 138 {
 139         debugmsg("add constructor from archive_node", LOGLEVEL_CONSTRUCT);
 140 }
 141
 142 /** Unarchive the object. */
 143 ex add::unarchive(const archive_node &n, const lst &sym_lst)
 144 {
 145         return (new add(n, sym_lst))->setflag(status_flags::dynallocated);
 146 }
 147
 148 /** Archive the object. */
 149 void add::archive(archive_node &n) const
 150 {
 151         inherited::archive(n);
 152 }
 153
 154 //////////
 155 // functions overriding virtual functions from bases classes
 156 //////////
 157
 158 // public
 159
 160 basic * add::duplicate() const
 161 {
 162         debugmsg("add duplicate",LOGLEVEL_DUPLICATE);
 163         return new add(*this);
 164 }
 165
 166 void add::print(std::ostream & os, unsigned upper_precedence) const
 167 {
 168         debugmsg("add print",LOGLEVEL_PRINT);
 169         if (precedence<=upper_precedence) os << "(";
 170         numeric coeff;
 171         bool first = true;
 172         // first print the overall numeric coefficient, if present:
 173         if (!overall_coeff.is_zero()) {
 174                 os << overall_coeff;
 175                 first = false;
 176         }
 177         // then proceed with the remaining factors:
 178         for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
 179                 coeff = ex_to_numeric(cit->coeff);
 180                 if (!first) {
 181                         if (coeff.csgn()==-1) os << '-'; else os << '+';
 182                 } else {
 183                         if (coeff.csgn()==-1) os << '-';
 184                         first = false;
 185                 }
 186                 if (!coeff.is_equal(_num1()) &&
 187                     !coeff.is_equal(_num_1())) {
 188                         if (coeff.is_rational()) {
 189                                 if (coeff.is_negative())
 190                                         os << -coeff;
 191                                 else
 192                                         os << coeff;
 193                         } else {
 194                                 if (coeff.csgn()==-1)
 195                                         (-coeff).print(os, precedence);
 196                                 else
 197                                         coeff.print(os, precedence);
 198                         }
 199                         os << '*';
 200                 }
 201                 cit->rest.print(os, precedence);
 202         }
 203         if (precedence<=upper_precedence) os << ")";
 204 }
 205
 206 void add::printraw(std::ostream & os) const
 207 {
 208         debugmsg("add printraw",LOGLEVEL_PRINT);
 209
 210         os << "+(";
 211         for (epvector::const_iterator it=seq.begin(); it!=seq.end(); ++it) {
 212                 os << "(";
 213                 (*it).rest.bp->printraw(os);
 214                 os << ",";
 215                 (*it).coeff.bp->printraw(os);
 216                 os << "),";
 217         }
 218         os << ",hash=" << hashvalue << ",flags=" << flags;
 219         os << ")";
 220 }
 221
 222 void add::printcsrc(std::ostream & os, unsigned type, unsigned upper_precedence) const
 223 {
 224         debugmsg("add print csrc", LOGLEVEL_PRINT);
 225         if (precedence <= upper_precedence)
 226                 os << "(";
 227
 228         // Print arguments, separated by "+"
 229         epvector::const_iterator it = seq.begin();
 230         epvector::const_iterator itend = seq.end();
 231         while (it != itend) {
 232
 233                 // If the coefficient is -1, it is replaced by a single minus sign
 234                 if (it->coeff.compare(_num1()) == 0) {
 235                         it->rest.bp->printcsrc(os, type, precedence);
 236                 } else if (it->coeff.compare(_num_1()) == 0) {
 237                         os << "-";
 238                         it->rest.bp->printcsrc(os, type, precedence);
 239                 } else if (ex_to_numeric(it->coeff).numer().compare(_num1()) == 0) {
 240                         it->rest.bp->printcsrc(os, type, precedence);
 241                         os << "/";
 242                         ex_to_numeric(it->coeff).denom().printcsrc(os, type, precedence);
 243                 } else if (ex_to_numeric(it->coeff).numer().compare(_num_1()) == 0) {
 244                         os << "-";
 245                         it->rest.bp->printcsrc(os, type, precedence);
 246                         os << "/";
 247                         ex_to_numeric(it->coeff).denom().printcsrc(os, type, precedence);
 248                 } else {
 249                         it->coeff.bp->printcsrc(os, type, precedence);
 250                         os << "*";
 251                         it->rest.bp->printcsrc(os, type, precedence);
 252                 }
 253
 254                 // Separator is "+", except if the following expression would have a leading minus sign
 255                 it++;
 256                 if (it != itend && !(it->coeff.compare(_num0()) < 0 || (it->coeff.compare(_num1()) == 0 && is_ex_exactly_of_type(it->rest, numeric) && it->rest.compare(_num0()) < 0)))
 257                         os << "+";
 258         }
 259
 260         if (!overall_coeff.is_equal(_ex0())) {
 261                 if (overall_coeff.info(info_flags::positive)) os << '+';
 262                 overall_coeff.bp->printcsrc(os,type,precedence);
 263         }
 264
 265         if (precedence <= upper_precedence)
 266                 os << ")";
 267 }
 268
 269 bool add::info(unsigned inf) const
 270 {
 271         switch (inf) {
 272                 case info_flags::polynomial:
 273                 case info_flags::integer_polynomial:
 274                 case info_flags::cinteger_polynomial:
 275                 case info_flags::rational_polynomial:
 276                 case info_flags::crational_polynomial:
 277                 case info_flags::rational_function: {
 278                         for (epvector::const_iterator i=seq.begin(); i!=seq.end(); ++i) {
 279                                 if (!(recombine_pair_to_ex(*i).info(inf)))
 280                                         return false;
 281                         }
 282                         return overall_coeff.info(inf);
 283                 }
 284                 case info_flags::algebraic: {
 285                         for (epvector::const_iterator i=seq.begin(); i!=seq.end(); ++i) {
 286                                 if ((recombine_pair_to_ex(*i).info(inf)))
 287                                         return true;
 288                         }
 289                         return false;
 290                 }
 291         }
 292         return inherited::info(inf);
 293 }
 294
 295 int add::degree(const symbol & s) const
 296 {
 297         int deg = INT_MIN;
 298         if (!overall_coeff.is_equal(_ex0())) {
 299                 deg = 0;
 300         }
 301         int cur_deg;
 302         for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
 303                 cur_deg=(*cit).rest.degree(s);
 304                 if (cur_deg>deg) deg=cur_deg;
 305         }
 306         return deg;
 307 }
 308
 309 int add::ldegree(const symbol & s) const
 310 {
 311         int deg = INT_MAX;
 312         if (!overall_coeff.is_equal(_ex0())) {
 313                 deg = 0;
 314         }
 315         int cur_deg;
 316         for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
 317                 cur_deg = (*cit).rest.ldegree(s);
 318                 if (cur_deg<deg) deg=cur_deg;
 319         }
 320         return deg;
 321 }
 322
 323 ex add::coeff(const symbol & s, int n) const
 324 {
 325         epvector coeffseq;
 326         coeffseq.reserve(seq.size());
 327
 328         epvector::const_iterator it=seq.begin();
 329         while (it!=seq.end()) {
 330                 coeffseq.push_back(combine_ex_with_coeff_to_pair((*it).rest.coeff(s,n),
 331                                                                  (*it).coeff));
 332                 ++it;
 333         }
 334         if (n==0) {
 335                 return (new add(coeffseq,overall_coeff))->setflag(status_flags::dynallocated);
 336         }
 337         return (new add(coeffseq))->setflag(status_flags::dynallocated);
 338 }
 339
 340 ex add::eval(int level) const
 341 {
 342         // simplifications: +(;c) -> c
 343         //                  +(x;1) -> x
 344
 345         debugmsg("add eval",LOGLEVEL_MEMBER_FUNCTION);
 346
 347         epvector * evaled_seqp=evalchildren(level);
 348         if (evaled_seqp!=0) {
 349                 // do more evaluation later
 350                 return (new add(evaled_seqp,overall_coeff))->
 351                        setflag(status_flags::dynallocated);
 352         }
 353
 354 #ifdef DO_GINAC_ASSERT
 355         for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
 356                 GINAC_ASSERT(!is_ex_exactly_of_type((*cit).rest,add));
 357                 if (is_ex_exactly_of_type((*cit).rest,numeric)) {
 358                         dbgprint();
 359                 }
 360                 GINAC_ASSERT(!is_ex_exactly_of_type((*cit).rest,numeric));
 361         }
 362 #endif // def DO_GINAC_ASSERT
 363
 364         if (flags & status_flags::evaluated) {
 365                 GINAC_ASSERT(seq.size()>0);
 366                 GINAC_ASSERT((seq.size()>1)||!overall_coeff.is_equal(_ex0()));
 367                 return *this;
 368         }
 369
 370         int seq_size=seq.size();
 371         if (seq_size==0) {
 372                 // +(;c) -> c
 373                 return overall_coeff;
 374         } else if ((seq_size==1)&&overall_coeff.is_equal(_ex0())) {
 375                 // +(x;0) -> x
 376                 return recombine_pair_to_ex(*(seq.begin()));
 377         }
 378         return this->hold();
 379 }
 380
 381 exvector add::get_indices(void) const
 382 {
 383         // FIXME: all terms in the sum should have the same indices (compatible
 384         // tensors) however this is not checked, since there is no function yet
 385         // which compares indices (idxvector can be unsorted)
 386         if (seq.size()==0) {
 387                 return exvector();
 388         }
 389         return (seq.begin())->rest.get_indices();
 390 }
 391
 392 ex add::simplify_ncmul(const exvector & v) const
 393 {
 394         if (seq.size()==0) {
 395                 return inherited::simplify_ncmul(v);
 396         }
 397         return (*seq.begin()).rest.simplify_ncmul(v);
 398 }
 399
 400 // protected
 401
 402 /** Implementation of ex::diff() for a sum. It differentiates each term.
 403  *  @see ex::diff */
 404 ex add::derivative(const symbol & s) const
 405 {
 406         // D(a+b+c)=D(a)+D(b)+D(c)
 407         return (new add(diffchildren(s)))->setflag(status_flags::dynallocated);
 408 }
 409
 410 int add::compare_same_type(const basic & other) const
 411 {
 412         return inherited::compare_same_type(other);
 413 }
 414
 415 bool add::is_equal_same_type(const basic & other) const
 416 {
 417         return inherited::is_equal_same_type(other);
 418 }
 419
 420 unsigned add::return_type(void) const
 421 {
 422         if (seq.size()==0) {
 423                 return return_types::commutative;
 424         }
 425         return (*seq.begin()).rest.return_type();
 426 }
 427
 428 unsigned add::return_type_tinfo(void) const
 429 {
 430         if (seq.size()==0) {
 431                 return tinfo_key;
 432         }
 433         return (*seq.begin()).rest.return_type_tinfo();
 434 }
 435
 436 ex add::thisexpairseq(const epvector & v, const ex & oc) const
 437 {
 438         return (new add(v,oc))->setflag(status_flags::dynallocated);
 439 }
 440
 441 ex add::thisexpairseq(epvector * vp, const ex & oc) const
 442 {
 443         return (new add(vp,oc))->setflag(status_flags::dynallocated);
 444 }
 445
 446 expair add::split_ex_to_pair(const ex & e) const
 447 {
 448         if (is_ex_exactly_of_type(e,mul)) {
 449                 const mul & mulref=ex_to_mul(e);
 450                 ex numfactor=mulref.overall_coeff;
 451                 // mul * mulcopyp=static_cast<mul *>(mulref.duplicate());
 452                 mul * mulcopyp=new mul(mulref);
 453                 mulcopyp->overall_coeff=_ex1();
 454                 mulcopyp->clearflag(status_flags::evaluated);
 455                 mulcopyp->clearflag(status_flags::hash_calculated);
 456                 return expair(mulcopyp->setflag(status_flags::dynallocated),numfactor);
 457         }
 458         return expair(e,_ex1());
 459 }
 460
 461 expair add::combine_ex_with_coeff_to_pair(const ex & e,
 462                                                                                   const ex & c) const
 463 {
 464         GINAC_ASSERT(is_ex_exactly_of_type(c, numeric));
 465         ex one = _ex1();
 466         if (is_ex_exactly_of_type(e, mul)) {
 467                 const mul &mulref = ex_to_mul(e);
 468                 ex numfactor = mulref.overall_coeff;
 469                 mul *mulcopyp = new mul(mulref);
 470                 mulcopyp->overall_coeff = one;
 471                 mulcopyp->clearflag(status_flags::evaluated);
 472                 mulcopyp->clearflag(status_flags::hash_calculated);
 473                 mulcopyp->setflag(status_flags::dynallocated);
 474                 if (are_ex_trivially_equal(c, one)) {
 475                         return expair(*mulcopyp, numfactor);
 476                 } else if (are_ex_trivially_equal(numfactor, one)) {
 477                         return expair(*mulcopyp, c);
 478                 }
 479                 return expair(*mulcopyp, ex_to_numeric(numfactor).mul_dyn(ex_to_numeric(c)));
 480         } else if (is_ex_exactly_of_type(e, numeric)) {
 481                 if (are_ex_trivially_equal(c, one)) {
 482                         return expair(e, one);
 483                 }
 484                 return expair(ex_to_numeric(e).mul_dyn(ex_to_numeric(c)), one);
 485         }
 486         return expair(e, c);
 487 }
 488
 489 expair add::combine_pair_with_coeff_to_pair(const expair & p,
 490                                                                                         const ex & c) const
 491 {
 492         GINAC_ASSERT(is_ex_exactly_of_type(p.coeff,numeric));
 493         GINAC_ASSERT(is_ex_exactly_of_type(c,numeric));
 494
 495         if (is_ex_exactly_of_type(p.rest,numeric)) {
 496                 GINAC_ASSERT(ex_to_numeric(p.coeff).is_equal(_num1())); // should be normalized
 497                 return expair(ex_to_numeric(p.rest).mul_dyn(ex_to_numeric(c)),_ex1());
 498         }
 499
 500         return expair(p.rest,ex_to_numeric(p.coeff).mul_dyn(ex_to_numeric(c)));
 501 }
 502
 503 ex add::recombine_pair_to_ex(const expair & p) const
 504 {
 505         if (ex_to_numeric(p.coeff).is_equal(_num1()))
 506                 return p.rest;
 507         else
 508                 return p.rest*p.coeff;
 509 }
 510
 511 ex add::expand(unsigned options) const
 512 {
 513         if (flags & status_flags::expanded)
 514                 return *this;
 515
 516         epvector * vp = expandchildren(options);
 517         if (vp==0)
 518                 return *this;
 519
 520         return (new add(vp,overall_coeff))->setflag(status_flags::expanded | status_flags::dynallocated);
 521 }
 522
 523 //////////
 524 // static member variables
 525 //////////
 526
 527 // protected
 528
 529 unsigned add::precedence = 40;
 530
 531 #ifndef NO_NAMESPACE_GINAC
 532 } // namespace GiNaC
 533 #endif // ndef NO_NAMESPACE_GINAC