Added complex conjugation methods and GiNaC function "conjugate".
[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-2004 Johannes Gutenberg University Mainz, Germany
7  *
8  *  This program is free software; you can redistribute it and/or modify
9  *  it under the terms of the GNU General Public License as published by
10  *  the Free Software Foundation; either version 2 of the License, or
11  *  (at your option) any later version.
12  *
13  *  This program is distributed in the hope that it will be useful,
14  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
16  *  GNU General Public License for more details.
17  *
18  *  You should have received a copy of the GNU General Public License
19  *  along with this program; if not, write to the Free Software
20  *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
21  */
22
23 #include <vector>
24 #include <stdexcept>
25
26 #include "inifcns.h"
27 #include "ex.h"
28 #include "constant.h"
29 #include "lst.h"
30 #include "matrix.h"
31 #include "mul.h"
32 #include "power.h"
33 #include "operators.h"
34 #include "relational.h"
35 #include "pseries.h"
36 #include "symbol.h"
37 #include "symmetry.h"
38 #include "utils.h"
39
40 namespace GiNaC {
41
42 //////////
43 // complex conjugate
44 //////////
45
46 static ex conjugate_evalf(const ex & arg)
47 {
48         if (is_exactly_a<numeric>(arg)) {
49                 return ex_to<numeric>(arg).conjugate();
50         }
51         return conjugate(arg).hold();
52 }
53
54 static ex conjugate_eval(const ex & arg)
55 {
56         return arg.conjugate();
57 }
58
59 static void conjugate_print_latex(const ex & arg, const print_context & c)
60 {
61         c.s << "\bar{"; arg.print(c); c.s << "}";
62 }
63
64 static ex conjugate_conjugate(const ex & arg)
65 {
66         return arg;
67 }
68
69 REGISTER_FUNCTION(conjugate, eval_func(conjugate_eval).
70                        evalf_func(conjugate_evalf).
71                        print_func<print_latex>(conjugate_print_latex).
72                        conjugate_func(conjugate_conjugate));
73
74 //////////
75 // absolute value
76 //////////
77
78 static ex abs_evalf(const ex & arg)
79 {
80         if (is_exactly_a<numeric>(arg))
81                 return abs(ex_to<numeric>(arg));
82         
83         return abs(arg).hold();
84 }
85
86 static ex abs_eval(const ex & arg)
87 {
88         if (is_exactly_a<numeric>(arg))
89                 return abs(ex_to<numeric>(arg));
90         else
91                 return abs(arg).hold();
92 }
93
94 static void abs_print_latex(const ex & arg, const print_context & c)
95 {
96         c.s << "{|"; arg.print(c); c.s << "|}";
97 }
98
99 static void abs_print_csrc_float(const ex & arg, const print_context & c)
100 {
101         c.s << "fabs("; arg.print(c); c.s << ")";
102 }
103
104 static ex abs_conjugate(const ex & arg)
105 {
106         return abs(arg);
107 }
108
109 REGISTER_FUNCTION(abs, eval_func(abs_eval).
110                        evalf_func(abs_evalf).
111                        print_func<print_latex>(abs_print_latex).
112                        print_func<print_csrc_float>(abs_print_csrc_float).
113                        print_func<print_csrc_double>(abs_print_csrc_float).
114                        conjugate_func(abs_conjugate));
115
116
117 //////////
118 // Complex sign
119 //////////
120
121 static ex csgn_evalf(const ex & arg)
122 {
123         if (is_exactly_a<numeric>(arg))
124                 return csgn(ex_to<numeric>(arg));
125         
126         return csgn(arg).hold();
127 }
128
129 static ex csgn_eval(const ex & arg)
130 {
131         if (is_exactly_a<numeric>(arg))
132                 return csgn(ex_to<numeric>(arg));
133         
134         else if (is_exactly_a<mul>(arg) &&
135                  is_exactly_a<numeric>(arg.op(arg.nops()-1))) {
136                 numeric oc = ex_to<numeric>(arg.op(arg.nops()-1));
137                 if (oc.is_real()) {
138                         if (oc > 0)
139                                 // csgn(42*x) -> csgn(x)
140                                 return csgn(arg/oc).hold();
141                         else
142                                 // csgn(-42*x) -> -csgn(x)
143                                 return -csgn(arg/oc).hold();
144                 }
145                 if (oc.real().is_zero()) {
146                         if (oc.imag() > 0)
147                                 // csgn(42*I*x) -> csgn(I*x)
148                                 return csgn(I*arg/oc).hold();
149                         else
150                                 // csgn(-42*I*x) -> -csgn(I*x)
151                                 return -csgn(I*arg/oc).hold();
152                 }
153         }
154         
155         return csgn(arg).hold();
156 }
157
158 static ex csgn_series(const ex & arg,
159                       const relational & rel,
160                       int order,
161                       unsigned options)
162 {
163         const ex arg_pt = arg.subs(rel, subs_options::no_pattern);
164         if (arg_pt.info(info_flags::numeric)
165             && ex_to<numeric>(arg_pt).real().is_zero()
166             && !(options & series_options::suppress_branchcut))
167                 throw (std::domain_error("csgn_series(): on imaginary axis"));
168         
169         epvector seq;
170         seq.push_back(expair(csgn(arg_pt), _ex0));
171         return pseries(rel,seq);
172 }
173
174 static ex csgn_conjugate(const ex& arg)
175 {
176         return csgn(arg);
177 }
178
179 REGISTER_FUNCTION(csgn, eval_func(csgn_eval).
180                         evalf_func(csgn_evalf).
181                         series_func(csgn_series).
182                         conjugate_func(csgn_conjugate));
183
184
185 //////////
186 // Eta function: eta(x,y) == log(x*y) - log(x) - log(y).
187 // This function is closely related to the unwinding number K, sometimes found
188 // in modern literature: K(z) == (z-log(exp(z)))/(2*Pi*I).
189 //////////
190
191 static ex eta_evalf(const ex &x, const ex &y)
192 {
193         // It seems like we basically have to replicate the eval function here,
194         // since the expression might not be fully evaluated yet.
195         if (x.info(info_flags::positive) || y.info(info_flags::positive))
196                 return _ex0;
197
198         if (x.info(info_flags::numeric) &&      y.info(info_flags::numeric)) {
199                 const numeric nx = ex_to<numeric>(x);
200                 const numeric ny = ex_to<numeric>(y);
201                 const numeric nxy = ex_to<numeric>(x*y);
202                 int cut = 0;
203                 if (nx.is_real() && nx.is_negative())
204                         cut -= 4;
205                 if (ny.is_real() && ny.is_negative())
206                         cut -= 4;
207                 if (nxy.is_real() && nxy.is_negative())
208                         cut += 4;
209                 return evalf(I/4*Pi)*((csgn(-imag(nx))+1)*(csgn(-imag(ny))+1)*(csgn(imag(nxy))+1)-
210                                       (csgn(imag(nx))+1)*(csgn(imag(ny))+1)*(csgn(-imag(nxy))+1)+cut);
211         }
212
213         return eta(x,y).hold();
214 }
215
216 static ex eta_eval(const ex &x, const ex &y)
217 {
218         // trivial:  eta(x,c) -> 0  if c is real and positive
219         if (x.info(info_flags::positive) || y.info(info_flags::positive))
220                 return _ex0;
221
222         if (x.info(info_flags::numeric) &&      y.info(info_flags::numeric)) {
223                 // don't call eta_evalf here because it would call Pi.evalf()!
224                 const numeric nx = ex_to<numeric>(x);
225                 const numeric ny = ex_to<numeric>(y);
226                 const numeric nxy = ex_to<numeric>(x*y);
227                 int cut = 0;
228                 if (nx.is_real() && nx.is_negative())
229                         cut -= 4;
230                 if (ny.is_real() && ny.is_negative())
231                         cut -= 4;
232                 if (nxy.is_real() && nxy.is_negative())
233                         cut += 4;
234                 return (I/4)*Pi*((csgn(-imag(nx))+1)*(csgn(-imag(ny))+1)*(csgn(imag(nxy))+1)-
235                                  (csgn(imag(nx))+1)*(csgn(imag(ny))+1)*(csgn(-imag(nxy))+1)+cut);
236         }
237         
238         return eta(x,y).hold();
239 }
240
241 static ex eta_series(const ex & x, const ex & y,
242                      const relational & rel,
243                      int order,
244                      unsigned options)
245 {
246         const ex x_pt = x.subs(rel, subs_options::no_pattern);
247         const ex y_pt = y.subs(rel, subs_options::no_pattern);
248         if ((x_pt.info(info_flags::numeric) && x_pt.info(info_flags::negative)) ||
249             (y_pt.info(info_flags::numeric) && y_pt.info(info_flags::negative)) ||
250             ((x_pt*y_pt).info(info_flags::numeric) && (x_pt*y_pt).info(info_flags::negative)))
251                         throw (std::domain_error("eta_series(): on discontinuity"));
252         epvector seq;
253         seq.push_back(expair(eta(x_pt,y_pt), _ex0));
254         return pseries(rel,seq);
255 }
256
257 static ex eta_conjugate(const ex & x, const ex & y)
258 {
259         return -eta(x,y);
260 }
261
262 REGISTER_FUNCTION(eta, eval_func(eta_eval).
263                        evalf_func(eta_evalf).
264                        series_func(eta_series).
265                        latex_name("\\eta").
266                        set_symmetry(sy_symm(0, 1)).
267                        conjugate_func(eta_conjugate));
268
269
270 //////////
271 // dilogarithm
272 //////////
273
274 static ex Li2_evalf(const ex & x)
275 {
276         if (is_exactly_a<numeric>(x))
277                 return Li2(ex_to<numeric>(x));
278         
279         return Li2(x).hold();
280 }
281
282 static ex Li2_eval(const ex & x)
283 {
284         if (x.info(info_flags::numeric)) {
285                 // Li2(0) -> 0
286                 if (x.is_zero())
287                         return _ex0;
288                 // Li2(1) -> Pi^2/6
289                 if (x.is_equal(_ex1))
290                         return power(Pi,_ex2)/_ex6;
291                 // Li2(1/2) -> Pi^2/12 - log(2)^2/2
292                 if (x.is_equal(_ex1_2))
293                         return power(Pi,_ex2)/_ex12 + power(log(_ex2),_ex2)*_ex_1_2;
294                 // Li2(-1) -> -Pi^2/12
295                 if (x.is_equal(_ex_1))
296                         return -power(Pi,_ex2)/_ex12;
297                 // Li2(I) -> -Pi^2/48+Catalan*I
298                 if (x.is_equal(I))
299                         return power(Pi,_ex2)/_ex_48 + Catalan*I;
300                 // Li2(-I) -> -Pi^2/48-Catalan*I
301                 if (x.is_equal(-I))
302                         return power(Pi,_ex2)/_ex_48 - Catalan*I;
303                 // Li2(float)
304                 if (!x.info(info_flags::crational))
305                         return Li2(ex_to<numeric>(x));
306         }
307         
308         return Li2(x).hold();
309 }
310
311 static ex Li2_deriv(const ex & x, unsigned deriv_param)
312 {
313         GINAC_ASSERT(deriv_param==0);
314         
315         // d/dx Li2(x) -> -log(1-x)/x
316         return -log(_ex1-x)/x;
317 }
318
319 static ex Li2_series(const ex &x, const relational &rel, int order, unsigned options)
320 {
321         const ex x_pt = x.subs(rel, subs_options::no_pattern);
322         if (x_pt.info(info_flags::numeric)) {
323                 // First special case: x==0 (derivatives have poles)
324                 if (x_pt.is_zero()) {
325                         // method:
326                         // The problem is that in d/dx Li2(x==0) == -log(1-x)/x we cannot 
327                         // simply substitute x==0.  The limit, however, exists: it is 1.
328                         // We also know all higher derivatives' limits:
329                         // (d/dx)^n Li2(x) == n!/n^2.
330                         // So the primitive series expansion is
331                         // Li2(x==0) == x + x^2/4 + x^3/9 + ...
332                         // and so on.
333                         // We first construct such a primitive series expansion manually in
334                         // a dummy symbol s and then insert the argument's series expansion
335                         // for s.  Reexpanding the resulting series returns the desired
336                         // result.
337                         const symbol s;
338                         ex ser;
339                         // manually construct the primitive expansion
340                         for (int i=1; i<order; ++i)
341                                 ser += pow(s,i) / pow(numeric(i), _num2);
342                         // substitute the argument's series expansion
343                         ser = ser.subs(s==x.series(rel, order), subs_options::no_pattern);
344                         // maybe that was terminating, so add a proper order term
345                         epvector nseq;
346                         nseq.push_back(expair(Order(_ex1), order));
347                         ser += pseries(rel, nseq);
348                         // reexpanding it will collapse the series again
349                         return ser.series(rel, order);
350                         // NB: Of course, this still does not allow us to compute anything
351                         // like sin(Li2(x)).series(x==0,2), since then this code here is
352                         // not reached and the derivative of sin(Li2(x)) doesn't allow the
353                         // substitution x==0.  Probably limits *are* needed for the general
354                         // cases.  In case L'Hospital's rule is implemented for limits and
355                         // basic::series() takes care of this, this whole block is probably
356                         // obsolete!
357                 }
358                 // second special case: x==1 (branch point)
359                 if (x_pt.is_equal(_ex1)) {
360                         // method:
361                         // construct series manually in a dummy symbol s
362                         const symbol s;
363                         ex ser = zeta(_ex2);
364                         // manually construct the primitive expansion
365                         for (int i=1; i<order; ++i)
366                                 ser += pow(1-s,i) * (numeric(1,i)*(I*Pi+log(s-1)) - numeric(1,i*i));
367                         // substitute the argument's series expansion
368                         ser = ser.subs(s==x.series(rel, order), subs_options::no_pattern);
369                         // maybe that was terminating, so add a proper order term
370                         epvector nseq;
371                         nseq.push_back(expair(Order(_ex1), order));
372                         ser += pseries(rel, nseq);
373                         // reexpanding it will collapse the series again
374                         return ser.series(rel, order);
375                 }
376                 // third special case: x real, >=1 (branch cut)
377                 if (!(options & series_options::suppress_branchcut) &&
378                         ex_to<numeric>(x_pt).is_real() && ex_to<numeric>(x_pt)>1) {
379                         // method:
380                         // This is the branch cut: assemble the primitive series manually
381                         // and then add the corresponding complex step function.
382                         const symbol &s = ex_to<symbol>(rel.lhs());
383                         const ex point = rel.rhs();
384                         const symbol foo;
385                         epvector seq;
386                         // zeroth order term:
387                         seq.push_back(expair(Li2(x_pt), _ex0));
388                         // compute the intermediate terms:
389                         ex replarg = series(Li2(x), s==foo, order);
390                         for (size_t i=1; i<replarg.nops()-1; ++i)
391                                 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));
392                         // append an order term:
393                         seq.push_back(expair(Order(_ex1), replarg.nops()-1));
394                         return pseries(rel, seq);
395                 }
396         }
397         // all other cases should be safe, by now:
398         throw do_taylor();  // caught by function::series()
399 }
400
401 REGISTER_FUNCTION(Li2, eval_func(Li2_eval).
402                        evalf_func(Li2_evalf).
403                        derivative_func(Li2_deriv).
404                        series_func(Li2_series).
405                        latex_name("\\mbox{Li}_2"));
406
407 //////////
408 // trilogarithm
409 //////////
410
411 static ex Li3_eval(const ex & x)
412 {
413         if (x.is_zero())
414                 return x;
415         return Li3(x).hold();
416 }
417
418 REGISTER_FUNCTION(Li3, eval_func(Li3_eval).
419                        latex_name("\\mbox{Li}_3"));
420
421 //////////
422 // Derivatives of Riemann's Zeta-function  zetaderiv(0,x)==zeta(x)
423 //////////
424
425 static ex zetaderiv_eval(const ex & n, const ex & x)
426 {
427         if (n.info(info_flags::numeric)) {
428                 // zetaderiv(0,x) -> zeta(x)
429                 if (n.is_zero())
430                         return zeta(x);
431         }
432         
433         return zetaderiv(n, x).hold();
434 }
435
436 static ex zetaderiv_deriv(const ex & n, const ex & x, unsigned deriv_param)
437 {
438         GINAC_ASSERT(deriv_param<2);
439         
440         if (deriv_param==0) {
441                 // d/dn zeta(n,x)
442                 throw(std::logic_error("cannot diff zetaderiv(n,x) with respect to n"));
443         }
444         // d/dx psi(n,x)
445         return zetaderiv(n+1,x);
446 }
447
448 REGISTER_FUNCTION(zetaderiv, eval_func(zetaderiv_eval).
449                                  derivative_func(zetaderiv_deriv).
450                                  latex_name("\\zeta^\\prime"));
451
452 //////////
453 // factorial
454 //////////
455
456 static ex factorial_evalf(const ex & x)
457 {
458         return factorial(x).hold();
459 }
460
461 static ex factorial_eval(const ex & x)
462 {
463         if (is_exactly_a<numeric>(x))
464                 return factorial(ex_to<numeric>(x));
465         else
466                 return factorial(x).hold();
467 }
468
469 static ex factorial_conjugate(const ex & x)
470 {
471         return factorial(x);
472 }
473
474 REGISTER_FUNCTION(factorial, eval_func(factorial_eval).
475                              evalf_func(factorial_evalf).
476                              conjugate_func(factorial_conjugate));
477
478 //////////
479 // binomial
480 //////////
481
482 static ex binomial_evalf(const ex & x, const ex & y)
483 {
484         return binomial(x, y).hold();
485 }
486
487 static ex binomial_eval(const ex & x, const ex &y)
488 {
489         if (is_exactly_a<numeric>(x) && is_exactly_a<numeric>(y))
490                 return binomial(ex_to<numeric>(x), ex_to<numeric>(y));
491         else
492                 return binomial(x, y).hold();
493 }
494
495 // At the moment the numeric evaluation of a binomail function always
496 // gives a real number, but if this would be implemented using the gamma
497 // function, also complex conjugation should be changed (or rather, deleted).
498 static ex binomial_conjugate(const ex & x, const ex & y)
499 {
500         return binomial(x,y);
501 }
502
503 REGISTER_FUNCTION(binomial, eval_func(binomial_eval).
504                             evalf_func(binomial_evalf).
505                             conjugate_func(binomial_conjugate));
506
507 //////////
508 // Order term function (for truncated power series)
509 //////////
510
511 static ex Order_eval(const ex & x)
512 {
513         if (is_exactly_a<numeric>(x)) {
514                 // O(c) -> O(1) or 0
515                 if (!x.is_zero())
516                         return Order(_ex1).hold();
517                 else
518                         return _ex0;
519         } else if (is_exactly_a<mul>(x)) {
520                 const mul &m = ex_to<mul>(x);
521                 // O(c*expr) -> O(expr)
522                 if (is_exactly_a<numeric>(m.op(m.nops() - 1)))
523                         return Order(x / m.op(m.nops() - 1)).hold();
524         }
525         return Order(x).hold();
526 }
527
528 static ex Order_series(const ex & x, const relational & r, int order, unsigned options)
529 {
530         // Just wrap the function into a pseries object
531         epvector new_seq;
532         GINAC_ASSERT(is_a<symbol>(r.lhs()));
533         const symbol &s = ex_to<symbol>(r.lhs());
534         new_seq.push_back(expair(Order(_ex1), numeric(std::min(x.ldegree(s), order))));
535         return pseries(r, new_seq);
536 }
537
538 static ex Order_conjugate(const ex & x)
539 {
540         return Order(x);
541 }
542
543 // Differentiation is handled in function::derivative because of its special requirements
544
545 REGISTER_FUNCTION(Order, eval_func(Order_eval).
546                          series_func(Order_series).
547                          latex_name("\\mathcal{O}").
548                          conjugate_func(Order_conjugate));
549
550 //////////
551 // Solve linear system
552 //////////
553
554 ex lsolve(const ex &eqns, const ex &symbols, unsigned options)
555 {
556         // solve a system of linear equations
557         if (eqns.info(info_flags::relation_equal)) {
558                 if (!symbols.info(info_flags::symbol))
559                         throw(std::invalid_argument("lsolve(): 2nd argument must be a symbol"));
560                 const ex sol = lsolve(lst(eqns),lst(symbols));
561                 
562                 GINAC_ASSERT(sol.nops()==1);
563                 GINAC_ASSERT(is_exactly_a<relational>(sol.op(0)));
564                 
565                 return sol.op(0).op(1); // return rhs of first solution
566         }
567         
568         // syntax checks
569         if (!eqns.info(info_flags::list)) {
570                 throw(std::invalid_argument("lsolve(): 1st argument must be a list"));
571         }
572         for (size_t i=0; i<eqns.nops(); i++) {
573                 if (!eqns.op(i).info(info_flags::relation_equal)) {
574                         throw(std::invalid_argument("lsolve(): 1st argument must be a list of equations"));
575                 }
576         }
577         if (!symbols.info(info_flags::list)) {
578                 throw(std::invalid_argument("lsolve(): 2nd argument must be a list"));
579         }
580         for (size_t i=0; i<symbols.nops(); i++) {
581                 if (!symbols.op(i).info(info_flags::symbol)) {
582                         throw(std::invalid_argument("lsolve(): 2nd argument must be a list of symbols"));
583                 }
584         }
585         
586         // build matrix from equation system
587         matrix sys(eqns.nops(),symbols.nops());
588         matrix rhs(eqns.nops(),1);
589         matrix vars(symbols.nops(),1);
590         
591         for (size_t r=0; r<eqns.nops(); r++) {
592                 const ex eq = eqns.op(r).op(0)-eqns.op(r).op(1); // lhs-rhs==0
593                 ex linpart = eq;
594                 for (size_t c=0; c<symbols.nops(); c++) {
595                         const ex co = eq.coeff(ex_to<symbol>(symbols.op(c)),1);
596                         linpart -= co*symbols.op(c);
597                         sys(r,c) = co;
598                 }
599                 linpart = linpart.expand();
600                 rhs(r,0) = -linpart;
601         }
602         
603         // test if system is linear and fill vars matrix
604         for (size_t i=0; i<symbols.nops(); i++) {
605                 vars(i,0) = symbols.op(i);
606                 if (sys.has(symbols.op(i)))
607                         throw(std::logic_error("lsolve: system is not linear"));
608                 if (rhs.has(symbols.op(i)))
609                         throw(std::logic_error("lsolve: system is not linear"));
610         }
611         
612         matrix solution;
613         try {
614                 solution = sys.solve(vars,rhs,options);
615         } catch (const std::runtime_error & e) {
616                 // Probably singular matrix or otherwise overdetermined system:
617                 // It is consistent to return an empty list
618                 return lst();
619         }
620         GINAC_ASSERT(solution.cols()==1);
621         GINAC_ASSERT(solution.rows()==symbols.nops());
622         
623         // return list of equations of the form lst(var1==sol1,var2==sol2,...)
624         lst sollist;
625         for (size_t i=0; i<symbols.nops(); i++)
626                 sollist.append(symbols.op(i)==solution(i,0));
627         
628         return sollist;
629 }
630
631 /* Force inclusion of functions from inifcns_gamma and inifcns_zeta
632  * for static lib (so ginsh will see them). */
633 unsigned force_include_tgamma = tgamma_SERIAL::serial;
634 unsigned force_include_zeta1 = zeta1_SERIAL::serial;
635
636 } // namespace GiNaC