850ed04a5515e666f18ca3e42aa25a7957196bf7
[ginac.git] / 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 (arg.info(info_flags::negative) || (-arg).info(info_flags::nonnegative))
277                 return -arg;
278
279         if (is_ex_the_function(arg, abs))
280                 return arg;
281
282         if (is_ex_the_function(arg, exp))
283                 return exp(arg.op(0).real_part());
284
285         if (is_exactly_a<power>(arg)) {
286                 const ex& base = arg.op(0);
287                 const ex& exponent = arg.op(1);
288                 if (base.info(info_flags::positive) || exponent.info(info_flags::real))
289                         return pow(abs(base), exponent.real_part());
290         }
291
292         if (is_ex_the_function(arg, conjugate_function))
293                 return abs(arg.op(0));
294
295         if (is_ex_the_function(arg, step))
296                 return arg;
297
298         return abs(arg).hold();
299 }
300
301 static ex abs_expand(const ex & arg, unsigned options)
302 {
303         if ((options & expand_options::expand_transcendental)
304                 && is_exactly_a<mul>(arg)) {
305                 exvector prodseq;
306                 prodseq.reserve(arg.nops());
307                 for (const_iterator i = arg.begin(); i != arg.end(); ++i) {
308                         if (options & expand_options::expand_function_args)
309                                 prodseq.push_back(abs(i->expand(options)));
310                         else
311                                 prodseq.push_back(abs(*i));
312                 }
313                 return dynallocate<mul>(prodseq).setflag(status_flags::expanded);
314         }
315
316         if (options & expand_options::expand_function_args)
317                 return abs(arg.expand(options)).hold();
318         else
319                 return abs(arg).hold();
320 }
321
322 static ex abs_expl_derivative(const ex & arg, const symbol & s)
323 {
324         ex diff_arg = arg.diff(s);
325         return (diff_arg*arg.conjugate()+arg*diff_arg.conjugate())/2/abs(arg);
326 }
327
328 static void abs_print_latex(const ex & arg, const print_context & c)
329 {
330         c.s << "{|"; arg.print(c); c.s << "|}";
331 }
332
333 static void abs_print_csrc_float(const ex & arg, const print_context & c)
334 {
335         c.s << "fabs("; arg.print(c); c.s << ")";
336 }
337
338 static ex abs_conjugate(const ex & arg)
339 {
340         return abs(arg).hold();
341 }
342
343 static ex abs_real_part(const ex & arg)
344 {
345         return abs(arg).hold();
346 }
347
348 static ex abs_imag_part(const ex& arg)
349 {
350         return 0;
351 }
352
353 static ex abs_power(const ex & arg, const ex & exp)
354 {
355         if ((is_a<numeric>(exp) && ex_to<numeric>(exp).is_even()) || exp.info(info_flags::even)) {
356                 if (arg.info(info_flags::real) || arg.is_equal(arg.conjugate()))
357                         return power(arg, exp);
358                 else
359                         return power(arg, exp/2)*power(arg.conjugate(), exp/2);
360         } else
361                 return power(abs(arg), exp).hold();
362 }
363
364 bool abs_info(const ex & arg, unsigned inf)
365 {
366         switch (inf) {
367                 case info_flags::integer:
368                 case info_flags::even:
369                 case info_flags::odd:
370                 case info_flags::prime:
371                         return arg.info(inf);
372                 case info_flags::nonnegint:
373                         return arg.info(info_flags::integer);
374                 case info_flags::nonnegative:
375                 case info_flags::real:
376                         return true;
377                 case info_flags::negative:
378                         return false;
379                 case info_flags::positive:
380                         return arg.info(info_flags::positive) || arg.info(info_flags::negative);
381                 case info_flags::has_indices: {
382                         if (arg.info(info_flags::has_indices))
383                                 return true;
384                         else
385                                 return false;
386                 }
387         }
388         return false;
389 }
390
391 REGISTER_FUNCTION(abs, eval_func(abs_eval).
392                        evalf_func(abs_evalf).
393                        expand_func(abs_expand).
394                        expl_derivative_func(abs_expl_derivative).
395                        info_func(abs_info).
396                        print_func<print_latex>(abs_print_latex).
397                        print_func<print_csrc_float>(abs_print_csrc_float).
398                        print_func<print_csrc_double>(abs_print_csrc_float).
399                        conjugate_func(abs_conjugate).
400                        real_part_func(abs_real_part).
401                        imag_part_func(abs_imag_part).
402                        power_func(abs_power));
403
404 //////////
405 // Step function
406 //////////
407
408 static ex step_evalf(const ex & arg)
409 {
410         if (is_exactly_a<numeric>(arg))
411                 return step(ex_to<numeric>(arg));
412         
413         return step(arg).hold();
414 }
415
416 static ex step_eval(const ex & arg)
417 {
418         if (is_exactly_a<numeric>(arg))
419                 return step(ex_to<numeric>(arg));
420         
421         else if (is_exactly_a<mul>(arg) &&
422                  is_exactly_a<numeric>(arg.op(arg.nops()-1))) {
423                 numeric oc = ex_to<numeric>(arg.op(arg.nops()-1));
424                 if (oc.is_real()) {
425                         if (oc > 0)
426                                 // step(42*x) -> step(x)
427                                 return step(arg/oc).hold();
428                         else
429                                 // step(-42*x) -> step(-x)
430                                 return step(-arg/oc).hold();
431                 }
432                 if (oc.real().is_zero()) {
433                         if (oc.imag() > 0)
434                                 // step(42*I*x) -> step(I*x)
435                                 return step(I*arg/oc).hold();
436                         else
437                                 // step(-42*I*x) -> step(-I*x)
438                                 return step(-I*arg/oc).hold();
439                 }
440         }
441         
442         return step(arg).hold();
443 }
444
445 static ex step_series(const ex & arg,
446                       const relational & rel,
447                       int order,
448                       unsigned options)
449 {
450         const ex arg_pt = arg.subs(rel, subs_options::no_pattern);
451         if (arg_pt.info(info_flags::numeric)
452             && ex_to<numeric>(arg_pt).real().is_zero()
453             && !(options & series_options::suppress_branchcut))
454                 throw (std::domain_error("step_series(): on imaginary axis"));
455         
456         epvector seq { expair(step(arg_pt), _ex0) };
457         return pseries(rel, std::move(seq));
458 }
459
460 static ex step_conjugate(const ex& arg)
461 {
462         return step(arg).hold();
463 }
464
465 static ex step_real_part(const ex& arg)
466 {
467         return step(arg).hold();
468 }
469
470 static ex step_imag_part(const ex& arg)
471 {
472         return 0;
473 }
474
475 REGISTER_FUNCTION(step, eval_func(step_eval).
476                         evalf_func(step_evalf).
477                         series_func(step_series).
478                         conjugate_func(step_conjugate).
479                         real_part_func(step_real_part).
480                         imag_part_func(step_imag_part));
481
482 //////////
483 // Complex sign
484 //////////
485
486 static ex csgn_evalf(const ex & arg)
487 {
488         if (is_exactly_a<numeric>(arg))
489                 return csgn(ex_to<numeric>(arg));
490         
491         return csgn(arg).hold();
492 }
493
494 static ex csgn_eval(const ex & arg)
495 {
496         if (is_exactly_a<numeric>(arg))
497                 return csgn(ex_to<numeric>(arg));
498         
499         else if (is_exactly_a<mul>(arg) &&
500                  is_exactly_a<numeric>(arg.op(arg.nops()-1))) {
501                 numeric oc = ex_to<numeric>(arg.op(arg.nops()-1));
502                 if (oc.is_real()) {
503                         if (oc > 0)
504                                 // csgn(42*x) -> csgn(x)
505                                 return csgn(arg/oc).hold();
506                         else
507                                 // csgn(-42*x) -> -csgn(x)
508                                 return -csgn(arg/oc).hold();
509                 }
510                 if (oc.real().is_zero()) {
511                         if (oc.imag() > 0)
512                                 // csgn(42*I*x) -> csgn(I*x)
513                                 return csgn(I*arg/oc).hold();
514                         else
515                                 // csgn(-42*I*x) -> -csgn(I*x)
516                                 return -csgn(I*arg/oc).hold();
517                 }
518         }
519         
520         return csgn(arg).hold();
521 }
522
523 static ex csgn_series(const ex & arg,
524                       const relational & rel,
525                       int order,
526                       unsigned options)
527 {
528         const ex arg_pt = arg.subs(rel, subs_options::no_pattern);
529         if (arg_pt.info(info_flags::numeric)
530             && ex_to<numeric>(arg_pt).real().is_zero()
531             && !(options & series_options::suppress_branchcut))
532                 throw (std::domain_error("csgn_series(): on imaginary axis"));
533         
534         epvector seq { expair(csgn(arg_pt), _ex0) };
535         return pseries(rel, std::move(seq));
536 }
537
538 static ex csgn_conjugate(const ex& arg)
539 {
540         return csgn(arg).hold();
541 }
542
543 static ex csgn_real_part(const ex& arg)
544 {
545         return csgn(arg).hold();
546 }
547
548 static ex csgn_imag_part(const ex& arg)
549 {
550         return 0;
551 }
552
553 static ex csgn_power(const ex & arg, const ex & exp)
554 {
555         if (is_a<numeric>(exp) && exp.info(info_flags::positive) && ex_to<numeric>(exp).is_integer()) {
556                 if (ex_to<numeric>(exp).is_odd())
557                         return csgn(arg).hold();
558                 else
559                         return power(csgn(arg), _ex2).hold();
560         } else
561                 return power(csgn(arg), exp).hold();
562 }
563
564
565 REGISTER_FUNCTION(csgn, eval_func(csgn_eval).
566                         evalf_func(csgn_evalf).
567                         series_func(csgn_series).
568                         conjugate_func(csgn_conjugate).
569                         real_part_func(csgn_real_part).
570                         imag_part_func(csgn_imag_part).
571                         power_func(csgn_power));
572
573
574 //////////
575 // Eta function: eta(x,y) == log(x*y) - log(x) - log(y).
576 // This function is closely related to the unwinding number K, sometimes found
577 // in modern literature: K(z) == (z-log(exp(z)))/(2*Pi*I).
578 //////////
579
580 static ex eta_evalf(const ex &x, const ex &y)
581 {
582         // It seems like we basically have to replicate the eval function here,
583         // since the expression might not be fully evaluated yet.
584         if (x.info(info_flags::positive) || y.info(info_flags::positive))
585                 return _ex0;
586
587         if (x.info(info_flags::numeric) &&      y.info(info_flags::numeric)) {
588                 const numeric nx = ex_to<numeric>(x);
589                 const numeric ny = ex_to<numeric>(y);
590                 const numeric nxy = ex_to<numeric>(x*y);
591                 int cut = 0;
592                 if (nx.is_real() && nx.is_negative())
593                         cut -= 4;
594                 if (ny.is_real() && ny.is_negative())
595                         cut -= 4;
596                 if (nxy.is_real() && nxy.is_negative())
597                         cut += 4;
598                 return evalf(I/4*Pi)*((csgn(-imag(nx))+1)*(csgn(-imag(ny))+1)*(csgn(imag(nxy))+1)-
599                                       (csgn(imag(nx))+1)*(csgn(imag(ny))+1)*(csgn(-imag(nxy))+1)+cut);
600         }
601
602         return eta(x,y).hold();
603 }
604
605 static ex eta_eval(const ex &x, const ex &y)
606 {
607         // trivial:  eta(x,c) -> 0  if c is real and positive
608         if (x.info(info_flags::positive) || y.info(info_flags::positive))
609                 return _ex0;
610
611         if (x.info(info_flags::numeric) &&      y.info(info_flags::numeric)) {
612                 // don't call eta_evalf here because it would call Pi.evalf()!
613                 const numeric nx = ex_to<numeric>(x);
614                 const numeric ny = ex_to<numeric>(y);
615                 const numeric nxy = ex_to<numeric>(x*y);
616                 int cut = 0;
617                 if (nx.is_real() && nx.is_negative())
618                         cut -= 4;
619                 if (ny.is_real() && ny.is_negative())
620                         cut -= 4;
621                 if (nxy.is_real() && nxy.is_negative())
622                         cut += 4;
623                 return (I/4)*Pi*((csgn(-imag(nx))+1)*(csgn(-imag(ny))+1)*(csgn(imag(nxy))+1)-
624                                  (csgn(imag(nx))+1)*(csgn(imag(ny))+1)*(csgn(-imag(nxy))+1)+cut);
625         }
626         
627         return eta(x,y).hold();
628 }
629
630 static ex eta_series(const ex & x, const ex & y,
631                      const relational & rel,
632                      int order,
633                      unsigned options)
634 {
635         const ex x_pt = x.subs(rel, subs_options::no_pattern);
636         const ex y_pt = y.subs(rel, subs_options::no_pattern);
637         if ((x_pt.info(info_flags::numeric) && x_pt.info(info_flags::negative)) ||
638             (y_pt.info(info_flags::numeric) && y_pt.info(info_flags::negative)) ||
639             ((x_pt*y_pt).info(info_flags::numeric) && (x_pt*y_pt).info(info_flags::negative)))
640                         throw (std::domain_error("eta_series(): on discontinuity"));
641         epvector seq { expair(eta(x_pt,y_pt), _ex0) };
642         return pseries(rel, std::move(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 { expair(Order(_ex1), order) };
746                         ser += pseries(rel, std::move(nseq));
747                         // reexpanding it will collapse the series again
748                         return ser.series(rel, order);
749                         // NB: Of course, this still does not allow us to compute anything
750                         // like sin(Li2(x)).series(x==0,2), since then this code here is
751                         // not reached and the derivative of sin(Li2(x)) doesn't allow the
752                         // substitution x==0.  Probably limits *are* needed for the general
753                         // cases.  In case L'Hospital's rule is implemented for limits and
754                         // basic::series() takes care of this, this whole block is probably
755                         // obsolete!
756                 }
757                 // second special case: x==1 (branch point)
758                 if (x_pt.is_equal(_ex1)) {
759                         // method:
760                         // construct series manually in a dummy symbol s
761                         const symbol s;
762                         ex ser = zeta(_ex2);
763                         // manually construct the primitive expansion
764                         for (int i=1; i<order; ++i)
765                                 ser += pow(1-s,i) * (numeric(1,i)*(I*Pi+log(s-1)) - numeric(1,i*i));
766                         // substitute the argument's series expansion
767                         ser = ser.subs(s==x.series(rel, order), subs_options::no_pattern);
768                         // maybe that was terminating, so add a proper order term
769                         epvector nseq { expair(Order(_ex1), order) };
770                         ser += pseries(rel, std::move(nseq));
771                         // reexpanding it will collapse the series again
772                         return ser.series(rel, order);
773                 }
774                 // third special case: x real, >=1 (branch cut)
775                 if (!(options & series_options::suppress_branchcut) &&
776                         ex_to<numeric>(x_pt).is_real() && ex_to<numeric>(x_pt)>1) {
777                         // method:
778                         // This is the branch cut: assemble the primitive series manually
779                         // and then add the corresponding complex step function.
780                         const symbol &s = ex_to<symbol>(rel.lhs());
781                         const ex point = rel.rhs();
782                         const symbol foo;
783                         epvector seq;
784                         // zeroth order term:
785                         seq.push_back(expair(Li2(x_pt), _ex0));
786                         // compute the intermediate terms:
787                         ex replarg = series(Li2(x), s==foo, order);
788                         for (size_t i=1; i<replarg.nops()-1; ++i)
789                                 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));
790                         // append an order term:
791                         seq.push_back(expair(Order(_ex1), replarg.nops()-1));
792                         return pseries(rel, std::move(seq));
793                 }
794         }
795         // all other cases should be safe, by now:
796         throw do_taylor();  // caught by function::series()
797 }
798
799 static ex Li2_conjugate(const ex & x)
800 {
801         // conjugate(Li2(x))==Li2(conjugate(x)) unless on the branch cuts which
802         // run along the positive real axis beginning at 1.
803         if (x.info(info_flags::negative)) {
804                 return Li2(x).hold();
805         }
806         if (is_exactly_a<numeric>(x) &&
807             (!x.imag_part().is_zero() || x < *_num1_p)) {
808                 return Li2(x.conjugate());
809         }
810         return conjugate_function(Li2(x)).hold();
811 }
812
813 REGISTER_FUNCTION(Li2, eval_func(Li2_eval).
814                        evalf_func(Li2_evalf).
815                        derivative_func(Li2_deriv).
816                        series_func(Li2_series).
817                        conjugate_func(Li2_conjugate).
818                        latex_name("\\mathrm{Li}_2"));
819
820 //////////
821 // trilogarithm
822 //////////
823
824 static ex Li3_eval(const ex & x)
825 {
826         if (x.is_zero())
827                 return x;
828         return Li3(x).hold();
829 }
830
831 REGISTER_FUNCTION(Li3, eval_func(Li3_eval).
832                        latex_name("\\mathrm{Li}_3"));
833
834 //////////
835 // Derivatives of Riemann's Zeta-function  zetaderiv(0,x)==zeta(x)
836 //////////
837
838 static ex zetaderiv_eval(const ex & n, const ex & x)
839 {
840         if (n.info(info_flags::numeric)) {
841                 // zetaderiv(0,x) -> zeta(x)
842                 if (n.is_zero())
843                         return zeta(x).hold();
844         }
845         
846         return zetaderiv(n, x).hold();
847 }
848
849 static ex zetaderiv_deriv(const ex & n, const ex & x, unsigned deriv_param)
850 {
851         GINAC_ASSERT(deriv_param<2);
852         
853         if (deriv_param==0) {
854                 // d/dn zeta(n,x)
855                 throw(std::logic_error("cannot diff zetaderiv(n,x) with respect to n"));
856         }
857         // d/dx psi(n,x)
858         return zetaderiv(n+1,x);
859 }
860
861 REGISTER_FUNCTION(zetaderiv, eval_func(zetaderiv_eval).
862                                  derivative_func(zetaderiv_deriv).
863                                  latex_name("\\zeta^\\prime"));
864
865 //////////
866 // factorial
867 //////////
868
869 static ex factorial_evalf(const ex & x)
870 {
871         return factorial(x).hold();
872 }
873
874 static ex factorial_eval(const ex & x)
875 {
876         if (is_exactly_a<numeric>(x))
877                 return factorial(ex_to<numeric>(x));
878         else
879                 return factorial(x).hold();
880 }
881
882 static void factorial_print_dflt_latex(const ex & x, const print_context & c)
883 {
884         if (is_exactly_a<symbol>(x) ||
885             is_exactly_a<constant>(x) ||
886                 is_exactly_a<function>(x)) {
887                 x.print(c); c.s << "!";
888         } else {
889                 c.s << "("; x.print(c); c.s << ")!";
890         }
891 }
892
893 static ex factorial_conjugate(const ex & x)
894 {
895         return factorial(x).hold();
896 }
897
898 static ex factorial_real_part(const ex & x)
899 {
900         return factorial(x).hold();
901 }
902
903 static ex factorial_imag_part(const ex & x)
904 {
905         return 0;
906 }
907
908 REGISTER_FUNCTION(factorial, eval_func(factorial_eval).
909                              evalf_func(factorial_evalf).
910                              print_func<print_dflt>(factorial_print_dflt_latex).
911                              print_func<print_latex>(factorial_print_dflt_latex).
912                              conjugate_func(factorial_conjugate).
913                              real_part_func(factorial_real_part).
914                              imag_part_func(factorial_imag_part));
915
916 //////////
917 // binomial
918 //////////
919
920 static ex binomial_evalf(const ex & x, const ex & y)
921 {
922         return binomial(x, y).hold();
923 }
924
925 static ex binomial_sym(const ex & x, const numeric & y)
926 {
927         if (y.is_integer()) {
928                 if (y.is_nonneg_integer()) {
929                         const unsigned N = y.to_int();
930                         if (N == 0) return _ex1;
931                         if (N == 1) return x;
932                         ex t = x.expand();
933                         for (unsigned i = 2; i <= N; ++i)
934                                 t = (t * (x + i - y - 1)).expand() / i;
935                         return t;
936                 } else
937                         return _ex0;
938         }
939
940         return binomial(x, y).hold();
941 }
942
943 static ex binomial_eval(const ex & x, const ex &y)
944 {
945         if (is_exactly_a<numeric>(y)) {
946                 if (is_exactly_a<numeric>(x) && ex_to<numeric>(x).is_integer())
947                         return binomial(ex_to<numeric>(x), ex_to<numeric>(y));
948                 else
949                         return binomial_sym(x, ex_to<numeric>(y));
950         } else
951                 return binomial(x, y).hold();
952 }
953
954 // At the moment the numeric evaluation of a binomial function always
955 // gives a real number, but if this would be implemented using the gamma
956 // function, also complex conjugation should be changed (or rather, deleted).
957 static ex binomial_conjugate(const ex & x, const ex & y)
958 {
959         return binomial(x,y).hold();
960 }
961
962 static ex binomial_real_part(const ex & x, const ex & y)
963 {
964         return binomial(x,y).hold();
965 }
966
967 static ex binomial_imag_part(const ex & x, const ex & y)
968 {
969         return 0;
970 }
971
972 REGISTER_FUNCTION(binomial, eval_func(binomial_eval).
973                             evalf_func(binomial_evalf).
974                             conjugate_func(binomial_conjugate).
975                             real_part_func(binomial_real_part).
976                             imag_part_func(binomial_imag_part));
977
978 //////////
979 // Order term function (for truncated power series)
980 //////////
981
982 static ex Order_eval(const ex & x)
983 {
984         if (is_exactly_a<numeric>(x)) {
985                 // O(c) -> O(1) or 0
986                 if (!x.is_zero())
987                         return Order(_ex1).hold();
988                 else
989                         return _ex0;
990         } else if (is_exactly_a<mul>(x)) {
991                 const mul &m = ex_to<mul>(x);
992                 // O(c*expr) -> O(expr)
993                 if (is_exactly_a<numeric>(m.op(m.nops() - 1)))
994                         return Order(x / m.op(m.nops() - 1)).hold();
995         }
996         return Order(x).hold();
997 }
998
999 static ex Order_series(const ex & x, const relational & r, int order, unsigned options)
1000 {
1001         // Just wrap the function into a pseries object
1002         GINAC_ASSERT(is_a<symbol>(r.lhs()));
1003         const symbol &s = ex_to<symbol>(r.lhs());
1004         epvector new_seq { expair(Order(_ex1), numeric(std::min(x.ldegree(s), order))) };
1005         return pseries(r, std::move(new_seq));
1006 }
1007
1008 static ex Order_conjugate(const ex & x)
1009 {
1010         return Order(x).hold();
1011 }
1012
1013 static ex Order_real_part(const ex & x)
1014 {
1015         return Order(x).hold();
1016 }
1017
1018 static ex Order_imag_part(const ex & x)
1019 {
1020         if(x.info(info_flags::real))
1021                 return 0;
1022         return Order(x).hold();
1023 }
1024
1025 static ex Order_expl_derivative(const ex & arg, const symbol & s)
1026 {
1027         return Order(arg.diff(s));
1028 }
1029
1030 REGISTER_FUNCTION(Order, eval_func(Order_eval).
1031                          series_func(Order_series).
1032                          latex_name("\\mathcal{O}").
1033                          expl_derivative_func(Order_expl_derivative).
1034                          conjugate_func(Order_conjugate).
1035                          real_part_func(Order_real_part).
1036                          imag_part_func(Order_imag_part));
1037
1038 //////////
1039 // Solve linear system
1040 //////////
1041
1042 ex lsolve(const ex &eqns, const ex &symbols, unsigned options)
1043 {
1044         // solve a system of linear equations
1045         if (eqns.info(info_flags::relation_equal)) {
1046                 if (!symbols.info(info_flags::symbol))
1047                         throw(std::invalid_argument("lsolve(): 2nd argument must be a symbol"));
1048                 const ex sol = lsolve(lst{eqns}, lst{symbols});
1049                 
1050                 GINAC_ASSERT(sol.nops()==1);
1051                 GINAC_ASSERT(is_exactly_a<relational>(sol.op(0)));
1052                 
1053                 return sol.op(0).op(1); // return rhs of first solution
1054         }
1055         
1056         // syntax checks
1057         if (!eqns.info(info_flags::list)) {
1058                 throw(std::invalid_argument("lsolve(): 1st argument must be a list or an equation"));
1059         }
1060         for (size_t i=0; i<eqns.nops(); i++) {
1061                 if (!eqns.op(i).info(info_flags::relation_equal)) {
1062                         throw(std::invalid_argument("lsolve(): 1st argument must be a list of equations"));
1063                 }
1064         }
1065         if (!symbols.info(info_flags::list)) {
1066                 throw(std::invalid_argument("lsolve(): 2nd argument must be a list or a symbol"));
1067         }
1068         for (size_t i=0; i<symbols.nops(); i++) {
1069                 if (!symbols.op(i).info(info_flags::symbol)) {
1070                         throw(std::invalid_argument("lsolve(): 2nd argument must be a list of symbols"));
1071                 }
1072         }
1073         
1074         // build matrix from equation system
1075         matrix sys(eqns.nops(),symbols.nops());
1076         matrix rhs(eqns.nops(),1);
1077         matrix vars(symbols.nops(),1);
1078         
1079         for (size_t r=0; r<eqns.nops(); r++) {
1080                 const ex eq = eqns.op(r).op(0)-eqns.op(r).op(1); // lhs-rhs==0
1081                 ex linpart = eq;
1082                 for (size_t c=0; c<symbols.nops(); c++) {
1083                         const ex co = eq.coeff(ex_to<symbol>(symbols.op(c)),1);
1084                         linpart -= co*symbols.op(c);
1085                         sys(r,c) = co;
1086                 }
1087                 linpart = linpart.expand();
1088                 rhs(r,0) = -linpart;
1089         }
1090         
1091         // test if system is linear and fill vars matrix
1092         for (size_t i=0; i<symbols.nops(); i++) {
1093                 vars(i,0) = symbols.op(i);
1094                 if (sys.has(symbols.op(i)))
1095                         throw(std::logic_error("lsolve: system is not linear"));
1096                 if (rhs.has(symbols.op(i)))
1097                         throw(std::logic_error("lsolve: system is not linear"));
1098         }
1099         
1100         matrix solution;
1101         try {
1102                 solution = sys.solve(vars,rhs,options);
1103         } catch (const std::runtime_error & e) {
1104                 // Probably singular matrix or otherwise overdetermined system:
1105                 // It is consistent to return an empty list
1106                 return lst{};
1107         }
1108         GINAC_ASSERT(solution.cols()==1);
1109         GINAC_ASSERT(solution.rows()==symbols.nops());
1110         
1111         // return list of equations of the form lst{var1==sol1,var2==sol2,...}
1112         lst sollist;
1113         for (size_t i=0; i<symbols.nops(); i++)
1114                 sollist.append(symbols.op(i)==solution(i,0));
1115         
1116         return sollist;
1117 }
1118
1119 //////////
1120 // Find real root of f(x) numerically
1121 //////////
1122
1123 const numeric
1124 fsolve(const ex& f_in, const symbol& x, const numeric& x1, const numeric& x2)
1125 {
1126         if (!x1.is_real() || !x2.is_real()) {
1127                 throw std::runtime_error("fsolve(): interval not bounded by real numbers");
1128         }
1129         if (x1==x2) {
1130                 throw std::runtime_error("fsolve(): vanishing interval");
1131         }
1132         // xx[0] == left interval limit, xx[1] == right interval limit.
1133         // fx[0] == f(xx[0]), fx[1] == f(xx[1]).
1134         // We keep the root bracketed: xx[0]<xx[1] and fx[0]*fx[1]<0.
1135         numeric xx[2] = { x1<x2 ? x1 : x2,
1136                           x1<x2 ? x2 : x1 };
1137         ex f;
1138         if (is_a<relational>(f_in)) {
1139                 f = f_in.lhs()-f_in.rhs();
1140         } else {
1141                 f = f_in;
1142         }
1143         const ex fx_[2] = { f.subs(x==xx[0]).evalf(),
1144                             f.subs(x==xx[1]).evalf() };
1145         if (!is_a<numeric>(fx_[0]) || !is_a<numeric>(fx_[1])) {
1146                 throw std::runtime_error("fsolve(): function does not evaluate numerically");
1147         }
1148         numeric fx[2] = { ex_to<numeric>(fx_[0]),
1149                           ex_to<numeric>(fx_[1]) };
1150         if (!fx[0].is_real() || !fx[1].is_real()) {
1151                 throw std::runtime_error("fsolve(): function evaluates to complex values at interval boundaries");
1152         }
1153         if (fx[0]*fx[1]>=0) {
1154                 throw std::runtime_error("fsolve(): function does not change sign at interval boundaries");
1155         }
1156
1157         // The Newton-Raphson method has quadratic convergence!  Simply put, it
1158         // replaces x with x-f(x)/f'(x) at each step.  -f/f' is the delta:
1159         const ex ff = normal(-f/f.diff(x));
1160         int side = 0;  // Start at left interval limit.
1161         numeric xxprev;
1162         numeric fxprev;
1163         do {
1164                 xxprev = xx[side];
1165                 fxprev = fx[side];
1166                 ex dx_ = ff.subs(x == xx[side]).evalf();
1167                 if (!is_a<numeric>(dx_))
1168                         throw std::runtime_error("fsolve(): function derivative does not evaluate numerically");
1169                 xx[side] += ex_to<numeric>(dx_);
1170                 // Now check if Newton-Raphson method shot out of the interval 
1171                 bool bad_shot = (side == 0 && xx[0] < xxprev) || 
1172                                 (side == 1 && xx[1] > xxprev) || xx[0] > xx[1];
1173                 if (!bad_shot) {
1174                         // Compute f(x) only if new x is inside the interval.
1175                         // The function might be difficult to compute numerically
1176                         // or even ill defined outside the interval. Also it's
1177                         // a small optimization. 
1178                         ex f_x = f.subs(x == xx[side]).evalf();
1179                         if (!is_a<numeric>(f_x))
1180                                 throw std::runtime_error("fsolve(): function does not evaluate numerically");
1181                         fx[side] = ex_to<numeric>(f_x);
1182                 }
1183                 if (bad_shot) {
1184                         // Oops, Newton-Raphson method shot out of the interval.
1185                         // Restore, and try again with the other side instead!
1186                         xx[side] = xxprev;
1187                         fx[side] = fxprev;
1188                         side = !side;
1189                         xxprev = xx[side];
1190                         fxprev = fx[side];
1191
1192                         ex dx_ = ff.subs(x == xx[side]).evalf();
1193                         if (!is_a<numeric>(dx_))
1194                                 throw std::runtime_error("fsolve(): function derivative does not evaluate numerically [2]");
1195                         xx[side] += ex_to<numeric>(dx_);
1196
1197                         ex f_x = f.subs(x==xx[side]).evalf();
1198                         if (!is_a<numeric>(f_x))
1199                                 throw std::runtime_error("fsolve(): function does not evaluate numerically [2]");
1200                         fx[side] = ex_to<numeric>(f_x);
1201                 }
1202                 if ((fx[side]<0 && fx[!side]<0) || (fx[side]>0 && fx[!side]>0)) {
1203                         // Oops, the root isn't bracketed any more.
1204                         // Restore, and perform a bisection!
1205                         xx[side] = xxprev;
1206                         fx[side] = fxprev;
1207
1208                         // Ah, the bisection! Bisections converge linearly. Unfortunately,
1209                         // they occur pretty often when Newton-Raphson arrives at an x too
1210                         // close to the result on one side of the interval and
1211                         // f(x-f(x)/f'(x)) turns out to have the same sign as f(x) due to
1212                         // precision errors! Recall that this function does not have a
1213                         // precision goal as one of its arguments but instead relies on
1214                         // x converging to a fixed point. We speed up the (safe but slow)
1215                         // bisection method by mixing in a dash of the (unsafer but faster)
1216                         // secant method: Instead of splitting the interval at the
1217                         // arithmetic mean (bisection), we split it nearer to the root as
1218                         // determined by the secant between the values xx[0] and xx[1].
1219                         // Don't set the secant_weight to one because that could disturb
1220                         // the convergence in some corner cases!
1221                         constexpr double secant_weight = 0.984375;  // == 63/64 < 1
1222                         numeric xxmid = (1-secant_weight)*0.5*(xx[0]+xx[1])
1223                             + secant_weight*(xx[0]+fx[0]*(xx[0]-xx[1])/(fx[1]-fx[0]));
1224                         ex fxmid_ = f.subs(x == xxmid).evalf();
1225                         if (!is_a<numeric>(fxmid_))
1226                                 throw std::runtime_error("fsolve(): function does not evaluate numerically [3]");
1227                         numeric fxmid = ex_to<numeric>(fxmid_);
1228                         if (fxmid.is_zero()) {
1229                                 // Luck strikes...
1230                                 return xxmid;
1231                         }
1232                         if ((fxmid<0 && fx[side]>0) || (fxmid>0 && fx[side]<0)) {
1233                                 side = !side;
1234                         }
1235                         xxprev = xx[side];
1236                         fxprev = fx[side];
1237                         xx[side] = xxmid;
1238                         fx[side] = fxmid;
1239                 }
1240         } while (xxprev!=xx[side]);
1241         return xxprev;
1242 }
1243
1244
1245 /* Force inclusion of functions from inifcns_gamma and inifcns_zeta
1246  * for static lib (so ginsh will see them). */
1247 unsigned force_include_tgamma = tgamma_SERIAL::serial;
1248 unsigned force_include_zeta1 = zeta1_SERIAL::serial;
1249
1250 } // namespace GiNaC