15. Aritmetik med flyttal: problem och begränsningar

Flyttal representeras i datorhårdvara som bas 2 (binära) fraktioner. Till exempel har det decimala bråket 0,625 värdet 6/10 + 2/100 + 5/1000, och på samma sätt har det binära bråket 0,101 värdet 1/2 + 0/4 + 1/8. Dessa två bråk har identiska värden, den enda verkliga skillnaden är att det första är skrivet i bråknotation i bas 10 och det andra i bas 2.

Tyvärr kan de flesta decimala fraktioner inte representeras exakt som binära fraktioner. En följd av detta är att de decimala flyttal som du anger i allmänhet endast approximeras av de binära flyttal som faktiskt lagras i maskinen.

Problemet är lättare att förstå till en början i bas 10. Tänk på bråket 1/3. Du kan approximera det som ett bas 10-bråk:

0.3

eller, bättre,

0.33

eller, bättre,

0.333

och så vidare. Oavsett hur många siffror du är villig att skriva ner kommer resultatet aldrig att bli exakt 1/3, men det kommer att bli en allt bättre approximation av 1/3.

På samma sätt, oavsett hur många bas 2-siffror du är villig att använda, kan decimalvärdet 0,1 inte representeras exakt som ett bas 2-brott. I bas 2 är 1/10 det oändligt upprepade bråket

0.0001100110011001100110011001100110011001100110011...

Om du stannar vid ett ändligt antal bitar får du en approximation. På de flesta maskiner idag approximeras floats med hjälp av ett binärt bråk där täljaren använder de första 53 bitarna med början med den mest signifikanta biten och där nämnaren är en tvåpotens. I fallet 1/10 är det binära bråket 3602879701896397 / 2 ** 55, vilket är nära men inte exakt lika med det sanna värdet 1/10.

Många användare är inte medvetna om approximationen på grund av hur värdena visas. Python skriver bara ut en decimal approximation av det sanna decimalvärdet av den binära approximationen som lagras av maskinen. Om Python skulle skriva ut det sanna decimalvärdet för den binära approximation som lagrats för 0,1, skulle det på de flesta maskiner behöva visa:

>>> 0.1
0.1000000000000000055511151231257827021181583404541015625

Det är fler siffror än de flesta tycker är användbart, så Python håller antalet siffror hanterbart genom att visa ett avrundat värde istället:

>>> 1 / 10
0.1

Kom ihåg att även om det utskrivna resultatet ser ut som det exakta värdet av 1/10, är det faktiska lagrade värdet det närmaste representerbara binära bråket.

Intressant nog finns det många olika decimaltal som delar samma närmaste approximativa binära fraktion. Till exempel approximeras siffrorna 0.1 och 0.10000000000000001 och 0.1000000000000000055511151231257827021181583404541015625 alla av 3602879701896397 / 2 ** 55. Eftersom alla dessa decimalvärden har samma approximation kan vilket som helst av dem visas samtidigt som invariansen eval(repr(x)) == x bevaras.

Historiskt sett har Pythons prompt och den inbyggda funktionen repr() valt den med 17 signifikanta siffror, 0.10000000000000001. Från och med Python 3.1 kan Python (på de flesta system) nu välja det kortaste av dessa och helt enkelt visa 0.1.

Observera att detta ligger i själva naturen för binärt flyttal: detta är inte en bugg i Python, och det är inte heller en bugg i din kod. Du kommer att se samma sak i alla språk som stöder din hårdvaras aritmetik med flyttal (även om vissa språk kanske inte visar skillnaden som standard eller i alla utdatalägen).

Om du vill ha en trevligare utmatning kan du använda strängformatering för att få fram ett begränsat antal signifikanta siffror:

>>> format(math.pi, '.12g')  # give 12 significant digits
'3.14159265359'

>>> format(math.pi, '.2f')   # give 2 digits after the point
'3.14'

>>> repr(math.pi)
'3.141592653589793'

Det är viktigt att inse att detta i själva verket är en illusion: du avrundar helt enkelt displayen av det verkliga maskinvärdet.

En illusion kan ge upphov till en annan. Eftersom 0,1 till exempel inte är exakt 1/10, kan summan av tre värden av 0,1 inte heller ge exakt 0,3:

>>> 0.1 + 0.1 + 0.1 == 0.3
False

Eftersom 0,1 inte kan komma närmare det exakta värdet av 1/10 och 0,3 inte kan komma närmare det exakta värdet av 3/10, kan förrundning med funktionen round() inte heller hjälpa:

>>> round(0.1, 1) + round(0.1, 1) + round(0.1, 1) == round(0.3, 1)
False

Även om talen inte kan göras närmare sina avsedda exakta värden, kan funktionen math.isclose() vara användbar för att jämföra inexakta värden:

>>> math.isclose(0.1 + 0.1 + 0.1, 0.3)
True

Alternativt kan funktionen round() användas för att jämföra grova approximationer:

>>> round(math.pi, ndigits=2) == round(22 / 7, ndigits=2)
True

Binär aritmetik med flyttal innehåller många överraskningar som denna. Problemet med ”0,1” förklaras i detalj nedan, i avsnittet ”Representation Error”. Se Exempel på flyttalsproblem för en trevlig sammanfattning av hur binär flyttalsaritmetik fungerar och vilka typer av problem som är vanliga i praktiken. Se även The Perils of Floating Point för en mer fullständig redogörelse för andra vanliga överraskningar.

Som det står i slutet: ”Det finns inga enkla svar.” Var ändå inte onödigt försiktig med flyttal! Felen i Pythons flyttalsoperationer ärvs från flyttalshårdvaran, och på de flesta maskiner är de i storleksordningen högst 1 del i 2**53 per operation. Det är mer än tillräckligt för de flesta uppgifter, men du måste komma ihåg att det inte är decimalaritmetik och att varje flyttalsoperation kan drabbas av ett nytt avrundningsfel.

Även om patologiska fall finns, för de flesta vardaglig användning av flyttalsaritmetik du kommer att se det resultat du förväntar dig i slutändan om du helt enkelt avrunda visningen av dina slutresultat till antalet decimaler du förväntar dig. str() oftast räcker, och för finare kontroll se str.format() metodens format specificerare i Format String Syntax.

För användningsfall som kräver exakt decimal representation, försök använda modulen decimal som implementerar decimal aritmetik som är lämplig för redovisningstillämpningar och högprecisionstillämpningar.

En annan form av exakt aritmetik stöds av modulen fractions som implementerar aritmetik baserad på rationella tal (så att tal som 1/3 kan representeras exakt).

Om du är en flitig användare av flyttalsoperationer bör du ta en titt på NumPy-paketet och många andra paket för matematiska och statistiska operationer som tillhandahålls av SciPy-projektet. Se <https://scipy.org>.

Python tillhandahåller verktyg som kan hjälpa till vid de sällsynta tillfällen då du verkligen vill veta det exakta värdet av en float. Metoden float.as_integer_ratio() uttrycker värdet av en float som ett bråk:

>>> x = 3.14159
>>> x.as_integer_ratio()
(3537115888337719, 1125899906842624)

Eftersom förhållandet är exakt kan det användas för att förlustfritt återskapa originalvärdet:

>>> x == 3537115888337719 / 1125899906842624
True

Metoden float.hex() uttrycker en float i hexadecimal (bas 16), vilket återigen ger det exakta värde som lagras av din dator:

>>> x.hex()
'0x1.921f9f01b866ep+1'

Denna exakta hexadecimala representation kan användas för att rekonstruera floatvärdet exakt:

>>> x == float.fromhex('0x1.921f9f01b866ep+1')
True

Eftersom representationen är exakt är den användbar för att på ett tillförlitligt sätt överföra värden mellan olika versioner av Python (plattformsoberoende) och utbyta data med andra språk som stöder samma format (t.ex. Java och C99).

Ett annat användbart verktyg är funktionen sum() som hjälper till att minska precisionsförlusten vid summering. Den använder utökad precision för mellanliggande avrundningssteg när värden läggs till i en löpande totalsumma. Det kan göra skillnad för den totala noggrannheten så att felen inte ackumuleras till den punkt där de påverkar den slutliga totalsumman:

>>> 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 == 1.0
False
>>> sum([0.1] * 10) == 1.0
True

math.fsum() går längre och spårar alla ”förlorade siffror” när värden läggs till en löpande summa så att resultatet bara har en enda avrundning. Detta är långsammare än sum() men kommer att vara mer exakt i ovanliga fall där stora ingångsvärden mestadels upphäver varandra och ger en slutlig summa nära noll:

>>> arr = [-0.10430216751806065, -266310978.67179024, 143401161448607.16,
...        -143401161400469.7, 266262841.31058735, -0.003244936839808227]
>>> float(sum(map(Fraction, arr)))   # Exact summation with single rounding
8.042173697819788e-13
>>> math.fsum(arr)                   # Single rounding
8.042173697819788e-13
>>> sum(arr)                         # Multiple roundings in extended precision
8.042178034628478e-13
>>> total = 0.0
>>> for x in arr:
...     total += x                   # Multiple roundings in standard precision
...
>>> total                            # Straight addition has no correct digits!
-0.0051575902860057365

15.1. Representation Fel

I detta avsnitt förklaras exemplet ”0,1” i detalj och det visas hur du själv kan utföra en exakt analys av fall som detta. Grundläggande kännedom om binär flyttalsrepresentation förutsätts.

Representation error hänvisar till det faktum att vissa (de flesta, faktiskt) decimalfraktioner inte kan representeras exakt som binära (bas 2) fraktioner. Detta är den främsta anledningen till att Python (eller Perl, C, C++, Java, Fortran och många andra) ofta inte visar det exakta decimaltal du förväntar dig.

Varför är det så? 1/10 är inte exakt representerbar som en binär fraktion. Sedan åtminstone år 2000 använder nästan alla maskiner IEEE 754 binär aritmetik med flyttal och nästan alla plattformar mappar Python-flottal till IEEE 754 binär64-värden med ”dubbel precision”. IEEE 754 binary64-värden innehåller 53 bitars precision, så vid inmatning strävar datorn efter att konvertera 0,1 till det närmaste bråket den kan av formen J/2**N där J är ett heltal som innehåller exakt 53 bitar. Omskrivning

1 / 10 ~= J / (2**N)

som

J ~= 2**N / 10

och med tanke på att J har exakt 53 bitar (är >= 2**52 men < 2**53), är det bästa värdet för N 56:

>>> 2**52 <=  2**56 // 10  < 2**53
True

Det vill säga 56 är det enda värdet för N som gör att J har exakt 53 bitar. Det bästa möjliga värdet för J är då denna kvot avrundad:

>>> q, r = divmod(2**56, 10)
>>> r
6

Eftersom återstoden är mer än hälften av 10 erhålls den bästa approximationen genom att avrunda uppåt:

>>> q+1
7205759403792794

Därför är den bästa möjliga approximationen av 1/10 i IEEE 754 dubbel precision:

7205759403792794 / 2 ** 56

Genom att dividera både täljaren och nämnaren med två minskar bråket till:

3602879701896397 / 2 ** 55

Observera att eftersom vi avrundade uppåt är detta faktiskt lite större än 1/10; om vi inte hade avrundat uppåt skulle kvoten ha varit lite mindre än 1/10. Men det kan inte i något fall vara exakt 1/10!

Så datorn ”ser” aldrig 1/10: vad den ser är den exakta fraktionen som anges ovan, den bästa IEEE 754-dubbelapproximationen som den kan få:

>>> 0.1 * 2 ** 55
3602879701896397.0

Om vi multiplicerar bråket med 10**55 kan vi se värdet med 55 decimalers noggrannhet:

>>> 3602879701896397 * 10 ** 55 // 2 ** 55
1000000000000000055511151231257827021181583404541015625

vilket innebär att det exakta talet som lagras i datorn är lika med decimalvärdet 0,1000000000000000055511151231257827021181583404541015625. Istället för att visa hela decimalvärdet avrundar många språk (inklusive äldre versioner av Python) resultatet till 17 signifikanta siffror:

>>> format(0.1, '.17f')
'0.10000000000000001'

Modulerna fractions och decimal gör dessa beräkningar enkla:

>>> from decimal import Decimal
>>> from fractions import Fraction

>>> Fraction.from_float(0.1)
Fraction(3602879701896397, 36028797018963968)

>>> (0.1).as_integer_ratio()
(3602879701896397, 36028797018963968)

>>> Decimal.from_float(0.1)
Decimal('0.1000000000000000055511151231257827021181583404541015625')

>>> format(Decimal.from_float(0.1), '.17')
'0.10000000000000001'