trancendental function branch-point sanity

diff --git a/check/exam_pseries.cpp b/check/exam_pseries.cpp
index 54f4a80620f51f195e18099175358e6ee2ce0920..cb65a53361de1ce5fa8edfda060edfab572c1412 100644 (file)
--- a/check/exam_pseries.cpp
+++ b/check/exam_pseries.cpp
@@ -220,6 +220,73 @@ static unsigned exam_series9(void)
        return check_series(e,2,d,5);
 }
 
        return check_series(e,2,d,5);
 }
 

+// Series expansion of logarithms around branch points
+static unsigned exam_series10(void)
+{
+       unsigned result = 0;
+       ex e, d;
+       symbol a("a");
+       
+       e = log(x);
+       d = log(x);
+       result += check_series(e,0,d,5);
+       
+       e = log(3/x);
+       d = log(3)-log(x);
+       result += check_series(e,0,d,5);
+       
+       e = log(3*pow(x,2));
+       d = log(3)+2*log(x);
+       result += check_series(e,0,d,5);
+       
+       // These ones must not be expanded because it would result in a branch cut
+       // running in the wrong direction. (Other systems tend to get this wrong.)
+       e = log(-x);
+       d = e;
+       result += check_series(e,0,d,5);
+       
+       e = log(I*(x-123));
+       d = e;
+       result += check_series(e,123,d,5);
+       
+       e = log(a*x);
+       d = e;  // we don't know anything about a!
+       result += check_series(e,0,d,5);
+       
+       e = log((1-x)/x);
+       d = log(1-x) - (x-1) + pow(x-1,2)/2 - pow(x-1,3)/3 + Order(pow(x-1,4));
+       result += check_series(e,1,d,4);
+       
+       return result;
+}
+
+// Series expansion of other functions around branch points
+static unsigned exam_series11(void)
+{
+       unsigned result = 0;
+       ex e, d;
+       
+       // NB: Mma and Maple give different results, but they agree if one
+       // takes into account that by assumption |x|<1.
+       e = atan(x);
+       d = (I*log(2)/2-I*log(1+I*x)/2) + (x-I)/4 + I*pow(x-I,2)/16 + Order(pow(x-I,3));
+       result += check_series(e,I,d,3);
+       
+       // NB: here, at -I, Mathematica disagrees, but it is wrong -- they
+       // pick up a complex phase by incorrectly expanding logarithms.
+       e = atan(x);
+       d = (-I*log(2)/2+I*log(1-I*x)/2) + (x+I)/4 - I*pow(x+I,2)/16 + Order(pow(x+I,3));
+       result += check_series(e,-I,d,3);
+       
+       // This is basically the same as above, the branch point is at +/-1:
+       e = atanh(x);
+       d = (-log(2)/2+log(x+1)/2) + (x+1)/4 + pow(x+1,2)/16 + Order(pow(x+1,3));
+       result += check_series(e,-1,d,3);
+       
+       return result;
+}
+
+

 unsigned exam_pseries(void)
 {
        unsigned result = 0;
 unsigned exam_pseries(void)
 {
        unsigned result = 0;


@@ -236,6 +303,8 @@ unsigned exam_pseries(void)
        result += exam_series7();  cout << '.' << flush;
        result += exam_series8();  cout << '.' << flush;
        result += exam_series9();  cout << '.' << flush;
        result += exam_series7();  cout << '.' << flush;
        result += exam_series8();  cout << '.' << flush;
        result += exam_series9();  cout << '.' << flush;

+       result += exam_series10();  cout << '.' << flush;
+       result += exam_series11();  cout << '.' << flush;

        
        if (!result) {
                cout << " passed " << endl;
        
        if (!result) {
                cout << " passed " << endl;

diff --git a/ginac/inifcns.cpp b/ginac/inifcns.cpp
index f410e7aa6202213d585a55de36f522e87d4fc9a7..d5c90e1d59c7c0382bec755be540f38c8c305224 100644 (file)
--- a/ginac/inifcns.cpp
+++ b/ginac/inifcns.cpp
@@ -114,7 +114,8 @@ static ex csgn_series(const ex & arg,
 {
        const ex arg_pt = arg.subs(rel);
        if (arg_pt.info(info_flags::numeric)
 {
        const ex arg_pt = arg.subs(rel);
        if (arg_pt.info(info_flags::numeric)

-        && ex_to_numeric(arg_pt).real().is_zero())
+           && ex_to_numeric(arg_pt).real().is_zero()
+           && !(options & series_options::suppress_branchcut))

                throw (std::domain_error("csgn_series(): on imaginary axis"));
        
        epvector seq;
                throw (std::domain_error("csgn_series(): on imaginary axis"));
        
        epvector seq;

diff --git a/ginac/inifcns_trans.cpp b/ginac/inifcns_trans.cpp
index c7dbfbe6b0508240d0666b1770c18f782568ab0c..cceaad99937e84eb2b72741b4040cfa451eb8255 100644 (file)
--- a/ginac/inifcns_trans.cpp
+++ b/ginac/inifcns_trans.cpp
@@ -158,7 +158,7 @@ static ex log_series(const ex &arg,
        } catch (pole_error) {
                must_expand_arg = true;
        }
        } catch (pole_error) {
                must_expand_arg = true;
        }

-       // or we are at the branch cut anyways
+       // or we are at the branch point anyways

        if (arg_pt.is_zero())
                must_expand_arg = true;
        
        if (arg_pt.is_zero())
                must_expand_arg = true;
        


@@ -176,23 +176,31 @@ static ex log_series(const ex &arg,
                const ex point = rel.rhs();
                const int n = argser.ldegree(*s);
                epvector seq;
                const ex point = rel.rhs();
                const int n = argser.ldegree(*s);
                epvector seq;

-               seq.push_back(expair(n*log(*s-point), _ex0()));
+               // construct what we carelessly called the n*log(x) term above
+               ex coeff = argser.coeff(*s, n);
+               // expand the log, but only if coeff is real and > 0, since otherwise
+               // it would make the branch cut run into the wrong direction
+               if (coeff.info(info_flags::positive))
+                       seq.push_back(expair(n*log(*s-point)+log(coeff), _ex0()));
+               else
+                       seq.push_back(expair(log(coeff*pow(*s-point, n)), _ex0()));

                if (!argser.is_terminating() || argser.nops()!=1) {
                        // in this case n more terms are needed
                if (!argser.is_terminating() || argser.nops()!=1) {
                        // in this case n more terms are needed

-                       ex newarg = ex_to_pseries(arg.series(rel, order+n, options)).shift_exponents(-n).convert_to_poly(true);
+                       // (sadly, to generate them, we have to start from the beginning)
+                       ex newarg = ex_to_pseries((arg/coeff).series(rel, order+n, options)).shift_exponents(-n).convert_to_poly(true);

                        return pseries(rel, seq).add_series(ex_to_pseries(log(newarg).series(rel, order, options)));
                } else  // it was a monomial
                        return pseries(rel, seq);
        }
        if (!(options & series_options::suppress_branchcut) &&
                        return pseries(rel, seq).add_series(ex_to_pseries(log(newarg).series(rel, order, options)));
                } else  // it was a monomial
                        return pseries(rel, seq);
        }
        if (!(options & series_options::suppress_branchcut) &&

-                arg_pt.info(info_flags::negative)) {
+            arg_pt.info(info_flags::negative)) {

                // method:
                // This is the branch cut: assemble the primitive series manually and
                // then add the corresponding complex step function.
                const symbol *s = static_cast<symbol *>(rel.lhs().bp);
                const ex point = rel.rhs();
                const symbol foo;
                // method:
                // This is the branch cut: assemble the primitive series manually and
                // then add the corresponding complex step function.
                const symbol *s = static_cast<symbol *>(rel.lhs().bp);
                const ex point = rel.rhs();
                const symbol foo;

-               ex replarg = series(log(arg), *s==foo, order, false).subs(foo==point);
+               ex replarg = series(log(arg), *s==foo, order).subs(foo==point);

                epvector seq;
                seq.push_back(expair(-I*csgn(arg*I)*Pi, _ex0()));
                seq.push_back(expair(Order(_ex1()), order));
                epvector seq;
                seq.push_back(expair(-I*csgn(arg*I)*Pi, _ex0()));
                seq.push_back(expair(Order(_ex1()), order));


@@ -608,7 +616,7 @@ static ex atan_deriv(const ex & x, unsigned deriv_param)
        return power(_ex1()+power(x,_ex2()), _ex_1());
 }
 
        return power(_ex1()+power(x,_ex2()), _ex_1());
 }
 

-static ex atan_series(const ex &x,
+static ex atan_series(const ex &arg,

                       const relational &rel,
                       int order,
                       unsigned options)
                       const relational &rel,
                       int order,
                       unsigned options)


@@ -620,16 +628,35 @@ static ex atan_series(const ex &x,
        // one running from -I down the imaginary axis.  The points I and -I are
        // poles.
        // On the branch cuts and the poles series expand
        // one running from -I down the imaginary axis.  The points I and -I are
        // poles.
        // On the branch cuts and the poles series expand

-       //     log((1+I*x)/(1-I*x))/(2*I)
+       //     (log(1+I*x)-log(1-I*x))/(2*I)

        // instead.
        // instead.

-       // (The constant term on the cut itself could be made simpler.)
-       const ex x_pt = x.subs(rel);
-       if (!(I*x_pt).info(info_flags::real))
+       const ex arg_pt = arg.subs(rel);
+       if (!(I*arg_pt).info(info_flags::real))

                throw do_taylor();     // Re(x) != 0
                throw do_taylor();     // Re(x) != 0

-       if ((I*x_pt).info(info_flags::real) && abs(I*x_pt)<_ex1())
+       if ((I*arg_pt).info(info_flags::real) && abs(I*arg_pt)<_ex1())

                throw do_taylor();     // Re(x) == 0, but abs(x)<1
                throw do_taylor();     // Re(x) == 0, but abs(x)<1

-       // if we got here we have to care for cuts and poles
-       return (log((1+I*x)/(1-I*x))/(2*I)).series(rel, order, options);
+       // care for the poles, using the defining formula for atan()...
+       if (arg_pt.is_equal(I) || arg_pt.is_equal(-I))
+               return ((log(1+I*arg)-log(1-I*arg))/(2*I)).series(rel, order, options);
+       if (!(options & series_options::suppress_branchcut)) {
+               // method:
+               // This is the branch cut: assemble the primitive series manually and
+               // then add the corresponding complex step function.
+               const symbol *s = static_cast<symbol *>(rel.lhs().bp);
+               const ex point = rel.rhs();
+               const symbol foo;
+               ex replarg = series(atan(arg), *s==foo, order).subs(foo==point);
+               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();
+               else
+                       Order0correction += log((I*arg_pt+_ex1())/(I*arg_pt+_ex_1()))*I*_ex1_2();
+               epvector seq;
+               seq.push_back(expair(Order0correction, _ex0()));
+               seq.push_back(expair(Order(_ex1()), order));
+               return series(replarg - pseries(rel, seq), rel, order);
+       }
+       throw do_taylor();

 }
 
 REGISTER_FUNCTION(atan, eval_func(atan_eval).
 }
 
 REGISTER_FUNCTION(atan, eval_func(atan_eval).


@@ -979,27 +1006,47 @@ static ex atanh_deriv(const ex & x, unsigned deriv_param)
        return power(_ex1()-power(x,_ex2()),_ex_1());
 }
 
        return power(_ex1()-power(x,_ex2()),_ex_1());
 }
 

-static ex atanh_series(const ex &x,
+static ex atanh_series(const ex &arg,

                        const relational &rel,
                        int order,
                        unsigned options)
 {
        GINAC_ASSERT(is_ex_exactly_of_type(rel.lhs(),symbol));
        // method:
                        const relational &rel,
                        int order,
                        unsigned options)
 {
        GINAC_ASSERT(is_ex_exactly_of_type(rel.lhs(),symbol));
        // method:

-       // Taylor series where there is no pole or cut falls back to atan_deriv.
+       // Taylor series where there is no pole or cut falls back to atanh_deriv.

        // There are two branch cuts, one runnig from 1 up the real axis and one
        // one running from -1 down the real axis.  The points 1 and -1 are poles
        // On the branch cuts and the poles series expand
        // There are two branch cuts, one runnig from 1 up the real axis and one
        // one running from -1 down the real axis.  The points 1 and -1 are poles
        // On the branch cuts and the poles series expand

-       //     log((1+x)/(1-x))/(2*I)
+       //     (log(1+x)-log(1-x))/2

        // instead.
        // instead.

-       // (The constant term on the cut itself could be made simpler.)
-       const ex x_pt = x.subs(rel);
-       if (!(x_pt).info(info_flags::real))
+       const ex arg_pt = arg.subs(rel);
+       if (!(arg_pt).info(info_flags::real))

                throw do_taylor();     // Im(x) != 0
                throw do_taylor();     // Im(x) != 0

-       if ((x_pt).info(info_flags::real) && abs(x_pt)<_ex1())
+       if ((arg_pt).info(info_flags::real) && abs(arg_pt)<_ex1())

                throw do_taylor();     // Im(x) == 0, but abs(x)<1
                throw do_taylor();     // Im(x) == 0, but abs(x)<1

-       // if we got here we have to care for cuts and poles
-       return (log((1+x)/(1-x))/2).series(rel, order, options);
+       // care for the poles, using the defining formula for atanh()...
+       if (arg_pt.is_equal(_ex1()) || arg_pt.is_equal(_ex_1()))
+               return ((log(_ex1()+arg)-log(_ex1()-arg))*_ex1_2()).series(rel, order, options);
+       // ...and the branch cuts (the discontinuity at the cut being just I*Pi)
+       if (!(options & series_options::suppress_branchcut)) {
+               // method:
+               // This is the branch cut: assemble the primitive series manually and
+               // then add the corresponding complex step function.
+               const symbol *s = static_cast<symbol *>(rel.lhs().bp);
+               const ex point = rel.rhs();
+               const symbol foo;
+               ex replarg = series(atanh(arg), *s==foo, order).subs(foo==point);
+               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();
+               else
+                       Order0correction += log((arg_pt+_ex1())/(arg_pt+_ex_1()))*_ex_1_2();
+               epvector seq;
+               seq.push_back(expair(Order0correction, _ex0()));
+               seq.push_back(expair(Order(_ex1()), order));
+               return series(replarg - pseries(rel, seq), rel, order);
+       }
+       throw do_taylor();

 }
 
 REGISTER_FUNCTION(atanh, eval_func(atanh_eval).
 }
 
 REGISTER_FUNCTION(atanh, eval_func(atanh_eval).