On Wed, Feb 13, 2019 at 07:30:53PM +0100, Jakub Jelinek wrote: > On Wed, Feb 13, 2019 at 07:17:53PM +0100, Thomas Koenig wrote: > > > Bootstrapped/regtested on x86_64-linux and i686-linux (together with > > > Richard's patch), ok for trunk? > > > > A couple of points: > > > > with this patch, we will have an algorithm that will evaluate NORM2 > > in a different way than before, possibly leading to regressions > > in some cases where we previously generated good results and > > would return +/- Inf now. > > I'm not aware of such a case, do you have examples what do you have in mind? > The adjustments are only done by powers of two, the only reason I'm not > using just mpfr_set_exp is to deal with the case where there are extreme > differences in exponents, say 1.0e-8100L*1.0e-8100L+1.0e+16382L*1.0e+16382L, > where because 1.0e+16382L*1.0e+16382L would overflow the very tiny temporary > result needs to be divided by a very large power of 2 and so I want to let > mpfr do the rounding. > > > Also, with this patch, we would use a very different algorithm > > at run-time to compile-time, potentially leading to differences > > far exceeding a few ulps. I am not sure I like that idea. > > The runtime code is already so much different that I don't see how this > holds. And it is actually much worse, exactly because it doesn't scale just > by powers of two and thus every scaling has a rounding error. > > > I'll look at the literature and generate some test cases, and > > come back to you. > > > > From curiosity, did you chose this algorithm from literature? > > No. But, first of all, the patch should change nothing at all unless > either at least one extremely large number (whose **2 is not representable > in the kind or close to that) appears, or the array is so huge that the > result is attacking the maximum. And, for the cases where it does change > something, unless the differences in exponents between numbers are extreme, > there should be again no rounding difference, the exponent of the numbers > will simply change and nothing else. >
I think we can agree that gfortran's norm2 needs some work. To get Richard's emin/emax patch into the tree so adequate testing can be done prior to the 9.1 release, I'll suggest that we XFAIL norm2_3.f90 while we work out the details for norm2. Yesteraday, I looked at Jakub's patch, and thought the approach was fine. I haven't had time to look at the library side. Perhaps, we can adopt Jakub's approach for the library as well. -- Steve
