Algoritmo più veloce per identificare la x più piccola e più grande che crea l'equazione a doppia precisione x + a == b vero

Question

Nel contesto dell'analisi statica, sono interessato a determinare i valori di x nell'allora succursale del condizionale sotto:

double x; 
x = …; 
if (x + a == b) 
{ 
    … 

a e b può essere assunto come costanti doppia precisione (generalizzando alle espressioni arbitrarie è la parte più facile del problema), e il compilatore può presumere seguire IEEE 754 strettamente (FLT_EVAL_METHOD è 0). Si può presumere che la modalità di arrotondamento in fase di esecuzione sia prossima a pari.

Se il calcolo con i costi era economico, sarebbe semplice: i valori per x sarebbero i numeri a precisione doppia contenuti nell'intervallo razionale (b - a - 0.5 * ulp1 (b) ... b - a + 0.5 * ulp2 (b)). I limiti dovrebbero essere inclusi se b è pari, escluso se b è dispari, e ulp1 e ulp2 sono due definizioni leggermente diverse di "ULP" che possono essere prese identiche se non si preoccupa di perdere un po 'di precisione sulle potenze di due.

Sfortunatamente, il calcolo con i costi può essere costoso. Si consideri che un'altra possibilità è ottenere ciascuno dei limiti per dicotomia, in 64 aggiunte a doppia precisione (ciascuna operazione che decide un bit del risultato). 128 aggiunte in virgola mobile per ottenere i limiti inferiore e superiore potrebbero essere più veloci di qualsiasi soluzione basata sulla matematica.

Mi chiedo se c'è un modo per migliorare l'idea di "128 floating-point additions". In realtà, ho una mia soluzione che comporta cambiamenti della modalità di arrotondamento e chiamate nextafter, ma non vorrei ostacolare lo stile di nessuno e far sì che perdano una soluzione più elegante di quella che attualmente possiedo. Inoltre, non sono sicuro che la modifica della modalità di arrotondamento due volte sia effettivamente più economica di 64 aggiunte in virgola mobile.

Answer 1

È già dato una soluzione bella ed elegante nella sua domanda:

Se informatica con razionali era a buon mercato, sarebbe semplice: i valori per x sarebbero i numeri a precisione doppia contenute nel razionale Intervallo (b - a - 0,5 * ulp1 (b) ... b - a + 0,5 * ulp2 (b)). I limiti devono essere inclusi se b è anche, esclusa se B è dispari, e ulp1 e ulp2 sono due definizioni leggermente diverse di “ULP” che può essere preso identici se uno non dispiace perdere un po 'di precisione sui poteri di Due.

Quello che segue è uno schizzo semidiretto di una soluzione parziale al problema basato su questo paragrafo. Spero che avrò la possibilità di farcela presto. Per ottenere una soluzione reale, dovrai gestire i subnormali, gli zeri, i NaN e tutte le altre cose divertenti. Assumerò che a e b sono, ad esempio, tali che 1e-300 < |a| < 1e300 e 1e-300 < |b| < 1e300 in modo che non si verifichi alcuna follia in nessun punto.

Assenza di overflow e underflow, è possibile ottenere ulp1(b) da b - nextafter(b, -1.0/0.0). È possibile ottenere ulp2(b) da nextafter(b, 1.0/0.0) - b.

Se b/2 <= a <= 2b, il teorema di Sterbenz ti dice che b - a è esatto. Quindi corrisponde al valore double più vicino al limite inferiore e double al limite superiore è double più vicino.Prova questi valori, e i valori immediatamente prima e dopo, e scegli l'intervallo più ampio che funzioni.

Se b > 2a, b - a > b/2. Il valore calcolato di b - a è al massimo di mezzo ulp. Uno ulp1 è al massimo due ulp, come lo è uno ulp2, quindi l'intervallo razionale che hai dato è al massimo di due ulp. Scopri quale dei cinque valori più vicini a b-a funziona.

Se a > 2b, un numero di b-a è almeno pari a un massimo di b; se qualcosa funziona, scommetto che dovrà essere tra i tre valori più vicini a b-a. Immagino che il caso in cui a e b abbiano segni diversi che funzionano allo stesso modo.

Ho scritto una piccola pila di codice C++ che implementa questa idea. Non ha fallito il test fuzz casuale (in alcune gamme diverse) prima di annoiarmi ad aspettare. Eccolo:

void addeq_range(double a, double b, double &xlo, double &xhi) { 
    if (a != a) return; // empty interval 
    if (b != b) { 
    if (a-a != 0) { xlo = xhi = -a; return; } 
    else return; // empty interval 
    } 
    if (b-b != 0) { 
    // TODO: handle me. 
    } 

    // b is now guaranteed to be finite. 
    if (a-a != 0) return; // empty interval 

    if (b < 0) { 
    addeq_range(-a, -b, xlo, xhi); 
    xlo = -xlo; 
    xhi = -xhi; 
    return; 
    } 

    // b is now guaranteed to be zero or positive finite and a is finite. 
    if (a >= b/2 && a <= 2*b) { 
    double upulp = nextafter(b, 1.0/0.0) - b; 
    double downulp = b - nextafter(b, -1.0/0.0); 
    xlo = (b-a) - downulp/2; 
    xhi = (b-a) + upulp/2; 
    if (xlo + a == b) { 
     xlo = nextafter(xlo, -1.0/0.0); 
     if (xlo + a != b) xlo = nextafter(xlo, 1.0/0.0); 
    } else xlo = nextafter(xlo, 1.0/0.0); 
    if (xhi + a == b) { 
     xhi = nextafter(xhi, 1.0/0.0); 
     if (xhi + a != b) xhi = nextafter(xhi, -1.0/0.0); 
    } else xhi = nextafter(xhi, -1.0/0.0); 
    } else { 
    double xmid = b-a; 
    if (xmid + a < b) { 
     xhi = xlo = nextafter(xmid, 1.0/0.0); 
     if (xhi + a != b) xhi = xmid; 
    } else if (xmid + a == b) { 
     xlo = nextafter(xmid, -1.0/0.0); 
     xhi = nextafter(xmid, 1.0/0.0); 
     if (xlo + a != b) xlo = xmid; 
     if (xhi + a != b) xhi = xmid; 
    } else { 
     xlo = xhi = nextafter(xmid, -1.0/0.0); 
     if (xlo + a != b) xlo = xmid; 
    } 
    } 
}