Существует ли правильное выражение константы с точки зрения числа с плавающей запятой для его msb? - PullRequest
Новая
       50

Существует ли правильное выражение константы с точки зрения числа с плавающей запятой для его msb?

0 голосов
/ 03 декабря 2018

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

Для целей этого вопроса мы можем игнорировать:

  • Почти переполнение или почти-поток значений (они могут обрабатываться конечным числом приложений ?: для масштабирования).
  • Отрицательные входы (они могут обрабатываться аналогичным образом).
  • Реализации, не соответствующие Приложению-F(на самом деле с ними ничего полезного не получается).
  • Странность вокруг избыточной точности (float_t и double_t могут использоваться с FLT_EVAL_METHOD и другими макросами float.h для безопасной обработки).

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

Обратите внимание, что эта задача эквивалентна нахождению "эпсилона" дляконкретное значение, то есть nextafter(x,INF)-x (или эквивалент в float или long double), а результат только масштабируется на DBL_EPSILON (или эквивалент для типа).Решения, которые сочтут вполне приемлемыми, если они проще.

У меня есть предложенное решение, которое я публикую в качестве самостоятельного ответа, но я не уверен, что оно правильное.

Ответы [ 3 ]

0 голосов
/ 03 декабря 2018

Вот код для поиска ULP.Он был вдохновлен алгоритмом 3.5 в Точное суммирование с плавающей точкой Зигфрендом М. Рампом, Такеши Огита и Синьити Оиси (который вычисляет 2 ⌈log 2 | p | ⌉ ): 

double ULP(double q)
{
    // SmallestPositive is the smallest positive floating-point number.
    static const double SmallestPositive = DBL_EPSILON * DBL_MIN;

    /*  Scale is .75 ULP, so multiplying it by any significand in [1, 2) yields
        something in [.75 ULP, 1.5 ULP) (even with rounding).
    */
    static const double Scale = 0.75 * DBL_EPSILON;

    q = fabs(q);

    // Handle denormals, and get the lowest normal exponent as a bonus.
    if (q < 2*DBL_MIN)
        return SmallestPositive;

    /*  Subtract from q something more than .5 ULP but less than 1.5 ULP.  That
        must produce q - 1 ULP.  Then subtract that from q, and we get 1 ULP.

        The significand 1 is of particular interest.  We subtract .75 ULP from
        q, which is midway between the greatest two floating-point numbers less
        than q.  Since we round to even, the lesser one is selected, which is
        less than q by 1 ULP of q, although 2 ULP of itself.
    */
    return q - (q - q * Scale);
}

fabs и if можно заменить на ?:.

Для справки: 2 ⌈log 2 | p | ⌉ Алгоритм: 

q = p / FLT_EPSILON
L = |(q+p) - q|
if L = 0
    L = |p|
0 голосов
/ 03 декабря 2018

Если , вы можете принять формат двоичной64 IEEE 754 и семантику (и, в частности, правильное округление арифметических операций) и режим округления с округлением до четности, то это хороший факт, что длялюбое не слишком маленькое не слишком большое положительное конечное double значение x, следующее представимое значение по сравнению с x всегда задается как x / 0x1.fffffffffffffp-1 (где 0x1.fffffffffffffp-1 - это просто 1.0 - 0.5 * DBL_EPSILON, обозначенное какшестнадцатеричный литерал).

Таким образом, мы можем получить наиболее значимый бит, который вы запрашиваете просто: 

(x / 0x1.fffffffffffffp-1 - x) * 0x1.0p+52

И, конечно, есть аналогичные результаты для float, предполагая IEEE 754Формат и семантика binary32.

Фактически, единственное нормальное положительное значение, для которого это не удается, - DBL_MAX, где результат деления переполняется до бесконечности.

Чтобы показать, что трюк деленияработает, достаточно доказать это для x в диапазоне 1.0 <= x < 2.0;Легко показать, что для любого x в этом диапазоне значение x / 0x1.fffffffffffffp-1 - x (где / представляет математическое деление в этом случае) лежит в полуоткрытом интервале (2^-53, 2^52], и из этого следует, что при округлениисвязь до четного (или фактически любой режим округления до ближайшего), x / 0x1.fffffffffffffp-1 округляет до следующего представимого значения.

Аналогично, при тех же предположениях, x * 0x1.fffffffffffffp-1 всегда является следующимпредставительное значение по сравнению с x.

0 голосов
/ 03 декабря 2018

Для примера предположим, что типом является float, и пусть x будет вводом.Первоначально я напишу это как последовательность операторов для удобства чтения, но они могут быть переведены непосредственно в макросы, которые производят константные выражения. 

float y = x*(1+FLT_EPSILON)-x;
if (y/FLT_EPSILON > x) y/=2;

Если бы мы могли обеспечить округление, начальное значение yдолжно быть именно то, что мы хотим.Однако, если верхние два бита x равны 1 и установлены младшие биты, или если мы выполнили случай округления до четного, x*(1+FLT_EPSILON) может превысить x на 2 единицы в последнем месте, а не просто1. Я не верю, что возможны какие-либо другие случаи, и я полагаю, что вторая строка полностью отвечает за это.

Записано в виде макросов: 

#define PRE_ULP(x) ((x)*(1+FLT_EPSILON)-(x))
#define ULP(x) ((PRE_ULP(x)/FLT_EPSILON>(x) ? PRE_ULP(x)/2 : PRE_ULP(x))

#define MSB_VAL(x) (ULP(x)/FLT_EPSILON)
...