ginac/inifcns.cpp

   1 /** @file inifcns.cpp
   2  *
   3  *  Implementation of GiNaC's initially known functions. */
   4
   5 /*
   6  *  GiNaC Copyright (C) 1999-2015 Johannes Gutenberg University Mainz, Germany
   7  *
   8  *  This program is free software; you can redistribute it and/or modify
   9  *  it under the terms of the GNU General Public License as published by
  10  *  the Free Software Foundation; either version 2 of the License, or
  11  *  (at your option) any later version.
  12  *
  13  *  This program is distributed in the hope that it will be useful,
  14  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  15  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  16  *  GNU General Public License for more details.
  17  *
  18  *  You should have received a copy of the GNU General Public License
  19  *  along with this program; if not, write to the Free Software
  20  *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  21  */
  22
  23 #include "inifcns.h"
  24 #include "ex.h"
  25 #include "constant.h"
  26 #include "lst.h"
  27 #include "fderivative.h"
  28 #include "matrix.h"
  29 #include "mul.h"
  30 #include "power.h"
  31 #include "operators.h"
  32 #include "relational.h"
  33 #include "pseries.h"
  34 #include "symbol.h"
  35 #include "symmetry.h"
  36 #include "utils.h"
  37
  38 #include <stdexcept>
  39 #include <vector>
  40
  41 namespace GiNaC {
  42
  43 //////////
  44 // complex conjugate
  45 //////////
  46
  47 static ex conjugate_evalf(const ex & arg)
  48 {
  49         if (is_exactly_a<numeric>(arg)) {
  50                 return ex_to<numeric>(arg).conjugate();
  51         }
  52         return conjugate_function(arg).hold();
  53 }
  54
  55 static ex conjugate_eval(const ex & arg)
  56 {
  57         return arg.conjugate();
  58 }
  59
  60 static void conjugate_print_latex(const ex & arg, const print_context & c)
  61 {
  62         c.s << "\\bar{"; arg.print(c); c.s << "}";
  63 }
  64
  65 static ex conjugate_conjugate(const ex & arg)
  66 {
  67         return arg;
  68 }
  69
  70 // If x is real then U.diff(x)-I*V.diff(x) represents both conjugate(U+I*V).diff(x)
  71 // and conjugate((U+I*V).diff(x))
  72 static ex conjugate_expl_derivative(const ex & arg, const symbol & s)
  73 {
  74         if (s.info(info_flags::real))
  75                 return conjugate(arg.diff(s));
  76         else {
  77                 exvector vec_arg;
  78                 vec_arg.push_back(arg);
  79                 return fderivative(ex_to<function>(conjugate(arg)).get_serial(),0,vec_arg).hold()*arg.diff(s);
  80         }
  81 }
  82
  83 static ex conjugate_real_part(const ex & arg)
  84 {
  85         return arg.real_part();
  86 }
  87
  88 static ex conjugate_imag_part(const ex & arg)
  89 {
  90         return -arg.imag_part();
  91 }
  92
  93 static bool func_arg_info(const ex & arg, unsigned inf)
  94 {
  95         // for some functions we can return the info() of its argument
  96         // (think of conjugate())
  97         switch (inf) {
  98                 case info_flags::polynomial:
  99                 case info_flags::integer_polynomial:
 100                 case info_flags::cinteger_polynomial:
 101                 case info_flags::rational_polynomial:
 102                 case info_flags::real:
 103                 case info_flags::rational:
 104                 case info_flags::integer:
 105                 case info_flags::crational:
 106                 case info_flags::cinteger:
 107                 case info_flags::even:
 108                 case info_flags::odd:
 109                 case info_flags::prime:
 110                 case info_flags::crational_polynomial:
 111                 case info_flags::rational_function:
 112                 case info_flags::algebraic:
 113                 case info_flags::positive:
 114                 case info_flags::negative:
 115                 case info_flags::nonnegative:
 116                 case info_flags::posint:
 117                 case info_flags::negint:
 118                 case info_flags::nonnegint:
 119                 case info_flags::has_indices:
 120                         return arg.info(inf);
 121         }
 122         return false;
 123 }
 124
 125 static bool conjugate_info(const ex & arg, unsigned inf)
 126 {
 127         return func_arg_info(arg, inf);
 128 }
 129
 130 REGISTER_FUNCTION(conjugate_function, eval_func(conjugate_eval).
 131                                       evalf_func(conjugate_evalf).
 132                                       expl_derivative_func(conjugate_expl_derivative).
 133                                       info_func(conjugate_info).
 134                                       print_func<print_latex>(conjugate_print_latex).
 135                                       conjugate_func(conjugate_conjugate).
 136                                       real_part_func(conjugate_real_part).
 137                                       imag_part_func(conjugate_imag_part).
 138                                       set_name("conjugate","conjugate"));
 139
 140 //////////
 141 // real part
 142 //////////
 143
 144 static ex real_part_evalf(const ex & arg)
 145 {
 146         if (is_exactly_a<numeric>(arg)) {
 147                 return ex_to<numeric>(arg).real();
 148         }
 149         return real_part_function(arg).hold();
 150 }
 151
 152 static ex real_part_eval(const ex & arg)
 153 {
 154         return arg.real_part();
 155 }
 156
 157 static void real_part_print_latex(const ex & arg, const print_context & c)
 158 {
 159         c.s << "\\Re"; arg.print(c); c.s << "";
 160 }
 161
 162 static ex real_part_conjugate(const ex & arg)
 163 {
 164         return real_part_function(arg).hold();
 165 }
 166
 167 static ex real_part_real_part(const ex & arg)
 168 {
 169         return real_part_function(arg).hold();
 170 }
 171
 172 static ex real_part_imag_part(const ex & arg)
 173 {
 174         return 0;
 175 }
 176
 177 // If x is real then Re(e).diff(x) is equal to Re(e.diff(x))
 178 static ex real_part_expl_derivative(const ex & arg, const symbol & s)
 179 {
 180         if (s.info(info_flags::real))
 181                 return real_part_function(arg.diff(s));
 182         else {
 183                 exvector vec_arg;
 184                 vec_arg.push_back(arg);
 185                 return fderivative(ex_to<function>(real_part(arg)).get_serial(),0,vec_arg).hold()*arg.diff(s);
 186         }
 187 }
 188
 189 REGISTER_FUNCTION(real_part_function, eval_func(real_part_eval).
 190                                       evalf_func(real_part_evalf).
 191                                       expl_derivative_func(real_part_expl_derivative).
 192                                       print_func<print_latex>(real_part_print_latex).
 193                                       conjugate_func(real_part_conjugate).
 194                                       real_part_func(real_part_real_part).
 195                                       imag_part_func(real_part_imag_part).
 196                                       set_name("real_part","real_part"));
 197
 198 //////////
 199 // imag part
 200 //////////
 201
 202 static ex imag_part_evalf(const ex & arg)
 203 {
 204         if (is_exactly_a<numeric>(arg)) {
 205                 return ex_to<numeric>(arg).imag();
 206         }
 207         return imag_part_function(arg).hold();
 208 }
 209
 210 static ex imag_part_eval(const ex & arg)
 211 {
 212         return arg.imag_part();
 213 }
 214
 215 static void imag_part_print_latex(const ex & arg, const print_context & c)
 216 {
 217         c.s << "\\Im"; arg.print(c); c.s << "";
 218 }
 219
 220 static ex imag_part_conjugate(const ex & arg)
 221 {
 222         return imag_part_function(arg).hold();
 223 }
 224
 225 static ex imag_part_real_part(const ex & arg)
 226 {
 227         return imag_part_function(arg).hold();
 228 }
 229
 230 static ex imag_part_imag_part(const ex & arg)
 231 {
 232         return 0;
 233 }
 234
 235 // If x is real then Im(e).diff(x) is equal to Im(e.diff(x))
 236 static ex imag_part_expl_derivative(const ex & arg, const symbol & s)
 237 {
 238         if (s.info(info_flags::real))
 239                 return imag_part_function(arg.diff(s));
 240         else {
 241                 exvector vec_arg;
 242                 vec_arg.push_back(arg);
 243                 return fderivative(ex_to<function>(imag_part(arg)).get_serial(),0,vec_arg).hold()*arg.diff(s);
 244         }
 245 }
 246
 247 REGISTER_FUNCTION(imag_part_function, eval_func(imag_part_eval).
 248                                       evalf_func(imag_part_evalf).
 249                                       expl_derivative_func(imag_part_expl_derivative).
 250                                       print_func<print_latex>(imag_part_print_latex).
 251                                       conjugate_func(imag_part_conjugate).
 252                                       real_part_func(imag_part_real_part).
 253                                       imag_part_func(imag_part_imag_part).
 254                                       set_name("imag_part","imag_part"));
 255
 256 //////////
 257 // absolute value
 258 //////////
 259
 260 static ex abs_evalf(const ex & arg)
 261 {
 262         if (is_exactly_a<numeric>(arg))
 263                 return abs(ex_to<numeric>(arg));
 264
 265         return abs(arg).hold();
 266 }
 267
 268 static ex abs_eval(const ex & arg)
 269 {
 270         if (is_exactly_a<numeric>(arg))
 271                 return abs(ex_to<numeric>(arg));
 272
 273         if (arg.info(info_flags::nonnegative))
 274                 return arg;
 275
 276         if (is_ex_the_function(arg, abs))
 277                 return arg;
 278
 279         if (is_ex_the_function(arg, exp))
 280                 return exp(arg.op(0).real_part());
 281
 282         if (is_exactly_a<power>(arg)) {
 283                 const ex& base = arg.op(0);
 284                 const ex& exponent = arg.op(1);
 285                 if (base.info(info_flags::positive) || exponent.info(info_flags::real))
 286                         return pow(abs(base), exponent.real_part());
 287         }
 288
 289         if (is_ex_the_function(arg, conjugate_function))
 290                 return abs(arg.op(0));
 291
 292         if (is_ex_the_function(arg, step))
 293                 return arg;
 294
 295         return abs(arg).hold();
 296 }
 297
 298 static ex abs_expand(const ex & arg, unsigned options)
 299 {
 300         if ((options & expand_options::expand_transcendental)
 301                 && is_exactly_a<mul>(arg)) {
 302                 exvector prodseq;
 303                 prodseq.reserve(arg.nops());
 304                 for (const_iterator i = arg.begin(); i != arg.end(); ++i) {
 305                         if (options & expand_options::expand_function_args)
 306                                 prodseq.push_back(abs(i->expand(options)));
 307                         else
 308                                 prodseq.push_back(abs(*i));
 309                 }
 310                 return (new mul(prodseq))->setflag(status_flags::dynallocated | status_flags::expanded);
 311         }
 312
 313         if (options & expand_options::expand_function_args)
 314                 return abs(arg.expand(options)).hold();
 315         else
 316                 return abs(arg).hold();
 317 }
 318
 319 static ex abs_expl_derivative(const ex & arg, const symbol & s)
 320 {
 321         ex diff_arg = arg.diff(s);
 322         return (diff_arg*arg.conjugate()+arg*diff_arg.conjugate())/2/abs(arg);
 323 }
 324
 325 static void abs_print_latex(const ex & arg, const print_context & c)
 326 {
 327         c.s << "{|"; arg.print(c); c.s << "|}";
 328 }
 329
 330 static void abs_print_csrc_float(const ex & arg, const print_context & c)
 331 {
 332         c.s << "fabs("; arg.print(c); c.s << ")";
 333 }
 334
 335 static ex abs_conjugate(const ex & arg)
 336 {
 337         return abs(arg).hold();
 338 }
 339
 340 static ex abs_real_part(const ex & arg)
 341 {
 342         return abs(arg).hold();
 343 }
 344
 345 static ex abs_imag_part(const ex& arg)
 346 {
 347         return 0;
 348 }
 349
 350 static ex abs_power(const ex & arg, const ex & exp)
 351 {
 352         if ((is_a<numeric>(exp) && ex_to<numeric>(exp).is_even()) || exp.info(info_flags::even)) {
 353                 if (arg.info(info_flags::real) || arg.is_equal(arg.conjugate()))
 354                         return power(arg, exp);
 355                 else
 356                         return power(arg, exp/2)*power(arg.conjugate(), exp/2);
 357         } else
 358                 return power(abs(arg), exp).hold();
 359 }
 360
 361 bool abs_info(const ex & arg, unsigned inf)
 362 {
 363         switch (inf) {
 364                 case info_flags::integer:
 365                 case info_flags::even:
 366                 case info_flags::odd:
 367                 case info_flags::prime:
 368                         return arg.info(inf);
 369                 case info_flags::nonnegint:
 370                         return arg.info(info_flags::integer);
 371                 case info_flags::nonnegative:
 372                 case info_flags::real:
 373                         return true;
 374                 case info_flags::negative:
 375                         return false;
 376                 case info_flags::positive:
 377                         return arg.info(info_flags::positive) || arg.info(info_flags::negative);
 378                 case info_flags::has_indices: {
 379                         if (arg.info(info_flags::has_indices))
 380                                 return true;
 381                         else
 382                                 return false;
 383                 }
 384         }
 385         return false;
 386 }
 387
 388 REGISTER_FUNCTION(abs, eval_func(abs_eval).
 389                        evalf_func(abs_evalf).
 390                        expand_func(abs_expand).
 391                        expl_derivative_func(abs_expl_derivative).
 392                        info_func(abs_info).
 393                        print_func<print_latex>(abs_print_latex).
 394                        print_func<print_csrc_float>(abs_print_csrc_float).
 395                        print_func<print_csrc_double>(abs_print_csrc_float).
 396                        conjugate_func(abs_conjugate).
 397                        real_part_func(abs_real_part).
 398                        imag_part_func(abs_imag_part).
 399                        power_func(abs_power));
 400
 401 //////////
 402 // Step function
 403 //////////
 404
 405 static ex step_evalf(const ex & arg)
 406 {
 407         if (is_exactly_a<numeric>(arg))
 408                 return step(ex_to<numeric>(arg));
 409
 410         return step(arg).hold();
 411 }
 412
 413 static ex step_eval(const ex & arg)
 414 {
 415         if (is_exactly_a<numeric>(arg))
 416                 return step(ex_to<numeric>(arg));
 417
 418         else if (is_exactly_a<mul>(arg) &&
 419                  is_exactly_a<numeric>(arg.op(arg.nops()-1))) {
 420                 numeric oc = ex_to<numeric>(arg.op(arg.nops()-1));
 421                 if (oc.is_real()) {
 422                         if (oc > 0)
 423                                 // step(42*x) -> step(x)
 424                                 return step(arg/oc).hold();
 425                         else
 426                                 // step(-42*x) -> step(-x)
 427                                 return step(-arg/oc).hold();
 428                 }
 429                 if (oc.real().is_zero()) {
 430                         if (oc.imag() > 0)
 431                                 // step(42*I*x) -> step(I*x)
 432                                 return step(I*arg/oc).hold();
 433                         else
 434                                 // step(-42*I*x) -> step(-I*x)
 435                                 return step(-I*arg/oc).hold();
 436                 }
 437         }
 438
 439         return step(arg).hold();
 440 }
 441
 442 static ex step_series(const ex & arg,
 443                       const relational & rel,
 444                       int order,
 445                       unsigned options)
 446 {
 447         const ex arg_pt = arg.subs(rel, subs_options::no_pattern);
 448         if (arg_pt.info(info_flags::numeric)
 449             && ex_to<numeric>(arg_pt).real().is_zero()
 450             && !(options & series_options::suppress_branchcut))
 451                 throw (std::domain_error("step_series(): on imaginary axis"));
 452
 453         epvector seq;
 454         seq.push_back(expair(step(arg_pt), _ex0));
 455         return pseries(rel,seq);
 456 }
 457
 458 static ex step_conjugate(const ex& arg)
 459 {
 460         return step(arg).hold();
 461 }
 462
 463 static ex step_real_part(const ex& arg)
 464 {
 465         return step(arg).hold();
 466 }
 467
 468 static ex step_imag_part(const ex& arg)
 469 {
 470         return 0;
 471 }
 472
 473 REGISTER_FUNCTION(step, eval_func(step_eval).
 474                         evalf_func(step_evalf).
 475                         series_func(step_series).
 476                         conjugate_func(step_conjugate).
 477                         real_part_func(step_real_part).
 478                         imag_part_func(step_imag_part));
 479
 480 //////////
 481 // Complex sign
 482 //////////
 483
 484 static ex csgn_evalf(const ex & arg)
 485 {
 486         if (is_exactly_a<numeric>(arg))
 487                 return csgn(ex_to<numeric>(arg));
 488
 489         return csgn(arg).hold();
 490 }
 491
 492 static ex csgn_eval(const ex & arg)
 493 {
 494         if (is_exactly_a<numeric>(arg))
 495                 return csgn(ex_to<numeric>(arg));
 496
 497         else if (is_exactly_a<mul>(arg) &&
 498                  is_exactly_a<numeric>(arg.op(arg.nops()-1))) {
 499                 numeric oc = ex_to<numeric>(arg.op(arg.nops()-1));
 500                 if (oc.is_real()) {
 501                         if (oc > 0)
 502                                 // csgn(42*x) -> csgn(x)
 503                                 return csgn(arg/oc).hold();
 504                         else
 505                                 // csgn(-42*x) -> -csgn(x)
 506                                 return -csgn(arg/oc).hold();
 507                 }
 508                 if (oc.real().is_zero()) {
 509                         if (oc.imag() > 0)
 510                                 // csgn(42*I*x) -> csgn(I*x)
 511                                 return csgn(I*arg/oc).hold();
 512                         else
 513                                 // csgn(-42*I*x) -> -csgn(I*x)
 514                                 return -csgn(I*arg/oc).hold();
 515                 }
 516         }
 517
 518         return csgn(arg).hold();
 519 }
 520
 521 static ex csgn_series(const ex & arg,
 522                       const relational & rel,
 523                       int order,
 524                       unsigned options)
 525 {
 526         const ex arg_pt = arg.subs(rel, subs_options::no_pattern);
 527         if (arg_pt.info(info_flags::numeric)
 528             && ex_to<numeric>(arg_pt).real().is_zero()
 529             && !(options & series_options::suppress_branchcut))
 530                 throw (std::domain_error("csgn_series(): on imaginary axis"));
 531
 532         epvector seq;
 533         seq.push_back(expair(csgn(arg_pt), _ex0));
 534         return pseries(rel,seq);
 535 }
 536
 537 static ex csgn_conjugate(const ex& arg)
 538 {
 539         return csgn(arg).hold();
 540 }
 541
 542 static ex csgn_real_part(const ex& arg)
 543 {
 544         return csgn(arg).hold();
 545 }
 546
 547 static ex csgn_imag_part(const ex& arg)
 548 {
 549         return 0;
 550 }
 551
 552 static ex csgn_power(const ex & arg, const ex & exp)
 553 {
 554         if (is_a<numeric>(exp) && exp.info(info_flags::positive) && ex_to<numeric>(exp).is_integer()) {
 555                 if (ex_to<numeric>(exp).is_odd())
 556                         return csgn(arg).hold();
 557                 else
 558                         return power(csgn(arg), _ex2).hold();
 559         } else
 560                 return power(csgn(arg), exp).hold();
 561 }
 562
 563
 564 REGISTER_FUNCTION(csgn, eval_func(csgn_eval).
 565                         evalf_func(csgn_evalf).
 566                         series_func(csgn_series).
 567                         conjugate_func(csgn_conjugate).
 568                         real_part_func(csgn_real_part).
 569                         imag_part_func(csgn_imag_part).
 570                         power_func(csgn_power));
 571
 572
 573 //////////
 574 // Eta function: eta(x,y) == log(x*y) - log(x) - log(y).
 575 // This function is closely related to the unwinding number K, sometimes found
 576 // in modern literature: K(z) == (z-log(exp(z)))/(2*Pi*I).
 577 //////////
 578
 579 static ex eta_evalf(const ex &x, const ex &y)
 580 {
 581         // It seems like we basically have to replicate the eval function here,
 582         // since the expression might not be fully evaluated yet.
 583         if (x.info(info_flags::positive) || y.info(info_flags::positive))
 584                 return _ex0;
 585
 586         if (x.info(info_flags::numeric) &&      y.info(info_flags::numeric)) {
 587                 const numeric nx = ex_to<numeric>(x);
 588                 const numeric ny = ex_to<numeric>(y);
 589                 const numeric nxy = ex_to<numeric>(x*y);
 590                 int cut = 0;
 591                 if (nx.is_real() && nx.is_negative())
 592                         cut -= 4;
 593                 if (ny.is_real() && ny.is_negative())
 594                         cut -= 4;
 595                 if (nxy.is_real() && nxy.is_negative())
 596                         cut += 4;
 597                 return evalf(I/4*Pi)*((csgn(-imag(nx))+1)*(csgn(-imag(ny))+1)*(csgn(imag(nxy))+1)-
 598                                       (csgn(imag(nx))+1)*(csgn(imag(ny))+1)*(csgn(-imag(nxy))+1)+cut);
 599         }
 600
 601         return eta(x,y).hold();
 602 }
 603
 604 static ex eta_eval(const ex &x, const ex &y)
 605 {
 606         // trivial:  eta(x,c) -> 0  if c is real and positive
 607         if (x.info(info_flags::positive) || y.info(info_flags::positive))
 608                 return _ex0;
 609
 610         if (x.info(info_flags::numeric) &&      y.info(info_flags::numeric)) {
 611                 // don't call eta_evalf here because it would call Pi.evalf()!
 612                 const numeric nx = ex_to<numeric>(x);
 613                 const numeric ny = ex_to<numeric>(y);
 614                 const numeric nxy = ex_to<numeric>(x*y);
 615                 int cut = 0;
 616                 if (nx.is_real() && nx.is_negative())
 617                         cut -= 4;
 618                 if (ny.is_real() && ny.is_negative())
 619                         cut -= 4;
 620                 if (nxy.is_real() && nxy.is_negative())
 621                         cut += 4;
 622                 return (I/4)*Pi*((csgn(-imag(nx))+1)*(csgn(-imag(ny))+1)*(csgn(imag(nxy))+1)-
 623                                  (csgn(imag(nx))+1)*(csgn(imag(ny))+1)*(csgn(-imag(nxy))+1)+cut);
 624         }
 625
 626         return eta(x,y).hold();
 627 }
 628
 629 static ex eta_series(const ex & x, const ex & y,
 630                      const relational & rel,
 631                      int order,
 632                      unsigned options)
 633 {
 634         const ex x_pt = x.subs(rel, subs_options::no_pattern);
 635         const ex y_pt = y.subs(rel, subs_options::no_pattern);
 636         if ((x_pt.info(info_flags::numeric) && x_pt.info(info_flags::negative)) ||
 637             (y_pt.info(info_flags::numeric) && y_pt.info(info_flags::negative)) ||
 638             ((x_pt*y_pt).info(info_flags::numeric) && (x_pt*y_pt).info(info_flags::negative)))
 639                         throw (std::domain_error("eta_series(): on discontinuity"));
 640         epvector seq;
 641         seq.push_back(expair(eta(x_pt,y_pt), _ex0));
 642         return pseries(rel,seq);
 643 }
 644
 645 static ex eta_conjugate(const ex & x, const ex & y)
 646 {
 647         return -eta(x, y).hold();
 648 }
 649
 650 static ex eta_real_part(const ex & x, const ex & y)
 651 {
 652         return 0;
 653 }
 654
 655 static ex eta_imag_part(const ex & x, const ex & y)
 656 {
 657         return -I*eta(x, y).hold();
 658 }
 659
 660 REGISTER_FUNCTION(eta, eval_func(eta_eval).
 661                        evalf_func(eta_evalf).
 662                        series_func(eta_series).
 663                        latex_name("\\eta").
 664                        set_symmetry(sy_symm(0, 1)).
 665                        conjugate_func(eta_conjugate).
 666                        real_part_func(eta_real_part).
 667                        imag_part_func(eta_imag_part));
 668
 669
 670 //////////
 671 // dilogarithm
 672 //////////
 673
 674 static ex Li2_evalf(const ex & x)
 675 {
 676         if (is_exactly_a<numeric>(x))
 677                 return Li2(ex_to<numeric>(x));
 678
 679         return Li2(x).hold();
 680 }
 681
 682 static ex Li2_eval(const ex & x)
 683 {
 684         if (x.info(info_flags::numeric)) {
 685                 // Li2(0) -> 0
 686                 if (x.is_zero())
 687                         return _ex0;
 688                 // Li2(1) -> Pi^2/6
 689                 if (x.is_equal(_ex1))
 690                         return power(Pi,_ex2)/_ex6;
 691                 // Li2(1/2) -> Pi^2/12 - log(2)^2/2
 692                 if (x.is_equal(_ex1_2))
 693                         return power(Pi,_ex2)/_ex12 + power(log(_ex2),_ex2)*_ex_1_2;
 694                 // Li2(-1) -> -Pi^2/12
 695                 if (x.is_equal(_ex_1))
 696                         return -power(Pi,_ex2)/_ex12;
 697                 // Li2(I) -> -Pi^2/48+Catalan*I
 698                 if (x.is_equal(I))
 699                         return power(Pi,_ex2)/_ex_48 + Catalan*I;
 700                 // Li2(-I) -> -Pi^2/48-Catalan*I
 701                 if (x.is_equal(-I))
 702                         return power(Pi,_ex2)/_ex_48 - Catalan*I;
 703                 // Li2(float)
 704                 if (!x.info(info_flags::crational))
 705                         return Li2(ex_to<numeric>(x));
 706         }
 707
 708         return Li2(x).hold();
 709 }
 710
 711 static ex Li2_deriv(const ex & x, unsigned deriv_param)
 712 {
 713         GINAC_ASSERT(deriv_param==0);
 714
 715         // d/dx Li2(x) -> -log(1-x)/x
 716         return -log(_ex1-x)/x;
 717 }
 718
 719 static ex Li2_series(const ex &x, const relational &rel, int order, unsigned options)
 720 {
 721         const ex x_pt = x.subs(rel, subs_options::no_pattern);
 722         if (x_pt.info(info_flags::numeric)) {
 723                 // First special case: x==0 (derivatives have poles)
 724                 if (x_pt.is_zero()) {
 725                         // method:
 726                         // The problem is that in d/dx Li2(x==0) == -log(1-x)/x we cannot
 727                         // simply substitute x==0.  The limit, however, exists: it is 1.
 728                         // We also know all higher derivatives' limits:
 729                         // (d/dx)^n Li2(x) == n!/n^2.
 730                         // So the primitive series expansion is
 731                         // Li2(x==0) == x + x^2/4 + x^3/9 + ...
 732                         // and so on.
 733                         // We first construct such a primitive series expansion manually in
 734                         // a dummy symbol s and then insert the argument's series expansion
 735                         // for s.  Reexpanding the resulting series returns the desired
 736                         // result.
 737                         const symbol s;
 738                         ex ser;
 739                         // manually construct the primitive expansion
 740                         for (int i=1; i<order; ++i)
 741                                 ser += pow(s,i) / pow(numeric(i), *_num2_p);
 742                         // substitute the argument's series expansion
 743                         ser = ser.subs(s==x.series(rel, order), subs_options::no_pattern);
 744                         // maybe that was terminating, so add a proper order term
 745                         epvector nseq;
 746                         nseq.push_back(expair(Order(_ex1), order));
 747                         ser += pseries(rel, nseq);
 748                         // reexpanding it will collapse the series again
 749                         return ser.series(rel, order);
 750                         // NB: Of course, this still does not allow us to compute anything
 751                         // like sin(Li2(x)).series(x==0,2), since then this code here is
 752                         // not reached and the derivative of sin(Li2(x)) doesn't allow the
 753                         // substitution x==0.  Probably limits *are* needed for the general
 754                         // cases.  In case L'Hospital's rule is implemented for limits and
 755                         // basic::series() takes care of this, this whole block is probably
 756                         // obsolete!
 757                 }
 758                 // second special case: x==1 (branch point)
 759                 if (x_pt.is_equal(_ex1)) {
 760                         // method:
 761                         // construct series manually in a dummy symbol s
 762                         const symbol s;
 763                         ex ser = zeta(_ex2);
 764                         // manually construct the primitive expansion
 765                         for (int i=1; i<order; ++i)
 766                                 ser += pow(1-s,i) * (numeric(1,i)*(I*Pi+log(s-1)) - numeric(1,i*i));
 767                         // substitute the argument's series expansion
 768                         ser = ser.subs(s==x.series(rel, order), subs_options::no_pattern);
 769                         // maybe that was terminating, so add a proper order term
 770                         epvector nseq;
 771                         nseq.push_back(expair(Order(_ex1), order));
 772                         ser += pseries(rel, nseq);
 773                         // reexpanding it will collapse the series again
 774                         return ser.series(rel, order);
 775                 }
 776                 // third special case: x real, >=1 (branch cut)
 777                 if (!(options & series_options::suppress_branchcut) &&
 778                         ex_to<numeric>(x_pt).is_real() && ex_to<numeric>(x_pt)>1) {
 779                         // method:
 780                         // This is the branch cut: assemble the primitive series manually
 781                         // and then add the corresponding complex step function.
 782                         const symbol &s = ex_to<symbol>(rel.lhs());
 783                         const ex point = rel.rhs();
 784                         const symbol foo;
 785                         epvector seq;
 786                         // zeroth order term:
 787                         seq.push_back(expair(Li2(x_pt), _ex0));
 788                         // compute the intermediate terms:
 789                         ex replarg = series(Li2(x), s==foo, order);
 790                         for (size_t i=1; i<replarg.nops()-1; ++i)
 791                                 seq.push_back(expair((replarg.op(i)/power(s-foo,i)).series(foo==point,1,options).op(0).subs(foo==s, subs_options::no_pattern),i));
 792                         // append an order term:
 793                         seq.push_back(expair(Order(_ex1), replarg.nops()-1));
 794                         return pseries(rel, seq);
 795                 }
 796         }
 797         // all other cases should be safe, by now:
 798         throw do_taylor();  // caught by function::series()
 799 }
 800
 801 static ex Li2_conjugate(const ex & x)
 802 {
 803         // conjugate(Li2(x))==Li2(conjugate(x)) unless on the branch cuts which
 804         // run along the positive real axis beginning at 1.
 805         if (x.info(info_flags::negative)) {
 806                 return Li2(x).hold();
 807         }
 808         if (is_exactly_a<numeric>(x) &&
 809             (!x.imag_part().is_zero() || x < *_num1_p)) {
 810                 return Li2(x.conjugate());
 811         }
 812         return conjugate_function(Li2(x)).hold();
 813 }
 814
 815 REGISTER_FUNCTION(Li2, eval_func(Li2_eval).
 816                        evalf_func(Li2_evalf).
 817                        derivative_func(Li2_deriv).
 818                        series_func(Li2_series).
 819                        conjugate_func(Li2_conjugate).
 820                        latex_name("\\mathrm{Li}_2"));
 821
 822 //////////
 823 // trilogarithm
 824 //////////
 825
 826 static ex Li3_eval(const ex & x)
 827 {
 828         if (x.is_zero())
 829                 return x;
 830         return Li3(x).hold();
 831 }
 832
 833 REGISTER_FUNCTION(Li3, eval_func(Li3_eval).
 834                        latex_name("\\mathrm{Li}_3"));
 835
 836 //////////
 837 // Derivatives of Riemann's Zeta-function  zetaderiv(0,x)==zeta(x)
 838 //////////
 839
 840 static ex zetaderiv_eval(const ex & n, const ex & x)
 841 {
 842         if (n.info(info_flags::numeric)) {
 843                 // zetaderiv(0,x) -> zeta(x)
 844                 if (n.is_zero())
 845                         return zeta(x).hold();
 846         }
 847
 848         return zetaderiv(n, x).hold();
 849 }
 850
 851 static ex zetaderiv_deriv(const ex & n, const ex & x, unsigned deriv_param)
 852 {
 853         GINAC_ASSERT(deriv_param<2);
 854
 855         if (deriv_param==0) {
 856                 // d/dn zeta(n,x)
 857                 throw(std::logic_error("cannot diff zetaderiv(n,x) with respect to n"));
 858         }
 859         // d/dx psi(n,x)
 860         return zetaderiv(n+1,x);
 861 }
 862
 863 REGISTER_FUNCTION(zetaderiv, eval_func(zetaderiv_eval).
 864                                  derivative_func(zetaderiv_deriv).
 865                                  latex_name("\\zeta^\\prime"));
 866
 867 //////////
 868 // factorial
 869 //////////
 870
 871 static ex factorial_evalf(const ex & x)
 872 {
 873         return factorial(x).hold();
 874 }
 875
 876 static ex factorial_eval(const ex & x)
 877 {
 878         if (is_exactly_a<numeric>(x))
 879                 return factorial(ex_to<numeric>(x));
 880         else
 881                 return factorial(x).hold();
 882 }
 883
 884 static void factorial_print_dflt_latex(const ex & x, const print_context & c)
 885 {
 886         if (is_exactly_a<symbol>(x) ||
 887             is_exactly_a<constant>(x) ||
 888                 is_exactly_a<function>(x)) {
 889                 x.print(c); c.s << "!";
 890         } else {
 891                 c.s << "("; x.print(c); c.s << ")!";
 892         }
 893 }
 894
 895 static ex factorial_conjugate(const ex & x)
 896 {
 897         return factorial(x).hold();
 898 }
 899
 900 static ex factorial_real_part(const ex & x)
 901 {
 902         return factorial(x).hold();
 903 }
 904
 905 static ex factorial_imag_part(const ex & x)
 906 {
 907         return 0;
 908 }
 909
 910 REGISTER_FUNCTION(factorial, eval_func(factorial_eval).
 911                              evalf_func(factorial_evalf).
 912                              print_func<print_dflt>(factorial_print_dflt_latex).
 913                              print_func<print_latex>(factorial_print_dflt_latex).
 914                              conjugate_func(factorial_conjugate).
 915                              real_part_func(factorial_real_part).
 916                              imag_part_func(factorial_imag_part));
 917
 918 //////////
 919 // binomial
 920 //////////
 921
 922 static ex binomial_evalf(const ex & x, const ex & y)
 923 {
 924         return binomial(x, y).hold();
 925 }
 926
 927 static ex binomial_sym(const ex & x, const numeric & y)
 928 {
 929         if (y.is_integer()) {
 930                 if (y.is_nonneg_integer()) {
 931                         const unsigned N = y.to_int();
 932                         if (N == 0) return _ex1;
 933                         if (N == 1) return x;
 934                         ex t = x.expand();
 935                         for (unsigned i = 2; i <= N; ++i)
 936                                 t = (t * (x + i - y - 1)).expand() / i;
 937                         return t;
 938                 } else
 939                         return _ex0;
 940         }
 941
 942         return binomial(x, y).hold();
 943 }
 944
 945 static ex binomial_eval(const ex & x, const ex &y)
 946 {
 947         if (is_exactly_a<numeric>(y)) {
 948                 if (is_exactly_a<numeric>(x) && ex_to<numeric>(x).is_integer())
 949                         return binomial(ex_to<numeric>(x), ex_to<numeric>(y));
 950                 else
 951                         return binomial_sym(x, ex_to<numeric>(y));
 952         } else
 953                 return binomial(x, y).hold();
 954 }
 955
 956 // At the moment the numeric evaluation of a binomail function always
 957 // gives a real number, but if this would be implemented using the gamma
 958 // function, also complex conjugation should be changed (or rather, deleted).
 959 static ex binomial_conjugate(const ex & x, const ex & y)
 960 {
 961         return binomial(x,y).hold();
 962 }
 963
 964 static ex binomial_real_part(const ex & x, const ex & y)
 965 {
 966         return binomial(x,y).hold();
 967 }
 968
 969 static ex binomial_imag_part(const ex & x, const ex & y)
 970 {
 971         return 0;
 972 }
 973
 974 REGISTER_FUNCTION(binomial, eval_func(binomial_eval).
 975                             evalf_func(binomial_evalf).
 976                             conjugate_func(binomial_conjugate).
 977                             real_part_func(binomial_real_part).
 978                             imag_part_func(binomial_imag_part));
 979
 980 //////////
 981 // Order term function (for truncated power series)
 982 //////////
 983
 984 static ex Order_eval(const ex & x)
 985 {
 986         if (is_exactly_a<numeric>(x)) {
 987                 // O(c) -> O(1) or 0
 988                 if (!x.is_zero())
 989                         return Order(_ex1).hold();
 990                 else
 991                         return _ex0;
 992         } else if (is_exactly_a<mul>(x)) {
 993                 const mul &m = ex_to<mul>(x);
 994                 // O(c*expr) -> O(expr)
 995                 if (is_exactly_a<numeric>(m.op(m.nops() - 1)))
 996                         return Order(x / m.op(m.nops() - 1)).hold();
 997         }
 998         return Order(x).hold();
 999 }
1000
1001 static ex Order_series(const ex & x, const relational & r, int order, unsigned options)
1002 {
1003         // Just wrap the function into a pseries object
1004         epvector new_seq;
1005         GINAC_ASSERT(is_a<symbol>(r.lhs()));
1006         const symbol &s = ex_to<symbol>(r.lhs());
1007         new_seq.push_back(expair(Order(_ex1), numeric(std::min(x.ldegree(s), order))));
1008         return pseries(r, new_seq);
1009 }
1010
1011 static ex Order_conjugate(const ex & x)
1012 {
1013         return Order(x).hold();
1014 }
1015
1016 static ex Order_real_part(const ex & x)
1017 {
1018         return Order(x).hold();
1019 }
1020
1021 static ex Order_imag_part(const ex & x)
1022 {
1023         if(x.info(info_flags::real))
1024                 return 0;
1025         return Order(x).hold();
1026 }
1027
1028 static ex Order_expl_derivative(const ex & arg, const symbol & s)
1029 {
1030         return Order(arg.diff(s));
1031 }
1032
1033 REGISTER_FUNCTION(Order, eval_func(Order_eval).
1034                          series_func(Order_series).
1035                          latex_name("\\mathcal{O}").
1036                          expl_derivative_func(Order_expl_derivative).
1037                          conjugate_func(Order_conjugate).
1038                          real_part_func(Order_real_part).
1039                          imag_part_func(Order_imag_part));
1040
1041 //////////
1042 // Solve linear system
1043 //////////
1044
1045 ex lsolve(const ex &eqns, const ex &symbols, unsigned options)
1046 {
1047         // solve a system of linear equations
1048         if (eqns.info(info_flags::relation_equal)) {
1049                 if (!symbols.info(info_flags::symbol))
1050                         throw(std::invalid_argument("lsolve(): 2nd argument must be a symbol"));
1051                 const ex sol = lsolve(lst(eqns),lst(symbols));
1052
1053                 GINAC_ASSERT(sol.nops()==1);
1054                 GINAC_ASSERT(is_exactly_a<relational>(sol.op(0)));
1055
1056                 return sol.op(0).op(1); // return rhs of first solution
1057         }
1058
1059         // syntax checks
1060         if (!eqns.info(info_flags::list)) {
1061                 throw(std::invalid_argument("lsolve(): 1st argument must be a list or an equation"));
1062         }
1063         for (size_t i=0; i<eqns.nops(); i++) {
1064                 if (!eqns.op(i).info(info_flags::relation_equal)) {
1065                         throw(std::invalid_argument("lsolve(): 1st argument must be a list of equations"));
1066                 }
1067         }
1068         if (!symbols.info(info_flags::list)) {
1069                 throw(std::invalid_argument("lsolve(): 2nd argument must be a list or a symbol"));
1070         }
1071         for (size_t i=0; i<symbols.nops(); i++) {
1072                 if (!symbols.op(i).info(info_flags::symbol)) {
1073                         throw(std::invalid_argument("lsolve(): 2nd argument must be a list of symbols"));
1074                 }
1075         }
1076
1077         // build matrix from equation system
1078         matrix sys(eqns.nops(),symbols.nops());
1079         matrix rhs(eqns.nops(),1);
1080         matrix vars(symbols.nops(),1);
1081
1082         for (size_t r=0; r<eqns.nops(); r++) {
1083                 const ex eq = eqns.op(r).op(0)-eqns.op(r).op(1); // lhs-rhs==0
1084                 ex linpart = eq;
1085                 for (size_t c=0; c<symbols.nops(); c++) {
1086                         const ex co = eq.coeff(ex_to<symbol>(symbols.op(c)),1);
1087                         linpart -= co*symbols.op(c);
1088                         sys(r,c) = co;
1089                 }
1090                 linpart = linpart.expand();
1091                 rhs(r,0) = -linpart;
1092         }
1093
1094         // test if system is linear and fill vars matrix
1095         for (size_t i=0; i<symbols.nops(); i++) {
1096                 vars(i,0) = symbols.op(i);
1097                 if (sys.has(symbols.op(i)))
1098                         throw(std::logic_error("lsolve: system is not linear"));
1099                 if (rhs.has(symbols.op(i)))
1100                         throw(std::logic_error("lsolve: system is not linear"));
1101         }
1102
1103         matrix solution;
1104         try {
1105                 solution = sys.solve(vars,rhs,options);
1106         } catch (const std::runtime_error & e) {
1107                 // Probably singular matrix or otherwise overdetermined system:
1108                 // It is consistent to return an empty list
1109                 return lst();
1110         }
1111         GINAC_ASSERT(solution.cols()==1);
1112         GINAC_ASSERT(solution.rows()==symbols.nops());
1113
1114         // return list of equations of the form lst(var1==sol1,var2==sol2,...)
1115         lst sollist;
1116         for (size_t i=0; i<symbols.nops(); i++)
1117                 sollist.append(symbols.op(i)==solution(i,0));
1118
1119         return sollist;
1120 }
1121
1122 //////////
1123 // Find real root of f(x) numerically
1124 //////////
1125
1126 const numeric
1127 fsolve(const ex& f_in, const symbol& x, const numeric& x1, const numeric& x2)
1128 {
1129         if (!x1.is_real() || !x2.is_real()) {
1130                 throw std::runtime_error("fsolve(): interval not bounded by real numbers");
1131         }
1132         if (x1==x2) {
1133                 throw std::runtime_error("fsolve(): vanishing interval");
1134         }
1135         // xx[0] == left interval limit, xx[1] == right interval limit.
1136         // fx[0] == f(xx[0]), fx[1] == f(xx[1]).
1137         // We keep the root bracketed: xx[0]<xx[1] and fx[0]*fx[1]<0.
1138         numeric xx[2] = { x1<x2 ? x1 : x2,
1139                           x1<x2 ? x2 : x1 };
1140         ex f;
1141         if (is_a<relational>(f_in)) {
1142                 f = f_in.lhs()-f_in.rhs();
1143         } else {
1144                 f = f_in;
1145         }
1146         const ex fx_[2] = { f.subs(x==xx[0]).evalf(),
1147                             f.subs(x==xx[1]).evalf() };
1148         if (!is_a<numeric>(fx_[0]) || !is_a<numeric>(fx_[1])) {
1149                 throw std::runtime_error("fsolve(): function does not evaluate numerically");
1150         }
1151         numeric fx[2] = { ex_to<numeric>(fx_[0]),
1152                           ex_to<numeric>(fx_[1]) };
1153         if (!fx[0].is_real() || !fx[1].is_real()) {
1154                 throw std::runtime_error("fsolve(): function evaluates to complex values at interval boundaries");
1155         }
1156         if (fx[0]*fx[1]>=0) {
1157                 throw std::runtime_error("fsolve(): function does not change sign at interval boundaries");
1158         }
1159
1160         // The Newton-Raphson method has quadratic convergence!  Simply put, it
1161         // replaces x with x-f(x)/f'(x) at each step.  -f/f' is the delta:
1162         const ex ff = normal(-f/f.diff(x));
1163         int side = 0;  // Start at left interval limit.
1164         numeric xxprev;
1165         numeric fxprev;
1166         do {
1167                 xxprev = xx[side];
1168                 fxprev = fx[side];
1169                 ex dx_ = ff.subs(x == xx[side]).evalf();
1170                 if (!is_a<numeric>(dx_))
1171                         throw std::runtime_error("fsolve(): function derivative does not evaluate numerically");
1172                 xx[side] += ex_to<numeric>(dx_);
1173                 // Now check if Newton-Raphson method shot out of the interval
1174                 bool bad_shot = (side == 0 && xx[0] < xxprev) ||
1175                                 (side == 1 && xx[1] > xxprev) || xx[0] > xx[1];
1176                 if (!bad_shot) {
1177                         // Compute f(x) only if new x is inside the interval.
1178                         // The function might be difficult to compute numerically
1179                         // or even ill defined outside the interval. Also it's
1180                         // a small optimization.
1181                         ex f_x = f.subs(x == xx[side]).evalf();
1182                         if (!is_a<numeric>(f_x))
1183                                 throw std::runtime_error("fsolve(): function does not evaluate numerically");
1184                         fx[side] = ex_to<numeric>(f_x);
1185                 }
1186                 if (bad_shot) {
1187                         // Oops, Newton-Raphson method shot out of the interval.
1188                         // Restore, and try again with the other side instead!
1189                         xx[side] = xxprev;
1190                         fx[side] = fxprev;
1191                         side = !side;
1192                         xxprev = xx[side];
1193                         fxprev = fx[side];
1194
1195                         ex dx_ = ff.subs(x == xx[side]).evalf();
1196                         if (!is_a<numeric>(dx_))
1197                                 throw std::runtime_error("fsolve(): function derivative does not evaluate numerically [2]");
1198                         xx[side] += ex_to<numeric>(dx_);
1199
1200                         ex f_x = f.subs(x==xx[side]).evalf();
1201                         if (!is_a<numeric>(f_x))
1202                                 throw std::runtime_error("fsolve(): function does not evaluate numerically [2]");
1203                         fx[side] = ex_to<numeric>(f_x);
1204                 }
1205                 if ((fx[side]<0 && fx[!side]<0) || (fx[side]>0 && fx[!side]>0)) {
1206                         // Oops, the root isn't bracketed any more.
1207                         // Restore, and perform a bisection!
1208                         xx[side] = xxprev;
1209                         fx[side] = fxprev;
1210
1211                         // Ah, the bisection! Bisections converge linearly. Unfortunately,
1212                         // they occur pretty often when Newton-Raphson arrives at an x too
1213                         // close to the result on one side of the interval and
1214                         // f(x-f(x)/f'(x)) turns out to have the same sign as f(x) due to
1215                         // precision errors! Recall that this function does not have a
1216                         // precision goal as one of its arguments but instead relies on
1217                         // x converging to a fixed point. We speed up the (safe but slow)
1218                         // bisection method by mixing in a dash of the (unsafer but faster)
1219                         // secant method: Instead of splitting the interval at the
1220                         // arithmetic mean (bisection), we split it nearer to the root as
1221                         // determined by the secant between the values xx[0] and xx[1].
1222                         // Don't set the secant_weight to one because that could disturb
1223                         // the convergence in some corner cases!
1224                         static const double secant_weight = 0.984375;  // == 63/64 < 1
1225                         numeric xxmid = (1-secant_weight)*0.5*(xx[0]+xx[1])
1226                             + secant_weight*(xx[0]+fx[0]*(xx[0]-xx[1])/(fx[1]-fx[0]));
1227                         ex fxmid_ = f.subs(x == xxmid).evalf();
1228                         if (!is_a<numeric>(fxmid_))
1229                                 throw std::runtime_error("fsolve(): function does not evaluate numerically [3]");
1230                         numeric fxmid = ex_to<numeric>(fxmid_);
1231                         if (fxmid.is_zero()) {
1232                                 // Luck strikes...
1233                                 return xxmid;
1234                         }
1235                         if ((fxmid<0 && fx[side]>0) || (fxmid>0 && fx[side]<0)) {
1236                                 side = !side;
1237                         }
1238                         xxprev = xx[side];
1239                         fxprev = fx[side];
1240                         xx[side] = xxmid;
1241                         fx[side] = fxmid;
1242                 }
1243         } while (xxprev!=xx[side]);
1244         return xxprev;
1245 }
1246
1247
1248 /* Force inclusion of functions from inifcns_gamma and inifcns_zeta
1249  * for static lib (so ginsh will see them). */
1250 unsigned force_include_tgamma = tgamma_SERIAL::serial;
1251 unsigned force_include_zeta1 = zeta1_SERIAL::serial;
1252
1253 } // namespace GiNaC