Floating complex numbers
Hi, I'm currently fixing some bugs in inifcns_nstdsums.cpp and came across the following strange behaviour: The GiNaC manual says "Complex numbers are automatically converted to real numbers if the imaginary part becomes zero". However, this does not seem to be the case, as the following lines show: ex w1 = numeric(3)+I*numeric(2); ex w2 = numeric(2)+I*numeric(2); ex w3 = w1-w2; if ( (w3 == 1) && ex_to<numeric>(w3).is_real() ) std::cout << "o.k., equal to one and real" << std::endl; ex z1 = numeric(3)+I*numeric("2.1"); ex z2 = numeric(2)+I*numeric("2.1"); ex z3 = z1-z2; if ( (z3 == 1) && !ex_to<numeric>(z3).is_real() ) std::cout << "bug, equal to one, but not real" << std::endl; The bug seems to be in CLN, as it does not downgrade z3 to an instance of cln::cl_R_ring. Currently this leads to an error in the file inifcns_nstdsums.cpp. In this file a variable z is first checked to be equal to +1 or -1, and later on a comparison z>0 is carried out. The comparison leads to the exception numeric::operator>(): complex inequality, since for the computer z is equal to +/- 1, but not real, similar to the last line in the code fragment above. I can work around this bug in the file inifcns_nstdsums.cpp, but the clean solution would be to fix this in CLN. Any chances for that ? Best wishes, Stefan
Hi Stefan, On 01/08/2014 09:05 PM, Stefan Weinzierl wrote:
Any chances for that ?
I've had a quick look at the problem. The attached patch to CLN seems to fix your problem and not cause any regressions on either CLN or GiNaC. Please feel free to test it. I am, however, not ready to commit it yet. I first need to find an answer to a question I am having with regard to this kind of conversions. I'll keep you posted... Best wishes -richy. -- Richard B. Kreckel <http://www.ginac.de/~kreckel/>
Hi Richy, I would also propose the attached patch for the tutorial to make it clear to all users, what GiNaC and CLN do with floating point complex numbers. The proposed modification describes the philosophy of Bruno and you. Best wishes, Stefan
Hi, On Fri, Jan 17, 2014 at 10:50 AM, Stefan Weinzierl <stefanw@thep.physik.uni-mainz.de> wrote:
I would also propose the attached patch for the tutorial to make it clear to all users, what GiNaC and CLN do with floating point complex numbers.
I think the issue is a bit more generic. That is GiNaC (and CLN) converts exact numbers into inexact ones only upon the user request (i.e. one need to call evalf()). Likewise inexact numbers are never converted into exact ones. This behavior is already documented. Do we really need to document it once more? Best regards, Alexei
Hi, attached is a patch, which fixes several bugs in the file inifcns_nstdsums.cpp. The bugs have been reported to me by Christoph Meyer and Andreas von Manteuffel. In detail, the bugs and its fixes are: - The method Li_projection of the classical polylogs calls the cln-function log(1-x). We better avoid calling this function with x=1. - The method mLi_numeric calls G_numeric, G_numeric expects the signs of the imaginary parts to be set correctly, so we do this. - In the methods G3_evalf and G3_eval the user supplied signs of the imaginary parts are only relevant for x real and positive. A negative sign in the case of x real and negative is irrelevant (there is no branch cut) and led to wrong results. This is now avoided. - The method S_num within the Nielsen polylogs can lead to an infinite recursion if the argument is close to the sixth root of unity or its conjugate. In more detail, the region abs(value)<=1 && abs(value)>0.95 && abs(1-value)<=1 && abs(1-value)>0.95 used to be mapped infinitely times onto itself. This is now avoided. This however reveals the next problem: The numerical convergence in this region is very slow. Within the Nielsen polylogs there is no transformation available to improve the convergence. However we can use the (1-x)/(1+x) transformation within the harmonic polylogs. In order to avoid another infinite recursion I have inserted a few hold()'s in the method H_evalf, otherwise we would fall back immediately again to Nielsen polylogs. The hold()'s should have been there anyway. Not included in this patch are workarounds for bugs caused by first comparing a complex floating number against +/-1, and in the case of equality using the operator < later on. This is logically correct and better fixed in cln. Best wishes, Stefan
On 01/12/2014 10:10 AM, Stefan Weinzierl wrote:
attached is a patch, which fixes several bugs in the file inifcns_nstdsums.cpp.
I have applied your four patches. Thanks, Stefan. -richy. -- Richard B. Kreckel <http://www.ginac.de/~kreckel/>
Hi Stefan, I'm reverting the patch in CLN. Please see Bruno's explanation below. (His email hasn't appeared on the list. I have no idea why.) I realized that the patch in CLN is wrong while I was trying to adjust the documentation. Then Bruno's email came in and settled the case. Regarding the original problem in GiNaC: I'm not sure it should be fixed at all. Considering Bruno's argument, it is neither CLN nor GiNaC that is wrong, but the calling code of Li2. It should be more careful when some of the x are computed as +/-1+0.0i. In any case, Li2 shouldn't convert the argument to +/-1 automatically. Please consider rewriting this calling code. Best wishes -richy. -- Richard B. Kreckel <http://www.ginac.de/~kreckel/>
Hi Richy, fine with me. If you revert the patch in CLN, you can then apply the attached patch to inifcns_nstdsums.cpp, this fixes the issues there. Best wishes, Stefan
Hi Stefan, On 01/16/2014 11:42 PM, Stefan Weinzierl wrote:
fine with me. If you revert the patch in CLN, you can then apply the attached patch to inifcns_nstdsums.cpp, this fixes the issues there.
What precisely *is* the issue? -richy. -- Richard B. Kreckel <http://www.ginac.de/~kreckel/>
Hi Richy, the method Li_eval checks if all entries of x_ are equal to +1 or -1. In this case, the function Li is converted to a multiple zeta value. If all but the first entry of x_ are equal to +1 or -1, the function Li is converted to a harmonic polylog H. After checking that an entry of x_ is equal to +1 or -1, the operator > is used. This may lead to an "numeric::operator>(): complex inequality" exception in the case where the real part of an entry of x_ is equal to +/-1 (either exact or as a float, this does not matter) and the imaginary part is equal to zero as a float. The attached patch avoids the improper use of the operator > in this cases and sets the +/-1 arguments of H and zeta as real numbers. For zeta the method Li_eval is modified, for H the appropriate place is a modification of the method convert_parameter_Li_to_H. Best wishes, Stefan
Hi, Stefan, Thanks for an explanation. I've got a few objections regarding the patch. I think handling the complex numbers with the imaginary part being floating-point zero should be done in zeta2_evalf() rather than to Li_eval(). The point is that (multiple) zeta values can appear not only as a result of polylogs evaluation. Also this kind of syntax + for (lst::const_iterator itx = x.begin(); itx != x.end(); ++itx) newx.append( *itx != _ex_1 ? 1 : -1 ); often gives rise to obscure bugs (C++ gives you enough rope to hang yourself), please avoid that. Use the following instead: for (lst::const_iterator itx = x.begin(); itx != x.end(); ++itx) { newx.append(*itx != _ex_1 ? 1 : -1); } Also could you please make a few test cases so the bug won't be re-introduced in the future (by changes in other parts of GiNaC, or CLN, or whatever)? Thanks a lot, Alexei
Hi, On 01/17/2014 09:22 AM, Stefan Weinzierl wrote:
the method Li_eval checks if all entries of x_ are equal to +1 or -1. In this case, the function Li is converted to a multiple zeta value. If all but the first entry of x_ are equal to +1 or -1, the function Li is converted to a harmonic polylog H. After checking that an entry of x_ is equal to +1 or -1, the operator > is used. This may lead to an "numeric::operator>(): complex inequality" exception in the case where the real part of an entry of x_ is equal to +/-1 (either exact or as a float, this does not matter) and the imaginary part is equal to zero as a float. The attached patch avoids the improper use of the operator > in this cases and sets the +/-1 arguments of H and zeta as real numbers. For zeta the method Li_eval is modified, for H the appropriate place is a modification of the method convert_parameter_Li_to_H.
I see: If one of the x_ is ±1+0.0*I, then it should be assumed to be ±1. As discussed, this would go against the principles of CLN. But those principles are not GiNaC's which simplifies 2.0*x-x-x to 0, not to 0.0*x. This just makes me wonder if the elimination of a floating-point zero imaginary part would make a useful transformation of class numeric. -richy. -- Richard B. Kreckel <http://www.ginac.de/~kreckel/>
Hi, On Sun, 19 Jan 2014, Richard B. Kreckel wrote:
I see: If one of the x_ is ±1+0.0*I, then it should be assumed to be ±1. As discussed, this would go against the principles of CLN. But those principles are not GiNaC's which simplifies 2.0*x-x-x to 0, not to 0.0*x.
This just makes me wonder if the elimination of a floating-point zero imaginary part would make a useful transformation of class numeric.
I'm not sure if we should do this. We have in GiNaC also the opposite case, namely that a floating 0.0 is not transformed into 0, as the following code fragment shows: ex t1 = numeric("2.1"); ex t3 = t1-t1; if ( (t3 == 0) && !ex_to<numeric>(t3).is_rational() ) std::cout << "equal to zero, but not rational" << std::endl; The core of the problem are two different philosophies, what a floating point number should be. These two philosophies are: 1. A floating point number represents an intervall, whose size depends on Digits. 2. A floating point number is a rational number, which is the closest approximation within Digits to the number we want. The GinaC class numeric is basically a wrapper around CLN, so this is a CLN-issue. In my opinion CLN is not strict in implementing one or the other philosophy. This can be seen by looking at the GiNaC-operators is_real() and ==, which boil down to the CLN methods instanceof(..., cln::cl_R_ring) and equal(...,...). As we have discussed, 1+0.0*I is not considered to be real, and thus an example of philosophy 1. On the other hand, comparing a float to an exact number as in 1+0.0*I==1, the float is converted to a rational number and compared (the details can be found in cln-1.3.3/src/real/elem/cl_R_equal.cc). This is basically philosophy 2. However, I would not argue to change this. CLN is a well-established library, and changing one of these operators might break other third-party programs. Besides, the manual of CLN is quite clear what these operators do. For GiNaC the documentation is not as precise, but this can be improved. Therefore my suggestions are: - a small correction to the manual (taking Alexei's comment into account) - a patch to inifcns_nstdsums.cpp (in order to please Alexei I have inserted some extra {...}, however the modification should be in Li_eval, there one knows that the arguments are actually +/-1, in zeta2_eval the only sensible thing to do would be to call hold(), which is not what we want.) A testcase is probably not so helpful, as it can only check for the re-introduction of the bug in the same place. It is more likely that the same bug (or should I say the same feature) in different disguise will hit us somewhere else. I have some hope that we can close this thread. Best wishes, Stefan
participants (3)
-
Alexei Sheplyakov
-
Richard B. Kreckel
-
Stefan Weinzierl