Quadruple tests failing with gfortran 4.7

[Juergen Reuter]

The following tests have screwed up numerics for quadruple precision and 64bit with gfortran 4.7. With 32bit they work! ep_1, circe1_4/5/6/7/8/9, beam_events_1, circe1_photons_4/5.

[Juergen Reuter]

So this is not really just numerical fluctuations for the test ep_1:

152,162c152,162
<    1       1000  2.1381104E+06  5.02E+05   23.47    7.42*   0.85
<    2       1000  2.5346387E+06  1.62E+05    6.40    2.02*   6.77
<    3       1000  2.3170573E+06  8.34E+04    3.60    1.14*  13.21
<    4       1000  2.3382587E+06  7.16E+04    3.06    0.97*  21.35
<    5       1000  2.2856695E+06  7.89E+04    3.45    1.09    9.17
< |-----------------------------------------------------------------------------|
<    5       5000  2.3293453E+06  4.30E+04    1.85    1.30    9.17    0.52   5
< |-----------------------------------------------------------------------------|
<    6       1000  2.4680568E+06  7.99E+04    3.24    1.02*   8.66
<    7       1000  2.3201634E+06  7.46E+04    3.22    1.02*   8.14
<    8       1000  2.3838999E+06  9.20E+04    3.86    1.22    5.73
---
>    1       1000  5.2715269E+06  1.58E+06   29.97    9.48*   0.55
>    2       1000  3.6650640E+07  1.37E+07   37.43   11.84    0.41
>    3       1000  2.8941880E+07  7.81E+06   26.97    8.53*   0.55
>    4       1000  4.2866974E+07  1.19E+07   27.76    8.78    0.75
>    5       1000  2.4872965E+07  2.67E+06   10.73    3.39*   1.88
> |-----------------------------------------------------------------------------|
>    5       5000  1.1542817E+07  1.32E+06   11.48    8.11    1.88   13.99   5
> |-----------------------------------------------------------------------------|
>    6       1000  2.6983659E+07  3.11E+06   11.51    3.64    1.45
>    7       1000  2.7432331E+07  2.51E+06    9.13    2.89*   1.47
>    8       1000  3.1026378E+07  4.38E+06   14.13    4.47    1.15
164c164
<    8       3000  2.3876957E+06  4.69E+04    1.96    1.08    5.73    0.92   3
---
>    8       3000  2.7878403E+07  1.78E+06    6.39    3.50    1.15    0.32   3
199,209c199,209
<    1       1000  2.4793091E+06  4.31E+05   17.39    5.50*   1.41
<    2       1000  2.3512041E+06  1.22E+05    5.18    1.64*   7.94
<    3       1000  2.3182192E+06  8.54E+04    3.68    1.16*  10.22
<    4       1000  2.3068350E+06  1.03E+05    4.47    1.42    5.23
<    5       1000  2.2911108E+06  8.73E+04    3.81    1.21*  13.66
< |-----------------------------------------------------------------------------|
<    5       5000  2.3147028E+06  4.79E+04    2.07    1.46   13.66    0.08   5
< |-----------------------------------------------------------------------------|
<    6       1000  2.3778027E+06  9.24E+04    3.88    1.23   10.18
<    7       1000  2.2492098E+06  8.91E+04    3.96    1.25    9.41
<    8       1000  2.2857040E+06  9.22E+04    4.04    1.28    6.95
---
>    1       1000  5.3119477E+06  1.69E+06   31.78   10.05*   0.60
>    2       1000  9.1371280E+06  1.51E+06   16.51    5.22*   1.14
>    3       1000  2.5726797E+07  7.14E+06   27.75    8.78    0.73
>    4       1000  1.8681845E+07  1.81E+06    9.70    3.07*   2.03
>    5       1000  4.2978316E+07  2.28E+07   53.08   16.78    0.37
> |-----------------------------------------------------------------------------|
>    5       5000  1.0890634E+07  9.46E+05    8.69    6.15    0.37    9.27   5
> |-----------------------------------------------------------------------------|
>    6       1000  2.2433460E+07  2.62E+06   11.68    3.69*   0.19
>    7       1000  2.7322556E+07  2.91E+06   10.65    3.37*   0.22
>    8       1000  1.9324002E+07  1.62E+06    8.39    2.65*   0.16
211c211
<    8       3000  2.3028833E+06  5.26E+04    2.29    1.25    6.95    0.53   3
---
>    8       3000  2.1493793E+07  1.25E+06    5.80    3.18    0.16    2.96   3

[Juergen Reuter]

Whoop, whoop: jenkins with 4.7.3 says: PASS. DESY computers 64bit with gfortran 4.7.1 say: FAIL. DESY computers 32bit with gfortran 4.7.1 say: PASS. WTF!?

[Juergen Reuter]

Uff: so here is what I learned: I thought it should have something to do with the failing testsuite of the library mpc v1.0.2 (a prerequisite of gcc) on the DESY cluster, such that I had to fall back on v1.0.1. The failing tests deal with roundings and exceptional values of atan and atanh function implementations. All compiler 4.7.x give quadruple test failures on the DESY machines. Then I thought it could be the difference between the Intel i5-2500 of my desktop and the Intel Xeon E5-2440 of the cluster computers. So I compiled exclusively on the i5: same result: tests are failing. I'm stymied: is it then one of the system libraries (libm.so.6 or libc.so.6 ?).

[Juergen Reuter]

I do understand more now: the culprit is the ISR structure function, I checked that neither the implementation of the dilogarithm (which only appears in the 3rd order ISR SF) nor WK's definition of the log_prec function (Taylor expansion of the logarithm close to x=1) are responsible. However, I believe it is the map_power_1 function and its inverse. Before there was a catch on too small arguments which is gone now. Further investigating ...

[Juergen Reuter]

Introducing the old rb > tiny(1._default)**eps barrier does not change anything for quadruple precision with 4.9.0 on MAC OS X, a slight numerical fluctuation for double precision, but completely screws up quadruple precision for 4.7.1 on DESY machines.

[Juergen Reuter]

I believe now that this is an incarnation of #514. The deviations in the tests on the DESY machines I think I can trace back to power law mapping for the ISR structure function, which involves an (1-x)^{eps. Here, eps is between 10}-2 and 10^{-3 while (1-x) can become as small as 10}-10 to 10^{-17. The 10byte precision we defined for quadruple gives deviations between my Laptop and the DESY machines for (1-x)=O(10}-17) already in the second digit. In the MC sampling this then sums up to the fluctuations leading to the failures of the test in which ISR is switched on. It seems that when I switch to true real(16) precision on the DESY machines, I get the same values. Here are the results: MAC OS X 64bit 4.7.1 quadruple according to kinds:

 ******************************************
 single =            4
 double =            8
 quad   =           10
 ******************************************
 alpha =    7.29927007299270072990E-0003
 log   =    13.7937490661827614124      
 eps   =    6.17741861636862526793E-0002
 tiny(1_double) =    2.2250738585072014E-308
 tiny(1_defaut) =    3.36210314311209350626E-4932
 tiny(1_quadru) =    3.36210314311209350626E-4932
 tiny(1_double)**eps =   9.88570351238654583308E-0020
 tiny(1_defaut)**eps =   2.30273948926999148916E-0305
 tiny(1_quadru)**eps =    2.30273948926999148916E-0305
 
 a          =    9.99999999999999999997E-0018
 1/eps      =    16.1879914913042858120      
 a**(1/eps) =    6.37007649851328957159E-0276

On the DESY machine, 64bit with gfortran 4.7.1, 10byte:

 ******************************************
 single =            4
 double =            8
 quad   =           10
 ******************************************
 alpha =    7.29927007299270072990E-0003
 log   =    13.7937490661827614124      
 eps   =    6.17741861636862526793E-0002
 tiny(1_double) =    2.2250738585072014E-308
 tiny(1_defaut) =    3.36210314311209350626E-4932
 tiny(1_quadru) =    3.36210314311209350626E-4932
 tiny(1_double)**eps =  0.958085168108412732555      
 tiny(1_defaut)**eps =  0.958085168108412732555      
 tiny(1_quadru)**eps =   0.958085168108412732555      
 
 a           =    9.99999999999999999997E-0018
 a16         =    9.99999999999999999999999999999999992E-0018
 1/eps       =    16.1879914913042858120      
 a**(1/eps)  =    6.11374201115035738970E-0276
 a16**(1/eps)=    6.37007649851328957003759085721633792E-0276

Maybe we should discuss about this!

[kilian]

Note: after r5549, beam_events_1 will probably work, but _2 and _3 will fail (they call ISR).

So is it the real(10) implementation in 4.7 which is different from 4.8 (and from MAC), or do the 4.8/4.9 versions resolve quadruple into real(16)?

I'm sorry that the ISR implementation depends on numerical details. Unfortunately, a naive cutoff does actually cut off events if we do not follow 1-x down to really tiny values. If I would resolve it like CIRCE and treat events below a cutoff specially (e.g., as a delta peak), we would not be able to map the structure function efficiently. And it appears as if the double-prec results are still ok, and not invalidated by (true) quad prec results ...?

In other words, the results agree within the Monte Carlo accuracy in all cases, as they should, and the numerical noise just lets the result fluctuate within errors?

[Juergen Reuter]

Finally resolved in r5554. Culprit indeed was numerical noise in the power mapping of the endpoint peaks in the ISR with quadruple precision. Since the quadruple precision does not change the results there (normally), we safeguard against this noise by explicitly enforcing double precision in the mapping routines. Happily closing.