Commit fc2526f4 authored by Yannick Ulrich's avatar Yannick Ulrich

GiNaC seems to have problems with some gpls, breaking the test

parent d6db09df
......@@ -701,7 +701,7 @@ CONTAINS
integer,parameter :: nieps = 30
integer,parameter :: ncmpl = 30
integer,parameter :: perweight(4) = (/ 100, 500, 5000, 15000 /)
integer,parameter :: perweight(4) = (/ 100, 500, 30000, 50000 /)
real(kind=prec), parameter :: rrange = 1.5
type(inum), dimension(nzero+nieps+ncmpl) :: basis
type(inum), dimension(size(perweight)) :: args
......@@ -745,7 +745,14 @@ CONTAINS
ans(1) = cmplx(geval([args(1:w),ione], w+1),kind=prec)
#endif
ans(2) = G(args(1:w),ione)
if ((abs(ans(1)) .gt. 1.e10).or.(abs(ans(2)) .gt. 1.e10)) then
print*," can't deal with",args," (seed was",oldseed,')'
print*,"GiNaC : ",ans(1)
print*,"handyG: ",ans(2)
cycle
endif
write(9) w,i,oldseed, abs(ans(1)-ans(2))
flush(9)
if(abs(ans(1)-ans(2)) > maxd) maxd = abs(ans(1)-ans(2))
if(abs(ans(1)-ans(2)) > tol) goto 123
enddo
......
......@@ -23,10 +23,15 @@ complex_t geval_(inum_t * z, int* n) {
s.append(z->i0);
z++;
}
ex ans = G(w,s,z->c.r).evalf();
return {
.r = ex_to<numeric>(evalf(real_part(ans))).to_double(),
.i = ex_to<numeric>(evalf(imag_part(ans))).to_double()
};
try {
ex ans = G(w,s,z->c.r).evalf();
return {
.r = ex_to<numeric>(evalf(real_part(ans))).to_double(),
.i = ex_to<numeric>(evalf(imag_part(ans))).to_double()
};
} catch (...) {
std::cout << "Caught!!\n";
return {.r = 1.e11, .i = 1.e11};
}
}
  • The problem that GiNaC throws is

    terminate called after throwing an instance of 'std::runtime_error'
      what():  map_trafo_H_1mx: cannot handle weights equal -1!

    at

    else if (parameter.op(0) == -1) {
        throw std::runtime_error("map_trafo_H_1mx: cannot handle weights equal -1!");
    }

    (taken from inifcns_nstdsums.cpp:2892)

    An example for one of these GPLs is

    G(x,y,x;1) \supset {\rm Li}_{1,2}\Big(\frac xy,1\Big) = H_{1,0,1}\Big(\frac xy\Big) = H_{1,0,1}(z)

    with x=\tfrac1{100}+\tfrac{\rm i}8, y=-\tfrac18, and z=-\tfrac2{25}-{\rm i}.

    This gets (internally) transformed to1

    H_{1,0,1}(z) \mapsto H_{-1,0,-1},(-1/z)

    to make sure that that \Re z>0

    ex res = 1; 
    
    // ensure that the realpart of the argument is positive
    if (cln::realpart(x) < 0) {
        x = -x;
        for (std::size_t i = 0; i < m.nops(); i++) {
            if (m.op(i) != 0) {
                m.let_op(i) = -m.op(i);
                res *= -1;
            }
        }
    }

    (taken from inifcns_nstdsums.cpp:3324)

    GiNaC then wants to speed up the evaluation of the HPL with arguments in the region 0.95\le|z|\le2.0 and

    \Big| \frac{1-x}{1+x} \Big| \nless 0.9

    (which is implemented by checking whether |x-9.53| \nleq 9.47) by performing the mapping

    H_{\vec w}(x) \mapsto H_{\vec w'}(1-x)
    // check transformations for 0.95 <= |x| < 2.0
    
    // |(1-x)/(1+x)| < 0.9 -> circular area with center=9.53+0i and radius=9.47
    if (cln::abs(x-9.53) <= 9.47) {
        // x -> (1-x)/(1+x)
        map_trafo_H_1mxt1px trafo;
        res *= trafo(H(m, xtemp).hold());
    } else {
        // x -> 1-x
        if (has_minus_one) {
            map_trafo_H_convert_to_Li filter;
            return filter(H(m, numeric(x)).hold()).evalf();
        }
        map_trafo_H_1mx trafo;
        res *= trafo(H(m, xtemp).hold());
    }

    (taken from inifcns_nstdsums.cpp:3301)

    This is what eventually causes the std::runtime_error. This should not actually happen because has_minus_one is checked.

    However, it is my belief that the transformation \Re z>0 is done without recomputing has_minus_one. This could be remedied by

    diff --git a/ginac/inifcns_nstdsums.cpp b/ginac/inifcns_nstdsums.cpp
    index f040e8ad..c08cd32b 100644
    --- a/ginac/inifcns_nstdsums.cpp
    +++ b/ginac/inifcns_nstdsums.cpp
    @@ -3300,6 +3300,9 @@ static ex H_evalf(const ex& x1, const ex& x2)
     		
     		// check transformations for 0.95 <= |x| < 2.0
     		
    +		for (auto const& i : m) {
    +			if (i == -1) has_minus_one = true;
    +		}
     		// |(1-x)/(1+x)| < 0.9 -> circular area with center=9.53+0i and radius=9.47
     		if (cln::abs(x-9.53) <= 9.47) {
     			// x -> (1-x)/(1+x)

    which is a rather ugly fix (due to my inability to write coherent C++ code)

    1. By which I mean that while evalf(H({1,0,1},z)); crashes, evalf(H({-1,0,-1},1/x)); does not.

    Edited by Yannick Ulrich
  • This was properly fixed in GiNaC 1.7.8

    commit 0eaae44cd9eb9fa987bb9cbd4341b0f4c8d2f495
    Author: Stefan Weinzierl <weinzierl@uni-mainz.de>
    Date:   Mon Oct 7 20:32:01 2019 +0200
    
        Fix bug in H_evalf: Flag has_minus_one is now computed where it is needed.
        
        This bug has been reported and fixed by Yannick Ulrich <yannick.ulrich@psi.ch>.
    
    diff --git a/ginac/inifcns_nstdsums.cpp b/ginac/inifcns_nstdsums.cpp
    index f040e8ad..c20526ba 100644
    --- a/ginac/inifcns_nstdsums.cpp
    +++ b/ginac/inifcns_nstdsums.cpp
    @@ -3223,7 +3223,6 @@ static ex H_evalf(const ex& x1, const ex& x2)
     			return filter(H(x1, xtemp).hold()).subs(xtemp==x2).evalf();
     		}
     		// ... and expand parameter notation
    -		bool has_minus_one = false;
     		lst m;
     		for (const auto & it : morg) {
     			if (it > 1) {
    @@ -3236,7 +3235,6 @@ static ex H_evalf(const ex& x1, const ex& x2)
     					m.append(0);
     				}
     				m.append(-1);
    -				has_minus_one = true;
     			} else {
     				m.append(it);
     			}
    @@ -3297,7 +3295,14 @@ static ex H_evalf(const ex& x1, const ex& x2)
     			}
     			return res.subs(xtemp == numeric(x)).evalf();
     		}
    -		
    +
    +		// check for letters (-1)
    +		bool has_minus_one = false;
    +		for (const auto & it : m) {
    +			if (it == -1)
    +				has_minus_one = true;
    +		}
    +
     		// check transformations for 0.95 <= |x| < 2.0
     		
     		// |(1-x)/(1+x)| < 0.9 -> circular area with center=9.53+0i and radius=9.47
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment