ginac/pseries.cpp

   1 /** @file pseries.cpp
   2  *
   3  *  Implementation of class for extended truncated power series and
   4  *  methods for series expansion. */
   5
   6 /*
   7  *  GiNaC Copyright (C) 1999-2001 Johannes Gutenberg University Mainz, Germany
   8  *
   9  *  This program is free software; you can redistribute it and/or modify
  10  *  it under the terms of the GNU General Public License as published by
  11  *  the Free Software Foundation; either version 2 of the License, or
  12  *  (at your option) any later version.
  13  *
  14  *  This program is distributed in the hope that it will be useful,
  15  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  17  *  GNU General Public License for more details.
  18  *
  19  *  You should have received a copy of the GNU General Public License
  20  *  along with this program; if not, write to the Free Software
  21  *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  22  */
  23
  24 #include <stdexcept>
  25
  26 #include "pseries.h"
  27 #include "add.h"
  28 #include "inifcns.h"
  29 #include "lst.h"
  30 #include "mul.h"
  31 #include "power.h"
  32 #include "relational.h"
  33 #include "symbol.h"
  34 #include "archive.h"
  35 #include "utils.h"
  36 #include "debugmsg.h"
  37
  38 namespace GiNaC {
  39
  40 GINAC_IMPLEMENT_REGISTERED_CLASS(pseries, basic)
  41
  42 /*
  43  *  Default ctor, dtor, copy ctor, assignment operator and helpers
  44  */
  45
  46 pseries::pseries() : basic(TINFO_pseries)
  47 {
  48         debugmsg("pseries default ctor", LOGLEVEL_CONSTRUCT);
  49 }
  50
  51 void pseries::copy(const pseries &other)
  52 {
  53         inherited::copy(other);
  54         seq = other.seq;
  55         var = other.var;
  56         point = other.point;
  57 }
  58
  59 void pseries::destroy(bool call_parent)
  60 {
  61         if (call_parent)
  62                 inherited::destroy(call_parent);
  63 }
  64
  65
  66 /*
  67  *  Other ctors
  68  */
  69
  70 /** Construct pseries from a vector of coefficients and powers.
  71  *  expair.rest holds the coefficient, expair.coeff holds the power.
  72  *  The powers must be integers (positive or negative) and in ascending order;
  73  *  the last coefficient can be Order(_ex1()) to represent a truncated,
  74  *  non-terminating series.
  75  *
  76  *  @param rel_  expansion variable and point (must hold a relational)
  77  *  @param ops_  vector of {coefficient, power} pairs (coefficient must not be zero)
  78  *  @return newly constructed pseries */
  79 pseries::pseries(const ex &rel_, const epvector &ops_) : basic(TINFO_pseries), seq(ops_)
  80 {
  81         debugmsg("pseries ctor from ex,epvector", LOGLEVEL_CONSTRUCT);
  82         GINAC_ASSERT(is_ex_exactly_of_type(rel_, relational));
  83         GINAC_ASSERT(is_ex_exactly_of_type(rel_.lhs(),symbol));
  84         point = rel_.rhs();
  85         var = *static_cast<symbol *>(rel_.lhs().bp);
  86 }
  87
  88
  89 /*
  90  *  Archiving
  91  */
  92
  93 /** Construct object from archive_node. */
  94 pseries::pseries(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst)
  95 {
  96         debugmsg("pseries ctor from archive_node", LOGLEVEL_CONSTRUCT);
  97         for (unsigned int i=0; true; ++i) {
  98                 ex rest;
  99                 ex coeff;
 100                 if (n.find_ex("coeff", rest, sym_lst, i) && n.find_ex("power", coeff, sym_lst, i))
 101                         seq.push_back(expair(rest, coeff));
 102                 else
 103                         break;
 104         }
 105         n.find_ex("var", var, sym_lst);
 106         n.find_ex("point", point, sym_lst);
 107 }
 108
 109 /** Unarchive the object. */
 110 ex pseries::unarchive(const archive_node &n, const lst &sym_lst)
 111 {
 112         return (new pseries(n, sym_lst))->setflag(status_flags::dynallocated);
 113 }
 114
 115 /** Archive the object. */
 116 void pseries::archive(archive_node &n) const
 117 {
 118         inherited::archive(n);
 119         epvector::const_iterator i = seq.begin(), iend = seq.end();
 120         while (i != iend) {
 121                 n.add_ex("coeff", i->rest);
 122                 n.add_ex("power", i->coeff);
 123                 ++i;
 124         }
 125         n.add_ex("var", var);
 126         n.add_ex("point", point);
 127 }
 128
 129 //////////
 130 // functions overriding virtual functions from bases classes
 131 //////////
 132
 133 void pseries::print(std::ostream &os, unsigned upper_precedence) const
 134 {
 135         debugmsg("pseries print", LOGLEVEL_PRINT);
 136         if (precedence<=upper_precedence) os << "(";
 137         // objects of type pseries must not have any zero entries, so the
 138         // trivial (zero) pseries needs a special treatment here:
 139         if (seq.size()==0)
 140                 os << '0';
 141         for (epvector::const_iterator i=seq.begin(); i!=seq.end(); ++i) {
 142                 // print a sign, if needed
 143                 if (i!=seq.begin())
 144                         os << '+';
 145                 if (!is_order_function(i->rest)) {
 146                         // print 'rest', i.e. the expansion coefficient
 147                         if (i->rest.info(info_flags::numeric) &&
 148                                 i->rest.info(info_flags::positive)) {
 149                                 os << i->rest;
 150                         } else
 151                                 os << "(" << i->rest << ')';
 152                         // print 'coeff', something like (x-1)^42
 153                         if (!i->coeff.is_zero()) {
 154                                 os << '*';
 155                                 if (!point.is_zero())
 156                                         os << '(' << var-point << ')';
 157                                 else
 158                                         os << var;
 159                                 if (i->coeff.compare(_ex1())) {
 160                                         os << '^';
 161                                         if (i->coeff.info(info_flags::negative))
 162                                                 os << '(' << i->coeff << ')';
 163                                         else
 164                                                 os << i->coeff;
 165                                 }
 166                         }
 167                 } else {
 168                         os << Order(power(var-point,i->coeff));
 169                 }
 170         }
 171         if (precedence<=upper_precedence) os << ")";
 172 }
 173
 174
 175 void pseries::printraw(std::ostream &os) const
 176 {
 177         debugmsg("pseries printraw", LOGLEVEL_PRINT);
 178         os << class_name() << "(" << var << ";" << point << ";";
 179         for (epvector::const_iterator i=seq.begin(); i!=seq.end(); ++i)
 180                 os << "(" << (*i).rest << "," << (*i).coeff << "),";
 181         os << ")";
 182 }
 183
 184
 185 void pseries::printtree(std::ostream & os, unsigned indent) const
 186 {
 187         debugmsg("pseries printtree",LOGLEVEL_PRINT);
 188         os << std::string(indent,' ') << class_name()
 189            << ", hash=" << hashvalue
 190            << " (0x" << std::hex << hashvalue << std::dec << ")"
 191            << ", flags=" << flags << std::endl;
 192         for (unsigned i=0; i<seq.size(); ++i) {
 193                 seq[i].rest.printtree(os,indent+delta_indent);
 194                 seq[i].coeff.printtree(os,indent+delta_indent);
 195                 if (i!=seq.size()-1)
 196                         os << std::string(indent+delta_indent,' ') << "-----" << std::endl;
 197         }
 198         var.printtree(os, indent+delta_indent);
 199         point.printtree(os, indent+delta_indent);
 200 }
 201
 202 int pseries::compare_same_type(const basic & other) const
 203 {
 204         GINAC_ASSERT(is_of_type(other, pseries));
 205         const pseries &o = static_cast<const pseries &>(other);
 206
 207         // first compare the lengths of the series...
 208         if (seq.size()>o.seq.size())
 209                 return 1;
 210         if (seq.size()<o.seq.size())
 211                 return -1;
 212
 213         // ...then the expansion point...
 214         int cmpval = var.compare(o.var);
 215         if (cmpval)
 216                 return cmpval;
 217         cmpval = point.compare(o.point);
 218         if (cmpval)
 219                 return cmpval;
 220
 221         // ...and if that failed the individual elements
 222         epvector::const_iterator it = seq.begin(), o_it = o.seq.begin();
 223         while (it!=seq.end() && o_it!=o.seq.end()) {
 224                 cmpval = it->compare(*o_it);
 225                 if (cmpval)
 226                         return cmpval;
 227                 ++it;
 228                 ++o_it;
 229         }
 230
 231         // so they are equal.
 232         return 0;
 233 }
 234
 235 /** Return the number of operands including a possible order term. */
 236 unsigned pseries::nops(void) const
 237 {
 238         return seq.size();
 239 }
 240
 241
 242 /** Return the ith term in the series when represented as a sum. */
 243 ex pseries::op(int i) const
 244 {
 245         if (i < 0 || unsigned(i) >= seq.size())
 246                 throw (std::out_of_range("op() out of range"));
 247         return seq[i].rest * power(var - point, seq[i].coeff);
 248 }
 249
 250
 251 ex &pseries::let_op(int i)
 252 {
 253         throw (std::logic_error("let_op not defined for pseries"));
 254 }
 255
 256
 257 /** Return degree of highest power of the series.  This is usually the exponent
 258  *  of the Order term.  If s is not the expansion variable of the series, the
 259  *  series is examined termwise. */
 260 int pseries::degree(const symbol &s) const
 261 {
 262         if (var.is_equal(s)) {
 263                 // Return last exponent
 264                 if (seq.size())
 265                         return ex_to_numeric((*(seq.end() - 1)).coeff).to_int();
 266                 else
 267                         return 0;
 268         } else {
 269                 epvector::const_iterator it = seq.begin(), itend = seq.end();
 270                 if (it == itend)
 271                         return 0;
 272                 int max_pow = INT_MIN;
 273                 while (it != itend) {
 274                         int pow = it->rest.degree(s);
 275                         if (pow > max_pow)
 276                                 max_pow = pow;
 277                         ++it;
 278                 }
 279                 return max_pow;
 280         }
 281 }
 282
 283 /** Return degree of lowest power of the series.  This is usually the exponent
 284  *  of the leading term.  If s is not the expansion variable of the series, the
 285  *  series is examined termwise.  If s is the expansion variable but the
 286  *  expansion point is not zero the series is not expanded to find the degree.
 287  *  I.e.: (1-x) + (1-x)^2 + Order((1-x)^3) has ldegree(x) 1, not 0. */
 288 int pseries::ldegree(const symbol &s) const
 289 {
 290         if (var.is_equal(s)) {
 291                 // Return first exponent
 292                 if (seq.size())
 293                         return ex_to_numeric((*(seq.begin())).coeff).to_int();
 294                 else
 295                         return 0;
 296         } else {
 297                 epvector::const_iterator it = seq.begin(), itend = seq.end();
 298                 if (it == itend)
 299                         return 0;
 300                 int min_pow = INT_MAX;
 301                 while (it != itend) {
 302                         int pow = it->rest.ldegree(s);
 303                         if (pow < min_pow)
 304                                 min_pow = pow;
 305                         ++it;
 306                 }
 307                 return min_pow;
 308         }
 309 }
 310
 311 /** Return coefficient of degree n in power series if s is the expansion
 312  *  variable.  If the expansion point is nonzero, by definition the n=1
 313  *  coefficient in s of a+b*(s-z)+c*(s-z)^2+Order((s-z)^3) is b (assuming
 314  *  the expansion took place in the s in the first place).
 315  *  If s is not the expansion variable, an attempt is made to convert the
 316  *  series to a polynomial and return the corresponding coefficient from
 317  *  there. */
 318 ex pseries::coeff(const symbol &s, int n) const
 319 {
 320         if (var.is_equal(s)) {
 321                 if (seq.size() == 0)
 322                         return _ex0();
 323
 324                 // Binary search in sequence for given power
 325                 numeric looking_for = numeric(n);
 326                 int lo = 0, hi = seq.size() - 1;
 327                 while (lo <= hi) {
 328                         int mid = (lo + hi) / 2;
 329                         GINAC_ASSERT(is_ex_exactly_of_type(seq[mid].coeff, numeric));
 330                         int cmp = ex_to_numeric(seq[mid].coeff).compare(looking_for);
 331                         switch (cmp) {
 332                                 case -1:
 333                                         lo = mid + 1;
 334                                         break;
 335                                 case 0:
 336                                         return seq[mid].rest;
 337                                 case 1:
 338                                         hi = mid - 1;
 339                                         break;
 340                                 default:
 341                                         throw(std::logic_error("pseries::coeff: compare() didn't return -1, 0 or 1"));
 342                         }
 343                 }
 344                 return _ex0();
 345         } else
 346                 return convert_to_poly().coeff(s, n);
 347 }
 348
 349 /** Does nothing. */
 350 ex pseries::collect(const symbol &s) const
 351 {
 352         return *this;
 353 }
 354
 355
 356 /** Evaluate coefficients. */
 357 ex pseries::eval(int level) const
 358 {
 359         if (level == 1)
 360                 return this->hold();
 361
 362         if (level == -max_recursion_level)
 363                 throw (std::runtime_error("pseries::eval(): recursion limit exceeded"));
 364
 365         // Construct a new series with evaluated coefficients
 366         epvector new_seq;
 367         new_seq.reserve(seq.size());
 368         epvector::const_iterator it = seq.begin(), itend = seq.end();
 369         while (it != itend) {
 370                 new_seq.push_back(expair(it->rest.eval(level-1), it->coeff));
 371                 ++it;
 372         }
 373         return (new pseries(relational(var,point), new_seq))->setflag(status_flags::dynallocated | status_flags::evaluated);
 374 }
 375
 376
 377 /** Evaluate coefficients numerically. */
 378 ex pseries::evalf(int level) const
 379 {
 380         if (level == 1)
 381                 return *this;
 382
 383         if (level == -max_recursion_level)
 384                 throw (std::runtime_error("pseries::evalf(): recursion limit exceeded"));
 385
 386         // Construct a new series with evaluated coefficients
 387         epvector new_seq;
 388         new_seq.reserve(seq.size());
 389         epvector::const_iterator it = seq.begin(), itend = seq.end();
 390         while (it != itend) {
 391                 new_seq.push_back(expair(it->rest.evalf(level-1), it->coeff));
 392                 ++it;
 393         }
 394         return (new pseries(relational(var,point), new_seq))->setflag(status_flags::dynallocated | status_flags::evaluated);
 395 }
 396
 397
 398 ex pseries::subs(const lst & ls, const lst & lr) const
 399 {
 400         // If expansion variable is being substituted, convert the series to a
 401         // polynomial and do the substitution there because the result might
 402         // no longer be a power series
 403         if (ls.has(var))
 404                 return convert_to_poly(true).subs(ls, lr);
 405
 406         // Otherwise construct a new series with substituted coefficients and
 407         // expansion point
 408         epvector newseq;
 409         newseq.reserve(seq.size());
 410         epvector::const_iterator it = seq.begin(), itend = seq.end();
 411         while (it != itend) {
 412                 newseq.push_back(expair(it->rest.subs(ls, lr), it->coeff));
 413                 ++it;
 414         }
 415         return (new pseries(relational(var,point.subs(ls, lr)), newseq))->setflag(status_flags::dynallocated);
 416 }
 417
 418
 419 /** Implementation of ex::expand() for a power series.  It expands all the
 420  *  terms individually and returns the resulting series as a new pseries. */
 421 ex pseries::expand(unsigned options) const
 422 {
 423         epvector newseq;
 424         for (epvector::const_iterator i=seq.begin(); i!=seq.end(); ++i) {
 425                 ex restexp = i->rest.expand();
 426                 if (!restexp.is_zero())
 427                         newseq.push_back(expair(restexp, i->coeff));
 428         }
 429         return (new pseries(relational(var,point), newseq))
 430                 ->setflag(status_flags::dynallocated | status_flags::expanded);
 431 }
 432
 433
 434 /** Implementation of ex::diff() for a power series.  It treats the series as a
 435  *  polynomial.
 436  *  @see ex::diff */
 437 ex pseries::derivative(const symbol & s) const
 438 {
 439         if (s == var) {
 440                 epvector new_seq;
 441                 epvector::const_iterator it = seq.begin(), itend = seq.end();
 442
 443                 // FIXME: coeff might depend on var
 444                 while (it != itend) {
 445                         if (is_order_function(it->rest)) {
 446                                 new_seq.push_back(expair(it->rest, it->coeff - 1));
 447                         } else {
 448                                 ex c = it->rest * it->coeff;
 449                                 if (!c.is_zero())
 450                                         new_seq.push_back(expair(c, it->coeff - 1));
 451                         }
 452                         ++it;
 453                 }
 454                 return pseries(relational(var,point), new_seq);
 455         } else {
 456                 return *this;
 457         }
 458 }
 459
 460
 461 /** Convert a pseries object to an ordinary polynomial.
 462  *
 463  *  @param no_order flag: discard higher order terms */
 464 ex pseries::convert_to_poly(bool no_order) const
 465 {
 466         ex e;
 467         epvector::const_iterator it = seq.begin(), itend = seq.end();
 468
 469         while (it != itend) {
 470                 if (is_order_function(it->rest)) {
 471                         if (!no_order)
 472                                 e += Order(power(var - point, it->coeff));
 473                 } else
 474                         e += it->rest * power(var - point, it->coeff);
 475                 ++it;
 476         }
 477         return e;
 478 }
 479
 480
 481 /** Returns true if there is no order term, i.e. the series terminates and
 482  *  false otherwise. */
 483 bool pseries::is_terminating(void) const
 484 {
 485         return seq.size() == 0 || !is_order_function((seq.end()-1)->rest);
 486 }
 487
 488
 489 /*
 490  *  Implementations of series expansion
 491  */
 492
 493 /** Default implementation of ex::series(). This performs Taylor expansion.
 494  *  @see ex::series */
 495 ex basic::series(const relational & r, int order, unsigned options) const
 496 {
 497         epvector seq;
 498         numeric fac(1);
 499         ex deriv = *this;
 500         ex coeff = deriv.subs(r);
 501         const symbol *s = static_cast<symbol *>(r.lhs().bp);
 502
 503         if (!coeff.is_zero())
 504                 seq.push_back(expair(coeff, numeric(0)));
 505
 506         int n;
 507         for (n=1; n<order; ++n) {
 508                 fac = fac.mul(numeric(n));
 509                 deriv = deriv.diff(*s).expand();
 510                 if (deriv.is_zero()) {
 511                         // Series terminates
 512                         return pseries(r, seq);
 513                 }
 514                 coeff = deriv.subs(r);
 515                 if (!coeff.is_zero())
 516                         seq.push_back(expair(fac.inverse() * coeff, numeric(n)));
 517         }
 518
 519         // Higher-order terms, if present
 520         deriv = deriv.diff(*s);
 521         if (!deriv.expand().is_zero())
 522                 seq.push_back(expair(Order(_ex1()), numeric(n)));
 523         return pseries(r, seq);
 524 }
 525
 526
 527 /** Implementation of ex::series() for symbols.
 528  *  @see ex::series */
 529 ex symbol::series(const relational & r, int order, unsigned options) const
 530 {
 531         epvector seq;
 532         const ex point = r.rhs();
 533         GINAC_ASSERT(is_ex_exactly_of_type(r.lhs(),symbol));
 534         const symbol *s = static_cast<symbol *>(r.lhs().bp);
 535
 536         if (this->is_equal(*s)) {
 537                 if (order > 0 && !point.is_zero())
 538                         seq.push_back(expair(point, _ex0()));
 539                 if (order > 1)
 540                         seq.push_back(expair(_ex1(), _ex1()));
 541                 else
 542                         seq.push_back(expair(Order(_ex1()), numeric(order)));
 543         } else
 544                 seq.push_back(expair(*this, _ex0()));
 545         return pseries(r, seq);
 546 }
 547
 548
 549 /** Add one series object to another, producing a pseries object that
 550  *  represents the sum.
 551  *
 552  *  @param other  pseries object to add with
 553  *  @return the sum as a pseries */
 554 ex pseries::add_series(const pseries &other) const
 555 {
 556         // Adding two series with different variables or expansion points
 557         // results in an empty (constant) series
 558         if (!is_compatible_to(other)) {
 559                 epvector nul;
 560                 nul.push_back(expair(Order(_ex1()), _ex0()));
 561                 return pseries(relational(var,point), nul);
 562         }
 563
 564         // Series addition
 565         epvector new_seq;
 566         epvector::const_iterator a = seq.begin();
 567         epvector::const_iterator b = other.seq.begin();
 568         epvector::const_iterator a_end = seq.end();
 569         epvector::const_iterator b_end = other.seq.end();
 570         int pow_a = INT_MAX, pow_b = INT_MAX;
 571         for (;;) {
 572                 // If a is empty, fill up with elements from b and stop
 573                 if (a == a_end) {
 574                         while (b != b_end) {
 575                                 new_seq.push_back(*b);
 576                                 ++b;
 577                         }
 578                         break;
 579                 } else
 580                         pow_a = ex_to_numeric((*a).coeff).to_int();
 581
 582                 // If b is empty, fill up with elements from a and stop
 583                 if (b == b_end) {
 584                         while (a != a_end) {
 585                                 new_seq.push_back(*a);
 586                                 ++a;
 587                         }
 588                         break;
 589                 } else
 590                         pow_b = ex_to_numeric((*b).coeff).to_int();
 591
 592                 // a and b are non-empty, compare powers
 593                 if (pow_a < pow_b) {
 594                         // a has lesser power, get coefficient from a
 595                         new_seq.push_back(*a);
 596                         if (is_order_function((*a).rest))
 597                                 break;
 598                         ++a;
 599                 } else if (pow_b < pow_a) {
 600                         // b has lesser power, get coefficient from b
 601                         new_seq.push_back(*b);
 602                         if (is_order_function((*b).rest))
 603                                 break;
 604                         ++b;
 605                 } else {
 606                         // Add coefficient of a and b
 607                         if (is_order_function((*a).rest) || is_order_function((*b).rest)) {
 608                                 new_seq.push_back(expair(Order(_ex1()), (*a).coeff));
 609                                 break;  // Order term ends the sequence
 610                         } else {
 611                                 ex sum = (*a).rest + (*b).rest;
 612                                 if (!(sum.is_zero()))
 613                                         new_seq.push_back(expair(sum, numeric(pow_a)));
 614                                 ++a;
 615                                 ++b;
 616                         }
 617                 }
 618         }
 619         return pseries(relational(var,point), new_seq);
 620 }
 621
 622
 623 /** Implementation of ex::series() for sums. This performs series addition when
 624  *  adding pseries objects.
 625  *  @see ex::series */
 626 ex add::series(const relational & r, int order, unsigned options) const
 627 {
 628         ex acc; // Series accumulator
 629
 630         // Get first term from overall_coeff
 631         acc = overall_coeff.series(r, order, options);
 632
 633         // Add remaining terms
 634         epvector::const_iterator it = seq.begin();
 635         epvector::const_iterator itend = seq.end();
 636         for (; it!=itend; ++it) {
 637                 ex op;
 638                 if (is_ex_exactly_of_type(it->rest, pseries))
 639                         op = it->rest;
 640                 else
 641                         op = it->rest.series(r, order, options);
 642                 if (!it->coeff.is_equal(_ex1()))
 643                         op = ex_to_pseries(op).mul_const(ex_to_numeric(it->coeff));
 644
 645                 // Series addition
 646                 acc = ex_to_pseries(acc).add_series(ex_to_pseries(op));
 647         }
 648         return acc;
 649 }
 650
 651
 652 /** Multiply a pseries object with a numeric constant, producing a pseries
 653  *  object that represents the product.
 654  *
 655  *  @param other  constant to multiply with
 656  *  @return the product as a pseries */
 657 ex pseries::mul_const(const numeric &other) const
 658 {
 659         epvector new_seq;
 660         new_seq.reserve(seq.size());
 661
 662         epvector::const_iterator it = seq.begin(), itend = seq.end();
 663         while (it != itend) {
 664                 if (!is_order_function(it->rest))
 665                         new_seq.push_back(expair(it->rest * other, it->coeff));
 666                 else
 667                         new_seq.push_back(*it);
 668                 ++it;
 669         }
 670         return pseries(relational(var,point), new_seq);
 671 }
 672
 673
 674 /** Multiply one pseries object to another, producing a pseries object that
 675  *  represents the product.
 676  *
 677  *  @param other  pseries object to multiply with
 678  *  @return the product as a pseries */
 679 ex pseries::mul_series(const pseries &other) const
 680 {
 681         // Multiplying two series with different variables or expansion points
 682         // results in an empty (constant) series
 683         if (!is_compatible_to(other)) {
 684                 epvector nul;
 685                 nul.push_back(expair(Order(_ex1()), _ex0()));
 686                 return pseries(relational(var,point), nul);
 687         }
 688
 689         // Series multiplication
 690         epvector new_seq;
 691
 692         const symbol *s = static_cast<symbol *>(var.bp);
 693         int a_max = degree(*s);
 694         int b_max = other.degree(*s);
 695         int a_min = ldegree(*s);
 696         int b_min = other.ldegree(*s);
 697         int cdeg_min = a_min + b_min;
 698         int cdeg_max = a_max + b_max;
 699
 700         int higher_order_a = INT_MAX;
 701         int higher_order_b = INT_MAX;
 702         if (is_order_function(coeff(*s, a_max)))
 703                 higher_order_a = a_max + b_min;
 704         if (is_order_function(other.coeff(*s, b_max)))
 705                 higher_order_b = b_max + a_min;
 706         int higher_order_c = std::min(higher_order_a, higher_order_b);
 707         if (cdeg_max >= higher_order_c)
 708                 cdeg_max = higher_order_c - 1;
 709
 710         for (int cdeg=cdeg_min; cdeg<=cdeg_max; ++cdeg) {
 711                 ex co = _ex0();
 712                 // c(i)=a(0)b(i)+...+a(i)b(0)
 713                 for (int i=a_min; cdeg-i>=b_min; ++i) {
 714                         ex a_coeff = coeff(*s, i);
 715                         ex b_coeff = other.coeff(*s, cdeg-i);
 716                         if (!is_order_function(a_coeff) && !is_order_function(b_coeff))
 717                                 co += a_coeff * b_coeff;
 718                 }
 719                 if (!co.is_zero())
 720                         new_seq.push_back(expair(co, numeric(cdeg)));
 721         }
 722         if (higher_order_c < INT_MAX)
 723                 new_seq.push_back(expair(Order(_ex1()), numeric(higher_order_c)));
 724         return pseries(relational(var,point), new_seq);
 725 }
 726
 727
 728 /** Implementation of ex::series() for product. This performs series
 729  *  multiplication when multiplying series.
 730  *  @see ex::series */
 731 ex mul::series(const relational & r, int order, unsigned options) const
 732 {
 733         ex acc; // Series accumulator
 734
 735         // Get first term from overall_coeff
 736         acc = overall_coeff.series(r, order, options);
 737
 738         // Multiply with remaining terms
 739         epvector::const_iterator it = seq.begin();
 740         epvector::const_iterator itend = seq.end();
 741         for (; it!=itend; ++it) {
 742                 ex op = it->rest;
 743                 if (op.info(info_flags::numeric)) {
 744                         // series * const (special case, faster)
 745                         ex f = power(op, it->coeff);
 746                         acc = ex_to_pseries(acc).mul_const(ex_to_numeric(f));
 747                         continue;
 748                 } else if (!is_ex_exactly_of_type(op, pseries))
 749                         op = op.series(r, order, options);
 750                 if (!it->coeff.is_equal(_ex1()))
 751                         op = ex_to_pseries(op).power_const(ex_to_numeric(it->coeff), order);
 752
 753                 // Series multiplication
 754                 acc = ex_to_pseries(acc).mul_series(ex_to_pseries(op));
 755         }
 756         return acc;
 757 }
 758
 759
 760 /** Compute the p-th power of a series.
 761  *
 762  *  @param p  power to compute
 763  *  @param deg  truncation order of series calculation */
 764 ex pseries::power_const(const numeric &p, int deg) const
 765 {
 766         // method:
 767         // let A(x) be this series and for the time being let it start with a
 768         // constant (later we'll generalize):
 769         //     A(x) = a_0 + a_1*x + a_2*x^2 + ...
 770         // We want to compute
 771         //     C(x) = A(x)^p
 772         //     C(x) = c_0 + c_1*x + c_2*x^2 + ...
 773         // Taking the derivative on both sides and multiplying with A(x) one
 774         // immediately arrives at
 775         //     C'(x)*A(x) = p*C(x)*A'(x)
 776         // Multiplying this out and comparing coefficients we get the recurrence
 777         // formula
 778         //     c_i = (i*p*a_i*c_0 + ((i-1)*p-1)*a_{i-1}*c_1 + ...
 779         //                    ... + (p-(i-1))*a_1*c_{i-1})/(a_0*i)
 780         // which can easily be solved given the starting value c_0 = (a_0)^p.
 781         // For the more general case where the leading coefficient of A(x) is not
 782         // a constant, just consider A2(x) = A(x)*x^m, with some integer m and
 783         // repeat the above derivation.  The leading power of C2(x) = A2(x)^2 is
 784         // then of course x^(p*m) but the recurrence formula still holds.
 785
 786         if (seq.size()==0) {
 787                 // as a spacial case, handle the empty (zero) series honoring the
 788                 // usual power laws such as implemented in power::eval()
 789                 if (p.real().is_zero())
 790                         throw (std::domain_error("pseries::power_const(): pow(0,I) is undefined"));
 791                 else if (p.real().is_negative())
 792                         throw (pole_error("pseries::power_const(): division by zero",1));
 793                 else
 794                         return *this;
 795         }
 796
 797         const symbol *s = static_cast<symbol *>(var.bp);
 798         int ldeg = ldegree(*s);
 799
 800         // Compute coefficients of the powered series
 801         exvector co;
 802         co.reserve(deg);
 803         co.push_back(power(coeff(*s, ldeg), p));
 804         bool all_sums_zero = true;
 805         for (int i=1; i<deg; ++i) {
 806                 ex sum = _ex0();
 807                 for (int j=1; j<=i; ++j) {
 808                         ex c = coeff(*s, j + ldeg);
 809                         if (is_order_function(c)) {
 810                                 co.push_back(Order(_ex1()));
 811                                 break;
 812                         } else
 813                                 sum += (p * j - (i - j)) * co[i - j] * c;
 814                 }
 815                 if (!sum.is_zero())
 816                         all_sums_zero = false;
 817                 co.push_back(sum / coeff(*s, ldeg) / numeric(i));
 818         }
 819
 820         // Construct new series (of non-zero coefficients)
 821         epvector new_seq;
 822         bool higher_order = false;
 823         for (int i=0; i<deg; ++i) {
 824                 if (!co[i].is_zero())
 825                         new_seq.push_back(expair(co[i], numeric(i) + p * ldeg));
 826                 if (is_order_function(co[i])) {
 827                         higher_order = true;
 828                         break;
 829                 }
 830         }
 831         if (!higher_order && !all_sums_zero)
 832                 new_seq.push_back(expair(Order(_ex1()), numeric(deg) + p * ldeg));
 833         return pseries(relational(var,point), new_seq);
 834 }
 835
 836
 837 /** Return a new pseries object with the powers shifted by deg. */
 838 pseries pseries::shift_exponents(int deg) const
 839 {
 840         epvector newseq(seq);
 841         for (epvector::iterator i=newseq.begin(); i!=newseq.end(); ++i)
 842                 i->coeff = i->coeff + deg;
 843         return pseries(relational(var, point), newseq);
 844 }
 845
 846
 847 /** Implementation of ex::series() for powers. This performs Laurent expansion
 848  *  of reciprocals of series at singularities.
 849  *  @see ex::series */
 850 ex power::series(const relational & r, int order, unsigned options) const
 851 {
 852         ex e;
 853         if (!is_ex_exactly_of_type(basis, pseries)) {
 854                 // Basis is not a series, may there be a singularity?
 855                 bool must_expand_basis = false;
 856                 try {
 857                         basis.subs(r);
 858                 } catch (pole_error) {
 859                         must_expand_basis = true;
 860                 }
 861
 862                 // Is the expression of type something^(-int)?
 863                 if (!must_expand_basis && !exponent.info(info_flags::negint))
 864                         return basic::series(r, order, options);
 865
 866                 // Is the expression of type 0^something?
 867                 if (!must_expand_basis && !basis.subs(r).is_zero())
 868                         return basic::series(r, order, options);
 869
 870                 // Singularity encountered, expand basis into series
 871                 e = basis.series(r, order, options);
 872         } else {
 873                 // Basis is a series
 874                 e = basis;
 875         }
 876
 877         // Power e
 878         return ex_to_pseries(e).power_const(ex_to_numeric(exponent), order);
 879 }
 880
 881
 882 /** Re-expansion of a pseries object. */
 883 ex pseries::series(const relational & r, int order, unsigned options) const
 884 {
 885         const ex p = r.rhs();
 886         GINAC_ASSERT(is_ex_exactly_of_type(r.lhs(),symbol));
 887         const symbol *s = static_cast<symbol *>(r.lhs().bp);
 888
 889         if (var.is_equal(*s) && point.is_equal(p)) {
 890                 if (order > degree(*s))
 891                         return *this;
 892                 else {
 893                         epvector new_seq;
 894                         epvector::const_iterator it = seq.begin(), itend = seq.end();
 895                         while (it != itend) {
 896                                 int o = ex_to_numeric(it->coeff).to_int();
 897                                 if (o >= order) {
 898                                         new_seq.push_back(expair(Order(_ex1()), o));
 899                                         break;
 900                                 }
 901                                 new_seq.push_back(*it);
 902                                 ++it;
 903                         }
 904                         return pseries(r, new_seq);
 905                 }
 906         } else
 907                 return convert_to_poly().series(r, order, options);
 908 }
 909
 910
 911 /** Compute the truncated series expansion of an expression.
 912  *  This function returns an expression containing an object of class pseries
 913  *  to represent the series. If the series does not terminate within the given
 914  *  truncation order, the last term of the series will be an order term.
 915  *
 916  *  @param r  expansion relation, lhs holds variable and rhs holds point
 917  *  @param order  truncation order of series calculations
 918  *  @param options  of class series_options
 919  *  @return an expression holding a pseries object */
 920 ex ex::series(const ex & r, int order, unsigned options) const
 921 {
 922         GINAC_ASSERT(bp!=0);
 923         ex e;
 924         relational rel_;
 925
 926         if (is_ex_exactly_of_type(r,relational))
 927                 rel_ = ex_to_relational(r);
 928         else if (is_ex_exactly_of_type(r,symbol))
 929                 rel_ = relational(r,_ex0());
 930         else
 931                 throw (std::logic_error("ex::series(): expansion point has unknown type"));
 932
 933         try {
 934                 e = bp->series(rel_, order, options);
 935         } catch (std::exception &x) {
 936                 throw (std::logic_error(std::string("unable to compute series (") + x.what() + ")"));
 937         }
 938         return e;
 939 }
 940
 941 //////////
 942 // static member variables
 943 //////////
 944
 945 // protected
 946
 947 unsigned pseries::precedence = 38;  // for clarity just below add::precedence
 948
 949 } // namespace GiNaC