somma array di doppi con intervallo di grande valore: algoritmo corretto

Question

Ho un algoritmo in cui ho bisogno di sommare (un sacco di volte) numeri doppi che vanno nella e-40 alla e + 40.

Array Esempio (scaricato casualmente dall'applicazione reale):

-2.06991e-05 
7.58132e-06 
-3.91367e-06 
7.38921e-07 
-5.33143e-09 
-4.13195e-11 
4.01724e-14 
6.03221e-17 
-4.4202e-20 
6.58873 
-1.22257 
-0.0606178 
0.00036508 
2.67599e-07 
0 
-627.061 
-59.048 
5.92985 
0.0885884 
0.000276455 
-2.02579e-07 

va da sé l'Sono consapevole degli effetti arrotondamento questo causerà, sto cercando di tenerla sotto controllo: il risultato finale dovrebbe non avere alcuna informazione mancante nella parte frazionaria del doppio o, se non evitabile, il risultato dovrebbe essere almeno accurato di n cifre (con n definito). Il risultato finale richiede qualcosa come 5 cifre più esponente.

Dopo qualche pensiero decente, ho finito con seguente algoritmo:

• ordinare l'array in modo che il più grande valore assoluto viene prima, più vicina a zero scorso.
• Aggiungi tutto in un ciclo

L'idea è che in questo caso, qualsiasi annullamento di grandi valori (negativi e positivi) non avrà un impatto valori minori quest'ultimo. In breve:

• (10e40 - 10e40) + 1 = 1: risultato è come previsto
• (1 + 10e40) - 10e40 = 0: non buona

ho finito utilizzando std :: multiset (punto di riferimento sul mio PC ha dato velocità del 20% superiore con lunghi il doppio rispetto ai normali doppie - mi trovo benissimo con risoluzione doppie) con una funzione di ordinamento personalizzata utilizzando STD: FAB.

È ancora piuttosto lento (ci vogliono 5 secondi per fare tutto) e ho ancora la sensazione di "hai perso qualcosa nel tuo algo". Qualsiasi raccomandazione:

1. per l'ottimizzazione della velocità. C'è un modo migliore per ordinare i prodotti intermedi? L'ordinamento di un set di 40 risultati intermedi (in genere) richiede circa il 70% del tempo totale di esecuzione.
2. per problemi mancati. C'è la possibilità di perdere ancora i dati critici (uno che avrebbe dovuto essere nella parte frazionaria del risultato finale)?

---

Su un quadro più grande, io sono l'attuazione reale coefficiente classi polinomiali di variabile immaginaria pura (impedenze elettriche: Z (JW)). Z è un grande polinom che rappresenta un sistema definito dall'utente, con un esponente di coefficiente che va molto lontano.
Il "grande" deriva dall'aggiunta di cose come Zc1 = 1/jC1w a Zc2 = 1/jC2w:
Zc1 + Zc2 = (C1C2 (jw)^2 + 0 (jw))/(C1 + C2) (jw)

In questo caso, con C1 e C2 in nanofarad (10e-9), C1C2 è già in 10e-18 (e ha iniziato solo ...)

mia funzione di ordinamento utilizzare una distanza di Manhattan di variabili complesse (perché, i miei sono in puro reale o immaginario puro):

struct manhattan_complex_distance 
{ 
     bool operator() (std::complex<long double> a, std::complex<long double> b) 
     { 
      return std::fabs(std::real(a) + std::imag(a)) > std::fabs(std::real(b) + std::imag(b)); 
     } 
}; 

e il mio multi set in azione:

std:complex<long double> get_value(std::vector<std::complex<long double>>& frequency_vector) 
{ 
    //frequency_vector is precalculated once for all to have at index n the value (jw)^n. 
    std::multiset<std::complex<long double>, manhattan_distance> temp_list; 
    for (int i=0; i<m_coeficients.size(); ++i) 
    { 
     // element of :  ℝ   *   ℂ 
     temp_list.insert(m_coeficients[i] * frequency_vector[i]); 
    } 
    std::complex<long double> ret=0; 
    for (auto i:temp_list) 
    { 
     // it is VERY important to start adding the big values before adding the small ones. 
     // in informatics, 10^60 - 10^60 + 1 = 1; while 1 + 10^60 - 10^60 = 0. Of course you'd expected to get 1, not 0. 
     ret += i; 
    } 
    return ret; 
} 

Il il progetto I ha è abilitato per C++ 11 (principalmente per il miglioramento degli strumenti matematici di lib e complessi)

ps: ho rifattorizzato il codice per renderlo di facile lettura, in realtà tutti i complessi e nomi doppi lunghi sono template: posso cambiare il tipo polinomiale in poco tempo o utilizzare la classe per polinomiale regolare di ℝ

Answer 1

Come GuyGreer suggerito, è possibile utilizzare Kahan summation:

double sum = 0.0; 
double c = 0.0; 
for (double value : values) { 
    double y = value - c; 
    double t = sum + y; 
    c = (t - sum) - y; 
    sum = t; 
} 

EDIT: Si dovrebbe anche prendere in considerazione utilizzando Horner's method per valutare il polinomio.

double value = coeffs[degree]; 
for (auto i = degree; i-- > 0;) { 
    value *= x; 
    value += coeffs[i]; 
} 

Answer 2

L'ordinamento dei dati è sulla strada giusta. Ma sicuramente dovresti riassumere dalla più piccola grandezza alla più grande, non dalla più grande alla più piccola. Sommando dal più grande al più piccolo, quando si arriva al più piccolo, allineando il valore successivo con la somma corrente è suscettibile di causare la maggior parte o tutti i bit del valore successivo a "cadere alla fine". Sommando invece dal più piccolo al più grande, i valori più piccoli hanno la possibilità di accumulare una somma di dimensioni decenti, per la quale più bit entreranno nel più grande. Combinato con la somma di Kahan, dovrebbe produrre una somma abbastanza precisa.

Answer 3

Primo: fai in modo che la matematica tenga traccia dell'errore. Sostituisci i tuoi doppi con i tipi che riconoscono gli errori e quando aggiungi o moltiplichi insieme due doppi calcola anche l'errore massimo .

Questo è l'unico modo per garantire che il codice produca risultati accurati pur essendo ragionevolmente veloce.

In secondo luogo, non utilizzare multiset. I contenitori associativi non sono per l'ordinamento, sono per mantenendo una raccolta ordinata, pur essendo in grado di aggiungere o rimuovere in modo incrementale elementi da esso in modo efficiente.

La possibilità di aggiungere/rimuovere elementi in modo incrementale significa che è basata su nodi e basata su nodi significa che è lenta in generale.

Se si desidera semplicemente una raccolta ordinata, iniziare con vector quindi std::sort.

Quindi, per ridurre al minimo l'errore, mantenere un elenco di elementi positivi e negativi. Inizia con zero come somma. Ora scegli il più piccolo tra gli elementi positivi o negativi in ​​modo tale che il totale della somma e quell'elemento sia più vicino allo zero.

Farlo con elementi che calcolano i limiti di errore.

Al termine, determinare se si dispone di 5 cifre di precisione, oppure no.

Questi doppi di propagazione degli errori dovrebbero essere utilizzati idealmente all'inizio dell'algoritmo.