Dotaz: Vzdálenost dvou reálných kladných čísel v C

[Thunder.m]

Provádím různé optimalizace kódu, dostal jsem se do části která se provádí asi nejvícekrát ze všech a to výpočet absolutní hodnoty rozdílu dvou čísel.

V tuto chvíli používám pro data signed int (4B) ale rád bych použil pouze unsigned short (2B).

Už jsem dělal nějaké pokusy v rychlosti, jde o kritickou část kódu která se nedá jinak zoptimalizovat.

Jsou dvě možnosti, při zpracování dat:

• data přetypovat na signed int a poté provést funkci abs(x - y)
• použít podmínku

  if (x > y) {
  z = x - y
  } else {
  z = y - x
  }
  

• Třetí metoda je již jen modifikací druhé a to vytvořit makro

  #define ABS_DIFF(X,Y) ((X > Y) ? (X - Y) : (Y - X))
  

Bohužel C zpracuje nejrychleji první metodu s přetypováním, ale i tak je rozdíl zpomalení docela výrazný (cca 30%) proti původnímu abs bez přetypování, nevíte jak se dá rychleji ideálně nativní funkcí spočítat rozdíl dvou reálných celočíselných čísel?

---

---

Odpovědi

[Thunder.m]

Oprava, reálných kladných celých čísel :)

[chrono]

Je nejaký špeciálny dôvod, prečo potrebuješ práve 2 bajtový typ?

[Thunder.m]

No tak protože hodnoty nabívají hodnot od 0 - 50 000, chci menší datový typ, protože pokud mám data načtená v operační paměti a že jich je opravdu velké množství, zabírají zbytečně 2x tolik místa.

[Thunder.m]

nabývají :D

[chrono]

V takom prípade (a aj vo veľa iných prípadoch) je potom najlepšie použiť priamo funkciu abs, pretože ju kompilátor optimalizuje (platí to minimálne pre gcc a llvm).

Predpokladám, že to ABS makro by sa dalo, minimálne pre gcc, urobiť tak, aby to generovalo lepší kód, ale je otázne, či to má zmysel a či to bude v priemere rýchlejšie ako tá abs funkcia.

[Pavel Löbl]

V glibc to vypada neak tak:

/* Return the absolute value of I.  */
int
abs (int i)
{
  return i < 0 ? -i : i;
}

na tom mym srotu (i686) to pri pouziti makra a prelozeni s -O1 vychazi podobne jako volani abs, a mezi pouzitim short a int neni taktez temer rozdil.

Dost mozny, ze to brzdi ten short. Prace s typem, kterej nema velikost slova muze bejt asi drazsi.

[chrono]

Rýchlejšie bude niečo ako

int abs(int i)
{
    int t = i >> (32 - 1);
    return (i ^ t) - t;
}

(a presne to isté generuje gcc pre 32 bitový int)

[Thunder.m]

No a napadá tě podobná funkce pro rozdíl unsigned celých čísel? Bitové operace mi nikdy moc nešly :)

[chrono]

Nič lepšie ako tú tvoju prvú možnosť asi nevymyslíš. :) (pretože 16 bitov na to odpočítavanie stačiť nebude, takže aj tak nakoniec skončíš pri aspoň 32 bitovom type).

Možno by sa dali použiť MMX inštrukcie, ale je otázne, či by to bolo rýchlejšie.

[l4m4]

Nejspíš byly, protože operují na 4 argumentech současně a existuje přímo instrukce na abs vektoru 16bitových čísel se znaménkem, viz můj příspěvek níže.

[Sixkillers]

Tak tady je trochu upravena varianta.

unsigned short my_abs(int i)
{
    int t = i >> (sizeof(int) * CHAR_BIT - 1);

    /* patent free */
    return (i + t) ^ t;
    /*return (i ^ t) - t;*/
}

Jinak ta varianta s minusem je možná patentově chráněna (viz.) :)

Např.:

unsigned short x = 5;
unsigned short y = 20;
unsigned short z = my_abs(x - y);

Přetypování my_abs(x - y) typicky zabere jeden procesorový cyklus. Přetypování int na unsigned short je bez výkonové penalizace.

[Josef Kufner]

A co se tomu přetypování vyhnout třeba tím, že by se už na začátku počítalo se správným typem? Nevím co tam okolo máš, ale snad by to mohlo jít nějak přerovnat. Případně bych zkusil se té absolutní hodnoty zbavit (nějak jinak zaručit, že tam budou vždy jen kladná čísla nebo aby nevadilo, že bude záporné).

[Thunder.m]

To bohužel nejde, unsigned short má rozsah +-32 768, ale já potřebuju porovnávat čísla od 0 do 50 000.

Myslím že ani nejde posunout si nulu (předem data přesunout do spodní části), protože bych narazil na možné překročení limitu unsigned shortu, např abs(-15 000 - 32000) nelze.

[Josef Kufner]

A co nějak zaručit, že a bude vždy větší neź b, takže bys nemusel dělat abs(), ale jen bys odčítal.

[Sten]

Funkce abs dělá:

return i < 0 ? -i : i;

Všechny tři metody jsou ve výsledném kódu totožné.

Doporučuji nepoužívat short, operace s ním jsou dražší než s intem (resp. longem na 64bitech).

[l4m4]

Pokud ten program skutečně tráví nejvíc času výpočtem absolutních hodnot, asi to dělá uvnitř vnitřních smyček na spoustě argumentů, a v tom případě se podívej na pabsw (v gcc __builtin_ia32_pabsw), psubw a podobné.

[Thunder.m]

Ano přesně tak, vnitřní smyčky probíhající cca 20x, tj. podobný kód

function soucet_abs_hodnot(unsigned short *input1, unsigned short *input2, unsigned int size) {
int z = 0, i;
for (i = 0; i < size; i++) z+= abs(input1[i] - input2[i]);
return z;
}

Takhle přesně vypadá funkce kterou chci zoptimalizovat.

[l4m4]

V tom případě se použití vektorových instrukcí nejspíš vyplatí. Kromě ručního psaní assembleru znám dvě cesty. Jednak gcc builtiny, které provádějí konkrétní instrukce procesoru, takže jsou platformově závislé (tj. budeš potřebovat přinejmenším fallback pro ostatní architektury, případně i více implementací kritického kódu). Druhak knihovna ORC, která umožňuje napsat jakýsi pseudokód, který pak převede do konkrétních instrukcí sama (i run-time); vypadá nadějně, je ovšem zatím poměrně dost ve vývoji.

[dooku]

Zkoušel jsem několik metod:
1> jednu zde jmenovanou - pomocí bitového posunu a nonekvivalence
2> nejjednodušší možnou porovnávající číslo s nulou pomocí podmíněného operátoru
3> abs() z math.h
4> mou metodu - zařadil jsem ji jen pro porovnání, nepočítal jsem s tím že by měla závratné výsledky:

int abs(int i) {
 if (i & 0x80000000) return -i; // hex hodnota se samozřejmě liší podle použitých číselných typů
 return i;
}

Každou metodou bylo porovnáno 4 000 000 000 čísel, tedy kousek pod limitem unsigned int. Nepřetypovával jsem, zajímal jsem se pouze o efektivnost algoritmů výpočtu absolutních hodnot. Výsledky jsou následující:
1> 21.65 sekund
2> 23.27 sekund
3> 15.44 sekund
4> 20.31 sekund

Měření jsem několikrát opakoval a výsledky se liší jen minimálně. Opět jsem se přesvědčil, že pokoušet se překonat optimalizované funkce je marné - je to stejné i v případě jiných matematických funkcí. Takže myslím, že pokud nechcete zabředávat do assembleru - což je sice zajímavé, avšak zdlouhavé a v mnoha případech nepraktické - nemá smysl zkoušet jiné algoritmy. Možná vektorové instrukce, ale to je docela kanón na vrabce. Lepší už to prostě nebude :)

[Thunder.m]

Bohužel v tomhle máš opravdu pravdu a protože počítání absolutních hodnot program tráví 99% procesorového času je to bohužel hodně znat, takže mi asi nezbývá nic jiného, než data předávat jako float, nebo integer a použít interní abs funkci, než jakékoliv jiné metody.

Data ve formátu unsigned int

int sim = 0;
int i;
for (i = 1; i <= input2[0]; i++) sim+= abs((int) input1[i] - (int) input2[i]);
return sim;

real	0m16.898s
user	0m56.828s

#define MACRO_DIST(X,Y) ((X < Y) ? (Y - X) : (X - Y))
unsigned short sim = 0;

for (i = 1; i <= input2[0]; i++) sim += MACRO_DIST(input1[i],input2[i]);

return sim;

real	0m20.070s
user	1m18.761s

Data ve formátu float

float sim = 0;
int i;
for (i = 1; i <= input2[0]; i++) sim+= fabs(input1[i] - input2[i]);
return sim;

real	0m12.351s
user	0m33.758s

[l4m4]

Pokud tím program tráví 99% času, tak se možná na vektorové gcc builtiny podívej. Aritmetické operace na integerech i floatech lze provádět na blocích 4 hodnot. Při -O3 taky gcc zkouší vektorizovat sám, ale nevím, jak je v tom úspěšný.

[dustin]

Nebo to napsat pod OpenCL nebo CUDA + GPU, tam by těch paralelních výpočtů mohlo běžet o řád víc:)

[Thunder.m]

Pokud se to co píšeš týká jen paralelních výpočtů tak mám celý program psaný vícevláknově přes pthread, škálování s počtem CPU je téměř lineární, měl jsem možnost otestovat až 12 fyzických CPU. Otázka je jestli by pak takové úpravy přinesly ještě nárust výkonu.

[stativ]

Určitě, protože pokud bys mohl pomocí SIMD počítat 4 hodnoty na každém procesoru, na 12 CPU by tě to paralelně počítalo 48 hodnot. A to je slušné zrychlení.

[stativ]

Pokud by ses rozhodl pro použití SSE, doporučuji mrknout na stránky Intel AVX. Dá se tam stáhnout „Intel Intrinsic Guide,“ což je nejlepší referenční příručka pro SSE co znám.

[Thunder.m]

Tak to zní skvěle, výpočty provádím na procesorech architektury Core, s podporou starších CPU se nepočítá. Pokusím se tedy o tom něco najít na internetu.

[Šangala]

Neznám celý problém, ale než pthread může být na úrovni výpočtů efektivnější OpenMP a je možné automaticky použít tolik thredů, kolik je jader k dispozici, bez nějakého dalšího programovaní.
Pokud ovšem máte nárust výkonu lineární s násobkem jader, tak to asi nepomůže :)

[a]

Mozno by bolo dobre, keby Ste popisali nasledujuce veci:

- odkial pochadzaju vstupne data (premenne "input1" a "input2")

- ake (matematicke) vlastnosti maju vstupne data

- ako sa vstupne data vyvijaju v case

[Thunder.m]

input1 a input2 jsou indexy seřazeného pole, které mají hodnoty od 0 - 50000, jsou celočíselné a vždy kladné, za celý chod programu se nemění. Obš pole jsou až třetí dimenzí nadřazených polí, snad jsem to popsal dostatečně. Obě pole mají průměrnou velikost cca 200 ale mohou mít až velikost 25 000.

[a]

Ake IPC (min, max, avg) dosahuje ten loop, co pocita vzdialenost.

[Thunder.m]

Aha, tak na takové úrovni ještě v C nejsem, nevím co tím myslíš, ale pokud myslíš kolik tím tráví času tak je to drtivá většina, co se týká dostupných různých polí input1 a input2 pak je to:

input1: 18000x200 input2: 5000x150x200

První je pole dvourozměrné, druhé je pole třírozměrné.

V rámci zásadních optimalizací už neporovnávám navzájem celé hodnoty, ale pouze vzdálenosti seřazených indexů a k tomu slouží ta funkce o které tu mluvíme, spočítá mi celkový rozdíl vzádlenosti mezi indexy polí.

[Josef Kufner]

Ono tohle už není ani tak o C jako spíš o té matematice za tím. Třeba by se to dalo nějak upravit, že bys například rozlišil několik případů, ty vyřešil samostatně a efektivněji a nakonec dal výsledky dohromady. Nebo si někde přidal mezivýsledky a později ušetřil výpočty... Současný algoritmus zjevně narazil na hranice svých možností, takže pokud chceš víc, musíš to vzít odjinud.

[a]

Vzhladom na velkost dat, asi to nebude tak, ze cele data by sa generovali nanovo kazdu minutu. Z toho mi vyplyva, ze N-te spustenie programu dostane priblizne tie iste 'input1' a 'input2' ako (N+1) spustenie. Vznika otazka, kedze sa 'input1' a 'input2' prilis nemenia, naco je potrebne program spustat opakovane (nestaci iba raz spustit?).

Alebo to je tak, ze treba vykonat 18000*5000*150 iteracii, co by celkovo zabralo vela casu (den, tyzden)?

[Thunder.m]

Ano v tuto chvíli je opravdu potřeba tolik iterací, na testovaném stroji Intel Core i7 2600K to zvládne za desítky vteřin.

[dustin]

Už jsem to psal jednou, ale nebylo by to vážně zralé počítat na mnoha desítkách paralelních jednotek GPU?

[__dark__]

Příloha:

• main.cpp (4090 bytů)

Ahoj, tady máš základní kód, který by ti mohl pomoct se odrazit dál. Implementace používá buď čisté C, kde si můžeš upravit MY_ABS makro (pokud to pomůže), pak SSE2, kde není instrukce abs, takže je nutné udělat pár instrukcí navíc, a SSSE3, kde se použije instrukce pabsd.

Neprováděl jsem nějaké velké testy, ale zdá se, že to počítá správně. Přeložit by to mělo být možné snad čímkoliv pro x86/amd64, zkoušel jsem jen MSVC.

[Thunder.m]

Příloha:

• main.cpp (2886 bytů)

Verze 2, trik s unsigned saturation je mnohem lepší, teoreticky lze vystačit s SSE2.

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

// Config.
#define USE_SSE2
#define USE_SSSE3

// SSE2.
#if defined(USE_SSE2)
#include <emmintrin.h>
#endif // USE_SSE2

// SSSE3.
#if defined(USE_SSSE3)
#include <tmmintrin.h>
#endif // USE_SSE3

#define ABS_C(_Value_) abs(_Value_)

int sum_abs_u16(
  const unsigned short* input1,
  const unsigned short* input2,
  size_t size)
{
  size_t i = size;
  int z = 0;

#if defined(USE_SSE2)
  if (i >= 20)
  {
    // Align.
    while ((intptr_t)input1 & (intptr_t)0xF)
    {
      z += ABS_C((int)input1[0] - (int)input2[0]);

      if (--i == 0) return z;
      input1++;
      input2++;
    }

    // Counter.
    __m128i cn = _mm_setzero_si128();
    __m128i zn = _mm_setzero_si128();

    // Large loop.
    while (i >= 16)
    {
      __m128i r0, r1;
      __m128i r2, r3;
      __m128i r4, r5;
      
      // Load input1, aligned.
      r0 = _mm_load_si128(reinterpret_cast<const __m128i*>(input1 + 0));
      r3 = _mm_load_si128(reinterpret_cast<const __m128i*>(input1 + 8));

      // Load input2, unaligned.
      r2 = _mm_loadu_si128(reinterpret_cast<const __m128i*>(input2 + 0));
      r4 = _mm_loadu_si128(reinterpret_cast<const __m128i*>(input2 + 8));

      // Get absolute values.
      r1 = _mm_subs_epu16(r2, r0);
      r0 = _mm_subs_epu16(r0, r2);

      r5 = _mm_subs_epu16(r3, r4);
      r4 = _mm_subs_epu16(r4, r3);

      r0 = _mm_add_epi16(r0, r1);
      r4 = _mm_add_epi16(r4, r5);

      // Unpack to 32-bit and sum.
      r1 = _mm_unpackhi_epi16(r0, zn);
      r5 = _mm_unpackhi_epi16(r4, zn);

      r0 = _mm_unpacklo_epi16(r0, zn);
      r4 = _mm_unpacklo_epi16(r4, zn);

      r0 = _mm_add_epi32(r0, r1);
      r4 = _mm_add_epi32(r4, r5);

      cn = _mm_add_epi32(cn, r0);
      cn = _mm_add_epi32(cn, r4);

      i -= 16;
      input1 += 16;
      input2 += 16;
    }

#if defined(USE_SSSE3)
    cn = _mm_hadd_epi32(cn, cn);
    cn = _mm_hadd_epi32(cn, cn);

    z += _mm_cvtsi128_si32(cn);
#else
    cn = _mm_add_epi32(cn, _mm_shuffle_epi32(cn, _MM_SHUFFLE(2, 3, 0, 1)));
    cn = _mm_add_epi32(cn, _mm_shuffle_epi32(cn, _MM_SHUFFLE(0, 1, 3, 2)));
    z += _mm_cvtsi128_si32(cn);
#endif
  }
#endif // USE_SSE2

  // Small loop.
  while (i > 0)
  {
    z += ABS_C((int)input1[0] - (int)input2[0]);

    i--;
    input1++;
    input2++;
  }

  return z;
}

int main(int argc, char* argv[])
{
  const unsigned short input1[40] = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 9, 8, 7, 6, 5, 4, 3, 2, 1, 0 };
  const unsigned short input2[40] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 9, 8, 7, 6, 5, 4, 3, 2, 1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 0 };

  int sum = sum_abs_u16(input1, input2, 40);
  printf("Sum=%d\n", sum);

  return 0;
}

[Thunder.m]

Tak tohle vypadá opravdu pěkně, s C sice umím, ale assembler jsem už dlouho nepoužíval a tak mi bude chvíli trvat pochopit jak to funguje, jakmile to implementuji poskytnu výsledky rychlosti.

[chochi]

Ono taky dost zalezi na prekladaci a jeho konfiguraci.
Podle popisku mas Core i7, ktere jiz podporuje SSSE3 (neplest si s SSE3), ktere jiz na funkci na absolutni fondnotu intu (pabsd).
Takze pokud se rekne gcc, ze muze pouzit SSSE3 (gcc -O3 -mssse3), tak vygeneruje jeste o kousek lepsi kod (protoze ma funkci abs).
Napr. tento kousek

  // Small loop.
  while (i > 0)
  {
    z += ABS_C((int)input1[0] - (int)input2[0]);
 
    i--;
    input1++;
    input2++;
  }

Prelozi nejak takhle:

	movl	8(%ebp), %ebx
	xorl	%edx, %edx
	xorl	%eax, %eax
	pxor	%xmm4, %xmm4
	pxor	%xmm5, %xmm5
	.p2align 4,,7
	.p2align 3
.L7:
	movdqu	(%ebx,%eax), %xmm2
	movdqu	(%edi,%eax), %xmm3
	movdqa	%xmm2, %xmm0
	movdqa	%xmm3, %xmm1
	punpcklwd	%xmm5, %xmm0
	punpcklwd	%xmm5, %xmm1
	addl	$1, %edx
	psubd	%xmm1, %xmm0
	punpckhwd	%xmm5, %xmm2
	pabsd	%xmm0, %xmm0
	punpckhwd	%xmm5, %xmm3
	paddd	%xmm4, %xmm0
	psubd	%xmm3, %xmm2
	addl	$16, %eax
	cmpl	%edx, %ecx
	pabsd	%xmm2, %xmm4
	paddd	%xmm0, %xmm4
	ja	.L7

+ nejaka ta omacko okolo.
Takze nekdy si prekladac poradi sam. To se asi nejlepe zjisti tak, ze si nechas prelozit kod do ASM a podivas se, zda to dostatecne optimalizoval.

[__dark__]

Jo, benchmark by byl super :)

V tom kódu je použitých jen pár triků. Pro portabilitu jsem zvolil makro USE_SSE2, které si definuješ, pokud bude cílová architektura mít SSE2. Je možné udělat i runtime check, ale to bych moc komplikoval.

Dále jsem využil toho, že neznaménková saturace US(A - B) dá buď kladný rozdíl těchto dvou hodnot nebo nulu (saturuje se, pokud je B větší než A), no takže abych získal absolutní hodnotu rodílu, použiju vlastně US(A - B) + US(B - A), výsledek není nikdy větší než 16-bitové číslo, takže využiju instrukci pro práci s 16bitovými páry hodnot. Na toto stačí pouhé SSE2 a je to podle mě asi to nejoptimálnější řešení (teda nenapadlo mě nic lepšího).

Dále tyto hodnoty převedu na 32bitové čísla (unpack), protože je potřebuju sčítat a 16bitová hodnota už není dostačující. Registr CN funguje jako paralelní akumulátor 4 hodnot, které se po vykonání hlavní smyčky sečtou a přičtou k výsledné proměnné Z.

Vychází mi to 22 instrukcí (+ instrukce na chod cyklu) pro akumulaci 16 hodnot. Problém bude spíš sběrnice než rychlost tohoto algoritmu. Pokud máš velké sady dat, které se nikdy nevlezou do cache, je možné udělat ještě další experimenty, ale ty můžou v neoptimálním případě vést i ke zpomalení.

(Na těch 22 instrukcí bych to napsal v asm, nevím zda to překladač zvládne stejně)