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;
}

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


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


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


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


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


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

Склонен согласиться. На фоне данного трюка использующего особенности формата хранения fp64 в памяти, данный UB совсем безобиден.
В частности одна из причин UB — неоднозначность формата хранения, вспомним 8087 с его 80-битными флоатами, даже если в коде они определены как FP64…
Примерно как соблюдение правил движения, когда убегаешь от полиции после ограбления
вспомним 8087 с его 80-битными флоатами
Извините, то никаких 80-битных флоатов 8087 в памяти не размещает. И с ним указанный трюк (при использовании memcpy) будет отлично работать.

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

P.S. Чисто для тех, кто в танке: в стандарте существует такая вещь, как is_iec559 — именно для того, чтобы всё описанное в статетье не было “ограблением банка”.
8087 может работать с 80-битными флоатами в памяти и это даже будет is_iec559 тип. Но правда 754 вроде требует сейчас от компилятора гарантировать что double это именно 64-битный тип и ни один компилятор без явной просьбы использовать long double подобный код не генерировал
8087 может работать с 80-битными флоатами в памяти и это даже будет is_iec559 тип.
Это закрывается банальным static_assert(sizeof(double) == sizeof(std::uit64_t));.

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

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

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

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

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

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

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

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

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

Но забавно что как в силу специфики реализации железа как раз в C с этим подходом и возникают проблемы на старом железе (современные процессоры от концепции сопроцессоров давно отказались: либо плавучка встроена в основной процессор, либо её нет вообще и она эмулируется).
Ну memcpy то целиком на 8 байт вызывать не стоит, наверное… Мы же об эффективности говорим?

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

Если добавить -O1, то не будет и них.
Спасибо, покопался. -O3 дает возможность всем компиляторам убить memcpy() кроме… msvc. Может там еще чего надо, я хз, не умею в него.

godbolt.org/z/GGnxr1
У MSVC нужный ключ называется /O2

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

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

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

Вы не путаете reinterpret_cast и bit_cast?

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

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

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


Но да, оно по факту работает так же как с указателями, я ошибся )

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

Разработчики АЛУ в железе не могли обойтись без обработки NaN.

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

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

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

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

Преимущество этого трюка перед FCMP в том, что FPU остаются свободными, и могут параллельно со сравнением вычислять что-нибудь ещё.

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


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

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

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

Теперь ждём трюков как складывать длинные целые числа.
Это может быть действительно полезно, если мы работаем на архитектуре, где нет аппаратной поддержки плавающей точки, и нет нужды работать ненормальными числами. Но на современном компьютере… Не верю!
Смотрим во что компилируется обычное наивное сравнение: godbolt.org/z/s3hMjh
Видим инструкцию ucomisd xmm0, xmm1 на x86.
Читаем здесь: mudongliang.github.io/x86/html/file_module_x86_id_316.html
Latency не более 7 тактов.
Вы думаете это медленнее чем остальная портянка с зависимостями по данным будет?
Latency не более 7 тактов.
Вы думаете это медленнее чем остальная портянка с зависимостями по данным будет?

А почему нет? SAR+SHR+XOR+CMP — четыре простейшие операции, каждая выполняется за такт.
Ну вот и докажите бенчмарками.
Вы, кстати, лукавите насчёт 4 операций, потому что это надо сделать два раза, с каждым аргументом.
Вы, кстати, лукавите насчёт 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

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


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>
У меня вот так
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
У Вас АЛУ тупо забит делением.

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

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

Я сейчас прогнал вышезапощенный код на ещё одной машине (Xeon 5118 @ 2.30GHz, GCC 9.3.1) и получил ещё большую разницу (3.20 против 2.71). У вас есть для этого экспериментального наблюдения иное объяснение, чем «ucomisd не может выполняться одновременно с вычислениями»?

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

Вся 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
У вас не включены расширения -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

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

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

Да, я вижу 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>

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

Если поменяете
if (gig[i] < min)
min = gig[i];
на
min = std::min(min,gig[i] );
то без деления int64 проигрывает.

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


Без деления:


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>
Кошмар…
У меня
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

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

Странно, у меня на 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

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

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

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

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

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

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

Уверен, что это не единственный сценарий, где этот трюк может дать выигрыш; но сходу не могу придумать, где ещё в тесном цикле много сравнений.

То есть вы предполагаете, что, заняв 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);
Нету там аргумента.
И при чём здесь вообще legacy FPU?

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 раза больше инструкций.

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

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

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

Кстати в коде для 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.

Вы ЭТО называете бенчмарком?
Где прогрев, где статистика, где банально устранение лишних обращений в память?
И странная строчка, да.
А почему нет? SAR+SHR+XOR+CMP — четыре простейшие операции, каждая выполняется за такт.

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


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

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

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

Я не понимаю, о чём вы: to_int64 нужно применить к каждому из сравниваемых чисел, как это можно сделать за один вызов?
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);

У вас код неоптимальный. Если условие сработает, то 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 все равно далеко.

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

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

При равномерном распределение разница не принципиальна. Мне заметить не удалось. Это условие очень редко выполняется.
Only those users with full accounts are able to leave comments. Log in, please.