Какой цели служит таблица поиска в этом коде?

Question

Ниже я адаптировал код Уильяма Кахана и К.С. Нг (см. блок комментариев внизу), написанный в 1986 году, для получения аппроксимации 1 / sqrt(x), где x — число двойной точности с плавающей запятой IEEE-754. Именно этот код адаптировали Клив Молер и Грег Уолш, чтобы он стал « быстрым обратным квадратным корнем », прославившимся благодаря Quake III.

      #include <stdio.h>
int main() {
    double x = 3.1415926535;
    double y, z;
    unsigned long int x0, y0, k;
    unsigned long long int xi, yi;
    static int T2[64]= {
    0x1500, 0x2ef8, 0x4d67, 0x6b02, 0x87be, 0xa395, 0xbe7a, 0xd866,
    0xf14a, 0x1091b,0x11fcd,0x13552,0x14999,0x15c98,0x16e34,0x17e5f,
    0x18d03,0x19a01,0x1a545,0x1ae8a,0x1b5c4,0x1bb01,0x1bfde,0x1c28d,
    0x1c2de,0x1c0db,0x1ba73,0x1b11c,0x1a4b5,0x1953d,0x18266,0x16be0,
    0x1683e,0x179d8,0x18a4d,0x19992,0x1a789,0x1b445,0x1bf61,0x1c989,
    0x1d16d,0x1d77b,0x1dddf,0x1e2ad,0x1e5bf,0x1e6e8,0x1e654,0x1e3cd,
    0x1df2a,0x1d635,0x1cb16,0x1be2c,0x1ae4e,0x19bde,0x1868e,0x16e2e,
    0x1527f,0x1334a,0x11051,0xe951, 0xbe01, 0x8e0d, 0x5924, 0x1edd };
    /* Convert double to unsigned long long (64 bits wide)  */
    xi = *(unsigned long long int*)&x;
    /* Extract only the first half of the bits of the integer */
    x0 = (xi & 0xffffffff00000000) >> 32;
    /* Perform bit shift and subtract from an integer constant to get k */
    k = 0x5fe80000 - (x0 >> 1);
    /* T2 is indexed by mostly by exponent bits. 
       Only 5 highest bits from orig 52 in mantissa are used to index T2 */
    y0 = k - T2[63 & (k >> 14)];
    /* Pad with zeros for the LS 32 bits, convert back to long long */
    yi = ((unsigned long long int) y0 << 32);
    /* Get double from bits making up unsigned long long */ 
    y = *(double*)&yi;
    /* 1/sqrt(pi) ~ 0.564 */
    printf("%lf\n", y);
    return 0;
}

Похоже, он делает то же самое, что и версия Quake (за исключением шага Ньютона-Рафсона).

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

    x2 = number * 0.5F;
    y  = number;
    i  = * ( long * ) &y;                       // evil floating point bit level hacking
    i  = 0x5f3759df - ( i >> 1 );               // what the fuck? 
    y  = * ( float * ) &i;
    y  = y * ( threehalfs - ( x2 * y * y ) );   // 1st iteration
//  y  = y * ( threehalfs - ( x2 * y * y ) );   // 2nd iteration, this can be removed

    return y;
}

Мы начинаем с приведения битов числа double к целому числу, которое является приблизительным двоичным логарифмом числа double, с которого мы начали (метод Митчелла).

Мы сдвигаем эти биты вправо и вычитаем их из константы (для этого сравнения не имеет значения, что младшие 32 бита были отброшены), что является приблизительным делением в логарифмической области.

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

Наконец, мы приводим целые биты (плюс 32 0 для половины LS) как двойное число с плавающей запятой, что дает нам приблизительный антилогарифм.

Итак, я понимаю (или думаю, что понимаю) 3 из 4 шагов здесь, но третий шаг — что делает таблица поиска, как она индексируется по диапазону двойников и почему?

Answer 1

      xi = *(unsigned long long int*)&x;

Это переинтерпретирует биты, представляющие как , который должен быть 64-битным типом. Этот метод переинтерпретации не определен стандартом C, но поддерживается многими компиляторами, возможно, с учетом параметров командной строки.

Переосмысление как делает биты доступными для манипуляций.

      /* Extract only the first half of the bits of the integer */
x0 = (xi & 0xffffffff00000000) >> 32;

Как говорится, теперь содержит старшие 32 бита . Это включает бит знака (который предположительно равен 0, поскольку мы должны брать только обратный квадратный корень из положительного числа), 11 бит, кодирующих показатель степени, и 20 бит основного кодирования мантиссы.

      /* Perform bit shift and subtract from an integer constant to get k */
k = 0x5fe80000 - (x0 >> 1);

начинает приближение квадратного корня. Давайте посмотрим, почему: двоичное представление с плавающей запятой равно ± f •2 ^e для некоторого f в определенном формате с фиксированной точкой и целого числа e в определенном диапазоне. Для положительных чисел мы, конечно, имеем f •2 ^e без знака. Пусть q и r будут частным и остатком от деления e на 2, поэтому e = 2 q + r . Тогда sqrt(f•2^e) = sqrt(f•2^r)•2^q. В , мы сдвинули биты экспоненты вправо, настроив нас на наличие q в поле экспоненты.

Однако не совсем так, потому что необработанное поле экспоненты смещено. На самом деле он содержит e +1023. Сдвиг дает нам e / 2+511½, причем эта дробь смещается за пределы поля. Но мы с этим разберемся. В частности, смещение необходимо сбросить до 1023, но это встроено в константу 5FE80000₁₆ .

Затем изменяет это на приближение обратного квадратного корня, просто потому, что отрицание показателя степени принимает обратное значение. Однако значимость все еще остается проблемой.

Затем следует кодирование этой очень грубой аппроксимации обратного квадратного корня с добавлением 5FE80000₁₆ .

Напомним, что кодирование экспоненты должно иметь смещение 1023 от экспоненты, но мы разрезали его пополам и инвертировали, поэтому имеет смещение -511½, которое необходимо скорректировать до +1023, поэтому нам нужно добавить 1534½. Это должно войти в поле экспоненты, которое находится в битах с 62 по 52 в 64-битной кодировке, то есть в битах с 30 по 20 в этой 32-битной части от старшей половины. Итак, нам нужно добавить 1534½, начиная с бита 20, или 1534½•2 ²⁰ = 1 609 039 872 = 5FE80000₁₆ .

Та да! Биты экспоненты являются точно битами экспоненты 2 ^{− q} . А его мантиссы — это сдвинутый бит экспоненты плюс исходные мантиссы, сдвинутые вправо на одну позицию (одна из которых смещена и потеряна). (За исключением того, что смещенный бит экспоненты был инвертирован — к нему был добавлен этот бит в 8 _16. )

Теперь у нас есть правильная экспонента, но биты мантиссы в битах с 19 по 0 дурацкие. Они содержат один бит из экспоненты и несколько инвертированных и сдвинутых битов из кодировки мантиссы. Чтобы получить аппроксимацию обратного квадратного корня, нам пришлось бы их исправить. Мы могли бы взять сдвинутые биты и вычислить исходное значение мантиссы (с точностью до 19 битов, которые мы получили от него) и умножить на 2 ^r , и это дало бы нам часть исходного числа, которая не учитывалась. ибо в показателе, который мы подготовили до сих пор. Затем мы берем обратный квадратный корень из этой части и подставляем его вместо мантиссы, и все готово.

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

извлекает биты с 19 по 14, оставляя их в битах с 5 по 0. Это смещенный бит из экспоненты и первые пять битов мантиссы. Таким образом, это самые значащие биты исходного числа, которые еще не учтены в вычисленном нами поле экспоненты.

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

В материале, на который вы ссылаетесь, не говорится, как была рассчитана таблица. Поскольку каждая запись таблицы индексируется только 6 битами, она будет использоваться для многих значения . Запись могла быть вычислена путем вычисления обратного квадратного корня для наименьшего значения этой записи, которое будет использоваться для, обратного квадратного корня для наибольшего значения эта запись будет использоваться и принимать некоторое значение между этими двумя. Затем запись в таблице устанавливается так, чтобы при ее вычитании из , он отменяет шесть битов, используемых для индексации таблицы, а затем добавляет биты, которые кодируют целевое значение. (Подобно корректировке экспоненты, добавленной 511½, а также 1023, мы хотим, чтобы запись в таблице вычитала нежелательные и добавьте нужные биты.)

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

Обратите внимание, что код является неполным и не должен использоваться как обратный квадратный корень без дополнительной работы или предосторожности. Я ожидаю, что он неправильно обрабатывает субнормаль, бесконечность или NaN.

Добавка

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

      uint32_t T2_NJ[64] =
    {
        0x01289, 0x02efb, 0x04d6a, 0x06b05, 0x087c3, 0x0a39b, 0x0be82, 0x0d86e,
        0x0f153, 0x10927, 0x11fdb, 0x13563, 0x149af, 0x15cb1, 0x16e57, 0x17e91,
        0x18d4a, 0x19a6e, 0x1a5e7, 0x1af9d, 0x1b777, 0x1bd59, 0x1c123, 0x1c2b7,
        0x1c1ef, 0x1bea5, 0x1b8ae, 0x1afdc, 0x1a3fc, 0x194d5, 0x18229, 0x16bb4,
        0x16874, 0x17a1c, 0x18aa4, 0x19a00, 0x1a824, 0x1b501, 0x1c08b, 0x1cab1,
        0x1d365, 0x1da95, 0x1e02e, 0x1e41f, 0x1e652, 0x1e6b1, 0x1e525, 0x1e194,
        0x1dbe4, 0x1d3f8, 0x1c9b0, 0x1bcea, 0x1ad82, 0x19b51, 0x1862c, 0x16de4,
        0x15248, 0x1331f, 0x1102e, 0x0e933, 0x0bde5, 0x08df6, 0x0590c, 0x01ec8
    };