Undersökning av IEEE 754

Syftet med denna laboration är att Du skall lära Dig något om flyttalsaritmetik. Övningen kan ge intrycket av att allting alltid går fel. Som tur är ser praktiken ut på ett annat sätt, det mesta fungerar nästan jämt. Man måste dock lära sig att känna igen och undvika "farliga" situationer.

Tips: till några övningar har jag länkat en tips-sida. Denna innehåller lösningen (nästan i alla fall). Så om Du vill ha en större utmaning kan Du vänta med att studera tipsen. Observera att länkarna till tips-sidorna inte syns om Du skriver ut denna sida på papper.
En del övningar är frivilliga.

IEEE 754 är en standard som föreskriver hur flyttalsaritmetik med +, -, *, / och sqrt skall utföras. Den anger det interna formatet för tillgängliga flyttalstyper, avrundningsmoder etc. Standarden beskriver inte hur exp, sin etc. skall fungera. En person som har betytt mycket för att standarden blev av är William Kahan, professor vid UCB (Univ. of California at Berkeley; samma ställe som ansvarar för BSD Unix). Kahan är, för övrigt, hedersdoktor vid Chalmers och har t.ex. också gjort HP-räknedosornas integrationsrutin och lösaren av f(x) = 0 problem. Historien bakom standarden.

Testa hur Matlab uppför sig när man beräknar :  a=1.0/0.0, b=1.0/a, sin(a), sin(b), log(a), log(b) ? Vad händer när man beräknar c=0.0/0.0, d=1.0/c, sin(d) ? Slutsatser? Hur uppför sig Ditt favoritprogramspråk (om Du har något)? Får Du körningsavbrott eller inte? Observera decimalpunkterna. Om Du räknar med heltal (som vi numeriker inte är intresserade av ) får Du helt andra resultat (i Matlab är dock alla tal flyttal så decimalpunkterna är inte nödvändiga). Man kan också vara tvungen att lura kompilatorn så att den inte gnäller direkt (t.ex. genom, zero = 0.0; a = 1.0/zero;).

Byt utskriftsformat till format long e . Observera att detta inte ändrar antalet siffror vi räknar med, utan endast utskriften. Skriv   logspace(-301, -324, 24)'   (notera ' transponatet). Ser det inte lite konstigt ut? Slutsats? Förklaring? Ledning: byt till  format hex  och ge sedan logspace-kommandot.

Vi vill, med hjälp av Matlab, studera när fullständig utskiftning inträffar. Dvs givet ett tal x så vill vi hitta ett tal delta så att x + delta blir x men där x + ett lite större delta blir lite större än x. Ett sätt att göra detta är att bilda (x + delta) - x . Om detta blir exakt 0 (Matlab skriver 0 och inte bara ett litet tal) har delta skiftats ut fullständigt. För att vi skall veta att inget väsentligt större delta fungerar vill vi, för samma delta, att (x + 1.01 * delta) - x skall bli skilt från noll. Vi har då en hygglig approximation av det kritiska deltat. Bestäm de tre delta-värdena när x antar följande tre värden:  2, 2e20 respektive 2e-20. Använd format short e när Du räknar. Slutsats? Notera att 2e20 betyder 2 × 1020, så det är inte exponentialfunktionen, exp,  som är inblandad.

Flyttalsberäkningar följer inte de vanliga matematiska lagarna. Till exempel så behöver inte (a+b)+c vara lika med a+(b+c) (den associativa lagen gäller inte). Detta gör att man kan få olika resultat beroende på i vilken ordning beräkningarna utförs. Ibland är det viktigt att tänka på ordning, ibland inte.

Beräkna, i enkel precision,  s1 = 1.0 + 1/2.0 + 1/3.0 + ...+ 1/1000000.0 samt s2 = 1/1000000.0 + 1/999999.0 + ... + 1/2.0 + 1.0. Försök att förklara varför Du får olika värden på s1 och s2. Vilket värde ligger närmare det exakta? Varför?   Ledning: Om du använder Matlabs single-funktion, räcker det inte att bilda single av slutsumman som du får från beräkningen i dubbel precision, det ger ett alltför bra resultat. Gör så här i stället: s = single(0); bilda sedan summan med en for-loop där du accumulerar summan i s. Om du vill beräkna summan i något annat programspråk, se Tips nedan.

Tips.

Sammanställ Dina svar till ovanstående övningar på fil.

Frivillig

Uttrycken (x-1)^6 och 1-6*x+15*x^2-20*x^3+15*x^4-6*x^5+x^6 är ekvivalenta om man räknar exakt. På en dator kommer de att uppföra sig olika dock. Plotta (help plot) uttrycken ovan (i samma plot) i Matlab då x antar värdena 0.995:1e-4:1.005. Förklara varför kurvorna ser så olika ut. Vad händer om vi plottar (x+1)^6  (och motsvarande expanderade form) istället (för samma x)?

Intressant men frivillig läsning:

Vad avrundningsfel kan ställa till med i verkligheten.

Två länkar till,  Collection of Software Bugs och Some disasters attributable to bad numerical computing

Skrämmande resultat. cos och sin på en HP.

Källkoden till Suns rutiner för elementära funktioner.

En av många sidor om Intels divisionsbug.

The Risk Digest.