X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Finifcns_trans.cpp;h=872308b93f54e234a99ebf493a4416b4d6b7f1a8;hp=67a578dd72b4f1963b752cf2359ba8af9b5ce839;hb=4405b29465293f3b6ab37745ff601f519b0256e2;hpb=dbd9c306a74f1cb258c0d15a346b973b39deaad2 diff --git a/ginac/inifcns_trans.cpp b/ginac/inifcns_trans.cpp index 67a578dd..872308b9 100644 --- a/ginac/inifcns_trans.cpp +++ b/ginac/inifcns_trans.cpp @@ -4,7 +4,7 @@ * functions. */ /* - * GiNaC Copyright (C) 1999-2003 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,6 +29,7 @@ #include "constant.h" #include "numeric.h" #include "power.h" +#include "operators.h" #include "relational.h" #include "symbol.h" #include "pseries.h" @@ -54,6 +55,7 @@ static ex exp_eval(const ex & x) if (x.is_zero()) { return _ex1; } + // exp(n*Pi*I/2) -> {+1|+I|-1|-I} const ex TwoExOverPiI=(_ex2*x)/(Pi*I); if (TwoExOverPiI.info(info_flags::integer)) { @@ -67,11 +69,12 @@ static ex exp_eval(const ex & x) if (z.is_equal(_num3)) return ex(-I); } + // exp(log(x)) -> x if (is_ex_the_function(x, log)) return x.op(0); - // exp(float) + // exp(float) -> float if (x.info(info_flags::numeric) && !x.info(info_flags::crational)) return exp(ex_to(x)); @@ -116,13 +119,18 @@ static ex log_eval(const ex & x) return (Pi*I*_num1_2); if (x.is_equal(-I)) // log(-I) -> -Pi*I/2 return (Pi*I*_num_1_2); - // log(float) + + // log(float) -> float if (!x.info(info_flags::crational)) return log(ex_to(x)); } + // log(exp(t)) -> t (if -Pi < t.imag() <= Pi): if (is_ex_the_function(x, exp)) { const ex &t = x.op(0); + if (is_a(t) && t.info(info_flags::real)) { + return t; + } if (t.info(info_flags::numeric)) { const numeric &nt = ex_to(t); if (nt.is_real()) @@ -141,14 +149,6 @@ static ex log_deriv(const ex & x, unsigned deriv_param) return power(x, _ex_1); } -// This is a strange workaround for a compiliation problem with the try statement -// below. With -O1 the exception is not caucht properly as of GCC-2.95.2, at -// least on i386. Version 2.95.4 seems to have fixed this silly problem, though. -// Funnily, with a simple extern declaration here it mysteriously works again. -#if defined(__GNUC__) && (__GNUC__==2) -extern "C" int putchar(int); -#endif - static ex log_series(const ex &arg, const relational &rel, int order, @@ -159,7 +159,7 @@ static ex log_series(const ex &arg, bool must_expand_arg = false; // maybe substitution of rel into arg fails because of a pole try { - arg_pt = arg.subs(rel); + arg_pt = arg.subs(rel, subs_options::no_pattern); } catch (pole_error) { must_expand_arg = true; } @@ -201,6 +201,22 @@ static ex log_series(const ex &arg, // in this case n more (or less) terms are needed // (sadly, to generate them, we have to start from the beginning) const ex newarg = ex_to((arg/coeff).series(rel, order+n, options)).shift_exponents(-n).convert_to_poly(true); + if (n == 0 && coeff == 1) { + epvector epv; + ex acc = (new pseries(rel, epv))->setflag(status_flags::dynallocated); + epv.reserve(2); + epv.push_back(expair(-1, _ex0)); + epv.push_back(expair(Order(_ex1), order)); + ex rest = pseries(rel, epv).add_series(argser); + for (int i = order-1; i>0; --i) { + epvector cterm; + cterm.reserve(1); + cterm.push_back(expair(i%2 ? _ex1/i : _ex_1/i, _ex0)); + acc = pseries(rel, cterm).add_series(ex_to(acc)); + acc = (ex_to(rest)).mul_series(ex_to(acc)); + } + return acc; + } return pseries(rel, seq).add_series(ex_to(log(newarg).series(rel, order, options))); } else // it was a monomial return pseries(rel, seq); @@ -213,7 +229,7 @@ static ex log_series(const ex &arg, const symbol &s = ex_to(rel.lhs()); const ex &point = rel.rhs(); const symbol foo; - const ex replarg = series(log(arg), s==foo, order).subs(foo==point); + const ex replarg = series(log(arg), s==foo, order).subs(foo==point, subs_options::no_pattern); epvector seq; seq.push_back(expair(-I*csgn(arg*I)*Pi, _ex0)); seq.push_back(expair(Order(_ex1), order)); @@ -275,15 +291,18 @@ static ex sin_eval(const ex & x) if (z.is_equal(_num30)) // sin(Pi/2) -> 1 return sign; } - + if (is_exactly_a(x)) { const ex &t = x.op(0); + // sin(asin(x)) -> x if (is_ex_the_function(x, asin)) return t; + // sin(acos(x)) -> sqrt(1-x^2) if (is_ex_the_function(x, acos)) return sqrt(_ex1-power(t,_ex2)); + // sin(atan(x)) -> x/sqrt(1+x^2) if (is_ex_the_function(x, atan)) return t*power(_ex1+power(t,_ex2),_ex_1_2); @@ -292,6 +311,10 @@ static ex sin_eval(const ex & x) // sin(float) -> float if (x.info(info_flags::numeric) && !x.info(info_flags::crational)) return sin(ex_to(x)); + + // sin() is odd + if (x.info(info_flags::negative)) + return -sin(-x); return sin(x).hold(); } @@ -356,15 +379,18 @@ static ex cos_eval(const ex & x) if (z.is_equal(_num30)) // cos(Pi/2) -> 0 return _ex0; } - + if (is_exactly_a(x)) { const ex &t = x.op(0); + // cos(acos(x)) -> x if (is_ex_the_function(x, acos)) return t; + // cos(asin(x)) -> sqrt(1-x^2) if (is_ex_the_function(x, asin)) return sqrt(_ex1-power(t,_ex2)); + // cos(atan(x)) -> 1/sqrt(1+x^2) if (is_ex_the_function(x, atan)) return power(_ex1+power(t,_ex2),_ex_1_2); @@ -374,6 +400,10 @@ static ex cos_eval(const ex & x) if (x.info(info_flags::numeric) && !x.info(info_flags::crational)) return cos(ex_to(x)); + // cos() is even + if (x.info(info_flags::negative)) + return cos(-x); + return cos(x).hold(); } @@ -433,15 +463,18 @@ static ex tan_eval(const ex & x) if (z.is_equal(_num30)) // tan(Pi/2) -> infinity throw (pole_error("tan_eval(): simple pole",1)); } - + if (is_exactly_a(x)) { const ex &t = x.op(0); + // tan(atan(x)) -> x if (is_ex_the_function(x, atan)) return t; + // tan(asin(x)) -> x/sqrt(1+x^2) if (is_ex_the_function(x, asin)) return t*power(_ex1-power(t,_ex2),_ex_1_2); + // tan(acos(x)) -> sqrt(1-x^2)/x if (is_ex_the_function(x, acos)) return power(t,_ex_1)*sqrt(_ex1-power(t,_ex2)); @@ -452,6 +485,10 @@ static ex tan_eval(const ex & x) return tan(ex_to(x)); } + // tan() is odd + if (x.info(info_flags::negative)) + return -tan(-x); + return tan(x).hold(); } @@ -472,11 +509,11 @@ static ex tan_series(const ex &x, // method: // Taylor series where there is no pole falls back to tan_deriv. // On a pole simply expand sin(x)/cos(x). - const ex x_pt = x.subs(rel); + const ex x_pt = x.subs(rel, subs_options::no_pattern); if (!(2*x_pt/Pi).info(info_flags::odd)) throw do_taylor(); // caught by function::series() // if we got here we have to care for a simple pole - return (sin(x)/cos(x)).series(rel, order+2, options); + return (sin(x)/cos(x)).series(rel, order, options); } REGISTER_FUNCTION(tan, eval_func(tan_eval). @@ -500,24 +537,34 @@ static ex asin_evalf(const ex & x) static ex asin_eval(const ex & x) { if (x.info(info_flags::numeric)) { + // asin(0) -> 0 if (x.is_zero()) return x; + // asin(1/2) -> Pi/6 if (x.is_equal(_ex1_2)) return numeric(1,6)*Pi; + // asin(1) -> Pi/2 if (x.is_equal(_ex1)) return _num1_2*Pi; + // asin(-1/2) -> -Pi/6 if (x.is_equal(_ex_1_2)) return numeric(-1,6)*Pi; + // asin(-1) -> -Pi/2 if (x.is_equal(_ex_1)) return _num_1_2*Pi; + // asin(float) -> float if (!x.info(info_flags::crational)) return asin(ex_to(x)); + + // asin() is odd + if (x.info(info_flags::negative)) + return -asin(-x); } return asin(x).hold(); @@ -551,24 +598,34 @@ static ex acos_evalf(const ex & x) static ex acos_eval(const ex & x) { if (x.info(info_flags::numeric)) { + // acos(1) -> 0 if (x.is_equal(_ex1)) return _ex0; + // acos(1/2) -> Pi/3 if (x.is_equal(_ex1_2)) return _ex1_3*Pi; + // acos(0) -> Pi/2 if (x.is_zero()) return _ex1_2*Pi; + // acos(-1/2) -> 2/3*Pi if (x.is_equal(_ex_1_2)) return numeric(2,3)*Pi; + // acos(-1) -> Pi if (x.is_equal(_ex_1)) return Pi; + // acos(float) -> float if (!x.info(info_flags::crational)) return acos(ex_to(x)); + + // acos(-x) -> Pi-acos(x) + if (x.info(info_flags::negative)) + return Pi-acos(-x); } return acos(x).hold(); @@ -602,20 +659,29 @@ static ex atan_evalf(const ex & x) static ex atan_eval(const ex & x) { if (x.info(info_flags::numeric)) { + // atan(0) -> 0 if (x.is_zero()) return _ex0; + // atan(1) -> Pi/4 if (x.is_equal(_ex1)) return _ex1_4*Pi; + // atan(-1) -> -Pi/4 if (x.is_equal(_ex_1)) return _ex_1_4*Pi; + if (x.is_equal(I) || x.is_equal(-I)) throw (pole_error("atan_eval(): logarithmic pole",0)); + // atan(float) -> float if (!x.info(info_flags::crational)) return atan(ex_to(x)); + + // atan() is odd + if (x.info(info_flags::negative)) + return -atan(-x); } return atan(x).hold(); @@ -643,7 +709,7 @@ static ex atan_series(const ex &arg, // On the branch cuts and the poles series expand // (log(1+I*x)-log(1-I*x))/(2*I) // instead. - const ex arg_pt = arg.subs(rel); + const ex arg_pt = arg.subs(rel, subs_options::no_pattern); if (!(I*arg_pt).info(info_flags::real)) throw do_taylor(); // Re(x) != 0 if ((I*arg_pt).info(info_flags::real) && abs(I*arg_pt)<_ex1) @@ -658,7 +724,7 @@ static ex atan_series(const ex &arg, const symbol &s = ex_to(rel.lhs()); const ex &point = rel.rhs(); const symbol foo; - const ex replarg = series(atan(arg), s==foo, order).subs(foo==point); + const ex replarg = series(atan(arg), s==foo, order).subs(foo==point, subs_options::no_pattern); ex Order0correction = replarg.op(0)+csgn(arg)*Pi*_ex_1_2; if ((I*arg_pt)<_ex0) Order0correction += log((I*arg_pt+_ex_1)/(I*arg_pt+_ex1))*I*_ex_1_2; @@ -685,19 +751,79 @@ REGISTER_FUNCTION(atan, eval_func(atan_eval). static ex atan2_evalf(const ex &y, const ex &x) { if (is_exactly_a(y) && is_exactly_a(x)) - return atan2(ex_to(y), ex_to(x)); + return atan(ex_to(y), ex_to(x)); return atan2(y, x).hold(); } static ex atan2_eval(const ex & y, const ex & x) { - if (y.info(info_flags::numeric) && !y.info(info_flags::crational) && - x.info(info_flags::numeric) && !x.info(info_flags::crational)) { - return atan2_evalf(y,x); + if (y.info(info_flags::numeric) && x.info(info_flags::numeric)) { + + if (y.is_zero()) { + + // atan(0, 0) -> 0 + if (x.is_zero()) + return _ex0; + + // atan(0, x), x real and positive -> 0 + if (x.info(info_flags::positive)) + return _ex0; + + // atan(0, x), x real and negative -> -Pi + if (x.info(info_flags::negative)) + return _ex_1*Pi; + } + + if (x.is_zero()) { + + // atan(y, 0), y real and positive -> Pi/2 + if (y.info(info_flags::positive)) + return _ex1_2*Pi; + + // atan(y, 0), y real and negative -> -Pi/2 + if (y.info(info_flags::negative)) + return _ex_1_2*Pi; + } + + if (y.is_equal(x)) { + + // atan(y, y), y real and positive -> Pi/4 + if (y.info(info_flags::positive)) + return _ex1_4*Pi; + + // atan(y, y), y real and negative -> -3/4*Pi + if (y.info(info_flags::negative)) + return numeric(-3, 4)*Pi; + } + + if (y.is_equal(-x)) { + + // atan(y, -y), y real and positive -> 3*Pi/4 + if (y.info(info_flags::positive)) + return numeric(3, 4)*Pi; + + // atan(y, -y), y real and negative -> -Pi/4 + if (y.info(info_flags::negative)) + return _ex_1_4*Pi; + } + + // atan(float, float) -> float + if (!y.info(info_flags::crational) && !x.info(info_flags::crational)) + return atan(ex_to(y), ex_to(x)); + + // atan(real, real) -> atan(y/x) +/- Pi + if (y.info(info_flags::real) && x.info(info_flags::real)) { + if (x.info(info_flags::positive)) + return atan(y/x); + else if(y.info(info_flags::positive)) + return atan(y/x)+Pi; + else + return atan(y/x)-Pi; + } } - - return atan2(y,x).hold(); + + return atan2(y, x).hold(); } static ex atan2_deriv(const ex & y, const ex & x, unsigned deriv_param) @@ -731,10 +857,18 @@ static ex sinh_evalf(const ex & x) static ex sinh_eval(const ex & x) { if (x.info(info_flags::numeric)) { - if (x.is_zero()) // sinh(0) -> 0 + + // sinh(0) -> 0 + if (x.is_zero()) return _ex0; - if (!x.info(info_flags::crational)) // sinh(float) -> float + + // sinh(float) -> float + if (!x.info(info_flags::crational)) return sinh(ex_to(x)); + + // sinh() is odd + if (x.info(info_flags::negative)) + return -sinh(-x); } if ((x/Pi).info(info_flags::numeric) && @@ -743,12 +877,15 @@ static ex sinh_eval(const ex & x) if (is_exactly_a(x)) { const ex &t = x.op(0); + // sinh(asinh(x)) -> x if (is_ex_the_function(x, asinh)) return t; + // sinh(acosh(x)) -> sqrt(x-1) * sqrt(x+1) if (is_ex_the_function(x, acosh)) return sqrt(t-_ex1)*sqrt(t+_ex1); + // sinh(atanh(x)) -> x/sqrt(1-x^2) if (is_ex_the_function(x, atanh)) return t*power(_ex1-power(t,_ex2),_ex_1_2); @@ -785,10 +922,18 @@ static ex cosh_evalf(const ex & x) static ex cosh_eval(const ex & x) { if (x.info(info_flags::numeric)) { - if (x.is_zero()) // cosh(0) -> 1 + + // cosh(0) -> 1 + if (x.is_zero()) return _ex1; - if (!x.info(info_flags::crational)) // cosh(float) -> float + + // cosh(float) -> float + if (!x.info(info_flags::crational)) return cosh(ex_to(x)); + + // cosh() is even + if (x.info(info_flags::negative)) + return cosh(-x); } if ((x/Pi).info(info_flags::numeric) && @@ -797,12 +942,15 @@ static ex cosh_eval(const ex & x) if (is_exactly_a(x)) { const ex &t = x.op(0); + // cosh(acosh(x)) -> x if (is_ex_the_function(x, acosh)) return t; + // cosh(asinh(x)) -> sqrt(1+x^2) if (is_ex_the_function(x, asinh)) return sqrt(_ex1+power(t,_ex2)); + // cosh(atanh(x)) -> 1/sqrt(1-x^2) if (is_ex_the_function(x, atanh)) return power(_ex1-power(t,_ex2),_ex_1_2); @@ -839,10 +987,18 @@ static ex tanh_evalf(const ex & x) static ex tanh_eval(const ex & x) { if (x.info(info_flags::numeric)) { - if (x.is_zero()) // tanh(0) -> 0 + + // tanh(0) -> 0 + if (x.is_zero()) return _ex0; - if (!x.info(info_flags::crational)) // tanh(float) -> float + + // tanh(float) -> float + if (!x.info(info_flags::crational)) return tanh(ex_to(x)); + + // tanh() is odd + if (x.info(info_flags::negative)) + return -tanh(-x); } if ((x/Pi).info(info_flags::numeric) && @@ -851,12 +1007,15 @@ static ex tanh_eval(const ex & x) if (is_exactly_a(x)) { const ex &t = x.op(0); + // tanh(atanh(x)) -> x if (is_ex_the_function(x, atanh)) return t; + // tanh(asinh(x)) -> x/sqrt(1+x^2) if (is_ex_the_function(x, asinh)) return t*power(_ex1+power(t,_ex2),_ex_1_2); + // tanh(acosh(x)) -> sqrt(x-1)*sqrt(x+1)/x if (is_ex_the_function(x, acosh)) return sqrt(t-_ex1)*sqrt(t+_ex1)*power(t,_ex_1); @@ -882,11 +1041,11 @@ static ex tanh_series(const ex &x, // method: // Taylor series where there is no pole falls back to tanh_deriv. // On a pole simply expand sinh(x)/cosh(x). - const ex x_pt = x.subs(rel); + const ex x_pt = x.subs(rel, subs_options::no_pattern); if (!(2*I*x_pt/Pi).info(info_flags::odd)) throw do_taylor(); // caught by function::series() // if we got here we have to care for a simple pole - return (sinh(x)/cosh(x)).series(rel, order+2, options); + return (sinh(x)/cosh(x)).series(rel, order, options); } REGISTER_FUNCTION(tanh, eval_func(tanh_eval). @@ -910,12 +1069,18 @@ static ex asinh_evalf(const ex & x) static ex asinh_eval(const ex & x) { if (x.info(info_flags::numeric)) { + // asinh(0) -> 0 if (x.is_zero()) return _ex0; + // asinh(float) -> float if (!x.info(info_flags::crational)) return asinh(ex_to(x)); + + // asinh() is odd + if (x.info(info_flags::negative)) + return -asinh(-x); } return asinh(x).hold(); @@ -948,18 +1113,26 @@ static ex acosh_evalf(const ex & x) static ex acosh_eval(const ex & x) { if (x.info(info_flags::numeric)) { + // acosh(0) -> Pi*I/2 if (x.is_zero()) return Pi*I*numeric(1,2); + // acosh(1) -> 0 if (x.is_equal(_ex1)) return _ex0; + // acosh(-1) -> Pi*I if (x.is_equal(_ex_1)) return Pi*I; + // acosh(float) -> float if (!x.info(info_flags::crational)) return acosh(ex_to(x)); + + // acosh(-x) -> Pi*I-acosh(x) + if (x.info(info_flags::negative)) + return Pi*I-acosh(-x); } return acosh(x).hold(); @@ -992,15 +1165,22 @@ static ex atanh_evalf(const ex & x) static ex atanh_eval(const ex & x) { if (x.info(info_flags::numeric)) { + // atanh(0) -> 0 if (x.is_zero()) return _ex0; + // atanh({+|-}1) -> throw if (x.is_equal(_ex1) || x.is_equal(_ex_1)) throw (pole_error("atanh_eval(): logarithmic pole",0)); + // atanh(float) -> float if (!x.info(info_flags::crational)) return atanh(ex_to(x)); + + // atanh() is odd + if (x.info(info_flags::negative)) + return -atanh(-x); } return atanh(x).hold(); @@ -1027,7 +1207,7 @@ static ex atanh_series(const ex &arg, // On the branch cuts and the poles series expand // (log(1+x)-log(1-x))/2 // instead. - const ex arg_pt = arg.subs(rel); + const ex arg_pt = arg.subs(rel, subs_options::no_pattern); if (!(arg_pt).info(info_flags::real)) throw do_taylor(); // Im(x) != 0 if ((arg_pt).info(info_flags::real) && abs(arg_pt)<_ex1) @@ -1043,7 +1223,7 @@ static ex atanh_series(const ex &arg, const symbol &s = ex_to(rel.lhs()); const ex &point = rel.rhs(); const symbol foo; - const ex replarg = series(atanh(arg), s==foo, order).subs(foo==point); + const ex replarg = series(atanh(arg), s==foo, order).subs(foo==point, subs_options::no_pattern); ex Order0correction = replarg.op(0)+csgn(I*arg)*Pi*I*_ex1_2; if (arg_pt<_ex0) Order0correction += log((arg_pt+_ex_1)/(arg_pt+_ex1))*_ex1_2;