X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Finifcns.cpp;h=decef1a8806abd633420ecf5798b46ff7f2bc135;hp=7081eb5b333ebe22f7fd365e80a0a513d60e0603;hb=9606fc35b00ef890f632f98eddcca26631917e10;hpb=9961dfdd7a383f09b1040e195a8817a85d945e58 diff --git a/ginac/inifcns.cpp b/ginac/inifcns.cpp index 7081eb5b..decef1a8 100644 --- a/ginac/inifcns.cpp +++ b/ginac/inifcns.cpp @@ -3,7 +3,7 @@ * Implementation of GiNaC's initially known functions. */ /* - * GiNaC Copyright (C) 1999-2001 Johannes Gutenberg University Mainz, Germany + * GiNaC Copyright (C) 1999-2004 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,39 +29,90 @@ #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" namespace GiNaC { +////////// +// complex conjugate +////////// + +static ex conjugate_evalf(const ex & arg) +{ + if (is_exactly_a(arg)) { + return ex_to(arg).conjugate(); + } + return conjugate_function(arg).hold(); +} + +static ex conjugate_eval(const ex & arg) +{ + return arg.conjugate(); +} + +static void conjugate_print_latex(const ex & arg, const print_context & c) +{ + c.s << "\\bar{"; arg.print(c); c.s << "}"; +} + +static ex conjugate_conjugate(const ex & arg) +{ + return arg; +} + +REGISTER_FUNCTION(conjugate_function, eval_func(conjugate_eval). + evalf_func(conjugate_evalf). + print_func(conjugate_print_latex). + conjugate_func(conjugate_conjugate). + set_name("conjugate","conjugate")); + ////////// // absolute value ////////// static ex abs_evalf(const ex & arg) { - BEGIN_TYPECHECK - TYPECHECK(arg,numeric) - END_TYPECHECK(abs(arg)) + if (is_exactly_a(arg)) + return abs(ex_to(arg)); - return abs(ex_to_numeric(arg)); + return abs(arg).hold(); } static ex abs_eval(const ex & arg) { - if (is_ex_exactly_of_type(arg, numeric)) - return abs(ex_to_numeric(arg)); + 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 << ")"; +} + +static ex abs_conjugate(const ex & arg) +{ + return abs(arg); +} + 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). + conjugate_func(abs_conjugate)); ////////// @@ -70,21 +121,20 @@ REGISTER_FUNCTION(abs, eval_func(abs_eval). static ex csgn_evalf(const ex & arg) { - BEGIN_TYPECHECK - TYPECHECK(arg,numeric) - END_TYPECHECK(csgn(arg)) + if (is_exactly_a(arg)) + return csgn(ex_to(arg)); - return csgn(ex_to_numeric(arg)); + return csgn(arg).hold(); } static ex csgn_eval(const ex & arg) { - if (is_ex_exactly_of_type(arg, numeric)) - return csgn(ex_to_numeric(arg)); + if (is_exactly_a(arg)) + return csgn(ex_to(arg)); - else if (is_ex_of_type(arg, mul) && - is_ex_of_type(arg.op(arg.nops()-1),numeric)) { - numeric oc = ex_to_numeric(arg.op(arg.nops()-1)); + 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) @@ -111,75 +161,111 @@ static ex csgn_series(const ex & arg, int order, unsigned options) { - const ex arg_pt = arg.subs(rel); + const ex arg_pt = arg.subs(rel, subs_options::no_pattern); if (arg_pt.info(info_flags::numeric) - && ex_to_numeric(arg_pt).real().is_zero() + && 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())); + seq.push_back(expair(csgn(arg_pt), _ex0)); return pseries(rel,seq); } +static ex csgn_conjugate(const ex& arg) +{ + return csgn(arg); +} + REGISTER_FUNCTION(csgn, eval_func(csgn_eval). evalf_func(csgn_evalf). - series_func(csgn_series)); + series_func(csgn_series). + conjugate_func(csgn_conjugate)); ////////// -// Eta function: log(x*y) == log(x) + log(y) + eta(x,y). +// 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) +static ex eta_evalf(const ex &x, const ex &y) { - BEGIN_TYPECHECK - TYPECHECK(x,numeric) - TYPECHECK(y,numeric) - END_TYPECHECK(eta(x,y)) - - numeric xim = imag(ex_to_numeric(x)); - numeric yim = imag(ex_to_numeric(y)); - numeric xyim = imag(ex_to_numeric(x*y)); - return evalf(I/4*Pi)*((csgn(-xim)+1)*(csgn(-yim)+1)*(csgn(xyim)+1)-(csgn(xim)+1)*(csgn(yim)+1)*(csgn(-xyim)+1)); + // 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) +static ex eta_eval(const ex &x, const ex &y) { - if (is_ex_exactly_of_type(x, numeric) && - is_ex_exactly_of_type(y, numeric)) { + // 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()! - numeric xim = imag(ex_to_numeric(x)); - numeric yim = imag(ex_to_numeric(y)); - numeric xyim = imag(ex_to_numeric(x*y)); - return (I/4)*Pi*((csgn(-xim)+1)*(csgn(-yim)+1)*(csgn(xyim)+1)-(csgn(xim)+1)*(csgn(yim)+1)*(csgn(-xyim)+1)); + 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 & arg1, - const ex & arg2, +static ex eta_series(const ex & x, const ex & y, const relational & rel, int order, unsigned options) { - const ex arg1_pt = arg1.subs(rel); - const ex arg2_pt = arg2.subs(rel); - if (ex_to_numeric(arg1_pt).imag().is_zero() || - ex_to_numeric(arg2_pt).imag().is_zero() || - ex_to_numeric(arg1_pt*arg2_pt).imag().is_zero()) { - throw (std::domain_error("eta_series(): on discontinuity")); - } + 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(arg1_pt,arg2_pt), _ex0())); + seq.push_back(expair(eta(x_pt,y_pt), _ex0)); return pseries(rel,seq); } +static ex eta_conjugate(const ex & x, const ex & y) +{ + return -eta(x,y); +} + REGISTER_FUNCTION(eta, eval_func(eta_eval). evalf_func(eta_evalf). series_func(eta_series). - latex_name("\\eta")); + latex_name("\\eta"). + set_symmetry(sy_symm(0, 1)). + conjugate_func(eta_conjugate)); ////////// @@ -188,11 +274,10 @@ REGISTER_FUNCTION(eta, eval_func(eta_eval). static ex Li2_evalf(const ex & x) { - BEGIN_TYPECHECK - TYPECHECK(x,numeric) - END_TYPECHECK(Li2(x)) + if (is_exactly_a(x)) + return Li2(ex_to(x)); - return Li2(ex_to_numeric(x)); // -> numeric Li2(numeric) + return Li2(x).hold(); } static ex Li2_eval(const ex & x) @@ -200,25 +285,25 @@ static ex Li2_eval(const ex & x) if (x.info(info_flags::numeric)) { // Li2(0) -> 0 if (x.is_zero()) - return _ex0(); + return _ex0; // Li2(1) -> Pi^2/6 - if (x.is_equal(_ex1())) - return power(Pi,_ex2())/_ex6(); + 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(); + 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(); + 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; + 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 power(Pi,_ex2)/_ex_48 - Catalan*I; // Li2(float) if (!x.info(info_flags::crational)) - return Li2_evalf(x); + return Li2(ex_to(x)); } return Li2(x).hold(); @@ -229,12 +314,12 @@ 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; + return -log(_ex1-x)/x; } static ex Li2_series(const ex &x, const relational &rel, int order, unsigned options) { - const ex x_pt = x.subs(rel); + 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()) { @@ -254,12 +339,12 @@ static ex Li2_series(const ex &x, const relational &rel, int order, unsigned opt ex ser; // manually construct the primitive expansion for (int i=1; i=1 (branch cut) if (!(options & series_options::suppress_branchcut) && - ex_to_numeric(x_pt).is_real() && ex_to_numeric(x_pt)>1) { + 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 = static_cast(rel.lhs().bp); + 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())); + seq.push_back(expair(Li2(x_pt), _ex0)); // compute the intermediate terms: - ex replarg = series(Li2(x), *s==foo, order); - for (unsigned 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 ////////// @@ -345,14 +461,20 @@ static ex factorial_evalf(const ex & x) static ex factorial_eval(const ex & x) { - if (is_ex_exactly_of_type(x, numeric)) - return factorial(ex_to_numeric(x)); + if (is_exactly_a(x)) + return factorial(ex_to(x)); else return factorial(x).hold(); } +static ex factorial_conjugate(const ex & x) +{ + return factorial(x); +} + REGISTER_FUNCTION(factorial, eval_func(factorial_eval). - evalf_func(factorial_evalf)); + evalf_func(factorial_evalf). + conjugate_func(factorial_conjugate)); ////////// // binomial @@ -363,16 +485,46 @@ static ex binomial_evalf(const ex & x, const ex & y) return binomial(x, y).hold(); } +static ex binomial_sym(const ex & x, const numeric & y) +{ + if (y.is_integer()) { + if (y.is_nonneg_integer()) { + const unsigned N = y.to_int(); + if (N == 0) return _num0; + if (N == 1) return x; + ex t = x.expand(); + for (unsigned i = 2; i <= N; ++i) + t = (t * (x + i - y - 1)).expand() / i; + return t; + } else + return _num0; + } + + 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 + if (is_exactly_a(y)) { + if (is_exactly_a(x) && ex_to(x).is_integer()) + return binomial(ex_to(x), ex_to(y)); + else + return binomial_sym(x, ex_to(y)); + } else return binomial(x, y).hold(); } +// At the moment the numeric evaluation of a binomail function always +// gives a real number, but if this would be implemented using the gamma +// function, also complex conjugation should be changed (or rather, deleted). +static ex binomial_conjugate(const ex & x, const ex & y) +{ + return binomial(x,y); +} + REGISTER_FUNCTION(binomial, eval_func(binomial_eval). - evalf_func(binomial_evalf)); + evalf_func(binomial_evalf). + conjugate_func(binomial_conjugate)); ////////// // Order term function (for truncated power series) @@ -380,17 +532,17 @@ REGISTER_FUNCTION(binomial, eval_func(binomial_eval). static ex Order_eval(const ex & x) { - if (is_ex_exactly_of_type(x, numeric)) { + if (is_exactly_a(x)) { // O(c) -> O(1) or 0 if (!x.is_zero()) - return Order(_ex1()).hold(); + return Order(_ex1).hold(); else - return _ex0(); - } else if (is_ex_exactly_of_type(x, mul)) { - mul *m = static_cast(x.bp); + return _ex0; + } else if (is_exactly_a(x)) { + const mul &m = ex_to(x); // O(c*expr) -> O(expr) - if (is_ex_exactly_of_type(m->op(m->nops() - 1), numeric)) - return Order(x / m->op(m->nops() - 1)).hold(); + if (is_exactly_a(m.op(m.nops() - 1))) + return Order(x / m.op(m.nops() - 1)).hold(); } return Order(x).hold(); } @@ -399,49 +551,38 @@ static ex Order_series(const ex & x, const relational & r, int order, unsigned o { // 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(std::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); } +static ex Order_conjugate(const ex & x) +{ + return Order(x); +} + // Differentiation is handled in function::derivative because of its special requirements REGISTER_FUNCTION(Order, eval_func(Order_eval). series_func(Order_series). - latex_name("\\mathcal{O}")); - -////////// -// 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)); + latex_name("\\mathcal{O}"). + conjugate_func(Order_conjugate)); ////////// // Solve linear system ////////// -ex lsolve(const ex &eqns, const ex &symbols) +ex lsolve(const ex &eqns, const ex &symbols, unsigned options) { // 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)); + const ex sol = lsolve(lst(eqns),lst(symbols)); GINAC_ASSERT(sol.nops()==1); - GINAC_ASSERT(is_ex_exactly_of_type(sol.op(0),relational)); + GINAC_ASSERT(is_exactly_a(sol.op(0))); return sol.op(0).op(1); // return rhs of first solution } @@ -450,7 +591,7 @@ ex lsolve(const ex &eqns, const ex &symbols) if (!eqns.info(info_flags::list)) { throw(std::invalid_argument("lsolve(): 1st argument must be a list")); } - for (unsigned i=0; i(symbols.op(c)),1); linpart -= co*symbols.op(c); - sys.set(r,c,co); + sys(r,c) = co; } linpart = linpart.expand(); - rhs.set(r,0,-linpart); + rhs(r,0) = -linpart; } // test if system is linear and fill vars matrix - for (unsigned i=0; i