X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Finifcns.cpp;h=443022a556da6c9f48e062e2381b4be05b0e0ae3;hp=5802325f944345c448a1bd34c3268d81c6126703;hb=def23d34c68a383ce3d7da0227b984c8291a3bf9;hpb=82c60cc9f90e6045e75e4aa4c48b29bcd7c5331f diff --git a/ginac/inifcns.cpp b/ginac/inifcns.cpp index 5802325f..443022a5 100644 --- a/ginac/inifcns.cpp +++ b/ginac/inifcns.cpp @@ -3,7 +3,7 @@ * Implementation of GiNaC's initially known functions. */ /* - * GiNaC Copyright (C) 1999-2000 Johannes Gutenberg University Mainz, Germany + * GiNaC Copyright (C) 1999-2003 Johannes Gutenberg University Mainz, Germany * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -29,183 +29,330 @@ #include "lst.h" #include "matrix.h" #include "mul.h" -#include "ncmul.h" -#include "numeric.h" #include "power.h" +#include "operators.h" #include "relational.h" #include "pseries.h" #include "symbol.h" +#include "symmetry.h" #include "utils.h" -#ifndef NO_NAMESPACE_GINAC namespace GiNaC { -#endif // ndef NO_NAMESPACE_GINAC ////////// // absolute value ////////// -static ex abs_evalf(const ex & x) +static ex abs_evalf(const ex & arg) { - BEGIN_TYPECHECK - TYPECHECK(x,numeric) - END_TYPECHECK(abs(x)) - - return abs(ex_to_numeric(x)); + if (is_exactly_a(arg)) + return abs(ex_to(arg)); + + return abs(arg).hold(); } -static ex abs_eval(const ex & x) +static ex abs_eval(const ex & arg) { - if (is_ex_exactly_of_type(x, numeric)) - return abs(ex_to_numeric(x)); - else - return abs(x).hold(); + if (is_exactly_a(arg)) + return abs(ex_to(arg)); + else + return abs(arg).hold(); +} + +static void abs_print_latex(const ex & arg, const print_context & c) +{ + c.s << "{|"; arg.print(c); c.s << "|}"; +} + +static void abs_print_csrc_float(const ex & arg, const print_context & c) +{ + c.s << "fabs("; arg.print(c); c.s << ")"; } REGISTER_FUNCTION(abs, eval_func(abs_eval). - evalf_func(abs_evalf)); + evalf_func(abs_evalf). + print_func(abs_print_latex). + print_func(abs_print_csrc_float). + print_func(abs_print_csrc_float)); ////////// // Complex sign ////////// -static ex csgn_evalf(const ex & x) +static ex csgn_evalf(const ex & arg) { - BEGIN_TYPECHECK - TYPECHECK(x,numeric) - END_TYPECHECK(csgn(x)) - - return csgn(ex_to_numeric(x)); + if (is_exactly_a(arg)) + return csgn(ex_to(arg)); + + return csgn(arg).hold(); } -static ex csgn_eval(const ex & x) +static ex csgn_eval(const ex & arg) { - if (is_ex_exactly_of_type(x, numeric)) - return csgn(ex_to_numeric(x)); - - else if (is_ex_exactly_of_type(x, mul)) { - numeric oc = ex_to_numeric(x.op(x.nops()-1)); - if (oc.is_real()) { - if (oc > 0) - // csgn(42*x) -> csgn(x) - return csgn(x/oc).hold(); - else - // csgn(-42*x) -> -csgn(x) - return -csgn(x/oc).hold(); - } - if (oc.real().is_zero()) { - if (oc.imag() > 0) - // csgn(42*I*x) -> csgn(I*x) - return csgn(I*x/oc).hold(); - else - // csgn(-42*I*x) -> -csgn(I*x) - return -csgn(I*x/oc).hold(); - } + if (is_exactly_a(arg)) + return csgn(ex_to(arg)); + + else if (is_exactly_a(arg) && + is_exactly_a(arg.op(arg.nops()-1))) { + numeric oc = ex_to(arg.op(arg.nops()-1)); + if (oc.is_real()) { + if (oc > 0) + // csgn(42*x) -> csgn(x) + return csgn(arg/oc).hold(); + else + // csgn(-42*x) -> -csgn(x) + return -csgn(arg/oc).hold(); + } + if (oc.real().is_zero()) { + if (oc.imag() > 0) + // csgn(42*I*x) -> csgn(I*x) + return csgn(I*arg/oc).hold(); + else + // csgn(-42*I*x) -> -csgn(I*x) + return -csgn(I*arg/oc).hold(); + } } - - return csgn(x).hold(); + + return csgn(arg).hold(); } -static ex csgn_series(const ex & x, const relational & rel, int order) +static ex csgn_series(const ex & arg, + const relational & rel, + int order, + unsigned options) { - const ex x_pt = x.subs(rel); - if (x_pt.info(info_flags::numeric)) { - if (ex_to_numeric(x_pt).real().is_zero()) - throw (std::domain_error("csgn_series(): on imaginary axis")); - epvector seq; - seq.push_back(expair(csgn(x_pt), _ex0())); - return pseries(rel,seq); - } - epvector seq; - seq.push_back(expair(csgn(x_pt), _ex0())); - return pseries(rel,seq); + const ex arg_pt = arg.subs(rel, subs_options::no_pattern); + if (arg_pt.info(info_flags::numeric) + && ex_to(arg_pt).real().is_zero() + && !(options & series_options::suppress_branchcut)) + throw (std::domain_error("csgn_series(): on imaginary axis")); + + epvector seq; + seq.push_back(expair(csgn(arg_pt), _ex0)); + return pseries(rel,seq); } REGISTER_FUNCTION(csgn, eval_func(csgn_eval). evalf_func(csgn_evalf). series_func(csgn_series)); + +////////// +// Eta function: eta(x,y) == log(x*y) - log(x) - log(y). +// This function is closely related to the unwinding number K, sometimes found +// in modern literature: K(z) == (z-log(exp(z)))/(2*Pi*I). +////////// + +static ex eta_evalf(const ex &x, const ex &y) +{ + // It seems like we basically have to replicate the eval function here, + // since the expression might not be fully evaluated yet. + if (x.info(info_flags::positive) || y.info(info_flags::positive)) + return _ex0; + + if (x.info(info_flags::numeric) && y.info(info_flags::numeric)) { + const numeric nx = ex_to(x); + const numeric ny = ex_to(y); + const numeric nxy = ex_to(x*y); + int cut = 0; + if (nx.is_real() && nx.is_negative()) + cut -= 4; + if (ny.is_real() && ny.is_negative()) + cut -= 4; + if (nxy.is_real() && nxy.is_negative()) + cut += 4; + return evalf(I/4*Pi)*((csgn(-imag(nx))+1)*(csgn(-imag(ny))+1)*(csgn(imag(nxy))+1)- + (csgn(imag(nx))+1)*(csgn(imag(ny))+1)*(csgn(-imag(nxy))+1)+cut); + } + + return eta(x,y).hold(); +} + +static ex eta_eval(const ex &x, const ex &y) +{ + // trivial: eta(x,c) -> 0 if c is real and positive + if (x.info(info_flags::positive) || y.info(info_flags::positive)) + return _ex0; + + if (x.info(info_flags::numeric) && y.info(info_flags::numeric)) { + // don't call eta_evalf here because it would call Pi.evalf()! + const numeric nx = ex_to(x); + const numeric ny = ex_to(y); + const numeric nxy = ex_to(x*y); + int cut = 0; + if (nx.is_real() && nx.is_negative()) + cut -= 4; + if (ny.is_real() && ny.is_negative()) + cut -= 4; + if (nxy.is_real() && nxy.is_negative()) + cut += 4; + return (I/4)*Pi*((csgn(-imag(nx))+1)*(csgn(-imag(ny))+1)*(csgn(imag(nxy))+1)- + (csgn(imag(nx))+1)*(csgn(imag(ny))+1)*(csgn(-imag(nxy))+1)+cut); + } + + return eta(x,y).hold(); +} + +static ex eta_series(const ex & x, const ex & y, + const relational & rel, + int order, + unsigned options) +{ + const ex x_pt = x.subs(rel, subs_options::no_pattern); + const ex y_pt = y.subs(rel, subs_options::no_pattern); + if ((x_pt.info(info_flags::numeric) && x_pt.info(info_flags::negative)) || + (y_pt.info(info_flags::numeric) && y_pt.info(info_flags::negative)) || + ((x_pt*y_pt).info(info_flags::numeric) && (x_pt*y_pt).info(info_flags::negative))) + throw (std::domain_error("eta_series(): on discontinuity")); + epvector seq; + seq.push_back(expair(eta(x_pt,y_pt), _ex0)); + return pseries(rel,seq); +} + +REGISTER_FUNCTION(eta, eval_func(eta_eval). + evalf_func(eta_evalf). + series_func(eta_series). + latex_name("\\eta"). + set_symmetry(sy_symm(0, 1))); + + ////////// // dilogarithm ////////// +static ex Li2_evalf(const ex & x) +{ + if (is_exactly_a(x)) + return Li2(ex_to(x)); + + return Li2(x).hold(); +} + static ex Li2_eval(const ex & x) { - // Li2(0) -> 0 - if (x.is_zero()) - return x; - // Li2(1) -> Pi^2/6 - if (x.is_equal(_ex1())) - return power(Pi,_ex2())/_ex6(); - // Li2(1/2) -> Pi^2/12 - log(2)^2/2 - if (x.is_equal(_ex1_2())) - return power(Pi,_ex2())/_ex12() + power(log(_ex2()),_ex2())*_ex_1_2(); - // Li2(-1) -> -Pi^2/12 - if (x.is_equal(_ex_1())) - return -power(Pi,_ex2())/_ex12(); - // Li2(I) -> -Pi^2/48+Catalan*I - if (x.is_equal(I)) - return power(Pi,_ex2())/_ex_48() + Catalan*I; - // Li2(-I) -> -Pi^2/48-Catalan*I - if (x.is_equal(-I)) - return power(Pi,_ex2())/_ex_48() - Catalan*I; - return Li2(x).hold(); + if (x.info(info_flags::numeric)) { + // Li2(0) -> 0 + if (x.is_zero()) + return _ex0; + // Li2(1) -> Pi^2/6 + if (x.is_equal(_ex1)) + return power(Pi,_ex2)/_ex6; + // Li2(1/2) -> Pi^2/12 - log(2)^2/2 + if (x.is_equal(_ex1_2)) + return power(Pi,_ex2)/_ex12 + power(log(_ex2),_ex2)*_ex_1_2; + // Li2(-1) -> -Pi^2/12 + if (x.is_equal(_ex_1)) + return -power(Pi,_ex2)/_ex12; + // Li2(I) -> -Pi^2/48+Catalan*I + if (x.is_equal(I)) + return power(Pi,_ex2)/_ex_48 + Catalan*I; + // Li2(-I) -> -Pi^2/48-Catalan*I + if (x.is_equal(-I)) + return power(Pi,_ex2)/_ex_48 - Catalan*I; + // Li2(float) + if (!x.info(info_flags::crational)) + return Li2(ex_to(x)); + } + + return Li2(x).hold(); } static ex Li2_deriv(const ex & x, unsigned deriv_param) { - GINAC_ASSERT(deriv_param==0); - - // d/dx Li2(x) -> -log(1-x)/x - return -log(1-x)/x; + GINAC_ASSERT(deriv_param==0); + + // d/dx Li2(x) -> -log(1-x)/x + return -log(_ex1-x)/x; } -static ex Li2_series(const ex &x, const relational &rel, int order) +static ex Li2_series(const ex &x, const relational &rel, int order, unsigned options) { - const ex x_pt = x.subs(rel); - if (!x_pt.is_zero() && !x_pt.is_equal(_ex1())) - throw do_taylor(); // caught by function::series() - // First case: x==0 (derivatives have poles) - if (x_pt.is_zero()) { - // method: - // The problem is that in d/dx Li2(x==0) == -log(1-x)/x we cannot - // simply substitute x==0. The limit, however, exists: it is 1. We - // also know all higher derivatives' limits: (d/dx)^n Li2(x) == n!/n^2. - // So the primitive series expansion is Li2(x==0) == x + x^2/4 + x^3/9 - // and so on. - // We first construct such a primitive series expansion manually in - // a dummy symbol s and then insert the argument's series expansion - // for s. Reexpanding the resulting series returns the desired result. - const symbol s; - ex ser; - // construct manually the primitive expansion - for (int i=1; i=1 (branch cut) - return pseries(); - // TODO: Li2_series should do something around branch point? - // Careful: may involve logs! + const ex x_pt = x.subs(rel, subs_options::no_pattern); + if (x_pt.info(info_flags::numeric)) { + // First special case: x==0 (derivatives have poles) + if (x_pt.is_zero()) { + // method: + // The problem is that in d/dx Li2(x==0) == -log(1-x)/x we cannot + // simply substitute x==0. The limit, however, exists: it is 1. + // We also know all higher derivatives' limits: + // (d/dx)^n Li2(x) == n!/n^2. + // So the primitive series expansion is + // Li2(x==0) == x + x^2/4 + x^3/9 + ... + // and so on. + // We first construct such a primitive series expansion manually in + // a dummy symbol s and then insert the argument's series expansion + // for s. Reexpanding the resulting series returns the desired + // result. + const symbol s; + ex ser; + // manually construct the primitive expansion + for (int i=1; i=1 (branch cut) + if (!(options & series_options::suppress_branchcut) && + ex_to(x_pt).is_real() && ex_to(x_pt)>1) { + // method: + // This is the branch cut: assemble the primitive series manually + // and then add the corresponding complex step function. + const symbol &s = ex_to(rel.lhs()); + const ex point = rel.rhs(); + const symbol foo; + epvector seq; + // zeroth order term: + seq.push_back(expair(Li2(x_pt), _ex0)); + // compute the intermediate terms: + ex replarg = series(Li2(x), s==foo, order); + for (size_t i=1; i zeta(x) + if (n.is_zero()) + return zeta(x); + } + + return zetaderiv(n, x).hold(); +} + +static ex zetaderiv_deriv(const ex & n, const ex & x, unsigned deriv_param) +{ + GINAC_ASSERT(deriv_param<2); + + if (deriv_param==0) { + // d/dn zeta(n,x) + throw(std::logic_error("cannot diff zetaderiv(n,x) with respect to n")); + } + // d/dx psi(n,x) + return zetaderiv(n+1,x); +} + +REGISTER_FUNCTION(zetaderiv, eval_func(zetaderiv_eval). + derivative_func(zetaderiv_deriv). + latex_name("\\zeta^\\prime")); ////////// // factorial @@ -226,15 +405,15 @@ REGISTER_FUNCTION(Li3, eval_func(Li3_eval)); static ex factorial_evalf(const ex & x) { - return factorial(x).hold(); + return factorial(x).hold(); } static ex factorial_eval(const ex & x) { - if (is_ex_exactly_of_type(x, numeric)) - return factorial(ex_to_numeric(x)); - else - return factorial(x).hold(); + if (is_exactly_a(x)) + return factorial(ex_to(x)); + else + return factorial(x).hold(); } REGISTER_FUNCTION(factorial, eval_func(factorial_eval). @@ -246,15 +425,15 @@ REGISTER_FUNCTION(factorial, eval_func(factorial_eval). static ex binomial_evalf(const ex & x, const ex & y) { - return binomial(x, y).hold(); + return binomial(x, y).hold(); } static ex binomial_eval(const ex & x, const ex &y) { - if (is_ex_exactly_of_type(x, numeric) && is_ex_exactly_of_type(y, numeric)) - return binomial(ex_to_numeric(x), ex_to_numeric(y)); - else - return binomial(x, y).hold(); + if (is_exactly_a(x) && is_exactly_a(y)) + return binomial(ex_to(x), ex_to(y)); + else + return binomial(x, y).hold(); } REGISTER_FUNCTION(binomial, eval_func(binomial_eval). @@ -266,169 +445,121 @@ REGISTER_FUNCTION(binomial, eval_func(binomial_eval). static ex Order_eval(const ex & x) { - if (is_ex_exactly_of_type(x, numeric)) { - - // O(c)=O(1) - return Order(_ex1()).hold(); - - } else if (is_ex_exactly_of_type(x, mul)) { - - mul *m = static_cast(x.bp); - if (is_ex_exactly_of_type(m->op(m->nops() - 1), numeric)) { - - // O(c*expr)=O(expr) - return Order(x / m->op(m->nops() - 1)).hold(); - } + if (is_exactly_a(x)) { + // O(c) -> O(1) or 0 + if (!x.is_zero()) + return Order(_ex1).hold(); + else + return _ex0; + } else if (is_exactly_a(x)) { + const mul &m = ex_to(x); + // O(c*expr) -> O(expr) + if (is_exactly_a(m.op(m.nops() - 1))) + return Order(x / m.op(m.nops() - 1)).hold(); } return Order(x).hold(); } -static ex Order_series(const ex & x, const relational & r, int order) +static ex Order_series(const ex & x, const relational & r, int order, unsigned options) { // Just wrap the function into a pseries object epvector new_seq; - GINAC_ASSERT(is_ex_exactly_of_type(r.lhs(),symbol)); - const symbol *s = static_cast(r.lhs().bp); - new_seq.push_back(expair(Order(_ex1()), numeric(min(x.ldegree(*s), order)))); + GINAC_ASSERT(is_a(r.lhs())); + const symbol &s = ex_to(r.lhs()); + new_seq.push_back(expair(Order(_ex1), numeric(std::min(x.ldegree(s), order)))); return pseries(r, new_seq); } // Differentiation is handled in function::derivative because of its special requirements REGISTER_FUNCTION(Order, eval_func(Order_eval). - series_func(Order_series)); - -////////// -// Inert partial differentiation operator -////////// - -static ex Derivative_eval(const ex & f, const ex & l) -{ - if (!is_ex_exactly_of_type(f, function)) { - throw(std::invalid_argument("Derivative(): 1st argument must be a function")); - } - if (!is_ex_exactly_of_type(l, lst)) { - throw(std::invalid_argument("Derivative(): 2nd argument must be a list")); - } - return Derivative(f, l).hold(); -} - -REGISTER_FUNCTION(Derivative, eval_func(Derivative_eval)); + series_func(Order_series). + latex_name("\\mathcal{O}")); ////////// // Solve linear system ////////// -ex lsolve(const ex &eqns, const ex &symbols) -{ - // solve a system of linear equations - if (eqns.info(info_flags::relation_equal)) { - if (!symbols.info(info_flags::symbol)) { - throw(std::invalid_argument("lsolve: 2nd argument must be a symbol")); - } - ex sol=lsolve(lst(eqns),lst(symbols)); - - GINAC_ASSERT(sol.nops()==1); - GINAC_ASSERT(is_ex_exactly_of_type(sol.op(0),relational)); - - return sol.op(0).op(1); // return rhs of first solution - } - - // syntax checks - if (!eqns.info(info_flags::list)) { - throw(std::invalid_argument("lsolve: 1st argument must be a list")); - } - for (unsigned i=0; i(sol.op(0))); + + return sol.op(0).op(1); // return rhs of first solution + } + + // syntax checks + if (!eqns.info(info_flags::list)) { + throw(std::invalid_argument("lsolve(): 1st argument must be a list")); + } + for (size_t i=0; i(symbols.op(c)),1); + linpart -= co*symbols.op(c); + sys(r,c) = co; + } + linpart = linpart.expand(); + rhs(r,0) = -linpart; + } + + // test if system is linear and fill vars matrix + for (size_t i=0; i