#ifndef around namespace GiNaC { }
[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 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 "ncmul.h"
33 #include "numeric.h"
34 #include "power.h"
35 #include "relational.h"
36 #include "series.h"
37 #include "symbol.h"
38
39 #ifndef NO_GINAC_NAMESPACE
40 namespace GiNaC {
41 #endif // ndef NO_GINAC_NAMESPACE
42
43 //////////
44 // dilogarithm
45 //////////
46
47 static ex Li2_eval(ex const & x)
48 {
49     if (x.is_zero())
50         return x;
51     if (x.is_equal(exONE()))
52         return power(Pi, 2) / 6;
53     if (x.is_equal(exMINUSONE()))
54         return -power(Pi, 2) / 12;
55     return Li2(x).hold();
56 }
57
58 REGISTER_FUNCTION(Li2, Li2_eval, NULL, NULL, NULL);
59
60 //////////
61 // trilogarithm
62 //////////
63
64 static ex Li3_eval(ex const & x)
65 {
66     if (x.is_zero())
67         return x;
68     return Li3(x).hold();
69 }
70
71 REGISTER_FUNCTION(Li3, Li3_eval, NULL, NULL, NULL);
72
73 //////////
74 // factorial
75 //////////
76
77 static ex factorial_evalf(ex const & x)
78 {
79     return factorial(x).hold();
80 }
81
82 static ex factorial_eval(ex const & x)
83 {
84     if (is_ex_exactly_of_type(x, numeric))
85         return factorial(ex_to_numeric(x));
86     else
87         return factorial(x).hold();
88 }
89
90 REGISTER_FUNCTION(factorial, factorial_eval, factorial_evalf, NULL, NULL);
91
92 //////////
93 // binomial
94 //////////
95
96 static ex binomial_evalf(ex const & x, ex const & y)
97 {
98     return binomial(x, y).hold();
99 }
100
101 static ex binomial_eval(ex const & x, ex const &y)
102 {
103     if (is_ex_exactly_of_type(x, numeric) && is_ex_exactly_of_type(y, numeric))
104         return binomial(ex_to_numeric(x), ex_to_numeric(y));
105     else
106         return binomial(x, y).hold();
107 }
108
109 REGISTER_FUNCTION(binomial, binomial_eval, binomial_evalf, NULL, NULL);
110
111 //////////
112 // Order term function (for truncated power series)
113 //////////
114
115 static ex Order_eval(ex const & x)
116 {
117         if (is_ex_exactly_of_type(x, numeric)) {
118
119                 // O(c)=O(1)
120                 return Order(exONE()).hold();
121
122         } else if (is_ex_exactly_of_type(x, mul)) {
123
124                 mul *m = static_cast<mul *>(x.bp);
125                 if (is_ex_exactly_of_type(m->op(m->nops() - 1), numeric)) {
126
127                         // O(c*expr)=O(expr)
128                         return Order(x / m->op(m->nops() - 1)).hold();
129                 }
130         }
131         return Order(x).hold();
132 }
133
134 static ex Order_series(ex const & x, symbol const & s, ex const & point, int order)
135 {
136         // Just wrap the function into a series object
137         epvector new_seq;
138         new_seq.push_back(expair(Order(exONE()), numeric(min(x.ldegree(s), order))));
139         return series(s, point, new_seq);
140 }
141
142 REGISTER_FUNCTION(Order, Order_eval, NULL, NULL, Order_series);
143
144 /** linear solve. */
145 ex lsolve(ex const &eqns, ex const &symbols)
146 {
147     // solve a system of linear equations
148     if (eqns.info(info_flags::relation_equal)) {
149         if (!symbols.info(info_flags::symbol)) {
150             throw(std::invalid_argument("lsolve: 2nd argument must be a symbol"));
151         }
152         ex sol=lsolve(lst(eqns),lst(symbols));
153         
154         GINAC_ASSERT(sol.nops()==1);
155         GINAC_ASSERT(is_ex_exactly_of_type(sol.op(0),relational));
156         
157         return sol.op(0).op(1); // return rhs of first solution
158     }
159     
160     // syntax checks
161     if (!eqns.info(info_flags::list)) {
162         throw(std::invalid_argument("lsolve: 1st argument must be a list"));
163     }
164     for (int i=0; i<eqns.nops(); i++) {
165         if (!eqns.op(i).info(info_flags::relation_equal)) {
166             throw(std::invalid_argument("lsolve: 1st argument must be a list of equations"));
167         }
168     }
169     if (!symbols.info(info_flags::list)) {
170         throw(std::invalid_argument("lsolve: 2nd argument must be a list"));
171     }
172     for (int i=0; i<symbols.nops(); i++) {
173         if (!symbols.op(i).info(info_flags::symbol)) {
174             throw(std::invalid_argument("lsolve: 2nd argument must be a list of symbols"));
175         }
176     }
177     
178     // build matrix from equation system
179     matrix sys(eqns.nops(),symbols.nops());
180     matrix rhs(eqns.nops(),1);
181     matrix vars(symbols.nops(),1);
182
183     for (int r=0; r<eqns.nops(); r++) {
184         ex eq=eqns.op(r).op(0)-eqns.op(r).op(1); // lhs-rhs==0
185         ex linpart=eq;
186         for (int c=0; c<symbols.nops(); c++) {
187             ex co=eq.coeff(ex_to_symbol(symbols.op(c)),1);
188             linpart -= co*symbols.op(c);
189             sys.set(r,c,co);
190         }
191         linpart=linpart.expand();
192         rhs.set(r,0,-linpart);
193     }
194     
195     // test if system is linear and fill vars matrix
196     for (int i=0; i<symbols.nops(); i++) {
197         vars.set(i,0,symbols.op(i));
198         if (sys.has(symbols.op(i))) {
199             throw(std::logic_error("lsolve: system is not linear"));
200         }
201         if (rhs.has(symbols.op(i))) {
202             throw(std::logic_error("lsolve: system is not linear"));
203         }
204     }
205     
206     //matrix solution=sys.solve(rhs);
207     matrix solution;
208     try {
209         solution=sys.fraction_free_elim(vars,rhs);
210     } catch (runtime_error const & e) {
211         // probably singular matrix (or other error)
212         // return empty solution list
213         // cerr << e.what() << endl;
214         return lst();
215     }
216     
217     // return a list of equations
218     if (solution.cols()!=1) {
219         throw(std::runtime_error("lsolve: strange number of columns returned from matrix::solve"));
220     }
221     if (solution.rows()!=symbols.nops()) {
222         cout << "symbols.nops()=" << symbols.nops() << endl;
223         cout << "solution.rows()=" << solution.rows() << endl;
224         throw(std::runtime_error("lsolve: strange number of rows returned from matrix::solve"));
225     }
226     
227     // return list of the form lst(var1==sol1,var2==sol2,...)
228     lst sollist;
229     for (int i=0; i<symbols.nops(); i++) {
230         sollist.append(symbols.op(i)==solution(i,0));
231     }
232     
233     return sollist;
234 }
235
236 /** non-commutative power. */
237 ex ncpower(ex const &basis, unsigned exponent)
238 {
239     if (exponent==0) {
240         return exONE();
241     }
242
243     exvector v;
244     v.reserve(exponent);
245     for (unsigned i=0; i<exponent; ++i) {
246         v.push_back(basis);
247     }
248
249     return ncmul(v,1);
250 }
251
252 /** Force inclusion of functions from initcns_gamma and inifcns_zeta
253  *  for static lib (so ginsh will see them). */
254 unsigned force_include_gamma = function_index_gamma;
255 unsigned force_include_zeta = function_index_zeta;
256
257 #ifndef NO_GINAC_NAMESPACE
258 } // namespace GiNaC
259 #endif // ndef NO_GINAC_NAMESPACE