Снижение прочности при делении с плавающей запятой вручную

Question

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

Наши профи дали нам следующую подсказку:

Подсказка: узнайте, как числа с плавающей запятой одинарной точности IEEE 754 представлены в памяти.

Вот фрагмент кода: (a имеет типdouble*)

      for (int i = 0; i < N; ++i) {
    a[i] += i / 5.3;   
}

Сначала я попытался изучить вывод компилятора для этого, снятого на godbolt. Я попытался скомпилировать его без какой-либо оптимизации: (примечание: я скопировал только соответствующую часть цикла for)

          mov     eax, DWORD PTR [rbp-4]
    cdqe
    lea     rdx, [0+rax*8]
    mov     rax, QWORD PTR [rbp-16]
    add     rax, rdx
    movsd   xmm1, QWORD PTR [rax]
    cvtsi2sd        xmm0, DWORD PTR [rbp-4]    //division relevant
    movsd   xmm2, QWORD PTR .LC0[rip]          //division relevant
    divsd   xmm0, xmm2                         //division relevant
    mov     eax, DWORD PTR [rbp-4]
    cdqe
    lea     rdx, [0+rax*8]
    mov     rax, QWORD PTR [rbp-16]
    add     rax, rdx
    addsd   xmm0, xmm1
    movsd   QWORD PTR [rax], xmm0

и с-O3:

      .L2:
        pshufd  xmm0, xmm2, 238             //division relevant
        cvtdq2pd        xmm1, xmm2          //division relevant
        movupd  xmm6, XMMWORD PTR [rax]
        add     rax, 32
        cvtdq2pd        xmm0, xmm0          //division relevant
        divpd   xmm1, xmm3                  //division relevant
        movupd  xmm5, XMMWORD PTR [rax-16]
        paddd   xmm2, xmm4
        divpd   xmm0, xmm3                  //division relevant
        addpd   xmm1, xmm6
        movups  XMMWORD PTR [rax-32], xmm1
        addpd   xmm0, xmm5
        movups  XMMWORD PTR [rax-16], xmm0
        cmp     rax, rbp
        jne     .L2

Я прокомментировал часть деления ассемблерного кода. Но этот вывод не помогает мне понять, как применить уменьшение прочности фрагмента. (Возможно, слишком много оптимизаций, чтобы полностью понять вывод)

Во-вторых, я попытался понять битовое представление части с плавающей запятой.5.3. Который:

      0   |10000001|01010011001100110011010
Sign|Exponent|Mantissa

Но и это мне не помогает.

Answer 1

ПРЕДОСТЕРЕЖЕНИЕ: через несколько дней я понял, что этот ответ неверен, поскольку он игнорирует последствия потери значимости (до субнормального или до нуля) при вычислении . В этом случае умножение результата на степень двойки является «точным», но не дает такого результата, как деление большего целого числа на было бы.

необходимо вычислять только для нечетных значений . Для четных значений, вы можете просто умножить на 2,0 значение (i/2)/5,3, которое уже было вычислено ранее в цикле.

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

Один из способов добиться этого состоит в том, чтобы перебирать все нечетные числа, меньшие , и после вычислениядля обработки индекса, чтобы также обрабатывать все индексы формы.

      if (N > 0) {
  a[0] += 0.0; // this is needed for strict IEEE 754 compliance lol
  for (int o = 1; o < N; o += 2) {
    double d = o / 5.3;
    int i = o;
    do { 
      a[i] += d;
      i += i;
      d += d;
    } while (i < N);
  }
}

Примечание: здесь не используется предоставленная подсказка «Запросить, как числа с плавающей запятой одинарной точности IEEE 754 представлены в памяти». Я думаю, что довольно хорошо знаю, как числа с плавающей запятой одинарной точности представлены в памяти, но я не понимаю, насколько это актуально, тем более что в коде нет значений одинарной точности или вычислений, которые нужно оптимизировать. Я думаю, что в способе выражения проблемы есть ошибка, но все же вышеизложенное технически является частичным ответом на поставленный вопрос.

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

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

Answer 2

Если мы примем определение Википедии , что

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

тогда мы можем применить здесь снижение силы, преобразовав дорогостоящее деление с плавающей запятой в умножение с плавающей запятой плюс два умножения-сложения с плавающей запятой (FMA). При условии, чтосопоставлен с IEEE-754, режим округления по умолчанию для вычислений с плавающей запятой — округление до ближайшего или четного, и чтоявляется 32-битным типом, мы можем доказать правильность преобразования с помощью простой исчерпывающей проверки:

      #include <stdio.h>
#include <stdlib.h>
#include <limits.h>
#include <math.h>

int main (void)
{
    const double rcp_5p3 = 1.0 / 5.3; // 0x1.826a439f656f2p-3
    int i = INT_MAX;
    do {
        double ref = i / 5.3;
        double res = fma (fma (-5.3, i * rcp_5p3, i), rcp_5p3, i * rcp_5p3);
        if (res != ref) {
            printf ("error: i=%2d  res=%23.13a  ref=%23.13a\n", i, res, ref);
            return EXIT_FAILURE;
        }
        i--;
    } while (i >= 0);
    return EXIT_SUCCESS;
}

Большинство современных экземпляров распространенных архитектур процессоров, таких как x86-64 и ARM64, имеют аппаратную поддержку FMA, поэтому ее можно напрямую сопоставить с соответствующей аппаратной инструкцией. Это должно быть подтверждено просмотром разборки сгенерированного двоичного файла. Там, где аппаратная поддержка FMA отсутствует, преобразование, очевидно, не должно применяться, так как программные реализацииработают медленно и иногда функционально некорректно.

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

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

В комментариях user139746 указывает, что существует альтернативный алгоритм, потенциально снижающий силу деления с плавающей запятой с помощью постоянного делителя времени компиляции, путем предварительного вычисления обратной величины делителя с большей точностью, чем исходная, и, в частности, как двойной двойной. Для получения дополнительной информации см. Н. Бризбарр и Ж.-М. Мюллер, «Правильно округленное умножение на константу произвольной точности», IEEE Transactions on Computers, 57(2): 162-174, февраль 2008 г., в котором также содержится руководство по определению того, безопасно ли это преобразование для какой-либо конкретной константы. Поскольку данный случай прост, я снова использовал исчерпывающий тест, чтобы показать, что он безопасен. Там, где это применимо, это уменьшит деление до одного FMA плюс умножение:

      #include <stdio.h>
#include <stdlib.h>
#include <limits.h>
#include <mathimf.h>

int main (void)
{
    const double rcp_5p3_hi =  1.8867924528301888e-1; //  0x1.826a439f656f2p-3
    const double rcp_5p3_lo = -7.2921377017921457e-18;// -0x1.0d084b1883f6e0p-57
    int i = INT_MAX;
    do {
        double ref = i / 5.3;
        double res = fma (i, rcp_5p3_hi, i * rcp_5p3_lo);
        if (res != ref) {
            printf ("i=%2d  res=%23.13a  ref=%23.13a\n", i, res, ref);
            return EXIT_FAILURE;
        }
        i--;
    } while (i >= 0);
    return EXIT_SUCCESS;
}