NO

Author Topic: Floating point SSE2 error on x64?  (Read 403 times)

Offline Dig386

  • Member
  • *
  • Posts: 0
Floating point SSE2 error on x64?
« on: September 03, 2019, 07:41:27 pm »
Hello!
It seems that I've found a bug in Pelles C compiler.
The attached computer program gives different results at different compiler settings. The difference is very far beyond rounding errors and the code doesn't include precision-sensitive things such as numeric differentiation.  It very intensively uses floating point arithmetics.

The next compilation keys (enabled speed optimization for x64 platform with SSE) cause problems
pocc   x64_glitch.c /Tx64-coff /Ot /arch:SSE2
polink x64_glitch.obj /OUT:x64_glitch.exe

The next variant (without optimization) causes no problems.
pocc   x64_glitch.c /Tx64-coff
polink x64_glitch.obj /OUT:x64_noglitch.exe

The correctness/incorrectness was checked using other compilers (MinGW, MSVC, Open Watcom) and scientific literature (this piece of code reproduces a fragment of water density equation).

Offline TimoVJL

  • Global Moderator
  • Member
  • *****
  • Posts: 1882
Re: Floating point SSE2 error on x64?
« Reply #1 on: September 04, 2019, 12:37:27 pm »
optimizer bug in indexing :(
Code: [Select]
#include <stdio.h>
#include <math.h>

static double Acoeff[10][7] =
{
{29.492937,  -5.1985860,  6.8335354, -0.1564104,  -6.3972405, -3.9661401, -0.69048554},
{-132.13917,  7.7779182, -26.149751, -0.72546108, 26.409282,  15.453061,   2.7407416},
{274.64632,  -33.301902,  65.326396, -9.2734289, -47.740374, -29.142470,  -5.1028070},
{-360.93828, -16.254622, -26.181978,  4.3125840,  56.323130,  29.568796,   3.9636085},
{342.18431, -177.31074,   0,          0,          0,          0,           0},
{-244.50042, 127.48742,   0,          0,          0,          0,           0},
{155.18535,  137.46153,   0,          0,          0,          0,           0},
{5.9728487,  155.97836,   0,          0,          0,          0,           0},
{-410.30848, 337.31180, -137.46618,   6.7874983, 136.87317,   79.847970,  13.041253},
{-416.05860,-209.88866, -733.96848,  10.401717,  645.81880,  399.17570,   71.531353}
};

#define A(I,J) (Acoeff[(I-1)][(J-1)])

double test(double d, double T)
{
int i, j;
for (j = 1; j <= 2; j++) {
for (i = 1; i <= 8; i++) {
printf("\tA(%d,%d) = %f\n", i, j, A(i,j));
}
printf("\n");
}
return 0;
}

int main(void)
{
double d = 1000, T = 273.15;
printf("%g\n", test(d, T));
return 0;
}
test for poasm:
Code: [Select]
extern printf :proc
includelib msvcrt.lib

.rdata
Acoeff real8 29.492937,  -5.1985860,  6.8335354, -0.1564104,  -6.3972405, -3.9661401, -0.69048554,
-132.13917,  7.7779182, -26.149751, -0.72546108, 26.409282,  15.453061,   2.7407416,
274.64632,  -33.301902,  65.326396, -9.2734289, -47.740374, -29.142470,  -5.1028070,
-360.93828, -16.254622, -26.181978,  4.3125840,  56.323130,  29.568796,   3.9636085,
342.18431, -177.31074,   0,          0,          0,          0,           0,
-244.50042, 127.48742,   0,          0,          0,          0,           0,
155.18535,  137.46153,   0,          0,          0,          0,           0,
5.9728487,  155.97836,   0,          0,          0,          0,           0,
-410.30848, 337.31180, -137.46618,   6.7874983, 136.87317,   79.847970,  13.041253,
-416.05860,-209.88866, -733.96848,  10.401717,  645.81880,  399.17570,   71.531353

fmt1 db 9h, 'A(i,j) = %f',13,10,0
crlf db 13,10,0
@81     db "%g",13,10,0
@80     dq 4071126666666666H                            ; 0008 _ 273.15
@79     dq 408F400000000000H                            ; 0010 _ 1000.0

.code

test_loop    PROC
        push    rbx                                     ; 0000 _ 53
        push    rsi                                     ; 0001 _ 56
        push    rdi                                     ; 0002 _ 57
        push    r12                                     ; 0003 _ 41: 54
        sub     rsp, 40                                 ; 0005 _ 48: 83. EC, 28
        mov     ebx, 1                                  ; 0009 _ BB, 00000001
        lea     rsi, qword ptr [Acoeff]                 ; 000E _ 48: 8D. 35, 00000000(rel)
;lea     rdi, qword ptr [rsi-38H] ; fix base 1
?_001:  ;lea     rdi, qword ptr [rsi+38H]                ; 0015 _ 48: 8D. 7E, 38
lea     rdi, qword ptr [rsi-38H] ; fix base 1
        mov     r12d, 8                                 ; 0019 _ 41: BC, 00000008
?_002:  movsxd  rax, ebx                                ; 001F _ 48: 63. C3
        movsd   xmm1, qword ptr [rdi+rax*8+30H]         ; 0022 _ F2: 0F 10. 4C C7, 30
        lea     rcx, qword ptr [fmt1]                   ; 0028 _ 48: 8D. 0D, 00000000(rel)
        movq    rdx, xmm1                               ; 002F _ 66 48: 0F 7E. CA
        call    printf                                  ; 0034 _ E8, 00000000(rel)
        add     rdi, 56                                 ; 0039 _ 48: 83. C7, 38
        sub     r12d, 1                                 ; 003D _ 41: 83. EC, 01
        jnz     ?_002                                   ; 0041 _ 75, DC
        lea     rcx, qword ptr [crlf]                   ; 0043 _ 48: 8D. 0D, 00000000(rel)
        call    printf                                  ; 004A _ E8, 00000000(rel)
        add     ebx, 1                                  ; 004F _ 83. C3, 01
        cmp     ebx, 2                                  ; 0052 _ 83. FB, 02
        jle     ?_001                                   ; 0055 _ 7E, BE
        xorpd   xmm0, xmm0                              ; 0057 _ 66: 0F 57. C0
        add     rsp, 40                                 ; 005B _ 48: 83. C4, 28
        pop     r12                                     ; 005F _ 41: 5C
        pop     rdi                                     ; 0061 _ 5F
        pop     rsi                                     ; 0062 _ 5E
        pop     rbx                                     ; 0063 _ 5B
        ret                                             ; 0064 _ C3
test_loop    ENDP

mainCRTStartup proc
        sub     rsp, 40                                 ; 0070 _ 48: 83. EC, 28
        movsd   xmm1, qword ptr [@80]                   ; 0074 _ F2: 0F 10. 0D, 00000000(rel)
        movsd   xmm0, qword ptr [@79]                   ; 007C _ F2: 0F 10. 05, 00000000(rel)
        call    test_loop                               ; 0084 _ E8, FFFFFF77
        movapd  xmm1, xmm0                              ; 0089 _ 66: 0F 28. C8
        lea     rcx, qword ptr [@81]                          ; 008D _ 48: 8D. 0D, 00000000(rel)
        movq    rdx, xmm1                               ; 0094 _ 66 48: 0F 7E. CA
        call    printf                                  ; 0099 _ E8, 00000000(rel)
        xor     eax, eax                                ; 009E _ 31. C0
        add     rsp, 40                                 ; 00A0 _ 48: 83. C4, 28
        ret                                             ; 00A4 _ C3
mainCRTStartup endp
end
« Last Edit: September 04, 2019, 03:46:24 pm by TimoVJL »
May the source be with you