Come verificare le dipendenze dei galleggianti

Question

Voglio determinare (in C++) se un numero a virgola mobile è l'inverso moltiplicativo di un altro numero a virgola mobile. Il problema è che devo usare una terza variabile per farlo. Per esempio questo codice:

float x=5,y=0.2; 
if(x==(1/y)) cout<<"They are the multiplicative inverse of eachother"<<endl; 
else cout<<"They are NOT the multiplicative inverse of eachother"<<endl; 

uscita volontà: "non sono ...", che è sbagliato e questo codice:

float x=5,y=0.2,z; 
z=1/y; 
if(x==z) cout<<"They are the multiplicative inverse of eachother"<<endl; 
else cout<<"They are NOT the multiplicative inverse of eachother"<<endl; 

uscita volontà: "sono ...", che è di destra .
perché sta succedendo questo?

Answer 1

The Float Precision problema

Hai due problemi qui, ma entrambi provengono dalla stessa radice

Non si può paragonare con precisione galleggianti. Non puoi sottrarre o dividere in modo preciso. Non puoi contare esattamente per loro. Qualsiasi operazione con loro potrebbe (e quasi sempre fa) portare qualche errore nel risultato. Anche a=0.2f non è un'operazione precisa. Le ragioni più profonde di ciò sono spiegate molto bene dagli autori delle altre risposte qui. _{(I miei ringraziamenti e voti a loro per quello.)}

Ecco il primo e più semplice errore. Si dovrebbe mai, mai , mai, mai, MAI uso su di loro == o il suo equivalente in qualsiasi lingua.

Invece di a==b, utilizzare invece Abs(a-b)<HighestPossibleError.

---

Ma questo non è l'unico problema nel vostro compito.

Abs(1/y-x)<HighestPossibleError non funziona. Almeno, non funzionerà abbastanza spesso. Perché?

Prendiamo la coppia x = 1000 ey = 0,001. Prendiamo l'errore relativo "iniziale" di y per 10 ^-6.

_{(Errore relativo = errore/valore).}

Gli errori relativi dei valori si aggiungono a moltiplicazione e divisione.

1/y è circa 1000. Il suo errore relativo è lo stesso 10 ^-6. ("1" non ha errori)

Ciò rende l'errore assoluto = 1000 * 10 ^-6 = 0,001. Quando sottrassi x più tardi, quell'errore sarà tutto ciò che rimane. (Gli errori assoluti si aggiungono all'aggiunta e alla sottrazione e l'errore di x è trascurabilmente piccolo.) Sicuramente, non contano su così grandi errori, HighestPossibleError sarebbe sicuramente impostata più bassa e il programma sarebbe buttare fuori un buon paio di x, y

Così, i due successivi regola per le operazioni galleggiante: non cercate di dividere una maggiore valuer by lesser one e Dio ti salvi di sottrarre i valori vicini dopo.

Esistono due modi semplici per evitare questo problema.

• Fondando quello di x, y ha il maggior valore addominali e dividere 1 dalla maggiore uno e solo successivamente per sottrarre il minore uno.

• Se si desidera confrontare 1/y against x, mentre si sta lavorando ancora con le lettere, non i valori, e le operazioni non fare errori, si moltiplicano entrambi i lati del confronto da y e hai 1 against x*y. _{(Di solito dovresti controllare i segni in quell'operazione, ma qui usiamo valori abs, quindi, è pulito.)} Il confronto dei risultati non ha alcuna divisione.

In una via più breve:

1/y V x <=> y*(1/y) V x*y <=> 1 V x*y 

Sappiamo già che tale confronto come 1 against x*y dovrebbe essere fatto in modo da:

const float HighestPossibleError=1e-10; 
if(Abs(x*y-1.0)<HighestPossibleError){... 

Questo è tutto.

---

P.S. Se si ha realmente bisogno tutto su una riga, utilizzare:

if(Abs(x*y-1.0)<1e-10){... 

Ma è cattivo stile. Non lo consiglierei.

P.P.S. Nel secondo esempio il compilatore ottimizza il codice in modo tale da impostare z a 5 prima di eseguire qualsiasi codice. Quindi, il controllo di 5 contro 5 funziona anche per i galleggianti.

Answer 2

Dovrai definire con precisione cosa significhi per due approssimazioni inversioni moltiplicative. Altrimenti, non saprai che cosa dovresti testare.

0.2 non ha una rappresentazione binaria esatta. Se memorizzi numeri che non hanno una rappresentazione esatta con precisione limitata, non otterrai risposte che siano esattamente corrette.

Le stesse cose accade in decimale. Ad esempio, 1/3 non ha una rappresentazione decimale esatta. È possibile memorizzarlo come .333333. Ma poi hai un problema. Sono gli inversi moltiplicativi 3 e .333333? Se li si moltiplica, si ottiene .999999. Se vuoi che la risposta sia "sì" dovrai creare un test per invers moltiplicativi che non sia semplice come moltiplicare e testare per l'uguaglianza a 1.

La stessa cosa accade con il binario.

Answer 3

Il problema è che 0.2 non può essere rappresentato esattamente in binario, perché la sua espansione binario ha un numero infinito di cifre:

1/5: 0.0011001100110011001100110011001100110011... 

Questo è simile a come 1/3 non possono essere rappresentati esattamente in decimale. Poiché x viene memorizzato in un float che ha un numero finito di bit, queste cifre avranno tagliato ad un certo punto, per esempio:

x: 0.0011001100110011001100110011001 

Il problema deriva dal fatto che le CPU utilizzano spesso una precisione superiore internamente, in modo da quando 've appena calcolato 1/y, il risultato avrà più cifre, e quando si carica x per confrontarli, x otterrà esteso per abbinare la precisione interna della CPU.

1/y: 0.0011001100110011001100110011001100110011001100110011 
    x: 0.0011001100110011001100110011001000000000000000000000 

Così quando si effettua un confronto diretto bit per bit, sono diversi.

Nel vostro secondo esempio, invece, memorizzando il risultato in una variabile significa che diventa troncato prima di fare il confronto, quindi confrontandoli a questa precisione, sono uguali:

x: 0.0011001100110011001100110011001 
    z: 0.0011001100110011001100110011001 

Molti compilatori hanno interruttori voi può consentire di forzare i valori intermedi a essere troncati ad ogni passo per coerenza, tuttavia il consiglio abituale è di evitare di fare confronti diretti tra valori a virgola mobile e controllare invece se differiscono di meno rispetto ad un valore epsilon, che è ciò che Gangnus is suggesting.

Answer 4

Ciò che colpisce è che qualunque sia la regola dell'arrotondamento, ci si aspetta che il risultato delle due versioni sia lo stesso (due volte errato o due volte corretto)!

Molto probabilmente, nel primo caso una promozione con maggiore precisione nei registri FPU avviene quando si valuta x == 1/y, mentre z = 1/y memorizza realmente il risultato di precisione singola.

Altri contributori hanno spiegato perché 5 == 1/0,2 può fallire, non è necessario ripeterlo.

Answer 5

Le discussioni in altre risposte sono grandiose e quindi non ripeterò nessuno di loro, ma non c'è codice. Ecco un po 'di codice per controllare se in realtà un paio di carri dà esattamente 1.0 quando moltiplicati.

il codice effettua alcune assunzioni/asserzioni (che normalmente incontrati nella piattaforma x86):
- float 's sono 32 bit binario (alias single precision) IEEE-754
- sia int' s o long 's sono a 32-bit (I ha deciso di non fare affidamento sulla disponibilità di uint32_t)
- memcpy() copie galleggianti in INT/Longs tale che 8873283.0f diventa 0x4B076543 (cioè certa "endianness" è previsto)

Un'ipotesi in più è questo :
- riceve i float effettivi che * si moltiplicano (ad es.moltiplicazione dei galleggianti non utilizzare valori di precisione più elevati che l'hardware matematica/libreria può utilizzare internamente)

#include <stdio.h> 
#include <string.h> 
#include <limits.h> 
#include <assert.h> 

#define C_ASSERT(expr) extern char CAssertExtern[(expr)?1:-1] 

#if UINT_MAX >= 0xFFFFFFFF 
typedef unsigned int uint32; 
#else 
typedef unsigned long uint32; 
#endif 
typedef unsigned long long uint64; 

C_ASSERT(CHAR_BIT == 8); 
C_ASSERT(sizeof(uint32) == 4); 
C_ASSERT(sizeof(float) == 4); 

int ProductIsOne(float f1, float f2) 
{ 
    uint32 m1, m2; 
    int e1, e2, s1, s2; 
    int e; 
    uint64 m; 

    // Make sure floats are 32-bit IEE754 and 
    // reinterpreted as integers as we expect 
    { 
    static const float testf = 8873283.0f; 
    uint32 testi; 
    memcpy(&testi, &testf, sizeof(testf)); 
    assert(testi == 0x4B076543); 
    } 

    memcpy(&m1, &f1, sizeof(f1)); 
    s1 = m1 >= 0x80000000; 
    m1 &= 0x7FFFFFFF; 
    e1 = m1 >> 23; 
    m1 &= 0x7FFFFF; 
    if (e1 > 0) m1 |= 0x800000; 

    memcpy(&m2, &f2, sizeof(f2)); 
    s2 = m2 >= 0x80000000; 
    m2 &= 0x7FFFFFFF; 
    e2 = m2 >> 23; 
    m2 &= 0x7FFFFF; 
    if (e2 > 0) m2 |= 0x800000; 

    if (e1 == 0xFF || e2 == 0xFF || s1 != s2) // Inf, NaN, different signs 
    return 0; 

    m = (uint64)m1 * m2; 

    if (!m || (m & (m - 1))) // not a power of 2 
    return 0; 

    e = e1 + !e1 - 0x7F - 23 + e2 + !e2 - 0x7F - 23; 
    while (m > 1) m >>= 1, e++; 

    return e == 0; 
} 

const float testData[][2] = 
{ 
    { .1f, 10.0f }, 
    { 0.5f, 2.0f }, 
    { 0.25f, 2.0f }, 
    { 4.0f, 0.25f }, 
    { 0.33333333f, 3.0f }, 
    { 0.00000762939453125f, 131072.0f }, // 2^-17 * 2^17 
    { 1.26765060022822940E30f, 7.88860905221011805E-31f }, // 2^100 * 2^-100 
    { 5.87747175411143754E-39f, 1.70141183460469232E38f }, // 2^-127 (denormalized) * 2^127 
}; 

int main(void) 
{ 
    int i; 
    for (i = 0; i < sizeof(testData)/sizeof(testData[0]); i++) 
    printf("%g * %g %c= 1\n", 
      testData[i][0], testData[i][1], 
      "!="[ProductIsOne(testData[i][0], testData[i][1])]); 
    return 0; 
} 

uscita (vedi a ideone.com):

0.1 * 10 != 1 
0.5 * 2 == 1 
0.25 * 2 != 1 
4 * 0.25 == 1 
0.333333 * 3 != 1 
7.62939e-06 * 131072 == 1 
1.26765e+30 * 7.88861e-31 == 1 
5.87747e-39 * 1.70141e+38 == 1