[x86, fpu] rounding

lilia_rass

Возможно ли при каких-либо условиях равенство:
1/2 + A - A == 2^40
Если да, то при каких.
(Всё по стандарту binary floating point IEEE-какой-то там)

Andbar

операции честные или где-то могли перепутать тип данных?

lilia_rass

Обычные сложение и вычитание плавающих чисел.

Maurog

если речь про
IEEE-какой-то там
, то, думаю, что нет

lilia_rass

Похоже, что всё-таки возможно :)

Maurog

на какой странице стандарта ты это обнаружил?

lilia_rass

Ну это же как раз и будет ответом к задачке :)

Maurog

я был не прав
C:


  int A = 0;

  assert(1/2 + A - A == 2^40);

:grin:

lilia_rass

Замечательно :)
Пример для этого я писал как раз на С, но всё-таки там имеется в виду возведение в степень.

bleyman

Не знаю, что говорит стандарт, но кажется на практике нельзя.
Можно получить единицу:
>>> (2.0 ** 52 + 1).hex
'0x1.0000000000001p+52'
>>> 2.0 ** 52 + 1) - 0.5).hex
'0x1.0000000000000p+52'
типа вот он округлил к ближайшему числу с нулевым младшим битом.
>>> 2.0 ** 52 + 1) - 0.25).hex
'0x1.0000000000001p+52'
а вот — не округлил. Видимо, при округлении оно смотрит на ещё ровно один бит (что в общем-то правильно).
Не знаю, может быть можно как-то его обмануть, используя 80-битность чисел внутри сопроцессора и избирательно - fast fp model, но он же и их тоже наверное правильно округляет...

lilia_rass

Округляет он всегда правильно, важно - в какую сторону :o
> при округлении оно смотрит на ещё ровно один бит
В стандарте на самом деле точно написано, что именно происходит при округлении.

bleyman

Ну расскажи же уже, что же там написано!

lilia_rass

Ну, как я понимаю, решение такое. Допустим, A=2^63 и точность мантиссы - одинарная. Тогда при сложении 1/2 и 2^63 в мантиссе не хватит бит для точного сложения, придётся округлять (будет прерывание PE). Для округления строятся два числа, влезающие в нужную точность, слева и справа от округляемого, и, в зависимости от параметров округления, округление идёт к одному из них. В данном случае данными числами будут 0 и 2^40 соответственно. При дефолтных параметрах округление идёт к ближайшему, поэтому 1/2 округлится до нуля и в результате выражения получим ноль. Но если выставить принудительное округление вверх (rounding control - up то 1/2 округлится до 2^40.
Пример кода под gcc/x86 :o 


#include <stdio.h>

#include <math.h>



void _fpu_inline(void);



main

{

    double A = pow( 2, 63  B = 0.5;



    printf("before inline:\n\tB = %g,\n\tA = %g,\n\tB + A - A = %g\n", B, A, B + A - A );



    _fpu_inline;



    printf("after inline:\n\tB = %g,\n\tA = %g,\n\tB + A - A = %g\n", B, A, B + A - A );

}



 .section    .text,"ax"



.globl  _fpu_inline

.type   _fpu_inline,@function

.align  16

_fpu_inline:

    finit

    fstcw  mem2

    orw    $0x800, mem2

    andw   $0xfcff, mem2

    fldcw  mem2



    ret



.section    .data



.align 64

mem2:           .4byte 0x0

Serab

Да, хитро, не допер (я все искал варианты с переполнением при сложении, на это rounding влияет, почему-то сразу отмел выбор между числами, отличающимися на 2^40. Все относительно :().
Вот, кому угодно, код для студии (а может и под gcc сработает).


#include <stdio.h>

#include <float.h>

#include <math.h>



int main

{

    _controlfp( _RC_UP, _MCW_RC );

    _controlfp( _PC_24, _MCW_PC );

    float two40 = pow( 2.0f, 40.0f );

    float A = pow( 2.0f, 63.0f );

    printf( "%f\n", 0.5f + A - A - two40 );



    return 0;

}

Maurog

Разжуйте, пожалуйста, поподробнее как вы пришли к выводу.
Для округления строятся два числа, влезающие в нужную точность, слева и справа от округляемого, и, в зависимости от параметров округления, округление идёт к одному из них. В данном случае данными числами будут 0 и 2^40 соответственно.
Я не смог найти в драфте каким образом 1\2 может округлиться до 2^40.
Несколько раз просмотрел документ http://www.validlab.com/754R/nonabelian.com/754/comments/Q75...
Буду рад, если вы приведете ссылки на этот документ.
Из этого документа я понял, что округление применяется к infinitely precise number чтобы получить binary floating point format, описано 5 способов (tiesToEven, tiesToAway, towardPositive, towardNegative, towardZero) и каждый из них выбирает closest to initial number (пункты 6.2.1 + 6.2.2). Ближайшим к 1\2 ну никак 2^40 у меня не получается (1\2 замечательно представляется в виде float-a).
Далее переходим к операциям (глава 7 и пункт 7.4.1) и понимаем, что даже 1\2 не надо пытаться округлять, потому что округление делается _после_ применения операции formatOf-addition(source1, source2 то есть как бы выполняется infinitely precise сложение, а затем врубают алгоритм округления, чтобы уложиться в бинарный формат.
В общем, магия какая-то :o

apl13

Короче, ХЗ, что там в стандартах, но FPU в интеле, походу, думает, что (2^63 + 1/2) при включенном округлении вверх = (2^63 + 2^40).

lilia_rass

Ну вроде никакой магии. Можно сначала выполнить infinitely precise add, но потом (1/2 + 2^63) всё равно не влезут в точность сингла. В этот момент мы округлим последний бит мантиссы до 1, что и будет 2^40.

Serab

при включенном округлении вверх
и еще если заставить его вычислять это в single, не в double и не в double extended.

lilia_rass

Для double и extended, думаю, можно аналогичные примеры придумать.

Serab

ну это пнятно :)

apl13

Oh, thank you, Cap. :p

Serab

Не, ну не достаточно просто считать в типе float, надо еще соответствующие флаги выставить. Это действительно так очевидно?

Maurog

вот теперь я понял
все же не 1\2 округляют, а округляют сумму
спасибо
Оставить комментарий
Имя или ник:
Комментарий: