From 951f99eefc0f5c89adf7d76fbdfd01b5a2a81f3e Mon Sep 17 00:00:00 2001 From: "Vladimir V. Kisil" Date: Wed, 23 Oct 2013 10:26:00 +0100 Subject: [PATCH] Expansion laws for exp, log and abs Expansion laws for exponent (exp(a+b)=exp(a) exp(b)), logarithm (log(ab)=log(a)log(b), if at least one a or b is positive) and absolute value (|zw|=|z| |w|). All expansions are performed if expand_options::expand_transcendental is given. Signed-off-by: Vladimir V. Kisil --- check/exam_inifcns.cpp | 67 ++++++++++++++++++++++++++++++++++++++++++- doc/tutorial/ginac.texi | 53 ++++++++++++++++++++++++++++++++++ ginac/flags.h | 3 +- ginac/inifcns.cpp | 22 ++++++++++++++ ginac/inifcns_trans.cpp | 76 +++++++++++++++++++++++++++++++++++++++++++++++++ 5 files changed, 219 insertions(+), 2 deletions(-) diff --git a/check/exam_inifcns.cpp b/check/exam_inifcns.cpp index 30d5b21..820c3e4 100644 --- a/check/exam_inifcns.cpp +++ b/check/exam_inifcns.cpp @@ -243,6 +243,10 @@ static unsigned inifcns_consist_abs() if (!abs(pow(x+I*y,a+I*b)).eval().is_equal(abs(pow(x+I*y,a+I*b)))) ++result; + // check expansion of abs + if (!abs(-7*z*a*p).expand(expand_options::expand_transcendental).is_equal(7*abs(z)*abs(a)*p)) + ++result; + if (!abs(z.conjugate()).eval().is_equal(abs(z))) ++result; @@ -261,6 +265,65 @@ static unsigned inifcns_consist_abs() return result; } +static unsigned inifcns_consist_exp() +{ + unsigned result = 0; + symbol a("a"), b("b"); + + if (!exp(a+b).expand(expand_options::expand_transcendental).is_equal(exp(a)*exp(b))) + ++result; + + // shall not be expanded since the arg is not add + if (!exp(pow(a+b,2)).expand(expand_options::expand_transcendental).is_equal(exp(pow(a+b,2)))) + ++result; + + // expand now + if (!exp(pow(a+b,2)).expand(expand_options::expand_function_args | expand_options::expand_transcendental) + .is_equal(exp(a*a)*exp(b*b)*exp(2*a*b))) + ++result; + + return result; +} + +static unsigned inifcns_consist_log() +{ + unsigned result = 0; + symbol z("a"), w("b"); + realsymbol a("a"), b("b"); + possymbol p("p"), q("q"); + + // do not expand + if (!log(z*w).expand(expand_options::expand_transcendental).is_equal(log(z*w))) + ++result; + + // do not expand + if (!log(a*b).expand(expand_options::expand_transcendental).is_equal(log(a*b))) + ++result; + + // shall expand + if (!log(p*q).expand(expand_options::expand_transcendental).is_equal(log(p) +log(q))) + ++result; + + // a bit more complicated + ex e1=log(-7*p*pow(q,3)*a*pow(b,2)*z*w).expand(expand_options::expand_transcendental); + ex e2=log(7)+log(p)+log(pow(q,3))+log(-z*a*w*pow(b,2)); + if (!e1.is_equal(e2)) + ++result; + + if (!ex(log(pow(p,a))).is_equal(a*log(p))) + ++result; + + // shall not do for non-real powers + if (ex(log(pow(p,z))).is_equal(z*log(p))) + ++result; + + // shall not do for non-positive basis + if (ex(log(pow(a,b))).is_equal(b*log(a))) + ++result; + + return result; +} + static unsigned inifcns_consist_various() { unsigned result = 0; @@ -270,7 +333,7 @@ static unsigned inifcns_consist_various() clog << "ERROR: binomial(n,0) != 1" << endl; ++result; } - + return result; } @@ -285,6 +348,8 @@ unsigned exam_inifcns() result += inifcns_consist_psi(); cout << '.' << flush; result += inifcns_consist_zeta(); cout << '.' << flush; result += inifcns_consist_abs(); cout << '.' << flush; + result += inifcns_consist_exp(); cout << '.' << flush; + result += inifcns_consist_log(); cout << '.' << flush; result += inifcns_consist_various(); cout << '.' << flush; return result; diff --git a/doc/tutorial/ginac.texi b/doc/tutorial/ginac.texi index 7324870..cac3cde 100644 --- a/doc/tutorial/ginac.texi +++ b/doc/tutorial/ginac.texi @@ -5926,6 +5926,59 @@ towards -infinity. These functions are continuous as the branch cut is approached coming around the finite endpoint of the cut in a counter clockwise direction. +@c +@subsection Expanding functions +@cindex expand trancedent functions +@cindex @code{expand_options::expand_transcendental} +@cindex @code{expand_options::expand_function_args} +GiNaC knows several expansion laws for trancedent functions, e.g. +@tex +$e^{a+b}=e^a e^b$, +$|zw|=|z|\cdot |w|$ +@end tex +@ifnottex +@command{exp(a+b)=exp(a) exp(b), |zw|=|z| |w|} +@end ifnottex +or +@tex +$\log(c*d)=\log(c)+\log(d)$, +@end tex +@ifnottex +@command{log(cd)=log(c)+log(d)} +@end ifnottex +(for positive +@tex +$c,\ d$ +@end tex +@ifnottex +@command{c, d} +@end ifnottex +). In order to use these rules you need to call @code{expand()} method +with the option @code{expand_options::expand_transcendental}. Another +relevant option is @code{expand_options::expand_function_args}. Their +usage and interaction can be seen from the following example: +@example +@{ + symbol x("x"), y("y"); + ex e=exp(pow(x+y,2)); + cout << e.expand() << endl; + // -> exp((x+y)^2) + cout << e.expand(expand_options::expand_transcendental) << endl; + // -> exp((x+y)^2) + cout << e.expand(expand_options::expand_function_args) << endl; + // -> exp(2*x*y+x^2+y^2) + cout << e.expand(expand_options::expand_function_args + | expand_options::expand_transcendental) << endl; + // -> exp(y^2)*exp(2*x*y)*exp(x^2) +@} +@end example +If both flags are set (as in the last call), then GiNaC tries to get +the maximal expansion. For example, for the exponent GiNaC firstly expands +the argument and then the function. For the logarithm and absolute value, +GiNaC uses the opposite order: firstly expands the function and then its +argument. Of course, a user can fine-tune this behaviour by sequential +calls of several @code{expand()} methods with desired flags. + @node Multiple polylogarithms, Complex expressions, Built-in functions, Methods and functions @c node-name, next, previous, up @subsection Multiple polylogarithms diff --git a/ginac/flags.h b/ginac/flags.h index 3d398cf..6da71b0 100644 --- a/ginac/flags.h +++ b/ginac/flags.h @@ -31,7 +31,8 @@ public: enum { expand_indexed = 0x0001, ///< expands (a+b).i to a.i+b.i expand_function_args = 0x0002, ///< expands the arguments of functions - expand_rename_idx = 0x0004 ///< used internally by mul::expand() + expand_rename_idx = 0x0004, ///< used internally by mul::expand() + expand_transcendental = 0x0008 ///< expands trancendental functions like log and exp }; }; diff --git a/ginac/inifcns.cpp b/ginac/inifcns.cpp index 1144ce1..153a3c1 100644 --- a/ginac/inifcns.cpp +++ b/ginac/inifcns.cpp @@ -254,6 +254,27 @@ static ex abs_eval(const ex & arg) return abs(arg).hold(); } +static ex abs_expand(const ex & arg, unsigned options) +{ + if ((options & expand_options::expand_transcendental) + && is_exactly_a(arg)) { + exvector prodseq; + prodseq.reserve(arg.nops()); + for (const_iterator i = arg.begin(); i != arg.end(); ++i) { + if (options & expand_options::expand_function_args) + prodseq.push_back(abs(i->expand(options))); + else + prodseq.push_back(abs(*i)); + } + return (new mul(prodseq))->setflag(status_flags::dynallocated | status_flags::expanded); + } + + if (options & expand_options::expand_function_args) + return abs(arg.expand(options)).hold(); + else + return abs(arg).hold(); +} + static void abs_print_latex(const ex & arg, const print_context & c) { c.s << "{|"; arg.print(c); c.s << "|}"; @@ -317,6 +338,7 @@ bool abs_info(const ex & arg, unsigned inf) REGISTER_FUNCTION(abs, eval_func(abs_eval). evalf_func(abs_evalf). + expand_func(abs_expand). info_func(abs_info). print_func(abs_print_latex). print_func(abs_print_csrc_float). diff --git a/ginac/inifcns_trans.cpp b/ginac/inifcns_trans.cpp index 10a3675..422fe39 100644 --- a/ginac/inifcns_trans.cpp +++ b/ginac/inifcns_trans.cpp @@ -24,6 +24,8 @@ #include "inifcns.h" #include "ex.h" #include "constant.h" +#include "add.h" +#include "mul.h" #include "numeric.h" #include "power.h" #include "operators.h" @@ -81,6 +83,27 @@ static ex exp_eval(const ex & x) return exp(x).hold(); } +static ex exp_expand(const ex & arg, unsigned options) +{ + ex exp_arg; + if (options & expand_options::expand_function_args) + exp_arg = arg.expand(options); + else + exp_arg=arg; + + if ((options & expand_options::expand_transcendental) + && is_exactly_a(exp_arg)) { + exvector prodseq; + prodseq.reserve(exp_arg.nops()); + for (const_iterator i = exp_arg.begin(); i != exp_arg.end(); ++i) + prodseq.push_back(exp(*i)); + + return (new mul(prodseq))->setflag(status_flags::dynallocated | status_flags::expanded); + } + + return exp(exp_arg).hold(); +} + static ex exp_deriv(const ex & x, unsigned deriv_param) { GINAC_ASSERT(deriv_param==0); @@ -107,6 +130,7 @@ static ex exp_conjugate(const ex & x) REGISTER_FUNCTION(exp, eval_func(exp_eval). evalf_func(exp_evalf). + expand_func(exp_expand). derivative_func(exp_deriv). real_part_func(exp_real_part). imag_part_func(exp_imag_part). @@ -150,6 +174,10 @@ static ex log_eval(const ex & x) if (t.info(info_flags::real)) return t; } + + // log(pow(a,b)) -> b*log(a) (if a>0 and b is real) + if (is_exactly_a(x) && x.op(0).info(info_flags::positive) && x.op(1).info(info_flags::real)) + return x.op(1)*log(x.op(0)); return log(x).hold(); } @@ -265,6 +293,53 @@ static ex log_imag_part(const ex & x) return atan2(GiNaC::imag_part(x), GiNaC::real_part(x)); } +static ex log_expand(const ex & arg, unsigned options) +{ + if ((options & expand_options::expand_transcendental) + && is_exactly_a(arg) && !arg.info(info_flags::indefinite)) { + exvector sumseq; + exvector prodseq; + sumseq.reserve(arg.nops()); + prodseq.reserve(arg.nops()); + bool possign=true; + + // searching for positive/negative factors + for (const_iterator i = arg.begin(); i != arg.end(); ++i) { + ex e; + if (options & expand_options::expand_function_args) + e=i->expand(options); + else + e=*i; + if (e.info(info_flags::positive)) + sumseq.push_back(log(e)); + else if (e.info(info_flags::negative)) { + sumseq.push_back(log(-e)); + possign = !possign; + } else + prodseq.push_back(e); + } + + if (sumseq.size() > 0) { + ex newarg; + if (options & expand_options::expand_function_args) + newarg=((possign?_ex1:_ex_1)*mul(prodseq)).expand(options); + else { + newarg=(possign?_ex1:_ex_1)*mul(prodseq); + ex_to(newarg).setflag(status_flags::purely_indefinite); + } + return add(sumseq)+log(newarg); + } else { + if (!(options & expand_options::expand_function_args)) + ex_to(arg).setflag(status_flags::purely_indefinite); + } + } + + if (options & expand_options::expand_function_args) + return log(arg.expand(options)).hold(); + else + return log(arg).hold(); +} + static ex log_conjugate(const ex & x) { // conjugate(log(x))==log(conjugate(x)) unless on the branch cut which @@ -281,6 +356,7 @@ static ex log_conjugate(const ex & x) REGISTER_FUNCTION(log, eval_func(log_eval). evalf_func(log_evalf). + expand_func(log_expand). derivative_func(log_deriv). series_func(log_series). real_part_func(log_real_part). -- 1.8.4.rc3