Perché GCC non ottimizza un * a * a * a * a * a a (a * a * a) * (a * a * a)?

Sto facendo un’ottimizzazione numerica su un’applicazione scientifica. Una cosa che ho notato è che GCC ottimizzerà il call pow(a,2) compilandolo in a*a , ma il call pow(a,6) non è ottimizzato e in realtà chiamerà la funzione di libreria pow , che rallenta notevolmente la prestazione. (Al contrario, Intel C ++ Compiler , eseguibile icc , eliminerà la chiamata alla libreria per pow(a,6) .)

Quello di cui sono curioso è che quando ho sostituito pow(a,6) con a*a*a*a*a*a usando GCC 4.5.1 e le opzioni ” -O3 -lm -funroll-loops -msse4 “, usa 5 istruzioni mulsd :

 movapd %xmm14, %xmm13 mulsd %xmm14, %xmm13 mulsd %xmm14, %xmm13 mulsd %xmm14, %xmm13 mulsd %xmm14, %xmm13 mulsd %xmm14, %xmm13 

mentre se scrivo (a*a*a)*(a*a*a) , produrrà

 movapd %xmm14, %xmm13 mulsd %xmm14, %xmm13 mulsd %xmm14, %xmm13 mulsd %xmm13, %xmm13 

che riduce il numero di istruzioni icc a 3. icc ha un comportamento simile.

Perché i compilatori non riconoscono questo trucco di ottimizzazione?

Perché il Floating Point Math non è associativo . Il modo in cui si raggruppano gli operandi in moltiplicazione in virgola mobile ha un effetto sull’accuratezza numerica della risposta.

Di conseguenza, molti compilatori sono molto prudenti nel riordinare i calcoli in virgola mobile a meno che non siano sicuri che la risposta rimarrà la stessa, oa meno che non dite loro che non vi importa dell’accuratezza numerica. Ad esempio: l’opzione -fassociative-math di gcc che consente a gcc di riassociare le operazioni in virgola mobile, o anche l’opzione -ffast-math che consente compromessi ancora più aggressivi di accuratezza e velocità.

Lambdageek sottolinea correttamente che poiché l’associatività non è valida per i numeri in virgola mobile, l’ottimizzazione di a*a*a*a*a*a a (a*a*a)*(a*a*a) può cambiare il valore. Questo è il motivo per cui non è consentito da C99 (a meno che non sia espressamente consentito dall’utente, tramite flag di compilazione o pragma). Generalmente, il presupposto è che il programmatore abbia scritto quello che ha fatto per una ragione, e il compilatore dovrebbe rispettarlo. Se vuoi (a*a*a)*(a*a*a) , scrivilo.

Questo può essere un dolore scrivere, però; perché il compilatore non può fare solo [ciò che consideri essere] la cosa giusta quando usi pow(a,6) ? Perché sarebbe la cosa sbagliata da fare. Su una piattaforma con una buona libreria matematica, pow(a,6) è significativamente più preciso di a*a*a*a*a*a o (a*a*a)*(a*a*a) . Solo per fornire alcuni dati, ho eseguito un piccolo esperimento sul mio Mac Pro, misurando il peggiore errore nella valutazione di ^ 6 per tutti i numeri a virgola mobile a precisione singola tra [1,2):

 worst relative error using powf(a, 6.f): 5.96e-08 worst relative error using (a*a*a)*(a*a*a): 2.94e-07 worst relative error using a*a*a*a*a*a: 2.58e-07 

L’uso di pow invece di un albero di moltiplicazione riduce l’errore associato a un fattore di 4 . I compilatori non dovrebbero (e generalmente non lo fanno) fare “ottimizzazioni” che aumentano l’errore a meno che l’utente non ne faccia una licenza (ad es. Via -ffast-math ).

Nota che GCC fornisce __builtin_powi(x,n) come alternativa a pow( ) , che dovrebbe generare un albero di moltiplicazione in linea. Usalo se vuoi scambiare la precisione con le prestazioni, ma non vuoi abilitare la matematica veloce.

Un altro caso simile: la maggior parte dei compilatori non ottimizzerà a + b + c + d a (a + b) + (c + d) (questa è un’ottimizzazione poiché la seconda espressione può essere meglio sottoposta a pipeline) e la valuta come data (es. come (((a + b) + c) + d) ). Anche questo è dovuto a casi angolari:

 float a = 1e35, b = 1e-5, c = -1e35, d = 1e-5; printf("%e %e\n", a + b + c + d, (a + b) + (c + d)); 

Questo produce 1.000000e-05 0.000000e+00

Fortran (progettato per il calcolo scientifico) ha un operatore di potenza integrato e, per quanto ne so, i compilatori di Fortran generalmente ottimizzeranno l’aumento dei poteri interi in modo simile a quello che descrivete. C / C ++ sfortunatamente non ha un power operator, solo la funzione libreria pow() . Ciò non impedisce ai compilatori intelligenti di trattare pow speciale e di calcolarlo in modo più rapido per casi speciali, ma sembra che lo facciano meno comunemente …

Alcuni anni fa stavo cercando di rendere più conveniente il calcolo dei poteri interi in modo ottimale, e ho scoperto quanto segue. È C ++, non C, e dipende comunque dal fatto che il compilatore sia piuttosto intelligente su come ottimizzare / integrare le cose. Comunque, spero che tu possa trovare utile nella pratica:

 template struct power_impl; template struct power_impl { template static T calc(const T &x) { if (N%2 == 0) return power_impl::calc(x*x); else if (N%3 == 0) return power_impl::calc(x*x*x); return power_impl::calc(x)*x; } }; template<> struct power_impl<0> { template static T calc(const T &) { return 1; } }; template inline T power(const T &x) { return power_impl::calc(x); } 

Chiarimento per i curiosi: questo non trova il modo ottimale per calcolare i poteri, ma dal momento che trovare la soluzione ottimale è un problema NP-completo e vale comunque la pena farlo per le piccole potenze (invece di usare il pow ), non c’è motivo di agitarsi con il dettaglio.

Quindi usalo solo come power<6>(a) .

Ciò semplifica la digitazione di potenze (non c’è bisogno di precisare 6 a se parens), e consente di avere questo tipo di ottimizzazione senza -ffast-math nel caso in cui si abbia qualcosa di preciso dipendente come una sum compensata (un esempio in cui l’ordine delle operazioni è essenziale).

Probabilmente si può anche dimenticare che questo è C ++ e basta usarlo nel programma C (se si compila con un compilatore C ++).

Spero che questo possa essere utile.

MODIFICARE:

Questo è ciò che ottengo dal mio compilatore:

Per a*a*a*a*a*a ,

  movapd %xmm1, %xmm0 mulsd %xmm1, %xmm0 mulsd %xmm1, %xmm0 mulsd %xmm1, %xmm0 mulsd %xmm1, %xmm0 mulsd %xmm1, %xmm0 

Per (a*a*a)*(a*a*a) ,

  movapd %xmm1, %xmm0 mulsd %xmm1, %xmm0 mulsd %xmm1, %xmm0 mulsd %xmm0, %xmm0 

Per il power<6>(a) ,

  mulsd %xmm0, %xmm0 movapd %xmm0, %xmm1 mulsd %xmm0, %xmm1 mulsd %xmm0, %xmm1 

Perché un numero a virgola mobile a 32 bit, come 1.024, non è 1.024. In un computer, 1.024 è un intervallo: da (1.024-e) a (1.024 + e), dove “e” rappresenta un errore. Alcune persone non si rendono conto di questo e credono anche che * in a * a sta per moltiplicazione di numeri arbitrari di precisione senza che vi siano errori associati a quei numeri. Il motivo per cui alcune persone non riescono a rendersi conto di ciò sono forse i calcoli matematici che hanno esercitato nelle scuole elementari: lavorare solo con numeri ideali senza errori, e credere che sia giusto ignorare semplicemente “e” mentre si esegue la moltiplicazione. Non vedono la “e” implicita in “float a = 1.2”, “a * a * a” e codici C simili.

Se la maggioranza dei programmatori riconosce (ed è in grado di eseguire) l’idea che l’espressione C a * a * a * a * a * a non sta effettivamente lavorando con numeri ideali, il compilatore GCC sarebbe quindi GRATUITO per ottimizzare “a * a * a * a * a * a “in dire” t = (a * a); t * t * t “che richiede un numero minore di moltiplicazioni. Ma sfortunatamente, il compilatore GCC non sa se il programmatore che scrive il codice pensa che “a” sia un numero con o senza un errore. E così GCC farà solo quello che sembra il codice sorgente – perché è ciò che GCC vede con il suo “occhio nudo”.

… una volta che sai che tipo di programmatore sei , puoi usare l’opzione “-ffast-math” per dire a GCC che “Ehi, GCC, so cosa sto facendo!”. Ciò consentirà a GCC di convertire un * a * a * a * a * a in un altro pezzo di testo – sembra diverso da un * a * a * a * a * a – ma calcola ancora un numero nell’intervallo di errore di a * a * a * a * a * a. Questo è OK, dal momento che sai già che stai lavorando con intervalli, non numeri ideali.

GCC in realtà ottimizza un * a * a * a * a * a a (a * a * a) * (a * a * a) quando a è un numero intero. Ho provato con questo comando:

 $ echo 'int f(int x) { return x*x*x*x*x*x; }' | gcc -o - -O2 -S -masm=intel -xc - 

Ci sono molte bandiere gcc ma niente di speciale. Significa: leggi da stdin; utilizzare il livello di ottimizzazione O2; elenco linguistico assembly di output anziché binario; la lista dovrebbe usare la syntax del linguaggio assembly assembly; l’input è in linguaggio C (di solito la lingua è dedotta dall’estensione del file di input, ma non c’è estensione del file durante la lettura da stdin); e scrivere su stdout.

Ecco la parte importante dell’output. L’ho annotato con alcuni commenti che indicano cosa sta succedendo nel linguaggio assembly:

  ; x is in edi to begin with. eax will be used as a temporary register. mov eax, edi ; temp1 = x imul eax, edi ; temp2 = x * temp1 imul eax, edi ; temp3 = x * temp2 imul eax, eax ; temp4 = temp3 * temp3 

Sto usando il sistema GCC su Linux Mint 16 Petra, un derivato di Ubuntu. Ecco la versione di gcc:

 $ gcc --version gcc (Ubuntu/Linaro 4.8.1-10ubuntu9) 4.8.1 

Come hanno notato altri poster, questa opzione non è ansible in virgola mobile, perché l’aritmetica in virgola mobile non è in realtà associativa.

Nessun poster ha menzionato la contrazione delle espressioni fluttuanti ancora (ISO C standard, 6.5p8 e 7.12.2). Se il FP_CONTRACT è impostato su ON , il compilatore può considerare un’espressione come a*a*a*a*a*a come singola operazione, come se fosse valutata esattamente con un singolo arrotondamento. Ad esempio, un compilatore può sostituirlo con una funzione di alimentazione interna che è sia più veloce che più accurata. Ciò è particolarmente interessante in quanto il comportamento è parzialmente controllato dal programmatore direttamente nel codice sorgente, mentre le opzioni del compilatore fornite dall’utente finale possono talvolta essere utilizzate in modo errato.

Lo stato predefinito del pragma FP_CONTRACT è definito dall’implementazione, in modo che un compilatore possa fare tali ottimizzazioni per impostazione predefinita. Quindi il codice portatile che deve seguire rigorosamente le regole IEEE 754 dovrebbe impostarlo esplicitamente su OFF .

Se un compilatore non supporta questo pragma, deve essere prudente evitando tale ottimizzazione, nel caso in cui lo sviluppatore abbia scelto di impostarlo su OFF .

GCC non supporta questo pragma, ma con le opzioni predefinite, si presuppone che sia ON ; quindi per gli obiettivi con una FMA hardware, se si desidera impedire la trasformazione a*b+c a fma (a, b, c), è necessario fornire un’opzione come -ffp-contract=off (per impostare esplicitamente il pragma a OFF ) o -std=c99 (per indicare a GCC di conformarsi ad una versione standard C, qui C99, quindi seguire il paragrafo precedente). In passato, quest’ultima opzione non impediva la trasformazione, il che significa che GCC non si conformava su questo punto: https://gcc.gnu.org/bugzilla/show_bug.cgi?id=37845

Come ha sottolineato Lambdageek, la moltiplicazione float non è associativa e puoi ottenere meno precisione, ma anche quando ottieni una maggiore precisione puoi argomentare contro l’ottimizzazione, perché vuoi un’applicazione deterministica. Ad esempio nel client / server di simulazione di gioco, in cui ogni cliente deve simulare lo stesso mondo, si desidera che i calcoli in virgola mobile siano deterministici.

Non mi sarei mai aspettato che questo caso fosse ottimizzato. Non può essere molto spesso dove un’espressione contiene sottoespressioni che possono essere raggruppate per rimuovere intere operazioni. Mi aspetterei che gli scrittori di compilatori investano il loro tempo in aree che sarebbero più probabili produrre miglioramenti evidenti, piuttosto che coprire un caso limite raramente riscontrato.

Sono stato sorpreso di apprendere dalle altre risposte che questa espressione poteva essere ottimizzata con gli appositi switch del compilatore. O l’ottimizzazione è banale, o è un caso limite di un’ottimizzazione molto più comune, oppure gli scrittori del compilatore erano estremamente accurati.

Non c’è niente di sbagliato nel fornire suggerimenti al compilatore come hai fatto qui. È una parte normale e attesa del processo di micro-ottimizzazione per riorganizzare le dichiarazioni e le espressioni per vedere quali differenze porteranno.

Mentre il compilatore può essere giustificato nel considerare le due espressioni per fornire risultati incoerenti (senza gli interruttori appropriati), non è necessario che tu sia vincolato da tale restrizione. La differenza sarà incredibilmente piccola – tant’è vero che se la differenza è importante per te, non dovresti utilizzare l’aritmetica in virgola mobile standard in primo luogo.

Le funzioni di libreria come “pow” sono solitamente elaborate con cura per ottenere il minimo errore ansible (nel caso generico). Solitamente si ottengono funzioni di approssimazione con spline (secondo il commento di Pascal l’implementazione più comune sembra essere l’ algoritmo Remez )

fondamentalmente la seguente operazione:

 pow(x,y); 

ha un errore intrinseco di circa la stessa grandezza dell’errore in ogni singola moltiplicazione o divisione .

Mentre la seguente operazione:

 float a=someValue; float b=a*a*a*a*a*a; 

ha un errore intrinseco maggiore di 5 volte l’errore di una singola moltiplicazione o divisione (perché si combinano 5 moltiplicazioni).

Il compilatore dovrebbe essere molto attento al tipo di ottimizzazione che sta facendo:

  1. se l’ottimizzazione di pow(a,6) su a*a*a*a*a*a può migliorare le prestazioni, ma ridurre drasticamente la precisione per i numeri in virgola mobile.
  2. se si ottimizza a*a*a*a*a*a a pow(a,6) si può effettivamente ridurre l’accuratezza perché “a” era un valore speciale che consente la moltiplicazione senza errore (una potenza di 2 o qualche piccolo numero intero)
  3. se si ottimizza la pow(a,6) in (a*a*a)*(a*a*a) o (a*a)*(a*a)*(a*a) può ancora esserci una perdita di precisione rispetto alla funzione pow .

In generale, sai che per i valori in virgola mobile arbitrari “pow” ha una precisione migliore di qualsiasi funzione tu possa eventualmente scrivere, ma in alcuni casi speciali moltiplicazioni multiple possono avere una precisione e prestazioni migliori, spetta allo sviluppatore scegliere ciò che è più appropriato, eventualmente commentando il codice in modo che nessun altro possa “ottimizzare” quel codice.

L’unica cosa che ha senso (l’opinione personale, e apparentemente una scelta in GCC senza una particolare ottimizzazione o flag di compilazione) da ottimizzare dovrebbe essere la sostituzione di “pow (a, 2)” con “a * a”. Questa sarebbe l’unica cosa sensata che un venditore di compilatori dovrebbe fare.

Ci sono già alcune buone risposte a questa domanda, ma per completezza volevo sottolineare che la sezione applicabile dello standard C è 5.1.2.2.3 / 15 (che è la stessa della sezione 1.9 / 9 nel C ++ 11 standard). Questa sezione afferma che gli operatori possono essere raggruppati solo se sono realmente associativi o commutativi.

gcc in realtà può fare questa ottimizzazione, anche per i numeri in virgola mobile. Per esempio,

 double foo(double a) { return a*a*a*a*a*a; } 

diventa

 foo(double): mulsd %xmm0, %xmm0 movapd %xmm0, %xmm1 mulsd %xmm0, %xmm1 mulsd %xmm1, %xmm0 ret 

con -O -funsafe-math-optimizations . Questo riordino viola tuttavia IEEE-754, quindi richiede la bandiera.

Gli interi firmati, come sottolineato da Peter Cordes in un commento, possono eseguire questa ottimizzazione senza -funsafe-math-optimizations poiché contiene esattamente quando non c’è overflow e se c’è overflow si ottiene un comportamento indefinito. Quindi ottieni

 foo(long): movq %rdi, %rax imulq %rdi, %rax imulq %rdi, %rax imulq %rax, %rax ret 

con solo -O . Per gli interi non firmati, è ancora più semplice dal momento che funzionano con i poteri mod di 2 e quindi possono essere riordinati liberamente anche in caso di overflow.