come dividere piccoli numeri di doppia precisione correttamente senza errori di precisione?

Question

Sto cercando di diagnosticare e correggere un bug, che si riduce a X/Y ottenendo un risultato instabile quando X e Y sono piccole:

In questo caso, sia cx e patharea aumentano senza intoppi. Il loro rapporto è un asintoto liscio ad alti numeri, ma irregolare per numeri "piccoli". L'ovvio primo pensiero è che stiamo raggiungendo il limite della precisione in virgola mobile, ma i numeri reali stessi non sono neanche lontanamente vicini. ActionScript "numero" tipi sono IEE 754 carri a doppia precisione, in modo da dovrebbero avere 15 cifre decimali di precisione (se ho letto bene).

Alcuni valori tipici del denominatore (patharea):

0.0000000002119123 
0.0000000002137313 
0.0000000002137313 
0.0000000002155502 
0.0000000002182787 
0.0000000002200977 
0.0000000002210072 

E numeratore (cx):

0.0000000922932995 
0.0000000930474444 
0.0000000930582124 
0.0000000938123574 
0.0000000950458711 
0.0000000958000159 
0.0000000962901528 
0.0000000970442977 
0.0000000977984426 

Ciascuno di questi cresce monotonicamente, ma il rapporto è caotico come visto sopra.

A numeri più grandi si assesta ad un'iperbole liscia.

Quindi, la mia domanda: qual è il modo corretto per gestire numeri molto piccoli quando è necessario dividerli uno con l'altro?

Ho pensato di moltiplicare numeratore e/o denominatore per 1000 in anticipo, ma non riuscivo a risolverlo.

Il codice effettivo in questione è la funzione recalculate()here. Calcola il baricentro di un poligono, ma quando il poligono è piccolo, il baricentro salta irregolare intorno al luogo, e può finire per una lunga distanza dal poligono. La serie di dati di cui sopra sono il risultato di spostare un nodo del poligono in una direzione costante (a mano, per questo motivo non è perfettamente liscia).

Questo è Adobe Flex 4.5.

Answer 1

Credo che il problema più probabile è causato dalla seguente riga nel codice:

sc = (lx*latp-lon*ly)*paint.map.scalefactor; 

Se il poligono è molto piccolo, quindi lx e lon sono quasi la stessa, come lo sono ly e latp. Sono entrambi molto grandi rispetto al risultato, quindi stai sottraendo due numeri che sono quasi uguali.

Per aggirare il problema, possiamo fare uso del fatto che:

x1*y2-x2*y1 = (x2+(x1-x2))*y2 - x2*(y2+(y1-y2)) 
      = x2*y2 + (x1-x2)*y2 - x2*y2 - x2*(y2-y1) 
      = (x1-x2)*y2 - x2*(y2-y1) 

Quindi, provate questo:

dlon = lx - lon 
dlat = ly - latp 
sc = (dlon*latp-lon*dlat)*paint.map.scalefactor; 

Il valore è matematicamente la stessa, ma i termini sono un ordine di grandezza minore, quindi l'errore dovrebbe essere anche un ordine di grandezza più piccolo.

Answer 2

Jeffrey Sax ha identificato correttamente il problema di base: perdita di precisione dalla combinazione di termini che sono (molto) più grandi del risultato finale. Il suggerito riscrittura elimina una parte del problema - a quanto pare sufficiente per il caso concreto, data la risposta felice.

È possibile che, tuttavia, se il poligono torna a essere (molto) più piccolo e/o più lontano dall'origine, l'imprecisione verrà nuovamente visualizzata. Nella formula riscritta i termini sono ancora un po 'più grandi della loro differenza.

Inoltre, nell'algoritmo esiste un altro problema di combinazione di numeri comparabili con segni diversi di &. I vari valori "sc" nei cicli successivi dell'iterazione sopra i bordi del poligono si combinano efficacemente in un numero finale che è (molto) più piccolo del singolo sc (i). (se hai un poligono convesso scoprirai che esiste una sequenza contigua di valori positivi e una sequenza contigua di valori negativi, nei poligoni non convessi i negativi e i positivi possono essere intrecciati).

L'algoritmo sta calcolando l'area del poligono aggiungendo aree di triangoli attraversate dai bordi e dall'origine, dove alcuni dei termini sono negativi (ogni volta che un bordo viene attraversato in senso orario, osservandolo da l'origine) e alcuni positivi (camminata antioraria sul bordo).

Si elimina TUTTI i problemi di perdita di precisione definendo l'origine in uno degli angoli del poligono, ad esempio (lx, ly) e quindi aggiungendo le superfici triangolari suddivise dai bordi e da quell'angolo (così: trasformando lon in (lon-lx) e latp in (latp-ly) - con il bonus aggiuntivo che è necessario elaborare due triangoli in meno, perché ovviamente i bordi che si collegano all'angolo di origine scelto producono zero superfici

Per l'area-parte è tutto.Per la parte centroide, ovviamente dovrai "trasformare indietro" il risultato nel frame originale, cioè aggiungere (lx, ly) alla fine