Быстрый floor()

Всем привет.
Собственно вот встала тут задача.
Надо реализовать функцию floor_i32() которая на входе принимает float, а на выходе выдаёт округлённое к минус бесконечности целое число (в целочисленном формате).

Пока сделал так:

Код (Text):

1. int floor_i32(float x)
2. {
3.   int i=x;
4.   if(i>x) i--;
5.   return i;
6. }

Но всё приходит к тому с чего началось.
А началось всё с того что floor() из RTL сильно медленный,
В приведённом выше коде проблема вся в приведении float к int.
Компилятор пихает функцию из RTL, а она медленная, т.к. содерхит дофига проверок, изменяет fpu control word плюс результат возвращает в __int64 - а это лишние вычисления.
Собственно и от floor() я отказался потому что она возвращает float, и вылазит то же приведение типов.
Вот возник вопрос - как сделать быстрое приведение float к int.
Теоретически мой сишный код написанный выше будет работать в любом случае независимо от режима округления fpu, т.е. меня устроит приведение float к int простыми fld/fistp, но как это компилятору объяснить...
Встроенный ассемблер не подходит по той причине что компилятор тогда функцию не инлайнит
Есть какие-нибудь трюки для приведения float к int, или для округления к минус бесконечности?

 

А так не покатит:

Код (Text):

1. int floor_i32(float x)
2. {
3.   void* tmp = (void*)&x; //возможно к void* и приводить не надо
4.   int* tmp2 = (int*)tmp; //нужно потестить
5.   int tmp3 = *tmp2;
6. //округляем руками tmp3, как нам надо
7.  
8.   return tmp3;
9. }

Пробовал, компилятор инлайнит и никаких преобразований не делает.

 

Это то же самое что и
i32=*(int *)&f32;
Только ты получиш не округлённое значение.
Мне для 1.7 надо получить 1, а твой код получит 3F800000h.
Парсить руками точно быстрее не будет.

 

сделай макросом

 

в С99 есть функции

Код (Text):

1. long int   lrint    (double x) ;
2. long int   lrintf   (float x) ;

http://mega-nerd.com/FPcast/
но в MSVS их нету.

 

Для чисел double:

Код (Text):

1. volatile double magic;
2. // вычисление magic
3. _control87( _PC_64 | _RC_DOWN, _MCW_PC | _MCW_RC );
4. volatile double x;
5. magic = 1.0;
6. while ( true )
7. {
8.   x = magic + 0.5;
9.   if ( x == magic )
10.     break;
11.   magic *= 2.0;
12. }
13. _control87( cw, 0xfffff );

Собственно сам floor():

Код (Text):

1. // на входе - _val типа double
2. static volatile double x;
3. x = _val + magic; // Теперь *(int*)( &x ) == floor( _val )
4. return *(int*)( &x );

 

И да, нужно перед вызовом этой функции выставить режим округления - какой выставишь, так и округлит.

 

это шутка?

 

Y_Mur
Это не шутка, а реальный код, используемый в одном моем проекте. Идея взята из книги Уоррена "Алгоритмические трюки для программистов". Использует особенности бинарного представления вещественных чисел.

 

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

Хмм... А что, одно сообщение к другому автоматически не подцепляется?

 

Y_Mur
почему шутка? вполне рабочий код

 

nobodyzzz
C99 не подходит. Поищу исходники реализации.
Atlantic
Я так понимаю что твой код это более общий вид следующего:

Код (Text):

1. int ifloor(float sf)
2. {
3.   int      i32;
4.   double df=6755399441055744.0; // 2^51+2^52
5.   df+=sf;
6.   i32=*(int *)&df;
7.   if(i32 > sf) i32--;
8.   return i32;
9. }

Здесь if(i32 > sf) i32--; стоит для того чтобы не изменять FPU Control Word - это очень длительная операция которая к тому же сбрасывает конвейер.
И это в принципе работает.
Я только понять не могу откуда берётся эта константа, и почему при её прибавлении в мантиссе оказывается округлённое число.
Может кто-нибудь объяснить?

Взято вот отсюда: http://www.df.lth.se/~john_e/gems/gem0042.html
PS: да, кстати, а где про это у Уоррена написано (в какой главе)? я там первым делом смотрел - но безуспешно...

 

cppasm
Первая единица в этой константе нужна чтобы гарантировать что порядок числа будет таким, чтобы бит с весом 2^0 находился в младшем бите мантисты. А вторая единица чтобы если число отрицательное, то не испортилась первая единица

 

Вроди понял.
Сочинил такое:

Код (Text):

1. int __inline f2int(float f32)
2. {
3.     int i32;
4.     float   tmp;
5.     tmp   = f32+12582912.0f;        // 2^22+2^23
6.     i32   = *(int *)&tmp & 0x3FFFFF;
7.     i32 <<= 10;
8.     i32 >>= 10;
9.     if(i32 > f32) i32--;
10.     return i32;
11. }

Диапазон значений -(2^21) ... (2^21)-1
Меня устраивает.
Если надо 32-битный результат - можно tmp сделать double и константу поменять.
Но вот возникли вопросы: закреплено ли в стандарте что компилятор Си обязан для вещественных чисел использовать формат IEEE?
Для другого формата хранения финт не пройдёт...

 

cppasm
Мне почему-то кажется что код где промежуточная переменная double будет работать быстрее чем код с двумя сдвигами.

 

Сейчас померяю.
С другой стороны Агнер Фог пишет что не рекомендуется смешивать вычисления float с double...
Хотя в варианте с float на 2 сдвига и одно И больше.

 

Померял.
Вариант с double быстрее.
Ещё чуть-чуть быстрее такой вариант:

Код (Text):

1. int f2int_f64r(float f32)
2. {
3.     int i32;
4.     double  f64,magic=6755399441055744.0;   // 2^51+2^52
5.     f64  =  f32+magic;
6.     i32  =  *(int *)&f64;
7.     f64 -=  magic;
8.     if(f64 > f32) i32--;
9.     return i32;
10. }

Но тут для нормальной работы требуется чтобы разрядность мантиссы была 53бита (это по умолчанию так, но можно изменить на 64).
Тогда прибавление 2^51+2^52 отсекает всю дробную часть.
Соответственно потом при вычитании этой же константы получаем исходное число, но без дробной части (округлённое).

Есть у кого ещё какие идеи как ускорить можно?

 

У меня есть такая версия:

Код (Text):

1. int ifl(float f32)
2. {
3.     if(-1.0<f32&&f32<0.0)//if(0xBF800000>*(int*)&f32&&*(int*)&f32>0x80000000)
4.         return -1;
5.     return (int)(f32+0x40000000)-0x40000000;
6. }

работает немного быстрее варианта с double, наверно потому что для приведения к int MSVS2005 использует _ftol2_sse, который использует cvttsd2si. Можно заменить условие тем что в комментариях, вроде еще быстрее работает, только это уже хак.

[add]посмотрел что за код со вторым вариантом генерится:

Код (Text):

1. 004011F4  mov         edx,dword ptr [esp+30h]
2. 004011F8  add         edx,7FFFFFFFh
3. 004011FE  cmp         edx,3F7FFFFEh
4. 00401204  ja          main+20Bh (40120Bh)
5. 00401206  or          eax,0FFFFFFFFh
6. 00401209  jmp         main+21Bh (40121Bh)
7. 0040120B  fld         dword ptr [esp+30h]
8. 0040120F  fadd        st,st(1)
9. 00401211  call        _ftol2_sse (401970h)
10. 00401216  sub         eax,40000000h
11. 0040121B

впечатляет

 

У м$явок есть идеи. в msvcrt

Код (Text):

1. oword_7C0412B0  xmmword 3FF00000000000003FF0000000000000h
2. oword_7C0412C0  xmmword 4330000000000000433h
3. oword_7C0412D0  xmmword 4330000000000000BFF0000000000000h
4. oword_7C0412E0  xmmword 80000000000000008000000000000000h
5. oword_7C0412F0  xmmword 7FFh       
6.  
7.         public floor
8. floor       proc near
9.  
10. var_10      = dword ptr -10h
11. var_C       = dword ptr -0Ch
12. var_8       = dword ptr -8
13. var_4       = dword ptr -4
14. arg_0       = qword ptr  4
15.  
16.         cmp dword_7C04B14C, 0
17.         jz  sub_7C0363B7
18.         sub esp, 8
19.         stmxcsr [esp+8+var_4]
20.         mov eax, [esp+8+var_4]
21.         and eax, 1F80h
22.         cmp eax, 1F80h
23.         jnz short loc_7C00990F
24.         fnstcw  word ptr [esp+8+var_8]
25.         mov ax, word ptr [esp+8+var_8]
26.         and ax, 7Fh
27.         cmp ax, 7Fh
28.  
29. loc_7C00990F:               ; CODE XREF:
30.         lea esp, [esp+8]
31.         jnz sub_7C0363B7
32.         movq    xmm0, [esp+arg_0]
33.         movapd  xmm2, oword ptr ds:oword_7C0412C0
34.         movapd  xmm1, xmm0
35.         movapd  xmm7, xmm0
36.         psrlq   xmm0, 34h
37.         movd    eax, xmm0
38.         andpd   xmm0, oword ptr ds:oword_7C0412F0
39.         psubd   xmm2, xmm0
40.         psrlq   xmm1, xmm2
41.         test    eax, 800h
42. ...и тому подобный мрак
43. floor       endp

вероятно нечто очень быстрое, что лишний CALL не жаль потратить. да уж..

зы:
хм... cppasm а вообще - я за вариант с float. ибо Intel C++ v8.0 (c опциями /O3 /Ot /G6 /arch:SSE) выдал вот такую фкусняшку под трете-пень.

Код (Text):

1.     PUBLIC _f2int
2. _f2int  PROC NEAR
3. ; parameter 1: 8 + esp
4. .B1.1:                          ; Preds .B1.0
5.         push      esi
6.         movss     xmm4, DWORD PTR [esp+8]
7.         movss     xmm3, DWORD PTR _2il0floatpacket.1
8.         addss     xmm3, xmm4
9.         movss     DWORD PTR [esp], xmm3
10.         mov       eax, DWORD PTR [esp]
11.         and       eax, 4194303
12.         shl       eax, 10
13.         sar       eax, 10
14.         lea       edx, DWORD PTR [eax-1]
15.         cvtsi2ss  xmm5, eax
16.         comiss    xmm4, xmm5
17.         cmovb     eax, edx
18.         pop       ecx
19.         ret
20.         ALIGN     4
21. _f2int ENDP

А на double он таково не выдает :'( (ток инструкшены под FPU )

ззы:
и такаяже картина c моим любимым GCC 3.4.6

Код (Text):

1.     .def    _f2int; .scl    2; 
2. _f2int:
3.     push    ecx
4.     fld DWORD PTR [esp+8]
5.     fld st(0)
6.     fadd    DWORD PTR LC0
7.     fstp    DWORD PTR [esp]
8.     mov eax, DWORD PTR [esp]
9.     sal eax, 10
10.     sar eax, 10
11.     cvtsi2ss    xmm0, eax
12.     lea edx, [eax-1]
13.     movss   DWORD PTR [esp], xmm0
14.     fld DWORD PTR [esp]
15.     fucomip st, st(1)
16.     fstp    st(0)
17.     cmova   eax, edx
18.     pop edx
19.     ret

 

Я проверял на Intel C++ 10
Разница между вариантами с float и double в пару тактов (1-2), проверял на 100000000 итерациях.
При включённом SSE вариант с float незначительно выигрывает.
Тут вся фигня в том, что при включённом SSE значительно (тактов 10-15) выигрывает обычное приведение к int.
Проблема вся в том что без SSE этот вариант отстаёт очень значительно.

На msvc 2005 и 2008 как с SSE, SSE2, так и без выигрывает вариант с double.

Вот и сиди выбирай :/

С вариантом Black_mirror сейчас буду разбираться...

PS: предпочтение отдаётся FPU, если код будет работать быстрее и на реализации с SSE - это просто ещё один плюс в его пользу, но не основной.

PPS: кстати вот замеры на 100000000 итераций
опции компилятора: icl /O3 /Ot /G6 /arch:SSE round.c

Код (Text):

1. Execution time (float)  ->      2103813975 clocks
2. Execution time (double) ->      2316319823 clocks
3. Execution time (c-cast) ->      1601596140 clocks

Не такое уж и огромное преимущество.