X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Finifcns.cpp;h=8103cd86bf3ec12ad02a2ebec48aff1b122b0c49;hp=20ea659ff52dee8286ce5ff09877f847a7c4a61e;hb=10365850aa3803337bfea1fc201b81b6752096d4;hpb=e710763e51b6fe11020bac880c44f426544471c2 diff --git a/ginac/inifcns.cpp b/ginac/inifcns.cpp index 20ea659f..8103cd86 100644 --- a/ginac/inifcns.cpp +++ b/ginac/inifcns.cpp @@ -3,7 +3,7 @@ * Implementation of GiNaC's initially known functions. */ /* - * GiNaC Copyright (C) 1999-2009 Johannes Gutenberg University Mainz, Germany + * GiNaC Copyright (C) 1999-2011 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 @@ -212,7 +212,7 @@ static void abs_print_csrc_float(const ex & arg, const print_context & c) static ex abs_conjugate(const ex & arg) { - return abs(arg); + return abs(arg).hold(); } static ex abs_real_part(const ex & arg) @@ -227,7 +227,8 @@ static ex abs_imag_part(const ex& arg) static ex abs_power(const ex & arg, const ex & exp) { - if (arg.is_equal(arg.conjugate()) && is_a(exp) && ex_to(exp).is_even()) + if (arg.is_equal(arg.conjugate()) && ((is_a(exp) && ex_to(exp).is_even()) + || exp.info(info_flags::even))) return power(arg, exp); else return power(abs(arg), exp).hold(); @@ -319,8 +320,8 @@ REGISTER_FUNCTION(step, eval_func(step_eval). evalf_func(step_evalf). series_func(step_series). conjugate_func(step_conjugate). - real_part_func(step_real_part). - imag_part_func(step_imag_part)); + real_part_func(step_real_part). + imag_part_func(step_imag_part)); ////////// // Complex sign @@ -398,7 +399,7 @@ static ex csgn_power(const ex & arg, const ex & exp) { if (is_a(exp) && exp.info(info_flags::positive) && ex_to(exp).is_integer()) { if (ex_to(exp).is_odd()) - return csgn(arg); + return csgn(arg).hold(); else return power(csgn(arg), _ex2).hold(); } else @@ -489,7 +490,7 @@ static ex eta_series(const ex & x, const ex & y, static ex eta_conjugate(const ex & x, const ex & y) { - return -eta(x, y); + return -eta(x, y).hold(); } static ex eta_real_part(const ex & x, const ex & y) @@ -643,10 +644,25 @@ static ex Li2_series(const ex &x, const relational &rel, int order, unsigned opt throw do_taylor(); // caught by function::series() } +static ex Li2_conjugate(const ex & x) +{ + // conjugate(Li2(x))==Li2(conjugate(x)) unless on the branch cuts which + // run along the positive real axis beginning at 1. + if (x.info(info_flags::negative)) { + return Li2(x).hold(); + } + if (is_exactly_a(x) && + (!x.imag_part().is_zero() || x < *_num1_p)) { + return Li2(x.conjugate()); + } + return conjugate_function(Li2(x)).hold(); +} + REGISTER_FUNCTION(Li2, eval_func(Li2_eval). evalf_func(Li2_evalf). derivative_func(Li2_deriv). series_func(Li2_series). + conjugate_func(Li2_conjugate). latex_name("\\mathrm{Li}_2")); ////////// @@ -672,7 +688,7 @@ static ex zetaderiv_eval(const ex & n, const ex & x) if (n.info(info_flags::numeric)) { // zetaderiv(0,x) -> zeta(x) if (n.is_zero()) - return zeta(x); + return zeta(x).hold(); } return zetaderiv(n, x).hold(); @@ -992,9 +1008,24 @@ fsolve(const ex& f_in, const symbol& x, const numeric& x1, const numeric& x2) do { xxprev = xx[side]; fxprev = fx[side]; - xx[side] += ex_to(ff.subs(x==xx[side]).evalf()); - fx[side] = ex_to(f.subs(x==xx[side]).evalf()); - if ((side==0 && xx[0]xxprev) || xx[0]>xx[1]) { + ex dx_ = ff.subs(x == xx[side]).evalf(); + if (!is_a(dx_)) + throw std::runtime_error("fsolve(): function derivative does not evaluate numerically"); + xx[side] += ex_to(dx_); + // Now check if Newton-Raphson method shot out of the interval + bool bad_shot = (side == 0 && xx[0] < xxprev) || + (side == 1 && xx[1] > xxprev) || xx[0] > xx[1]; + if (!bad_shot) { + // Compute f(x) only if new x is inside the interval. + // The function might be difficult to compute numerically + // or even ill defined outside the interval. Also it's + // a small optimization. + ex f_x = f.subs(x == xx[side]).evalf(); + if (!is_a(f_x)) + throw std::runtime_error("fsolve(): function does not evaluate numerically"); + fx[side] = ex_to(f_x); + } + if (bad_shot) { // Oops, Newton-Raphson method shot out of the interval. // Restore, and try again with the other side instead! xx[side] = xxprev; @@ -1002,8 +1033,16 @@ fsolve(const ex& f_in, const symbol& x, const numeric& x1, const numeric& x2) side = !side; xxprev = xx[side]; fxprev = fx[side]; - xx[side] += ex_to(ff.subs(x==xx[side]).evalf()); - fx[side] = ex_to(f.subs(x==xx[side]).evalf()); + + ex dx_ = ff.subs(x == xx[side]).evalf(); + if (!is_a(dx_)) + throw std::runtime_error("fsolve(): function derivative does not evaluate numerically [2]"); + xx[side] += ex_to(dx_); + + ex f_x = f.subs(x==xx[side]).evalf(); + if (!is_a(f_x)) + throw std::runtime_error("fsolve(): function does not evaluate numerically [2]"); + fx[side] = ex_to(f_x); } if ((fx[side]<0 && fx[!side]<0) || (fx[side]>0 && fx[!side]>0)) { // Oops, the root isn't bracketed any more. @@ -1027,7 +1066,10 @@ fsolve(const ex& f_in, const symbol& x, const numeric& x1, const numeric& x2) static const double secant_weight = 0.984375; // == 63/64 < 1 numeric xxmid = (1-secant_weight)*0.5*(xx[0]+xx[1]) + secant_weight*(xx[0]+fx[0]*(xx[0]-xx[1])/(fx[1]-fx[0])); - numeric fxmid = ex_to(f.subs(x==xxmid).evalf()); + ex fxmid_ = f.subs(x == xxmid).evalf(); + if (!is_a(fxmid_)) + throw std::runtime_error("fsolve(): function does not evaluate numerically [3]"); + numeric fxmid = ex_to(fxmid_); if (fxmid.is_zero()) { // Luck strikes... return xxmid;