Quale confronto in virgola mobile è più accurato e perché?

Question

Sto sperimentando diverse implementazioni del metodo Newton per il calcolo delle radici quadrate. Una decisione importante è quando terminare l'algoritmo.

Ovviamente non farà utilizzare la differenza assoluta tra y*y e x dove y è la stima corrente della radice quadrata di x, perché per grandi valori di x potrebbe non essere possibile rappresentare la sua radice quadrata con abbastanza precisione.

Quindi dovrei usare un criterio relativo. Ingenuamente avrei usato qualcosa del genere:

static int sqrt_good_enough(float x, float y) { 
    return fabsf(y*y - x)/x < EPS; 
} 

E sembra funzionare molto bene. Ma di recente ho iniziato a leggere Kernighan e di Plauger The Elements of Programming Style e danno un programma Fortran per lo stesso algoritmo nel capitolo 1, i cui criteri di terminazione, tradotto in C, sarebbe:

static int sqrt_good_enough(float x, float y) { 
    return fabsf(x/y - y) < EPS * y; 
} 

Entrambi sono matematicamente equivalenti , ma esiste un motivo per preferire un modulo rispetto all'altro?

Answer 1

I due in realtà non sono esattamente equivalenti matematicamente, a meno che non si scriva fabsf (y * y - x)/(y * y) < EPS per il primo. (scusate l'errore di battitura nel mio commento originale)

Ma penso che il punto chiave sia far sì che l'espressione qui corrisponda alla vostra formula per calcolare y nell'iterazione di Newton. Ad esempio se la tua formula y è y = (y + x/y)/2, dovresti usare lo stile di Kernighan e Plauger. Se è y = (y * y + x)/(2 * y) dovresti usare (y * y - x)/(y * y) < EPS.

In genere, il criterio di terminazione dovrebbe essere che abs (y (n + 1) - y (n)) sia sufficientemente piccolo (cioè più piccolo di y (n + 1) * EPS). Questo è il motivo per cui le due espressioni devono corrispondere. Se non corrispondono esattamente, è possibile che il test di terminazione decida che il residuo non è abbastanza piccolo mentre la differenza in y (n) è inferiore all'errore in virgola mobile, a causa del diverso ridimensionamento. Il risultato sarebbe un ciclo infinito perché y (n) ha smesso di cambiare e il criterio di terminazione non è mai stato soddisfatto.

Ad esempio, il seguente codice Matlab è esattamente lo stesso Newton risolutore come primo esempio, ma funziona sempre:

x = 6.800000000000002 
yprev = 0 
y = 2 
while abs(y*y - x) > eps*abs(y*y) 
    yprev = y; 
    y = 0.5*(y + x/y); 
end 

La versione C/C++ di esso ha lo stesso problema.

Answer 2

Non sono ancora equivalenti; quello inferiore è matematicamente equivalente a fabsf(y*y - x)/(y*y) < EPS. Il problema che vedo con il vostro è che se overflow y*y (probabilmente perché è FLT_MAX e viene scelto sfortunatamente), quindi la risoluzione potrebbe non verificarsi mai. La seguente interazione usa il doppio.

>>> import math 
>>> x = (2.0 - 2.0 ** -52) * 2.0 ** 1023 
>>> y = x/math.sqrt(x) 
>>> y * y - x 
inf 
>>> y == 0.5 * (y + x/y) 
True 

EDIT: come commento (ora cancellato) ha sottolineato, è anche un bene ad operazioni tra l'iterazione e il test di terminazione.

EDIT2: entrambi hanno probabilmente problemi con subnormale x. Il professionals normalizza x per evitare le complicazioni di entrambi gli estremi.