Войти
ПрограммированиеФорумОбщее

Трюки с float: быстрое вычисление логарифмов. (комментарии) (6 стр)

Страницы: 15 6 7 811 Следующая »
#75
7:10, 12 янв 2022

FordPerfect, вы реализовали поддержку денормализованных чисел? Это хорошо, и сам код смотрится эффектно. Дело в том, что для меня запись на C++ выглядит непривычно, я местами с трудом читаю ваш код, например смысл первого условия от меня ускользает. Для меня более понятен машинный код, особенно если там толковые комментарии.

Насчёт выбора коэффициентов. Как я понял вы меряете ошибку в единицах младшего разряда мантиссы? Дело в том, один-два младших разряда обычно содержат шум, вызванный ошибками округления при вычислениях, поэтому при выборе коэффициентов лучше пользоваться теоретической погрешностью (в предположении, что нет ошибок округления). Для достаточно гладкой и примерно постоянной функции сходимость аппроксимации во всех метриках почти одинаковая, поэтому мне, например, нравится среднеквадратическая относительная погрешность (для этого минимизируется соответствующий интеграл на промежутке аппроксимации). Она хороша тем, что после, например, умножения на x вычисленного значения многочлена относительная погрешность мало меняется. Грубо говоря, после умножения (деления) количество верных значащих цифр остаётся прежним.

#76
10:31, 12 янв 2022

Aenigma
> вы реализовали поддержку денормализованных чисел?
Да, раз они и так под if(нештатный_аргумент), так почему бы и не? Хвостовая рекурсия, поэтому надеюсь, что компилятор переделает её тупо на goto.

> например смысл первого условия от меня ускользает
Это (bits>>23)-1u>253u ?
Выражение (bits>>23) выделяет старшие 9 бит (знак + экспонента). Нормализованные положительные аргументы - в диапазоне 1..254 (0 - денормалы и ноль, 255 - бесконечность и NaN, 256..511 - отрицательные, включая отрицательный ноль). Соответственно, вычитание 1 переводит этот диапазон в 0..253 (числа и вычитание беззнаковые, поэтому (0u-1u)=232-1, т. е. больше 253u), который дальше и проверяется. Таким образом на верхнем уровне остаётся одна проверка а уже внутри (где производительность не очень важна) дальше разбираются отдельные случаи.

> Дело в том, один-два младших разряда обычно содержат шум, вызванный ошибками округления при вычислениях, поэтому при выборе коэффициентов лучше пользоваться теоретической погрешностью (в предположении, что нет ошибок округления).
Да как бы не совсем. Интересует именно точность реальной функции. Как минимум, пошевелить коэффициенты, чтобы они были точно представимыми в плавающей точке - вполне осмысленно.
Об этом пишет, например, Robin Green в свой достаточно известной статье:
https://basesandframes.files.wordpress.com/2016/05/fast-math-functions_p1.pdf
https://basesandframes.files.wordpress.com/2016/05/fast-math-functions_p2.pdf (Optimizing For Floating Point, стр. 9)

#77
10:44, 12 янв 2022

> (bits>>23)-1u>253u
Его можно заменить на lookup table на 512 значений, но не факт, что будет быстрее.

> нравится среднеквадратическая относительная погрешность
Среднеквадратичную погрешность бывает удобно считать аналитически.
Я так sin считал https://gamedev.ru/code/forum/?id=204567&page=6&m=3958963#m75 .

#78
16:02, 12 янв 2022

FordPerfect

Как минимум, пошевелить коэффициенты, чтобы они были точно представимыми в плавающей точке - вполне осмысленно.

В случае с одинарной точностью я соглашусь, что имеет смысл "пошевелить" коэффициенты, потому что вариантов не слишком много, можно наткнуться на хороший. А вот с двойной точностью я бы уже не стал так делать - определял бы теоретический оптимум. А чтобы коэффициенты были точно представлены, как вариант, "пошевелить" можно, варьируя критерий оптимизации. Хотя в случае аппроксимации дробно-рациональной функцией, как в методе Паде-Чебышёва, коэффициенты можно "пошевелить" другим способом: домножать числитель и знаменатель на одно и то же число, которое подбирается. Я, кстати, так делаю. В этом плане удобно.

Среднеквадратичную погрешность бывает удобно считать аналитически.

Да, хотя там тоже могут встретиться неберущиеся интегралы, и СЛАУ надо решить, т.е.придётся прибегать за помощью к вычислительным системам типа WolframAlfa. Кстати, определение коэффициентов разложений Чебышёва и Паде-Чебышёва делается ненамного сложнее.
Я имел в виду то, что оптимальность аппроксимации в смысле среквадратической относительной погрешности инвариантна относительно домножения на точную функцию. Это хорошо, потому что значения аппроксимирующего многочлена часто ещё домножаются на x или что-то ещё. Оптимальность в этом случае сохраняется. Хотя согласен, что это не столь важно.

#79
19:02, 12 янв 2022

> Хвостовая рекурсия, поэтому надеюсь, что компилятор переделает её тупо на goto.
Гоню. И это может помешать разворачиванию и инлайну.

#80
4:58, 13 янв 2022

Я точность плохо мерял.
Я сравнивал с библиотечной log2f (которая, вроде, даже не correctly rounded, а около 0.51 ULP), хотя осмысленнее сравнивать с точным значением, вместо которого в данном случае сойдёт взять double версию.
https://rextester.com/KEFMA90078

  Linf: abs=1.11e-07, rel=2.46e-07, 3.479 ulp
  RMS : abs= 2.1e-08, rel=4.72e-08, 0.581 ulp
#81
5:15, 13 янв 2022

> И это может помешать разворачиванию и инлайну.
Похоже да:
https://godbolt.org/z/7onrxvazj

#82
(Правка: 7:36) 7:29, 13 янв 2022

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

Я точность плохо мерял.

Да, двойной точности для анализа погрешности достаточно. Я для этих целей, когда строил аппроксимацию с двойной точностью, написал собственную простенькую функцию с использованием x87-команд. Они обеспечивают расширенную точность, которая сойдёт за точное значение при анализе точности double.

#83
(Правка: 8:04) 8:02, 13 янв 2022

Aenigma
> Может быть заранее формировать целую часть логарифма, которая в случае штатной ветки равна 0, а для денормализованных чисел -32, а потом две ветки сходятся.
Кажись так инлайнит:

float log2_v5(float x)
{
    float y=-127.0f;
    uint32_t bits,bias=0x00400000u;
    memcpy(&bits,&x,sizeof(float));
    if((bits>>23)-1u>253u)
    {
        if(!(x>0.0f)) return (x==0.0f?-INFINITY:nanf(""));
        if(x==+INFINITY) return +INFINITY;
        x*=4294967296.0f;
        memcpy(&bits,&x,sizeof(float));
        y-=32.0f;
    }
    bits+=bias;
    y+=float(int32_t(bits>>23));
    bits=((bits&0x007FFFFFu)|0x3F800000u)-bias;
    memcpy(&x,&bits,sizeof(float));
    x=(x-1.0f)/(x+1.0f);
    float x2=x*x;
    return y+x*(2.88539003e+00f
           +x2*(9.61808794e-01f
           +x2*(5.76110615e-01f
           +x2*(4.42826415e-01f))));
}

https://godbolt.org/z/3bsPoErW9

#84
(Правка: 13:27) 13:23, 13 янв 2022

FordPerfect, да, получилось очень неплохо. Компилятор в целом очень хорошо всё перевёл в машинный код, но некоторые моменты удивляют:

subl    $1, %edx

У меня даже Visual Studio в этом месте поставил команду dec (хотя в остальном скомпилировал плохо). Одно хорошо, что код на C++, даже если он плохо откомпилировался, будет работать почти на любой машине в любой системе, а код на ассемблере не везде можно вставить.

#85
(Правка: 19:09) 19:09, 13 янв 2022

Aenigma
Насколько я понимаю, dec меньше, но sub иногда быстрее.

The INC and DEC instructions modify only a subset of the bits in the flag register. This creates a dependence on all previous writes of the flag register. This is especially problematic when these instructions are on the critical path because they are used to change an address for a load on which many other instructions depend. Assembly/Compiler Coding Rule 33. (M impact, H generality) INC and DEC instructions should be replaced with ADD or SUB instructions, because ADD and SUB overwrite all flags, whereas INC and DEC do not, therefore creating false dependencies on earlier instructions that set the flags.

(см. https://stackoverflow.com/questions/19890157/gcc-doesnt-make-use-of-inc )

#86
(Правка: 21:43) 21:41, 13 янв 2022

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

#87
0:08, 14 янв 2022

Aenigma
Скорее всего практически заметной разницы нет. Но
> В данном случае следующая команда сравнения меняет все флаги, и тем самым разрывает зависимости по флагам.
вроде же важнее то, что результат (формально) зависит от предыдущей (флагоменяющей) команды (а "следующая команда" внутри if, куда, обычно, всё равно не заходит).

> На более старых процессорах dec была однозначно быстрее, потому что процессоры были не такие умные.
Это да. Впрочем, если проставить

-m32 -mtune=i386 -mavx2 -mfma -mfpmath=sse -msseregparm -O3

то будет dec: https://godbolt.org/z/jq8nsbbKv (но это довольно забавная комбинация опций).

#88
8:37, 14 янв 2022

Aenigma
> Одно хорошо, что код на C++, даже если он плохо откомпилировался, будет работать почти на любой машине в любой системе
Чтобы сломать мой код, достаточно плавучки, отличной от IEEE (напр. https://en.wikipedia.org/wiki/IBM_hexadecimal_floating-point ), или различающейся endianness у целых и плавучки (вроде встречается где-то).
Относительно редкие случаи, да.

#89
(Правка: 11:09) 9:05, 14 янв 2022

FordPerfect

вроде же важнее то, что результат (формально) зависит от предыдущей (флагоменяющей) команды (а "следующая команда" внутри if, куда, обычно, всё равно не заходит).

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

Выкладываю AVX-код для вычисления логарифма с одинарной точностью с поддержкой денормализованных чисел. Заодно немного "подкрутил" производительность по сравнению с предыдущим простым вариантом.

#define salc _emit 0xD6              // Определяем команду salc её машинным кодом
_declspec(naked) float _vectorcall ln1(float x)
{
  static const float ct[6] =         // Таблица констант
  {
    1.0f,                            // 1
    0.576395561f,                    // 5^5*b2
    0.961802196f,                    // 5^3*b1
    2.88539007f,                     // 5*b0
    0.439077714f,                    // 5^7*b3
    0.693147181f                     // ln 2
  };
  _asm
  {
    mov edx,-127                     // В edx - смещение порядка для нормализованных чисел float
    vmovd eax,xmm0                   // Прочитать в eax двоичное представление числа x
    cmp eax,0x7F800000               // Сравнить число x со значением Inf
    jge ln_exit                      // Если x=+Inf или +NaN, вернуть аргумент x без изменений
    ror eax,23                       // Циклически сдвинуть число так, чтобы порядок оказался в al
    movzx ecx, al                    // Получить в ecx порядок числа со смещением +127
    jecxz ln_denorm                  // Перейти, если x=0 или x денормализованное (т.е. близко к 0)
      ln_calc:                       // Точка входа для вычисления логарифма после нормализации числа 
    salc                             // al=-1, если установлен старший бит явной части мантиссы; иначе al=0
    sub cl,al                        // ecx=k+127, где k - добавка к двоичному логарифму (целая часть)
    add al,127                       // Изменить порядок x так, чтобы привести его к диапазону 0.75<=x<1.5
    add ecx,edx                      // ecx=k
    ror eax,9                        // Получить в eax двоичное представление числа и проверить его знак
    jnc ln_cont                      // Перейти, если x>0, т.е. аргумент корректный. Продолжить, если x<0
      ln_nan:                        // Точка входа для формирования результата NaN
    sbb eax,eax                      // eax=-1 (значение NaN) в случае x<0
      ln_inf:                        // Точка входа для формирования результата -Inf в случае x=0
    or eax,0xFF800000                // Сформировать в eax значение -Inf, если x=0; или NaN, если x<0
    vmovd xmm0,eax                   // Перенести результат из eax в xmm0
      ln_exit:                       // В xmm0 готов результат для нештатного значения x
    ret                              // Вернуться
      ln_denorm:                     // Точка входа для обработки денормализованного x, в т.ч. x=0
    btr eax,8                        // Прочитать в cf и сбросить знаковый бит числа
    bsr ecx,eax                      // Найти номер старшего единичного бита; установить zf, если x=0
    jz ln_inf                        // Если x=0, перейти для формирования результата -Inf
    jc ln_nan                        // Если x<0 (или -NaN), перейти для формирования результата NaN
    lea edx,[edx+2*ecx-31]           // Обновить смещение порядка для денормализованного числа
    neg ecx                          // 32-ecx - это необходимое количество сдвигов мантиссы влево
    shl eax,cl                       // Сформировать мантиссу нормализованного числа
    bt eax,31                        // Проверить старший бит явной части мантиссы
    jmp ln_calc                      // Число нормализовано. Перейти для вычисления логарифма
      ln_cont:                       // Продолжаем вычисления в штатном случае, когда x>0
    vmovd xmm0, eax                  // Перенести из eax в xmm0 результат или приведённое значение x
    mov eax, offset ct               // В eax адрес таблицы констант
    vaddss xmm1, xmm0, [eax]         // xmm1=x+1
    vsubss xmm0, xmm0, [eax]         // xmm0=x-1
    vcvtsi2ss xmm3, xmm3, ecx        // xmm3=k, это добавка к двоичному логарифму
    vdivss xmm0, xmm0, xmm1          // xmm0=(x-1)/(x+1)=t/5
    vmovss xmm1, [eax+16]            // Инициализировать сумму старшим коэффициентом: xmm1=5^7*b3
    vmulss xmm2, xmm0, xmm0          // xmm2=t^2/25 - подготовить аргумент для вычисления многочлена
    vfmadd213ss xmm1, xmm2, [eax+4]  // Выполнить первый шаг вычисления многочлена по схеме Горнера
    vfmadd213ss xmm1, xmm2, [eax+8]  // Выполнить второй шаг вычисления многочлена по схеме Горнера
    vfmadd213ss xmm1, xmm2, [eax+12] // Выполнить третий шаг вычисления многочлена по схеме Горнера
    vfmadd213ss xmm0, xmm1, xmm3     // xmm0=t*f(t)+k - получили готовый двоичный логарифм
    vmulss xmm0, xmm0, [eax+20]      // Перевести двоичный логарифм в натуральный
    ret                              // Вернуться
  }
}
Страницы: 15 6 7 811 Следующая »
ПрограммированиеФорумОбщее