К вопросу о точности вычислений с double

Realist

Предлагается заценить следующий код:
  

anar:~$ cat 2.cpp

#include <stdio.h>

int main

{

double a=0.795774999999999954614082753323600627481937408447265625;

double x=5.6967615525588122249445177658344618976116180419921875;

double b=0.5515889999999999682156470726113184355199337005615234375;

double y=8.2186925853987091716845725386519916355609893798828125;

double a1=a*x;

double a2=b*y;

double a3=a1-a2;

double res = a*x-b*y;

printf("res=%.60e \na1=%.60e\na2=%.60e\na3=%.60e\n",res,a1,a2,a3);

return 0;

}

anar:~$ g++ -O0 2.cpp

anar:~$ ./a.out

res=1.049507702965968292119214311242103576660156250000000000000000e-16

a1=4.533340424487488640181709342868998646736145019531250000000000e+00

a2=4.533340424487488640181709342868998646736145019531250000000000e+00

a3=0.000000000000000000000000000000000000000000000000000000000000e+00

anar:~$

yolki



res=-5.273559366969493567000000000000000000000e-16

a1=4.533340424487488400000000000000000000000e+00

a2=4.533340424487488928000000000000000000000e+00

a3=-5.273559366969493567000000000000000000000e-16






#include <stdio.h>

int main

{

long double a=0.795774999999999954614082753323600627481937408447265625;

long double x=5.6967615525588122249445177658344618976116180419921875;

long double b=0.5515889999999999682156470726113184355199337005615234375;

long double y=8.2186925853987091716845725386519916355609893798828125;

long double a1=a*x;

long double a2=b*y;

long double a3=a1-a2;

long double res = a*x-b*y;

printf("res=%.60Le \na1=%.60Le\na2=%.60Le\na3=%.60Le\n",res,a1,a2,a3);

return 0;

}

yolki

могу только предположить, что это флюктуации сопроцессора, который вычисляет всё в лонг-даблах.

yolki

ага, так оно и есть.
res вычисляется через long double, а a3 - через даблы. вот и получается то, что видим.
для гимнастики ума советую ещё подумать, почему такой код:


#include <stdio.h>

int main

{

 double one, eps;

 one=1.0;

 eps=1.0;

 while(one+eps>one)

 {

  eps/=2;

 }

 printf("%le\n",eps);

 return 0;

}

печатает ~1e-21, хотя вот такой код:


#include <stdio.h>

int main

{

 double one, eps, one_plus_eps;

 one=1.0;

 eps=1.0;

 one_plus_eps=one+eps;

 

 while(one_plus_eps>one)

 {

  eps/=2;

  one_plus_eps=one+eps;

 }

 printf("%le\n",eps);

 return 0;

}

печатает ~1e-16

vall

$./a.out
res=0.000000000000000000000000000000000000000000000000000000000000e+00
a1=4.533340424487488640181709342868998646736145019531250000000000e+00
a2=4.533340424487488640181709342868998646736145019531250000000000e+00
a3=0.000000000000000000000000000000000000000000000000000000000000e+00
/tmp
$gcc -v
Using built-in specs.
Target: x86_64-linux-gnu
Configured with: ../src/configure -v --enable-languages=c,c++,java,f95,objc,ada,treelang --prefix=/usr --enable-shared --with-system-zlib --libexecdir=/usr/lib --without-included-gettext --enable-threads=posix --enable-nls --program-suffix=-4.0 --enable-__cxa_atexit --enable-clocale=gnu --enable-libstdcxx-debug --enable-java-awt=gtk --enable-gtk-cairo --with-java-home=/usr/lib/jvm/java-1.4.2-gcj-4.0-1.4.2.0/jre --enable-mpfr --disable-werror --enable-checking=release x86_64-linux-gnu
Thread model: posix
gcc version 4.0.3 (Debian 4.0.3-1)

Marinavo_0507

        fldl    .LC0

        fstpl   -8(%ebp)        #  a

        fldl    .LC1

        fstpl   -16(%ebp)       #  x

        fldl    .LC2

        fstpl   -24(%ebp)       #  b

        fldl    .LC3

        fstpl   -32(%ebp)       #  y

        fldl    -8(%ebp)        #  a

        fmull   -16(%ebp)       #  x

        fstpl   -40(%ebp)       #  a1

        fldl    -24(%ebp)       #  b

        fmull   -32(%ebp)       #  y

        fstpl   -48(%ebp)       #  a2

        fldl    -40(%ebp)       #  a1

        fsubl   -48(%ebp)       #  a2

        fstpl   -56(%ebp)       #  a3

        fldl    -8(%ebp)        #  a

        fmull   -16(%ebp)       #  x

        fldl    -24(%ebp)       #  b

        fmull   -32(%ebp)       #  y

        fsubrp  %st, %st(1)

        fstpl   -64(%ebp)       #  res

        fldl    -56(%ebp)       #  a3

        fstpl   28(%esp)

        fldl    -48(%ebp)       #  a2

        fstpl   20(%esp)

        fldl    -40(%ebp)       #  a1

        fstpl   12(%esp)

        fldl    -64(%ebp)       #  res

        fstpl   4(%esp)

        movl    $.LC4, (%esp)

        call    printf

Регистры сопроцессора 80-битные, поэтому точность вычислений без промежуточного хранения выше.
Если включить какую-нибудь оптимизацию, то res посчитается на этапе компиляции, и будет нулевым.

vall

ты забыл к константам L приписать. они у тебя до дабла порезались.
если поставить их то
res=2.072994553792284477822249755263328552246093750000000000000000e-16
a1=4.533340424487488539567747736214187170844525098800659179687500e+00
a2=4.533340424487488332268292356985739388619549572467803955078125e+00
a3=2.072994553792284477822249755263328552246093750000000000000000e-16

Landstreicher

На x86_64 это дело считается через MMX/SSE:


        movlpd  -64(%rbp %xmm0        # a, a

        mulsd   -56(%rbp %xmm0        # x, tmp67

        movsd   %xmm0, -32(%rbp)        # tmp67, a1

        movlpd  -48(%rbp %xmm0        # b, b

        mulsd   -40(%rbp %xmm0        # y, tmp69

        movsd   %xmm0, -24(%rbp)        # tmp69, a2

        movlpd  -32(%rbp %xmm0        # a1, a1

        subsd   -24(%rbp %xmm0        # a2, tmp71

        movsd   %xmm0, -16(%rbp)        # tmp71, a3

        movlpd  -64(%rbp %xmm0        # a, a

        movsd   %xmm0, %xmm1    # a, D.2531

        mulsd   -56(%rbp %xmm1        # x, D.2531

        movlpd  -48(%rbp %xmm0        # b, b

        mulsd   -40(%rbp %xmm0        # y, D.2532

        movsd   %xmm1, %xmm2    # D.2531,

        subsd   %xmm0, %xmm2    # D.2532,

        movsd   %xmm2, %xmm0    #, tmp74

        movsd   %xmm0, -8(%rbp) # tmp74, res

Почему-то у них другая точность (по всей видимости, double)
Вообще, допустим я пишу программу, которая делает какие-то вычисления с числами с плавающей точкой. Требуется, чтобы на любой платформе (в данном случае --- x86_64, x86, разные режимы оптимизации) результат был абсолютно одинаков. Как этого достичь? Разве вычисления с floating point не регулируются стандартом IEEE?

Marinavo_0507

Требуется, чтобы на любой платформе (в данном случае --- x86_64, x86, разные режимы оптимизации) результат был абсолютно одинаков. Как этого достичь?
Используй свою реализацию арифметики с плавающей точкой через целочисленные операции.
Или не свою, их есть.

Landstreicher

> Используй свою реализацию арифметики с плавающей точкой через целочисленные операции.
> Или не свою, их есть.
Я правильно понимаю, что используя стандартный double (или long double) это не добиться?
Может быть можно компилятору сказать какие-нибудь ключи?

vall

насколько я понимаю компилятор делает только оптимизации которые не ухудшат точность.
и во всех вычислениях нужно считать какая погрешность может накопится к этой точке.
точного равенства вообще никогда тут ожидать нельзя.

Chupa

Требуется, чтобы на любой платформе (в данном случае --- x86_64, x86, разные режимы оптимизации) результат был абсолютно одинаков.
Для вычислений с плавающей точкой это очень странное требование, но можно попробовать всегда считать на пределе возможностей железки (80-битный long double в данном случае тогда компилятор не сможет уменьшить точность, а процессор - увеличить. Хотя даже тут могут возникнуть отклонения.

Sharp

Посмотри http://gcc.gnu.org/bugs.html#known раздел Non-bugs.
Там, похоже, именно про это и написано.

xronik111

О да, +1 http://gcc.gnu.org/PR323 -- the most often reported non-bug. Если не лень сходить по ссылке, то можно заценить количество duplicate'ов этого PR.
Оставить комментарий
Имя или ник:
Комментарий: