Как реализовать быстрый обратный sqrt без поведения undefined?

Из того, что я понял о правиле строгого aliasing, этот код для быстрый обратный квадратный корень приведет к поведению undefined в С++:

float Q_rsqrt( float number )
{
    long i;
    float x2, y;
    const float threehalfs = 1.5F;

    x2 = number * 0.5F;
    y  = number;
    i  = * ( long * ) &y; // type punning
    i  = 0x5f3759df - ( i >> 1 );
    y  = * ( float * ) &i;
    y  = y * ( threehalfs - ( x2 * y * y ) );

    return y;
}

Действительно ли этот код вызывает UB? Если да, то как это можно переопределить стандартным образом? Если нет, почему бы и нет?

Предположения: перед вызовом этой функции мы как-то проверили, что поплавки находятся в 32-битном формате IEEE 754, sizeof(long)==sizeof(float), а платформа малозначительна.

Ответы

Ответ 1

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

Ответ 2

Стандартным образом будет std::memcpy. Это должно быть достаточно стандартно в соответствии с вашими предположениями. Любой разумный компилятор превратит это в кучу перемещений регистра, если это возможно. Кроме того, мы можем также облегчить (или, по крайней мере, проверить) некоторые ваши сделанные предположения, используя С++ 11 static_assert и фиксированные ширины целочисленных типов из <cstdint>. В любом случае, Endianness не имеет значения, так как мы не работаем над любыми массивами здесь, если целочисленный тип малозначен, тип с плавающей запятой также.

float Q_rsqrt( float number )
{
    static_assert(std::numeric_limits<float>::is_iec559, 
                  "fast inverse square root requires IEEE-comliant 'float'");
    static_assert(sizeof(float)==sizeof(std::uint32_t), 
                  "fast inverse square root requires 'float' to be 32-bit");
    float x2 = number * 0.5F, y = number;
    std::uint32_t i;
    std::memcpy(&i, &y, sizeof(float));
    i  = 0x5f3759df - ( i >> 1 );
    std::memcpy(&y, &i, sizeof(float));
    return y * ( 1.5F - ( x2 * y * y ) );
}

Ответ 3

Я не думаю, что вы можете сделать это без UB в строгом стандартном смысле, просто потому, что он полагается на такие вещи, как конкретное представление значений float и long. Стандарт не указывает эти (специально) и дает реализациям свободу делать то, что они считают подходящим в этом отношении, в частности применяет ярлык "UB" ко всему поведению, которое будет зависеть от таких вещей.

Это не означает, что это будет тип "Heisenbug" или "носовые демоны" UB; скорее всего, реализация даст определения для этого поведения. Пока эти определения удовлетворяют предварительным условиям кода, это прекрасно. Предпосылками здесь являются:

Тип punning разрешен реализацией.
sizeof(float) == sizeof(long).
Внутреннее представление float таково, что если биты интерпретируются как long, 0x5f3759df и >> 1 имеют значение, необходимое для работы этого хака.

Однако невозможно перезаписать это полностью переносимым способом (т.е. без UB) - что бы вы сделали для платформ, которые используют, например, другую реализацию арифметики с плавающей запятой. Пока вы полагаетесь на специфику конкретной платформы, вы полагаетесь на поведение undefined в стандарте, но определенное платформой.

Ответ 4

Не. Алгоритм начинается с предположения о макете IEEE754, который не является допустимым стандартным допуском. Конкретная константа 0x5f3759df была бы совершенно неправильной для любого другого макета. Даже предположение о существовании такой константы является ошибочным. Я думаю, что он разбивается, как только вы, например, измените порядок полей внутри float.

Я также не уверен, что Q_rsqrt(-0.0f) работает правильно, даже если IEEE754. Оба нуля (положительные и отрицательные) должны быть равны и sqrt(x) является только undefined для x<-0.0f