How to become an author
Search
Write a publication
Pull to refresh

Comments 83

Приведенный код содержит UB. Вот так делать правильнее:


inline int64_t to_int64(double x) {
    int64_t a;
    memcpy(&a, &x, 8);
    uint64_t mask = (uint64_t)(a >> 63) >> 1;
    return a ^ mask;
}
Total votes 11: ↑11 and ↓0+11

Правильный ответ "по кочану! в стандарте так сказано!".
На самом же деле, есть паттерны UB, которые эксплуатируются в благих целях, и на которые компилятор смотрит благосклонно.
Например, чисто сишный приём преобразования типов через union.


Почему нельзя reinterpret_cast'ить типы?


1) Потому что выравнивание. Некоторые процессоры болезненно чувствительны к попытке прочитать невыравненные double и даже int.
Здесь неприменимо, так как у приёмника int64 выравнивание уж точно не сильнее, чем у источника double.


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


3) Потому что компилятор — психопат с мачете. Как увидит явное UB, так вместо того, чтобы писать вонинг, вставит код форматирования диска.
Покажите мне такие компиляторы.


4) Потому что психопат с мачете — ревьювер.
Ну, тут уж медицина бессильна.

Total votes 16: ↑12 and ↓4+8
Склонен согласиться. На фоне данного трюка использующего особенности формата хранения fp64 в памяти, данный UB совсем безобиден.
В частности одна из причин UB — неоднозначность формата хранения, вспомним 8087 с его 80-битными флоатами, даже если в коде они определены как FP64…
Примерно как соблюдение правил движения, когда убегаешь от полиции после ограбления
Total votes 6: ↑4 and ↓2+2
вспомним 8087 с его 80-битными флоатами
Извините, то никаких 80-битных флоатов 8087 в памяти не размещает. И с ним указанный трюк (при использовании memcpy) будет отлично работать.

Примерно как соблюдение правил движения, когда убегаешь от полиции после ограбления
Я принёс в банк чек с подписями и печатями и мне выдали денег! Я ограбил банк! Срочно убегать от полиции!

P.S. Чисто для тех, кто в танке: в стандарте существует такая вещь, как is_iec559 — именно для того, чтобы всё описанное в статетье не было “ограблением банка”.
This comment wasn’t rated yet0
8087 может работать с 80-битными флоатами в памяти и это даже будет is_iec559 тип. Но правда 754 вроде требует сейчас от компилятора гарантировать что double это именно 64-битный тип и ни один компилятор без явной просьбы использовать long double подобный код не генерировал
Total votes 2: ↑2 and ↓0+2
8087 может работать с 80-битными флоатами в памяти и это даже будет is_iec559 тип.
Это закрывается банальным static_assert(sizeof(double) == sizeof(std::uit64_t));.

ни один компилятор без явной просьбы использовать long double подобный код не генерировал
Это тоже правда, но это даже не главное. Главное, что вы можете всё это проверить во время компиляции… а вот UB — это таки UB.
Total votes 1: ↑1 and ↓0+1
Например, чисто сишный приём преобразования типов через union.
Например, чисто сишный приём преобразования типов через union.
Только это не «сишный приём». Стандарт это запрещает. gcc поддерживает, хотя и с оговорками, разработчики clang говорят что-то вроде “простейшие случаи такого рода мы отлавливаем и поддерживаем, но никаких попыток гарантировать его работоспособность во всех случаях мы не делаем”.

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

Что самое смешное — что для этого даже не нужны какие-то козни с компилятором. Например математические сопроцессоры (8087, Weitek и другие) работают параллельно основному процессору, но независимо от него. Соотвественно простейшая конструкция:
  union {
    double d;
    uint64_t i64;
  };
  d = sin(d);
  return i64;
без специальных усилий со стороны компилятора работать не будет: вызов sin запусит вычисление синуса (занимающее под сотню тактов) на сопроцессоре, после чего основной процессор прочитает то, что там было в памяти — и начнёт радостно использовать в вычислениях.

Как увидит явное UB, так вместо того, чтобы писать вонинг, вставит код форматирования диска.
Приводили много раз.

Конкретно с этой функцией вроде компиляторы ничего плохого не делают ещё. Но имеют право. Вот реальный пример (что самое смешное: dtoa как раз использует трюки, подобные описанному в статье… только код там чуток посложнее).
Total votes 5: ↑4 and ↓1+3
Стандарт это запрещает.

Сишный — разрешает. stackoverflow.com/questions/25664848/unions-and-type-punning

Конкретно с этой функцией вроде компиляторы ничего плохого не делают ещё. Но имеют право.

Имеют, но в обозримом будущем не будут — именно потому, что такой код используется сплошь и рядом. Обратите внимание, что в случае с dtoa эту (правомочную!) оптимизацию сочли багом в clang и в конечном счёте убрали.
This comment wasn’t rated yet0
Имеют, но в обозримом будущем не будут — именно потому, что такой код используется сплошь и рядом.
Те, кто “используют сплошь и рядом” могут включить -fno-strict-aliasing и избавиться от UB.

Сишный — разрешает.
Согласен. Я с C стараюсь не связываться, так что некоторые ньюансы не знаю.

Но забавно что как в силу специфики реализации железа как раз в C с этим подходом и возникают проблемы на старом железе (современные процессоры от концепции сопроцессоров давно отказались: либо плавучка встроена в основной процессор, либо её нет вообще и она эмулируется).
Total votes 2: ↑1 and ↓10
Ну memcpy то целиком на 8 байт вызывать не стоит, наверное… Мы же об эффективности говорим?
This comment wasn’t rated yet0

А никакого "вызова целиком" там и не происходит, всё memcpy упрощается до двух инструкций mov: https://godbolt.org/z/x83Mez

Total votes 2: ↑2 and ↓0+2
Если добавить -O1, то не будет и них.
Total votes 1: ↑1 and ↓0+1
Спасибо, покопался. -O3 дает возможность всем компиляторам убить memcpy() кроме… msvc. Может там еще чего надо, я хз, не умею в него.

godbolt.org/z/GGnxr1
This comment wasn’t rated yet0
У MSVC нужный ключ называется /O2
Total votes 1: ↑1 and ↓0+1

А у него нет опции O3, только O1 и O2...

Total votes 3: ↑3 and ↓0+3
И даже больше скажу я ловил креши из-за этого UB. На старых arm процессорах которые не умеют читать по не выровненному адресу…
Так что да memcpy оно более безопасно. В релизе (-O3) оно все равно упрощается до нескольких инструкций (если это безопасно)
Total votes 1: ↑1 and ↓0+1
Каст double к int64_t не может вызвать проблем с выравниванием.
Total votes 3: ↑2 and ↓1+1
Пример: читаем данные из файла. Сделали mmap. Адрес не выровнен и тип указателя из которого читаем вполне может быть int64_t…
Total votes 2: ↑1 and ↓10
Если бы адрес не был выровнен, то по этому адресу и double не мог бы лежать.
Total votes 2: ↑1 and ↓10
Что вам запрещает записать double в бинарный файл со смещением 123?
Total votes 1: ↑1 and ↓0+1
Ну вы ещё скажите, что void foo(double *bar) {...} не может разыменовывать *bar, а может только memcpy его в локальную переменную, потому что мало ли, вдруг bar прочитан из файла.
Total votes 4: ↑2 and ↓20
В какой части кода ошибка: в функции foo(double *bar) или в функции что ее вызывает?
И да в некоторых случаях foo(double *bar) может ожидать что адрес не выравненный.
Пример:
template<class Type, class Ptr>
Type readUnaligned(Ptr *ptr) {
Type rv;
memcpy(&rv, ptr, sizeof(Type));
return run
}
Total votes 2: ↑0 and ↓2-2
Об этом и речь: если после каста double к int64_t возникли проблемы с выравниванием, значит они и без каста бы были при обращении к исходному double. Значит, эти проблемы вызваны вовсе не кастом double к int64_t.
Total votes 6: ↑4 and ↓2+2
Об этом и было мое сообщение. В таких случаях безопаснее использовать memcpy а не каст. Оно решает проблему выравнивания.
Total votes 5: ↑1 and ↓4-3
Стандарт явно не описывает размеры float, double, long double. Есть только определение что точность long double выше double и точность double выше float.
И в совсем общем случае каст в int любого размера (как и memcpy) даст неопределенное поведение.
Поэтому «Каст double к int64_t не может вызвать проблем с выравниванием» справедливо только в частном случае когда размер double равен 8 байт.
Total votes 3: ↑3 and ↓0+3

В стандарте есть std::numeric_limits<double>::is_iec559 (и static_assert) для тех, кто хочет убедиться, что тип данных double соответствует IEEE 754.

Total votes 5: ↑5 and ↓0+5
Вы не путаете reinterpret_cast и bit_cast?

Если “играться с указателями” как в статье, то как раз reinterpret_cast отлично скомпилируется (и натворит делов).

Ипользовать нужно только bit_cast. Если у вас не C++20 компилятор, то можете своровать определение с cppreference (только в std его не суйте, а то конфликт будет после обновления компилятора).
Total votes 2: ↑2 and ↓0+2
Я вообще имел в виду reinterpret_cast с использованием ссылочных типов

int64_t a = reinterpret_cast<int64_t&>(x);


Но да, оно по факту работает так же как с указателями, я ошибся )
This comment wasn’t rated yet0

Забавный трюк. А есть сравнение производительности?
Есть подозрение, что разработчики АЛУ в железе сделали примерно то же самое.

This comment wasn’t rated yet0
Разработчики АЛУ в железе не могли обойтись без обработки NaN.
Total votes 3: ↑3 and ↓0+3

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

Total votes 2: ↑2 and ↓0+2

Вот только вычислительный код, когда натыкается на NaN, почему-то начинает работать на пару порядков медленнее.

Total votes 1: ↑1 and ↓0+1
Как корректнее делать инверсию знака: a = -a или a = 0 — a?
Корректна ли проверка перед делением на a != 0?
This comment wasn’t rated yet0
  1. результаты будут разными: если a==0., то -a вернёт отрицательный ноль, 0-a — положительный
  2. бесконечный результат может получаться и в других случаях, например 1e200/1e-200
This comment wasn’t rated yet0
а что происходит при проверке a != 0 и делении в случаях когда а равно +0 и -0?
This comment wasn’t rated yet0
Два ноля равны (и это единственная пара double, которые не совпадают побитово, но равны). При делении на ноль получается положительная или отрицательная бесконечность, в соответствии со XOR знаков делимого и ноля.
This comment wasn’t rated yet0

В ARM AArch64 есть инструкции специально для этого: FCMP и FCM*, которые исполняются за 2-3 такта: https://godbolt.org/z/cr4Mcq
На x86_64 — UCOMISD.

This comment wasn’t rated yet0
Преимущество этого трюка перед FCMP в том, что FPU остаются свободными, и могут параллельно со сравнением вычислять что-нибудь ещё.
Total votes 5: ↑0 and ↓5-5

Сомнительное преимущество:


  1. FPU большую часть времени простаивают, так как количество целочисленных инструкций заведомо больше. Если только это не в ручную написанный код.
  2. Пересылка данных между целочисленными регистрами и FP регистрами не бесплатно, на ARM'ах 2-3 такта.

Есть у вас результаты бенчмарок, показывающие преимущество вашего подхода?

Total votes 6: ↑6 and ↓0+6

Современные процессоры суперскалярны. Не знаю как в ARM, но в Intel и AMD имеется минимум по 2 FPU-модуля на ядро, отвечающих за арифметику. Пока один модуль осуществляет сравнение, второй — свободен.

Total votes 1: ↑1 and ↓0+1
Теперь ждём трюков как складывать длинные целые числа.
This comment wasn’t rated yet0
Это может быть действительно полезно, если мы работаем на архитектуре, где нет аппаратной поддержки плавающей точки, и нет нужды работать ненормальными числами. Но на современном компьютере… Не верю!
Смотрим во что компилируется обычное наивное сравнение: godbolt.org/z/s3hMjh
Видим инструкцию ucomisd xmm0, xmm1 на x86.
Читаем здесь: mudongliang.github.io/x86/html/file_module_x86_id_316.html
Latency не более 7 тактов.
Вы думаете это медленнее чем остальная портянка с зависимостями по данным будет?
Total votes 2: ↑2 and ↓0+2
Latency не более 7 тактов.
Вы думаете это медленнее чем остальная портянка с зависимостями по данным будет?

А почему нет? SAR+SHR+XOR+CMP — четыре простейшие операции, каждая выполняется за такт.
Total votes 1: ↑0 and ↓1-1
Ну вот и докажите бенчмарками.
Вы, кстати, лукавите насчёт 4 операций, потому что это надо сделать два раза, с каждым аргументом.
Total votes 3: ↑3 and ↓0+3
Вы, кстати, лукавите насчёт 4 операций, потому что это надо сделать два раза, с каждым аргументом.

Да, но между аргументами нет зависимостей по данным.

Ну вот и докажите бенчмарками.

Пожалуйте 
#include <stdbool.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>

inline int64_t to_int64(double x) {
	int64_t a = *(int64_t*)&x;
	uint64_t mask = (uint64_t)(a >> 63) >> 1;
	return a ^ mask;
}

inline bool is_smaller(double x1, double x2) {
	return to_int64(x1) < to_int64(x2);
}

int main() {
  long const size = 3*128*1024*1024;
  double *gig = malloc(size*sizeof(double));
  for (long i=0; i<size; i++)
    do {
      for (int j=0; j<8; j++)
        ((char*)(gig+i))[j] = rand();
    } while (gig[i]!=gig[i]);
  printf("3GB random doubles generated\n");

  clock_t begin = clock();
  double min = 1e200;
  double next = 1;
  for (long i=0; i<size; i++) {
    next = 1./(next+1);
    if (gig[i] < min)
      min = gig[i];
  }
  clock_t end = clock();
  printf("next: %f min: %f time: %f\n", next, min, (double)(end - begin) / CLOCKS_PER_SEC);

  begin = clock();
  min = 1e200;
  next = 1;
  for (long i=0; i<size; i++) {
    next = 1./(next+1);
    if (is_smaller(gig[i], min))
      min = gig[i];
  }
  end = clock();
  printf("next: %f min: %f time: %f\n", next, min, (double)(end - begin) / CLOCKS_PER_SEC);
}

GCC 7.5.0 -O3, E5-2643 @ 3.50GHz: 1.942552 против 1.926676
Total votes 1: ↑0 and ↓1-1

А нахрена там вот эта строчка?


next = 1./(next+1);

Львиная доля времени приходится именно на неё, и к тестируемой операции она не имеет никакого отношения.


Если её убрать, то у меня получается 0.333950 (native double) против 0.554674 (to_int64). Разница в полтора раза более чем существенная. Компилятор clang 10.0.1 -O3. Процессор Core i7-2600K.


P.S. А с этой строчкой — порядка 2.54 и 2.62.


Disassembly:


    1230:   f2 41 0f 10 44 df d8    movsd  xmm0,QWORD PTR [r15+rbx*8-0x28]
    1237:   f2 0f 5d c1             minsd  xmm0,xmm1
    123b:   f2 41 0f 10 4c df e0    movsd  xmm1,QWORD PTR [r15+rbx*8-0x20]
    1242:   f2 0f 5d c8             minsd  xmm1,xmm0
    1246:   f2 41 0f 10 44 df e8    movsd  xmm0,QWORD PTR [r15+rbx*8-0x18]
    124d:   f2 0f 5d c1             minsd  xmm0,xmm1
    1251:   f2 41 0f 10 4c df f0    movsd  xmm1,QWORD PTR [r15+rbx*8-0x10]
    1258:   f2 0f 5d c8             minsd  xmm1,xmm0
    125c:   f2 41 0f 10 44 df f8    movsd  xmm0,QWORD PTR [r15+rbx*8-0x8]
    1263:   f2 0f 5d c1             minsd  xmm0,xmm1
    1267:   f2 41 0f 10 0c df       movsd  xmm1,QWORD PTR [r15+rbx*8]
    126d:   f2 0f 5d c8             minsd  xmm1,xmm0
    1271:   48 83 c3 06             add    rbx,0x6
    1275:   48 81 fb 05 00 00 18    cmp    rbx,0x18000005
    127c:   75 b2                   jne    1230 <main+0xb0>

    12d0:   48 83 c3 02             add    rbx,0x2
    12d4:   66 0f 6f d1             movdqa xmm2,xmm1
    12d8:   48 81 fb 01 00 00 18    cmp    rbx,0x18000001
    12df:   74 69                   je     134a <main+0x1ca>
    12e1:   f3 41 0f 7e 44 df f8    movq   xmm0,QWORD PTR [r15+rbx*8-0x8]
    12e8:   f3 41 0f 7e 0c df       movq   xmm1,QWORD PTR [r15+rbx*8]
    12ee:   66 48 0f 7e c0          movq   rax,xmm0
    12f3:   48 89 c1                mov    rcx,rax
    12f6:   48 c1 f9 3f             sar    rcx,0x3f
    12fa:   48 d1 e9                shr    rcx,1
    12fd:   48 31 c1                xor    rcx,rax
    1300:   66 48 0f 7e d0          movq   rax,xmm2
    1305:   48 89 c2                mov    rdx,rax
    1308:   48 c1 fa 3f             sar    rdx,0x3f
    130c:   48 d1 ea                shr    rdx,1
    130f:   48 31 c2                xor    rdx,rax
    1312:   48 39 d1                cmp    rcx,rdx
    1315:   7c 04                   jl     131b <main+0x19b>
    1317:   66 0f 6f c2             movdqa xmm0,xmm2
    131b:   66 48 0f 7e c8          movq   rax,xmm1
    1320:   48 89 c1                mov    rcx,rax
    1323:   48 c1 f9 3f             sar    rcx,0x3f
    1327:   48 d1 e9                shr    rcx,1
    132a:   48 31 c1                xor    rcx,rax
    132d:   66 48 0f 7e c0          movq   rax,xmm0
    1332:   48 89 c2                mov    rdx,rax
    1335:   48 c1 fa 3f             sar    rdx,0x3f
    1339:   48 d1 ea                shr    rdx,1
    133c:   48 31 c2                xor    rdx,rax
    133f:   48 39 d1                cmp    rcx,rdx
    1342:   7c 8c                   jl     12d0 <main+0x150>
Total votes 3: ↑3 and ↓0+3
У меня вот так
g++ -O4 -march=native -mavx -funroll-loops test.cpp -o…
3GB random doubles generated
next: 6.180340e-01 min: -1.797693e+308 time: 1.986066
next: 6.180340e-01 min: -1.797693e+308 time: 1.993061

А без мин, только деление:
3GB random doubles generated
next: 6.180340e-01 min: 1.000000e+200 time: 1.909581
next: 6.180340e-01 min: 1.000000e+200 time: 1.910505

У Вас АЛУ тупо забит делением. Остальное — практически бесплатно.
Без деления:
3GB random doubles generated
next: 1.000000e+00 min: -1.797693e+308 time: 0.266345
next: 1.000000e+00 min: -1.797693e+308 time: 0.560713

То ест выигрыш — 0.0
Total votes 1: ↑1 and ↓0+1
У Вас АЛУ тупо забит делением.

И это намеренно: вполне реалистичный юзкейс, когда внутри итерации есть и вычисления с плавающей точкой, и сравнения. ucomisd не может выполняться одновременно с вычислениями, to_int64 — может.
Total votes 1: ↑0 and ↓1-1
Так все равно без фокусов — быстрее. И в большинстве случаев — намного быстрее (0.26 vs 0.65). На моей машине Ваш трюк проигрывает всегда. Мне кажется что модуль деления(DIV) и умножения (MAD) друг от друга не зависят. И планировшик задействует их параллельно, если может.
Total votes 1: ↑1 and ↓0+1
ucomisd не может выполняться одновременно с вычислениями,

Почему не может?

Total votes 1: ↑1 and ↓0+1
Я сейчас прогнал вышезапощенный код на ещё одной машине (Xeon 5118 @ 2.30GHz, GCC 9.3.1) и получил ещё большую разницу (3.20 против 2.71). У вас есть для этого экспериментального наблюдения иное объяснение, чем «ucomisd не может выполняться одновременно с вычислениями»?
This comment wasn’t rated yet0

А можно ассемблерный код?

This comment wasn’t rated yet0
Вся main() целиком 
00000000004010a0 <main>:
  4010a0:       41 57                   push   r15
  4010a2:       bf 00 00 00 c0          mov    edi,0xc0000000
  4010a7:       41 bf 08 00 00 c0       mov    r15d,0xc0000008
  4010ad:       41 56                   push   r14
  4010af:       41 55                   push   r13
  4010b1:       41 54                   push   r12
  4010b3:       55                      push   rbp
  4010b4:       53                      push   rbx
  4010b5:       48 83 ec 18             sub    rsp,0x18
  4010b9:       e8 c2 ff ff ff          call   401080 <malloc@plt>
  4010be:       49 89 c5                mov    r13,rax
  4010c1:       49 89 c6                mov    r14,rax
  4010c4:       48 8d 68 08             lea    rbp,[rax+0x8]
  4010c8:       49 01 c7                add    r15,rax
  4010cb:       49 89 c4                mov    r12,rax
  4010ce:       66 90                   xchg   ax,ax
  4010d0:       4c 89 e3                mov    rbx,r12
  4010d3:       0f 1f 44 00 00          nop    DWORD PTR [rax+rax*1+0x0]
  4010d8:       e8 b3 ff ff ff          call   401090 <rand@plt>
  4010dd:       48 83 c3 01             add    rbx,0x1
  4010e1:       88 43 ff                mov    BYTE PTR [rbx-0x1],al
  4010e4:       48 39 eb                cmp    rbx,rbp
  4010e7:       75 ef                   jne    4010d8 <main+0x38>
  4010e9:       f2 41 0f 10 04 24       movsd  xmm0,QWORD PTR [r12]
  4010ef:       66 0f 2e c0             ucomisd xmm0,xmm0
  4010f3:       7a db                   jp     4010d0 <main+0x30>
  4010f5:       75 d9                   jne    4010d0 <main+0x30>
  4010f7:       48 8d 6b 08             lea    rbp,[rbx+0x8]
  4010fb:       49 83 c4 08             add    r12,0x8
  4010ff:       4c 39 fd                cmp    rbp,r15
  401102:       75 cc                   jne    4010d0 <main+0x30>
  401104:       bf 10 20 40 00          mov    edi,0x402010
  401109:       bd 00 00 00 c0          mov    ebp,0xc0000000
  40110e:       e8 1d ff ff ff          call   401030 <puts@plt>
  401113:       4c 01 ed                add    rbp,r13
  401116:       e8 25 ff ff ff          call   401040 <clock@plt>
  40111b:       f2 0f 10 1d 25 0f 00    movsd  xmm3,QWORD PTR [rip+0xf25]
  401122:       00
  401123:       48 89 c3                mov    rbx,rax
  401126:       48 8b 05 23 0f 00 00    mov    rax,QWORD PTR [rip+0xf23]
  40112d:       66 0f 28 c3             movapd xmm0,xmm3
  401131:       66 48 0f 6e c8          movq   xmm1,rax
  401136:       66 2e 0f 1f 84 00 00    nop    WORD PTR cs:[rax+rax*1+0x0]
  40113d:       00 00 00
  401140:       f2 0f 58 c3             addsd  xmm0,xmm3
  401144:       66 0f 28 e3             movapd xmm4,xmm3
  401148:       f2 41 0f 10 55 00       movsd  xmm2,QWORD PTR [r13+0x0]
  40114e:       49 83 c5 08             add    r13,0x8
  401152:       f2 0f 5d d1             minsd  xmm2,xmm1
  401156:       f2 0f 5e e0             divsd  xmm4,xmm0
  40115a:       66 0f 28 ca             movapd xmm1,xmm2
  40115e:       66 0f 28 c4             movapd xmm0,xmm4
  401162:       4c 39 ed                cmp    rbp,r13
  401165:       75 d9                   jne    401140 <main+0xa0>
  401167:       f2 0f 11 64 24 08       movsd  QWORD PTR [rsp+0x8],xmm4
  40116d:       f2 0f 11 14 24          movsd  QWORD PTR [rsp],xmm2
  401172:       e8 c9 fe ff ff          call   401040 <clock@plt>
  401177:       f2 0f 10 0c 24          movsd  xmm1,QWORD PTR [rsp]
  40117c:       f2 0f 10 44 24 08       movsd  xmm0,QWORD PTR [rsp+0x8]
  401182:       66 0f ef d2             pxor   xmm2,xmm2
  401186:       48 29 d8                sub    rax,rbx
  401189:       bf 2d 20 40 00          mov    edi,0x40202d
  40118e:       f2 48 0f 2a d0          cvtsi2sd xmm2,rax
  401193:       b8 03 00 00 00          mov    eax,0x3
  401198:       f2 0f 5e 15 b8 0e 00    divsd  xmm2,QWORD PTR [rip+0xeb8]
  40119f:       00
  4011a0:       e8 ab fe ff ff          call   401050 <printf@plt>
  4011a5:       e8 96 fe ff ff          call   401040 <clock@plt>
  4011aa:       48 89 c3                mov    rbx,rax
  4011ad:       48 8b 05 94 0e 00 00    mov    rax,QWORD PTR [rip+0xe94]
  4011b4:       66 48 0f 6e d8          movq   xmm3,rax
  4011b9:       66 48 0f 6e c0          movq   xmm0,rax
  4011be:       48 8b 05 8b 0e 00 00    mov    rax,QWORD PTR [rip+0xe8b]
  4011c5:       66 48 0f 6e c8          movq   xmm1,rax
  4011ca:       66 0f 1f 44 00 00       nop    WORD PTR [rax+rax*1+0x0]
  4011d0:       f2 0f 58 c3             addsd  xmm0,xmm3
  4011d4:       66 0f 28 eb             movapd xmm5,xmm3
  4011d8:       f2 41 0f 10 16          movsd  xmm2,QWORD PTR [r14]
  4011dd:       66 48 0f 7e c9          movq   rcx,xmm1
  4011e2:       66 48 0f 7e d0          movq   rax,xmm2
  4011e7:       f2 0f 5e e8             divsd  xmm5,xmm0
  4011eb:       48 99                   cqo
  4011ed:       48 d1 ea                shr    rdx,1
  4011f0:       48 31 c2                xor    rdx,rax
  4011f3:       66 48 0f 7e c8          movq   rax,xmm1
  4011f8:       48 c1 f8 3f             sar    rax,0x3f
  4011fc:       48 d1 e8                shr    rax,1
  4011ff:       48 31 c8                xor    rax,rcx
  401202:       66 0f 28 c5             movapd xmm0,xmm5
  401206:       48 39 c2                cmp    rdx,rax
  401209:       7d 04                   jge    40120f <main+0x16f>
  40120b:       66 0f 28 ca             movapd xmm1,xmm2
  40120f:       49 83 c6 08             add    r14,0x8
  401213:       4c 39 f5                cmp    rbp,r14
  401216:       75 b8                   jne    4011d0 <main+0x130>
  401218:       f2 0f 11 4c 24 08       movsd  QWORD PTR [rsp+0x8],xmm1
  40121e:       f2 0f 11 04 24          movsd  QWORD PTR [rsp],xmm0
  401223:       e8 18 fe ff ff          call   401040 <clock@plt>
  401228:       f2 0f 10 4c 24 08       movsd  xmm1,QWORD PTR [rsp+0x8]
  40122e:       f2 0f 10 04 24          movsd  xmm0,QWORD PTR [rsp]
  401233:       66 0f ef d2             pxor   xmm2,xmm2
  401237:       48 29 d8                sub    rax,rbx
  40123a:       bf 2d 20 40 00          mov    edi,0x40202d
  40123f:       f2 48 0f 2a d0          cvtsi2sd xmm2,rax
  401244:       b8 03 00 00 00          mov    eax,0x3
  401249:       f2 0f 5e 15 07 0e 00    divsd  xmm2,QWORD PTR [rip+0xe07]
  401250:       00
  401251:       e8 fa fd ff ff          call   401050 <printf@plt>
  401256:       48 83 c4 18             add    rsp,0x18
  40125a:       31 c0                   xor    eax,eax
  40125c:       5b                      pop    rbx
  40125d:       5d                      pop    rbp
  40125e:       41 5c                   pop    r12
  401260:       41 5d                   pop    r13
  401262:       41 5e                   pop    r14
  401264:       41 5f                   pop    r15
  401266:       c3                      ret
This comment wasn’t rated yet0
У вас не включены расширения -march=native. Впрочем и у меня на старой i7-3770 CPU @ 3.40GHz есть небольшой выигрыш у Вашего кода при смешивании операций:
next: 6.180340e-01 min: -1.797693e+308 time: 2.448361
next: 6.180340e-01 min: -1.797693e+308 time: 2.381875

На более современных Core(TM) i5-8265U CPU @ 1.60GHz- проигрывает.
next: 6.180340e-01 min: -1.797693e+308 time: 1.987874
next: 6.180340e-01 min: -1.797693e+308 time: 2.020840

Вообщем, интересно, но выигрыш не стоит выделки…

Total votes 1: ↑1 and ↓0+1
ОК, значит от формулировки a-tk «это может быть действительно полезно, если мы работаем на архитектуре, где нет аппаратной поддержки плавающей точки» уже переходим к «это может быть действительно полезно на процессоре без AVX» :-)
Total votes 2: ↑1 and ↓10

Да, я вижу ucomisd в коде. Вот только он соответствует строчке:


   } while (gig[i]!=gig[i]);

Первый цикл же у вас выглядит так:


  401140:       f2 0f 58 c3             addsd  xmm0,xmm3
  401144:       66 0f 28 e3             movapd xmm4,xmm3
  401148:       f2 41 0f 10 55 00       movsd  xmm2,QWORD PTR [r13+0x0]
  40114e:       49 83 c5 08             add    r13,0x8
  401152:       f2 0f 5d d1             minsd  xmm2,xmm1
  401156:       f2 0f 5e e0             divsd  xmm4,xmm0
  40115a:       66 0f 28 ca             movapd xmm1,xmm2
  40115e:       66 0f 28 c4             movapd xmm0,xmm4
  401162:       4c 39 ed                cmp    rbp,r13
  401165:       75 d9                   jne    401140 <main+0xa0>

Второй так:


  4011d0:       f2 0f 58 c3             addsd  xmm0,xmm3
  4011d4:       66 0f 28 eb             movapd xmm5,xmm3
  4011d8:       f2 41 0f 10 16          movsd  xmm2,QWORD PTR [r14]
  4011dd:       66 48 0f 7e c9          movq   rcx,xmm1
  4011e2:       66 48 0f 7e d0          movq   rax,xmm2
  4011e7:       f2 0f 5e e8             divsd  xmm5,xmm0
  4011eb:       48 99                   cqo
  4011ed:       48 d1 ea                shr    rdx,1
  4011f0:       48 31 c2                xor    rdx,rax
  4011f3:       66 48 0f 7e c8          movq   rax,xmm1
  4011f8:       48 c1 f8 3f             sar    rax,0x3f
  4011fc:       48 d1 e8                shr    rax,1
  4011ff:       48 31 c8                xor    rax,rcx
  401202:       66 0f 28 c5             movapd xmm0,xmm5
  401206:       48 39 c2                cmp    rdx,rax
  401209:       7d 04                   jge    40120f <main+0x16f>
  40120b:       66 0f 28 ca             movapd xmm1,xmm2
  40120f:       49 83 c6 08             add    r14,0x8
  401213:       4c 39 f5                cmp    rbp,r14
  401216:       75 b8                   jne    4011d0 <main+0x130>
This comment wasn’t rated yet0

Сейчас прогнал на Xeon Gold 6226R (clang 10). Получил идентичное время 1.86 для обоих сценариев, причём это время не изменилось, когда я оставил только деление, а вычисление минимума закомментировал — это показывает, что деление вычисляется независимо от всего остального и не мешает процессору сравнивать числа, используя FPU. При комментировании вычисления деления вышло 0.42 и 0.38 соответственно: int64_t таки оказался быстрее, но чтобы выжать эту производительность, нужно сильно изгаляться с кодом.

Total votes 1: ↑1 and ↓0+1
Если поменяете
if (gig[i] < min)
min = gig[i];
на
min = std::min(min,gig[i] );
то без деления int64 проигрывает.
This comment wasn’t rated yet0

Ага, щас! Код для первого цикла не поменялся, зато второй стал сильнее быстрее.


Без деления:


FPU: next: 1.000000 min: -1.79769e+308 time: 0.440429
ALU: next: 1.000000 min: -1.79769e+308 time: 0.219319

inline int64_t to_int64(int64_t a) {
        uint64_t mask = (uint64_t)(a >> 63) >> 1;
        return a ^ mask;
}

...

  begin = clock();
  min = 1e200;
  int64_t min_double = *(int64_t*)&min;
  int64_t min64 = to_int64(min_double);
  next = 1;

  for (long i=0; i<size; i++) {
    // next = 1./(next+1);
    int64_t x = to_int64(gig64[i]);
    min64 = std::min(x, min64);
  }
  end = clock();

  min_double = to_int64(min64); // Этой же функцией можно обратно вернуть всё в double, там же простой xor
  min = *(double*)&min_double;

Таки в итоге соглашусь с tyomitch, что сравнение double-ов через ALU может осуществляться быстрее, чем через FPU, однако:


  1. За счёт SIMD можно делать множество сравнений за раз. Вместо vminsd юзаем vminpd и получаем ускорение в 4 раза (или меньше, если упрёмся по скорости в память).


  2. Деление и сравнение таки выполняются независимо друг от друга, на разных FPU-модулях.


  3. Чтобы ALU-код оказался эффективен, нужно сильно его вылизывать. А если же не выливать, то наивный FPU-код оказывается эффективнее.



Update. Таки нет, это не ALU оказался быстрее, это компилятор слишком умный соптимизировал второй цикл до следующего монстра (первый цикл же остался скалярным):


  401320:       c4 c1 7e 6f 5c df a0    vmovdqu ymm3,YMMWORD PTR [r15+rbx*8-0x60]
  401327:       c4 c1 7e 6f 64 df c0    vmovdqu ymm4,YMMWORD PTR [r15+rbx*8-0x40]
  40132e:       c4 c1 7e 6f 6c df e0    vmovdqu ymm5,YMMWORD PTR [r15+rbx*8-0x20]
  401335:       c4 c1 7e 6f 34 df       vmovdqu ymm6,YMMWORD PTR [r15+rbx*8]
  40133b:       62 f1 c5 28 72 e3 3f    vpsraq ymm7,ymm3,0x3f
  401342:       62 f1 bd 28 72 e4 3f    vpsraq ymm8,ymm4,0x3f
  401349:       62 f1 b5 28 72 e5 3f    vpsraq ymm9,ymm5,0x3f
  401350:       62 f1 ad 28 72 e6 3f    vpsraq ymm10,ymm6,0x3f
  401357:       c5 c5 73 d7 01          vpsrlq ymm7,ymm7,0x1
  40135c:       c4 c1 3d 73 d0 01       vpsrlq ymm8,ymm8,0x1
  401362:       c4 c1 35 73 d1 01       vpsrlq ymm9,ymm9,0x1
  401368:       c4 c1 2d 73 d2 01       vpsrlq ymm10,ymm10,0x1
  40136e:       c5 c5 ef db             vpxor  ymm3,ymm7,ymm3
  401372:       62 72 a5 28 39 db       vpminsq ymm11,ymm11,ymm3
  401378:       c5 bd ef dc             vpxor  ymm3,ymm8,ymm4
  40137c:       62 f2 fd 28 39 c3       vpminsq ymm0,ymm0,ymm3
  401382:       c5 b5 ef dd             vpxor  ymm3,ymm9,ymm5
  401386:       62 f2 f5 28 39 cb       vpminsq ymm1,ymm1,ymm3
  40138c:       c5 ad ef de             vpxor  ymm3,ymm10,ymm6
  401390:       62 f2 ed 28 39 d3       vpminsq ymm2,ymm2,ymm3
  401396:       48 83 c3 10             add    rbx,0x10
  40139a:       48 81 fb 0c 00 00 18    cmp    rbx,0x1800000c
  4013a1:       0f 85 79 ff ff ff       jne    401320 <main+0x1b0>
Total votes 2: ↑2 and ↓0+2
Кошмар…
У меня
std:min(double) time: 0.294996
if (gig[i] < min) time: 0.493518
is_smaller(gig[i], min64) time: 0.428310
x < min64 time: 0.501825
is_smaller(gig[i], min) time: 0.619796
min64 = std::min(x, min64) time: 0.411116

Резюме — А ну его. Пусть компилятор оптимизирует. :)
This comment wasn’t rated yet0

Странно, у меня на ARM64 этот код дает неправильный результат. Вот такой код работает:


  begin = clock();
  min = 1e200;
  int64_t min64 = to_int64(min);
  next = 1;
  long j = -1;

  for (long i=0; i<size; i++) {
    // next = 1./(next+1);
    int64_t x = to_int64(gig[i]);
    if (x < min64) {
      min64 = x;
      j = i;
    }
  }
  end = clock();

Результат:


3GB random doubles generated
next: 1.000000e+00 min: -1.797693e+308 time: 0.725408
next: 1.000000e+00 min: -1.797693e+308 time: 0.511060

Можно переименовать статью: ускорения поиска минимального числа в массиве вещественных чисел.
Но как только в цикле появится переход от целого к вещественному, то все ускорение пропадет.

This comment wasn’t rated yet0
Чтобы ALU-код оказался эффективен, нужно сильно его вылизывать. А если же не выливать, то наивный FPU-код оказывается эффективнее.

Ну я же показал ситуацию, когда наивный ALU-код обгоняет наивный FPU-код: с делениями и без AVX.

Резюме — А ну его. Пусть компилятор оптимизирует. :)

Согласен, что этот трюк ценнее разработчикам компиляторов, чем прикладным программистам. Тем более, что он относительно кроссплатформенный: применим всюду, где IEEE floats.

Можно переименовать статью: ускорения поиска минимального числа в массиве вещественных чисел.

Уверен, что это не единственный сценарий, где этот трюк может дать выигрыш; но сходу не могу придумать, где ещё в тесном цикле много сравнений.
Total votes 1: ↑1 and ↓0+1

То есть вы предполагаете, что, заняв FPU чем-то долгим, целочисленные операции могут выполняться быстрее? Звучит, в принципе, разумно, но сценарий крайне неудачный. Операция деления и операция сравнения выполняются на разных FPU-модулях (см. документацию Agner), поэтому друг другу они не мешают.


Более правильный сценарий. Взять процессор с Hyper-Threading или SMT, создать два потока и заставить их работать на одном физическом ядре, далее рассмотреть следующие сценарии:


  1. В одном потоке минимум считается через FPU, в другом — через целочисленные операции.
  2. В обоих потоках минимум считается через FPU.
  3. В обоих потоках минимум считается через ALU.

После чего просуммировать время выполнения и сделать выводы.
Также стоит оптимизировать ALU-код, чтобы to_int64 не вызывалась для обоих аргументов, например, так:


  int64_t min2 = to_int64(gig[0]);

  for (long i=0; i<size; i++) {
    next = 1./(next+1);
    int64_t x = to_int64(gig[i]);
    if (x < min2)
      min2 = x;
  }

  double min = to_double(min2);
Total votes 3: ↑3 and ↓0+3
Нету там аргумента.
И при чём здесь вообще legacy FPU?
Total votes 3: ↑1 and ↓2-1

AArch64, Neoverse N1


$ gcc --version
gcc (GCC) 7.3.1 20180712 (Red Hat 7.3.1-12)
$ gcc -O3 test.c
$ ./a.out
3GB random doubles generated
next: 6.180340e-01 min: -1.797693e+308 time: 2.417564
next: 6.180340e-01 min: -1.797693e+308 time: 2.417601
$ gcc -O3 test_nodiv.c
$ ./a.out
3GB random doubles generated
next: 1.000000e+00 min: -1.797693e+308 time: 0.725358
next: 1.000000e+00 min: -1.797693e+308 time: 1.451623

Получается, что нативная реализация в 2 раза быстрее.
Ассемблерный код:


.L6: // цикл с нативным сравнением
    ldr d0, [x21, x0, lsl 3]
    add x0, x0, 1
    fcmp    d0, d8
    fcsel   d8, d8, d0, pl
    cmp x0, x1
    bne .L6
...
.L8: // цикл с is_smaller
    ldr d0, [x21, x0, lsl 3]
    fmov    x1, d8
    add x0, x0, 1
    fmov    x3, d0
    asr x2, x1, 63
    eor x2, x1, x2, lsr 1
    asr x1, x3, 63
    eor x1, x3, x1, lsr 1
    cmp x2, x1
    fcsel   d8, d8, d0, le
    cmp x0, x4
    bne .L8

Видно, что в цикле с is_smaller нам приходится пересылать значения из FP регистров D0/D8 в целочисленные регистры X1/X3. Цикл с is_smaller содержит в 2 раза больше инструкций.

Total votes 1: ↑1 and ↓0+1
Видно, что в цикле с is_smaller нам приходится пересылать значения из FP регистров D0/D8 в целочисленные регистры X1/X3.

Это к GCC вопрос, зачем он вообще использует D0/D8 внутри цикла, и FCSEL вместо CSEL.
This comment wasn’t rated yet0

Потому что GCC предполагает, что для работы с double будут использоваться в основном FP инструкции, что в целом правда. Так как целочисленных регистров всегда не хватает, компилятор старается использовать все другие регистры.
Кстати в коде для x86_64, значения также пересылаются из FP регистров XMM в целочисленные регистры. И это код сгенерированный clang.

Total votes 1: ↑1 and ↓0+1
Кстати в коде для x86_64, значения также пересылаются из FP регистров XMM в целочисленные регистры. И это код сгенерированный clang.

Кстати да. Если читать сразу int64_t, а не перегонять double в int64_t, да ещё и применить оптимизацию по вычислению int64_t только один раз, то XMM-регистры не задействуются, и у меня выходит:


0.336116 FPU, 0.369728 ALU без деления в clang.
0.444854 FPU, 0.328019 ALU без деления в gcc.


Но с делением всё те же 2.52 и 2.62 (clang), 2.62 и 2.64 (gcc). Удивительно.


Clang почему-то не смог полностью избавиться от xmm-регистров, а gcc не смог задействовать minpd.

Total votes 2: ↑2 and ↓0+2
Вы ЭТО называете бенчмарком?
Где прогрев, где статистика, где банально устранение лишних обращений в память?
И странная строчка, да.
Total votes 4: ↑1 and ↓3-2
А почему нет? SAR+SHR+XOR+CMP — четыре простейшие операции, каждая выполняется за такт.

А остальная обвязка из 6–7 команд, получается, вообще за 0 тактов выполняется?


P.S. COMISD — 5 тактов в Zen 2, в современных процессорах Intel — порядка 2-3 тактов.

This comment wasn’t rated yet0
Какую «остальную обвязку из 6–7 команд» вы имеете в виду?
This comment wasn’t rated yet0
При фактическом использовании компилятор её инлайнит, и mov-ы исчезают вместе с setal и ret.
This comment wasn’t rated yet0

Допустим. Но компилятор все равно не понимает, что to_int64() нужно вызывать только один раз, а по не два раза на каждое сравнение.

This comment wasn’t rated yet0
Я не понимаю, о чём вы: to_int64 нужно применить к каждому из сравниваемых чисел, как это можно сделать за один вызов?
Total votes 2: ↑0 and ↓2-2
inline bool is_smaller(double x1, double x2) {
	return to_int64(x1) < to_int64(x2);
}
inline bool is_smaller(double x1, int64_t x2) {
	return to_int64(x1) < x2;
}
...
clock_t begin = clock();
  double min = 1e200;
  int64_t  min64 = to_int64(min);
  double next = 1;
  for (long i=0; i<size; i++) {
    next = 1./(next+1.0);
    if (is_smaller(gig[i], min64))
    {
	 min =    gig[i];
        min64 = to_int64(min);  
    }
  }
  clock_t end = clock();
  printf("next: %e min: %e time: %f\n", next, min, (double)(end - begin) / CLOCKS_PER_SEC);
Total votes 2: ↑2 and ↓0+2

У вас код неоптимальный. Если условие сработает, то to_int64 второй раз выполнится. Тогда уж так надо:


  begin = clock();
  min = 1e200;
  int64_t min64 = to_int64(min);
  next = 1;

  for (long i=0; i<size; i++) {
//    next = 1./(next+1);
    int64_t x = to_int64(gig[i]);
    if (x < min64)
    {
      min = gig[i];
      min64 = x;
    }
  }
  end = clock();
  printf("next: %f min: %f time: %f\n", next, min, (double)(end - begin) / CLOCKS_PER_SEC);

У меня это дало ускорение с 0.554674 до 0.474570, но до 0.333950 все равно далеко.

Total votes 1: ↑1 and ↓0+1
Да я знаю. Я просто для иллюстрации. Оптимальный min = std::min(min,gig[i] ); :) 0.268540

Если условие сработает, то to_int64 второй раз

При равномерном распределение разница не принципиальна. Мне заметить не удалось. Это условие очень редко выполняется.
This comment wasn’t rated yet0
Sign up to leave a comment.

Articles

Show the best of all time

Your account

Sections

Information

Services

Support
© 2006–2024, Habr