ginac/series.cpp

   1 /** @file series.cpp
   2  *
   3  *  Implementation of class for extended truncated power-series and
   4  *  methods for series expansion. */
   5
   6 #include "ginac.h"
   7
   8
   9 /*
  10  *  Default constructor, destructor, copy constructor, assignment operator and helpers
  11  */
  12
  13 series::series() : basic(TINFO_SERIES)
  14 {
  15     debugmsg("series default constructor", LOGLEVEL_CONSTRUCT);
  16 }
  17
  18 series::~series()
  19 {
  20     debugmsg("series destructor", LOGLEVEL_DESTRUCT);
  21     destroy(false);
  22 }
  23
  24 series::series(series const &other)
  25 {
  26     debugmsg("series copy constructor", LOGLEVEL_CONSTRUCT);
  27     copy(other);
  28 }
  29
  30 series const &series::operator=(series const & other)
  31 {
  32     debugmsg("series operator=", LOGLEVEL_ASSIGNMENT);
  33     if (this != &other) {
  34         destroy(true);
  35         copy(other);
  36     }
  37     return *this;
  38 }
  39
  40 void series::copy(series const &other)
  41 {
  42     inherited::copy(other);
  43     seq = other.seq;
  44     var = other.var;
  45     point = other.point;
  46 }
  47
  48 void series::destroy(bool call_parent)
  49 {
  50     if (call_parent)
  51         inherited::destroy(call_parent);
  52 }
  53
  54
  55 /*
  56  *  Other constructors
  57  */
  58
  59 /** Construct series from a vector of coefficients and powers.
  60  *  expair.rest holds the coefficient, expair.coeff holds the power.
  61  *  The powers must be integers (positive or negative) and in ascending order;
  62  *  the last coefficient can be Order(exONE()) to represent a truncated,
  63  *  non-terminating series.
  64  *
  65  *  @param var_  series variable (must hold a symbol)
  66  *  @param point_  expansion point
  67  *  @param ops_  vector of {coefficient, power} pairs (coefficient must not be zero)
  68  *  @return newly constructed series */
  69 series::series(ex const &var_, ex const &point_, epvector const &ops_)
  70     : basic(TINFO_SERIES), seq(ops_), var(var_), point(point_)
  71 {
  72     debugmsg("series constructor from ex,ex,epvector", LOGLEVEL_CONSTRUCT);
  73     ASSERT(is_ex_exactly_of_type(var_, symbol));
  74 }
  75
  76
  77 /*
  78  *  Functions overriding virtual functions from base classes
  79  */
  80
  81 basic *series::duplicate() const
  82 {
  83     debugmsg("series duplicate", LOGLEVEL_DUPLICATE);
  84     return new series(*this);
  85 }
  86
  87 // Highest degree of variable
  88 int series::degree(symbol const &s) const
  89 {
  90     if (var == s) {
  91         // Return last exponent
  92         if (seq.size())
  93             return ex_to_numeric((*(seq.end() - 1)).coeff).to_int();
  94         else
  95             return 0;
  96     } else {
  97         epvector::const_iterator it = seq.begin(), itend = seq.end();
  98         if (it == itend)
  99             return 0;
 100         int max_pow = INT_MIN;
 101         while (it != itend) {
 102             int pow = it->rest.degree(s);
 103             if (pow > max_pow)
 104                 max_pow = pow;
 105             it++;
 106         }
 107         return max_pow;
 108     }
 109 }
 110
 111 // Lowest degree of variable
 112 int series::ldegree(symbol const &s) const
 113 {
 114     if (var == s) {
 115         // Return first exponent
 116         if (seq.size())
 117             return ex_to_numeric((*(seq.begin())).coeff).to_int();
 118         else
 119             return 0;
 120     } else {
 121         epvector::const_iterator it = seq.begin(), itend = seq.end();
 122         if (it == itend)
 123             return 0;
 124         int min_pow = INT_MAX;
 125         while (it != itend) {
 126             int pow = it->rest.ldegree(s);
 127             if (pow < min_pow)
 128                 min_pow = pow;
 129             it++;
 130         }
 131         return min_pow;
 132     }
 133 }
 134
 135 // Coefficient of variable
 136 ex series::coeff(symbol const &s, int n) const
 137 {
 138     if (var == s) {
 139         epvector::const_iterator it = seq.begin(), itend = seq.end();
 140         while (it != itend) {
 141             int pow = ex_to_numeric(it->coeff).to_int();
 142             if (pow == n)
 143                 return it->rest;
 144             if (pow > n)
 145                 return exZERO();
 146             it++;
 147         }
 148         return exZERO();
 149     } else
 150         return convert_to_poly().coeff(s, n);
 151 }
 152
 153 ex series::eval(int level) const
 154 {
 155     if (level == 1)
 156         return this->hold();
 157
 158     // Construct a new series with evaluated coefficients
 159     epvector new_seq;
 160     new_seq.reserve(seq.size());
 161     epvector::const_iterator it = seq.begin(), itend = seq.end();
 162     while (it != itend) {
 163         new_seq.push_back(expair(it->rest.eval(level-1), it->coeff));
 164         it++;
 165     }
 166     return (new series(var, point, new_seq))->setflag(status_flags::dynallocated | status_flags::evaluated);
 167 }
 168
 169 ex series::evalf(int level) const
 170 {
 171     return convert_to_poly().evalf(level);
 172 }
 173
 174
 175 /*
 176  *  Construct expression (polynomial) out of series
 177  */
 178
 179 /** Convert a series object to an ordinary polynomial.
 180  *
 181  *  @param no_order flag: discard higher order terms */
 182 ex series::convert_to_poly(bool no_order) const
 183 {
 184     ex e;
 185     epvector::const_iterator it = seq.begin(), itend = seq.end();
 186
 187     while (it != itend) {
 188         if (is_order_function(it->rest)) {
 189             if (!no_order)
 190                 e += Order(power(var - point, it->coeff));
 191         } else
 192             e += it->rest * power(var - point, it->coeff);
 193         it++;
 194     }
 195     return e;
 196 }
 197
 198
 199 /*
 200  *  Implementation of series expansion
 201  */
 202
 203 /** Default implementation of ex::series(). This performs Taylor expansion.
 204  *  @see ex::series */
 205 ex basic::series(symbol const & s, ex const & point, int order) const
 206 {
 207     epvector seq;
 208     numeric fac(1);
 209     ex deriv = *this;
 210     ex coeff = deriv.subs(s == point);
 211     if (!coeff.is_zero())
 212         seq.push_back(expair(coeff, numeric(0)));
 213
 214     int n;
 215     for (n=1; n<order; n++) {
 216         fac = fac.mul(numeric(n));
 217         deriv = deriv.diff(s).expand();
 218         if (deriv.is_zero()) {
 219             // Series terminates
 220             return series::series(s, point, seq);
 221         }
 222         coeff = power(fac, -1) * deriv.subs(s == point);
 223         if (!coeff.is_zero())
 224             seq.push_back(expair(coeff, numeric(n)));
 225     }
 226
 227     // Higher-order terms, if present
 228     deriv = deriv.diff(s);
 229     if (!deriv.is_zero())
 230         seq.push_back(expair(Order(exONE()), numeric(n)));
 231     return series::series(s, point, seq);
 232 }
 233
 234
 235 /** Add one series object to another, producing a series object that represents
 236  *  the sum.
 237  *
 238  *  @param other  series object to add with
 239  *  @return the sum as a series */
 240 ex series::add_series(const series &other) const
 241 {
 242     // Adding two series with different variables or expansion points
 243     // results in an empty (constant) series
 244     if (!is_compatible_to(other)) {
 245         epvector nul;
 246         nul.push_back(expair(Order(exONE()), exZERO()));
 247         return series(var, point, nul);
 248     }
 249
 250     // Series addition
 251     epvector new_seq;
 252     epvector::const_iterator a = seq.begin();
 253     epvector::const_iterator b = other.seq.begin();
 254     epvector::const_iterator a_end = seq.end();
 255     epvector::const_iterator b_end = other.seq.end();
 256     int pow_a = INT_MAX, pow_b = INT_MAX;
 257     for (;;) {
 258         // If a is empty, fill up with elements from b and stop
 259         if (a == a_end) {
 260             while (b != b_end) {
 261                 new_seq.push_back(*b);
 262                 b++;
 263             }
 264             break;
 265         } else
 266             pow_a = ex_to_numeric((*a).coeff).to_int();
 267
 268         // If b is empty, fill up with elements from a and stop
 269         if (b == b_end) {
 270             while (a != a_end) {
 271                 new_seq.push_back(*a);
 272                 a++;
 273             }
 274             break;
 275         } else
 276             pow_b = ex_to_numeric((*b).coeff).to_int();
 277
 278         // a and b are non-empty, compare powers
 279         if (pow_a < pow_b) {
 280             // a has lesser power, get coefficient from a
 281             new_seq.push_back(*a);
 282             if (is_order_function((*a).rest))
 283                 break;
 284             a++;
 285         } else if (pow_b < pow_a) {
 286             // b has lesser power, get coefficient from b
 287             new_seq.push_back(*b);
 288             if (is_order_function((*b).rest))
 289                 break;
 290             b++;
 291         } else {
 292             // Add coefficient of a and b
 293             if (is_order_function((*a).rest) || is_order_function((*b).rest)) {
 294                 new_seq.push_back(expair(Order(exONE()), (*a).coeff));
 295                 break;  // Order term ends the sequence
 296             } else {
 297                 ex sum = (*a).rest + (*b).rest;
 298                 if (!(sum == exZERO()))
 299                     new_seq.push_back(expair(sum, numeric(pow_a)));
 300                 a++;
 301                 b++;
 302             }
 303         }
 304     }
 305     return series(var, point, new_seq);
 306 }
 307
 308
 309 /** Implementation of ex::series() for sums. This performs series addition when
 310  *  adding series objects.
 311  *  @see ex::series */
 312 /*
 313 ex add::series(symbol const & s, ex const & point, int order) const
 314 {
 315     ex acc; // Series accumulator
 316
 317     // Get first term
 318     epvector::const_iterator it = seq.begin();
 319     epvector::const_iterator itend = seq.end();
 320     if (it != itend) {
 321         if (is_ex_exactly_of_type(it->rest, series))
 322             acc = it->rest;
 323         else
 324             acc = it->rest.series(s, point, order);
 325         if (!it->coeff.is_equal(exONE()))
 326             acc = ex_to_series(acc).mul_const(ex_to_numeric(it->coeff));
 327         it++;
 328     }
 329
 330     // Add remaining terms
 331     for (; it!=itend; it++) {
 332         ex op;
 333         if (is_ex_exactly_of_type(it->rest, series))
 334             op = it->rest;
 335         else
 336             op = it->rest.series(s, point, order);
 337         if (!it->coeff.is_equal(exONE()))
 338             op = ex_to_series(op).mul_const(ex_to_numeric(it->coeff));
 339
 340         // Series addition
 341         acc = ex_to_series(acc).add_series(ex_to_series(op));
 342     }
 343     return acc;
 344 }
 345 */
 346 ex add::series(symbol const & s, ex const & point, int order) const
 347 {
 348     ex acc; // Series accumulator
 349
 350     // Get first term from overall_coeff
 351     acc = overall_coeff.series(s,point,order);
 352
 353     // Add remaining terms
 354     epvector::const_iterator it = seq.begin();
 355     epvector::const_iterator itend = seq.end();
 356     for (; it!=itend; it++) {
 357         ex op;
 358         if (is_ex_exactly_of_type(it->rest, series))
 359             op = it->rest;
 360         else
 361             op = it->rest.series(s, point, order);
 362         if (!it->coeff.is_equal(exONE()))
 363             op = ex_to_series(op).mul_const(ex_to_numeric(it->coeff));
 364
 365         // Series addition
 366         acc = ex_to_series(acc).add_series(ex_to_series(op));
 367     }
 368     return acc;
 369 }
 370
 371
 372 /** Multiply a series object with a numeric constant, producing a series object
 373  *  that represents the product.
 374  *
 375  *  @param other  constant to multiply with
 376  *  @return the product as a series */
 377 ex series::mul_const(const numeric &other) const
 378 {
 379     epvector new_seq;
 380     new_seq.reserve(seq.size());
 381
 382     epvector::const_iterator it = seq.begin(), itend = seq.end();
 383     while (it != itend) {
 384         if (!is_order_function(it->rest))
 385             new_seq.push_back(expair(it->rest * other, it->coeff));
 386         else
 387             new_seq.push_back(*it);
 388         it++;
 389     }
 390     return series(var, point, new_seq);
 391 }
 392
 393
 394 /** Multiply one series object to another, producing a series object that
 395  *  represents the product.
 396  *
 397  *  @param other  series object to multiply with
 398  *  @return the product as a series */
 399 ex series::mul_series(const series &other) const
 400 {
 401     // Multiplying two series with different variables or expansion points
 402     // results in an empty (constant) series
 403     if (!is_compatible_to(other)) {
 404         epvector nul;
 405         nul.push_back(expair(Order(exONE()), exZERO()));
 406         return series(var, point, nul);
 407     }
 408
 409     // Series multiplication
 410     epvector new_seq;
 411
 412     const symbol *s = static_cast<symbol *>(var.bp);
 413     int a_max = degree(*s);
 414     int b_max = other.degree(*s);
 415     int a_min = ldegree(*s);
 416     int b_min = other.ldegree(*s);
 417     int cdeg_min = a_min + b_min;
 418     int cdeg_max = a_max + b_max;
 419
 420     int higher_order_a = INT_MAX;
 421     int higher_order_b = INT_MAX;
 422     if (is_order_function(coeff(*s, a_max)))
 423         higher_order_a = a_max + b_min;
 424     if (is_order_function(other.coeff(*s, b_max)))
 425         higher_order_b = b_max + a_min;
 426     int higher_order_c = min(higher_order_a, higher_order_b);
 427     if (cdeg_max >= higher_order_c)
 428         cdeg_max = higher_order_c - 1;
 429
 430     for (int cdeg=cdeg_min; cdeg<=cdeg_max; cdeg++) {
 431         ex co = exZERO();
 432         // c(i)=a(0)b(i)+...+a(i)b(0)
 433         for (int i=a_min; cdeg-i>=b_min; i++) {
 434             ex a_coeff = coeff(*s, i);
 435             ex b_coeff = other.coeff(*s, cdeg-i);
 436             if (!is_order_function(a_coeff) && !is_order_function(b_coeff))
 437                 co += coeff(*s, i) * other.coeff(*s, cdeg-i);
 438         }
 439         if (!co.is_zero())
 440             new_seq.push_back(expair(co, numeric(cdeg)));
 441     }
 442     if (higher_order_c < INT_MAX)
 443         new_seq.push_back(expair(Order(exONE()), numeric(higher_order_c)));
 444     return series::series(var, point, new_seq);
 445 }
 446
 447
 448 /** Implementation of ex::series() for product. This performs series multiplication when multiplying series.
 449  *  @see ex::series */
 450 /*
 451 ex mul::series(symbol const & s, ex const & point, int order) const
 452 {
 453     ex acc; // Series accumulator
 454
 455     // Get first term
 456     epvector::const_iterator it = seq.begin();
 457     epvector::const_iterator itend = seq.end();
 458     if (it != itend) {
 459         if (is_ex_exactly_of_type(it->rest, series))
 460             acc = it->rest;
 461         else
 462             acc = it->rest.series(s, point, order);
 463         if (!it->coeff.is_equal(exONE()))
 464             acc = ex_to_series(acc).power_const(ex_to_numeric(it->coeff), order);
 465         it++;
 466     }
 467
 468     // Multiply with remaining terms
 469     for (; it!=itend; it++) {
 470         ex op = it->rest;
 471         if (op.info(info_flags::numeric)) {
 472             // series * const (special case, faster)
 473             ex f = power(op, it->coeff);
 474             acc = ex_to_series(acc).mul_const(ex_to_numeric(f));
 475             continue;
 476         } else if (!is_ex_exactly_of_type(op, series))
 477             op = op.series(s, point, order);
 478         if (!it->coeff.is_equal(exONE()))
 479             op = ex_to_series(op).power_const(ex_to_numeric(it->coeff), order);
 480
 481         // Series multiplication
 482         acc = ex_to_series(acc).mul_series(ex_to_series(op));
 483     }
 484     return acc;
 485 }
 486 */
 487
 488 ex mul::series(symbol const & s, ex const & point, int order) const
 489 {
 490     ex acc; // Series accumulator
 491
 492     // Get first term from overall_coeff
 493     acc = overall_coeff.series(s, point, order);
 494
 495     // Multiply with remaining terms
 496     epvector::const_iterator it = seq.begin();
 497     epvector::const_iterator itend = seq.end();
 498     for (; it!=itend; it++) {
 499         ex op = it->rest;
 500         if (op.info(info_flags::numeric)) {
 501             // series * const (special case, faster)
 502             ex f = power(op, it->coeff);
 503             acc = ex_to_series(acc).mul_const(ex_to_numeric(f));
 504             continue;
 505         } else if (!is_ex_exactly_of_type(op, series))
 506             op = op.series(s, point, order);
 507         if (!it->coeff.is_equal(exONE()))
 508             op = ex_to_series(op).power_const(ex_to_numeric(it->coeff), order);
 509
 510         // Series multiplication
 511         acc = ex_to_series(acc).mul_series(ex_to_series(op));
 512     }
 513     return acc;
 514 }
 515
 516
 517 /** Compute the p-th power of a series.
 518  *
 519  *  @param p  power to compute
 520  *  @param deg  truncation order of series calculation */
 521 ex series::power_const(const numeric &p, int deg) const
 522 {
 523     int i;
 524     const symbol *s = static_cast<symbol *>(var.bp);
 525     int ldeg = ldegree(*s);
 526
 527     // Calculate coefficients of powered series
 528     exvector co;
 529     co.reserve(deg);
 530     ex co0;
 531     co.push_back(co0 = power(coeff(*s, ldeg), p));
 532     bool all_sums_zero = true;
 533     for (i=1; i<deg; i++) {
 534         ex sum = exZERO();
 535         for (int j=1; j<=i; j++) {
 536             ex c = coeff(*s, j + ldeg);
 537             if (is_order_function(c)) {
 538                 co.push_back(Order(exONE()));
 539                 break;
 540             } else
 541                 sum += (p * j - (i - j)) * co[i - j] * c;
 542         }
 543         if (!sum.is_zero())
 544             all_sums_zero = false;
 545         co.push_back(co0 * sum / numeric(i));
 546     }
 547
 548     // Construct new series (of non-zero coefficients)
 549     epvector new_seq;
 550     bool higher_order = false;
 551     for (i=0; i<deg; i++) {
 552         if (!co[i].is_zero())
 553             new_seq.push_back(expair(co[i], numeric(i) + p * ldeg));
 554         if (is_order_function(co[i])) {
 555             higher_order = true;
 556             break;
 557         }
 558     }
 559     if (!higher_order && !all_sums_zero)
 560         new_seq.push_back(expair(Order(exONE()), numeric(deg) + p * ldeg));
 561     return series::series(var, point, new_seq);
 562 }
 563
 564
 565 /** Implementation of ex::series() for powers. This performs Laurent expansion
 566  *  of reciprocals of series at singularities.
 567  *  @see ex::series */
 568 ex power::series(symbol const & s, ex const & point, int order) const
 569 {
 570     ex e;
 571     if (!is_ex_exactly_of_type(basis, series)) {
 572         // Basis is not a series, may there be a singulary?
 573         if (!exponent.info(info_flags::negint))
 574             return basic::series(s, point, order);
 575
 576         // Expression is of type something^(-int), check for singularity
 577         if (!basis.subs(s == point).is_zero())
 578             return basic::series(s, point, order);
 579
 580         // Singularity encountered, expand basis into series
 581         e = basis.series(s, point, order);
 582     } else {
 583         // Basis is a series
 584         e = basis;
 585     }
 586
 587     // Power e
 588     return ex_to_series(e).power_const(ex_to_numeric(exponent), order);
 589 }
 590
 591
 592 /** Compute the truncated series expansion of an expression.
 593  *  This function returns an expression containing an object of class series to
 594  *  represent the series. If the series does not terminate within the given
 595  *  truncation order, the last term of the series will be an order term.
 596  *
 597  *  @param s  expansion variable
 598  *  @param point  expansion point
 599  *  @param order  truncation order of series calculations
 600  *  @return an expression holding a series object */
 601 ex ex::series(symbol const &s, ex const &point, int order) const
 602 {
 603     ASSERT(bp!=0);
 604     return bp->series(s, point, order);
 605 }
 606
 607
 608 // Global constants
 609 const series some_series;
 610 type_info const & typeid_series = typeid(some_series);