Számítás garantált pontossággal I.

A számítógépeken használt eljárásoknak vannak olyan problémái, amik extrém esetben az eredményeket teljesen használhatatlanná tehetik. Ezt a technika jelenlegi fejlettségi fokán már többé nem kell eltűrnünk. Kifinomult eljárások a számítási eredményeket egy matematikai bizonyítás bizonyosságával tudják előállítani. A dupla cikk első része a számábrázolási lehetőségekkel foglalkozik.

A szerencsétlenség a derült égből csapott le. A mérnökök hetek óta a lipcsei Kraftwerk Nord turbinái beállításával foglalkoztak. A gépek már gond nélkül működtek., amikor 1994. november elsején napközben egy hangos durranással a hajtómű széttört. „Egy alkatrész kirepült a csarnokból autónyi lyukat vágva a csarnok tetejébe, és 250 méterre landolt a vízművek területén”, ahogy a következő napi újsághír jelentette. A kiömlő olaj elöntötte a géptermet és tűz ütött ki; négy ember az életét vesztette, hatan súlyosan megsebesültek.
Bármi volt is a katasztrófa közvetlen kiváltó oka  egy alapos vizsgálati jelentés anyaghibát nevezett meg  egy rezonanciahatás biztos szerepet játszott. Még az a tekintélyes mozgási energia sem elég ekkora lendülettel való kivágódáshoz, ami egy ilyen forgó alkatrészben rejtőzik. Ha azonban a tengely fordulatszáma a mechanikai rendszer sajátfrekvenciája közelében van, aminél az hajlamos kilengeni, akkor ezek a kitérések az idő múltával pusztító amplitúdókig emelkedhetnek.
Emiatt a turbinák beindításakor a fordulatszámot lassan kell növelni, úgy hogy az lassan és egyenletesen melegedjen fel, mivel nagyobb hőmérsékleti ugrások az anyagot károsítják; másrészt bizonyos veszélyes fordulatszámokat  amik a rezonanciafrekvenciának megfelelnek  lehetőség szerint gyorsan meg kell haladni.
Hogyan lehet ezeket a sajátfrekvenciákat megtudni? Inkább nem próbálgatással: olyan baleseteket kockáztatnánk, mint amilyen a lipcsei eset is volt. A kicsinyített modelleken való kísérletezés is hamar elvi korlátokba ütközik. Mindenek előtt úgy szeretnénk a berendezést tervezni, hogy az bizonyos fordulatszám-tartományokban egyáltalán nem rendelkezik sajátfrekvenciával. A kész berendezés utólagos javítása nagyon költséges lenne. Tehát a részalkatrészeket kilengési tulajdonságaikkal együtt számítógéppel határozzák meg. Emiatt aztán azoknak a hardver és szoftver szakembereknek, akik a számításokat végzik, illetve akik ehhez segédeszközöket állítanak elő, komoly felelőssége van.
Hogyan bizonyosodhatunk meg arról, hogy a számítások a megkövetelt pontosságon belül vannak és korrektek? A nyugtalanító válasz: a jelenlegi módszerekkel egyáltalán nem. Általánosságban nincs olyan matematikai tétel, ami az eredmények korrektségét biztosítaná.
Hogy a fizikai viszonyokat milyen mértékben lehet matematikai (differenciál)egyenletekbe átírni, az ismeretes. Vannak olyan eljárások  főleg a végeselem módszerek, melyekkel ezek az egyenletek számítógéppel megoldható (algebrai) egyenletekké alakíthatók, éspedig úgy, hogy az ezek közelítése útján adódó hiba elfogadható korláton belül marad. A gépre csak egy nagy egyenletrendszer megoldása marad  időtrabló, de kevésbé látványos munka.
Azonban éppen ebben a lépésben adódhatnak durva hibák. Vannak olyan számítógéppel kapott eredmények, amelyeknek az első számjegye sem pontos  a korrekt programozás ellenére. A szoftver a dicsért verifikációs teszteket minden további nélkül teljesítette volna.
Min múlik ez akkor? Egy konvención, amely olyan általánosan elfogadott és  legalábbis az első pillantásra  olyan világos, hogy nehéz elképzelni, hogyan is mehetne másként: a lebegőpontos számokon. De van alternatíva. Az alapgondolat meghökkentően egyszerű: de a hatás csak alkalmas hardver és teljesen új koncepciójú szoftver egyesítésével mutatkozik meg. Ezekben a cikkekben ezért nem csak nagypontosságú számítási eljárásokról lesz szó, hanem egy ezeket megvalósító lapkáról és intervallum aritmetikáról is. Mindezekkel az összetevőkkel sikerülhet hosszú, komplikált számításokat korrekten végrehajtani  ugyanolyan biztonsággal, mint amit egy klasszikus matematikai bizonyítás ad.
Ez a cikk egy matematikus, Ulrich Kulisch és iskolája munkáit írja le. Amióta Kulisch 1966-ban professzor lett a Karlsruhei Egyetemen, mindkét alapötletet a kezdetektől a perfekt megoldásig fejlesztette  és teszi ma is. Az eredmények sokaságából itt csak egy kis részt lehet megmutatni.

A technika jelen állása: a lebegőpontos számábrázolás

Egy számot a számítógép például a következő alakban ábrázol: 1,726793381539 10-1, ami egy 1 és 10 közötti szám  az úgy nevezett mantissza  rögzített számú decimális számjeggyel a tizedesvessző mögött beszorozva 10 egy hatványával, hogy a tizedesvessző a megfelelő helyre kerüljön. Szokásos, hogy 16 decimális jegyet használnak. (Az, hogy a számítógép belül bináris számokat használ decimálisok helyett, a tárgyalásunk szempontjából jelentéktelen.) Ha egy számítás eredménye kiesik az 1 – 10 tartományból  felül vagy alul, akkor a tizedesvesszőt, úgy toljuk el, hogy a mantissza ismét 1 és 10 közé kerüljön, és az exponenst megfelelően korrigáljuk (normalizálás): innen származik a lebegőpontos (floating point) név.
A valós számok, amikkel az elmélet dolgozik, végtelen sok jeggyel rendelkeznek a tizedesvessző mögött. Egy számítógépen ábrázolt szám tehát mindig egy kerekített közelítése az érintett valós számnak. A lebegőpontos számábrázolás előnye, hogy az (általában) 16 értékes jegyet mindig hozni tudja, függetlenül attól, hogy milyen nagyságrendű a szám. Ha például a tizedes vessző előtt 6, utána 10 jegyet engednénk meg, akkor egy milliónál nagyobb számokkal nem is lehetne számolni, és 10-5 nagyságrendű számokra csak 6 számjegyet lehetne megadni, mert a tizedesvessző mögött mindig négy nullát kellene előbb írni. Bonyolult számításokban azonban szabályszerűen és elkerülhetetlenül előfordulnak nagyon nagy és nagyon kicsi számok is.
 Minden számítási lépésben történik kerekítés: két, a tizedes vessző mögött 16 jegyű szám szorzata 32 jegyet tartalmaz a tizedesvessző után. Ezekből az utolsó 16-ot elhagyjuk, vagy már ki se számítjuk, hogy így az eredmény beleférjen a rögzített számformátumba. Ezek amúgy is jelentőség nélküliek ha a két összeszorzandó, mint általában, már amúgy is kerekített számok. Ha két 1 és 10 közti számot összeadunk, és az eredmény 10 vagy nagyobb, akkor a normalizálás során a tizedesvessző egy jeggyel balra tolódik, és a most 17 tizedesjegyű számból az utolsó jegy kerekítéssel eltűnik.
Általában tehát minden számítási lépés kerekítési hibát is jelent. Ilyen kerekítési hibák aztán több számítási lépésen keresztül továbbadódnak, és az idő teltével gyarapodnak. De ez önmagában még nem lehet nagyon rossz. Tegyük fel, hogy van egy millió darab összeadandó lebegőpontos számunk, amelyek mind 1 és 10 között vannak, és sohasem felfelé-, mindig csak lefelé kerekítünk, úgy hogy az összes kerekítési hiba összeadódik. Ebben a kedvezőtlen esetben az utolsó decimális jegyben elkövetett hiba egy milliószor okozna tévedést. Az eredmény tehát nem 16, hanem csak 10 jegyre lenne pontos. Ezzel még mindig elégedett lehet egy felhasználó, főleg akkor, ha az adatai egyébként is csak 5 jegyre pontosak. Akkor van nagyobb baj, ha az összeadandó számok nagyon különböző nagyságrendűek (lásd a kiemelt példát az információvesztésről).
Már csak amiatt, hogy a két számot azonos exponens alakúra kell hozni, hogy a tizedesvessző a tizedesvesszőhöz igazodjon, a kisebb szám tizedes jegyeit jobbra kell tolni egy megfelelő számú hellyel. Balról nullákkal töltjük fel az üres helyeket, és a jobboldalon túlnyúló jegyeket kerekítjük. Ha az a + b eredményből a későbbi számítás során a ismét levonódik, akkor voltaképp b-t kellene kapnunk. De b-nek a kerekített jegyei visszavonhatatlanul elvesztek: az említett példánkban b eredeti 16 jegyéből csak 4 helyes maradt, ha pedig a és b több mint 16 tizedes távolságra lennének egymástól, akkor egy sem. Az a + b - a számítási művelet tehát lényegesen eltérő eredményt ad a számítógépen. Ezzel az összeadás egy lényeges tulajdonsága, a kommutativitás sérül meg.

1. Ábra. A pontosságbeli eltérés okozta számítási eredmény különbsége egy iterált függvény példáján: a sík minden pontjára az f(x,y) = (y-sign(x)|bx - c|1/2,a-x) függvényt értékeljük ki, majd az eredményre ismét az f függvényt, és így tovább (lásd Barry Martin: „Hüpfer”, Spektrum der Wissenschaft 11-1986, 6. o.). Az ábrák több millió ilyen függvényértéket mutatnak, egyszeres pontossággal (balra) és dupla pontossággal (jobbra), indulóérték szerint kiszínezve; a = 12, b = 3, c = 36. Bár eredetileg kicsi kerekítési hiba volt, a milliószám végrehajtott műveletek után komoly hatások mutatkoznak.

A probléma: értékes jegyek elvesztése

Másként kifejezve: amikor két közel azonos nagyságú számot kell egymásból kivonni, az eredmény annyi jegyet veszít a pontosságából, amennyi a vezető nullák számából adódik. Ez a pontos jegyek számának csökkenése a kivonás során, amelytől tartunk. Amikor a hibás számértéket a további számítások során felhasználjuk, például egy nagyobb számmal megszorozzuk, a levágási hiba tetszőlegesen nagy eltérést eredményezhet, és a teljes számítást használhatatlanná teheti.
A klasszikus numerikus matematika a levágási hibát elvi problémának tekinti, amit teljesen megoldani nem lehet, csak elkerülni, éspedig szabály szerint nem a számítás során, hanem korábban. Fáradságosan sok munkát kell befektetni, hogy a számítások rendjét úgy határozzuk meg, hogy a levágás egyáltalán ne léphessen fel. Elég buta dolog a példánkban először az a-hoz a b-t hozzáadni, majd ebből az a-t megint levonni. Ha ezt a két műveletet egyáltalán nem hajtanánk végre, csak az eredményként b-t adnánk meg, nem csak a pontosságot mentenénk meg, hanem a számítási idővel is takarékoskodnánk.
 Általában tehát olyan algoritmusokat kell keresni, amelyek a közel azonos nagyságú számok kivonása nélkül tudnak működni. Kell hogy legyenek ilyenek, ha az alapul szolgáló feladat „jól definiált”. A gyakorlat általában messze elmarad e nemes céloktól.

Lebegőpontos számok összeadása

Információvesztés és levágás

Hogy az a = 1,109632969473378 100 és a b = 6,748659387134009 10-12 számokat korrekt módon összeadjuk, a számítógép lebegőpontos aritmetikája a decimális helyeknek megfelelően egymás alá kell hogy írja a számokat:

	1,109632969473378
	+0,000000000006748659387134009
 
és aztán ahogy szoktuk, össze kell adni a két számot. Az eredménynek az eddigi 16-tal szemben 28 jegye van; mivel azonban a lebegőpontos számábrázolás csak 16 jegyet enged meg, az utolsó 12-t eldobjuk. 
Ezen kerekítési hiba miatt kerül ellentmondásba hosszú összeadási sor esetén eredményünk a kommutativitási szabállyal. Ha a következő elméletileg ekvivalens kifejezéseket kiszámítjuk egy hagyományos számítógépen:

1020 + 17   -10   + 130  - 1020   
1020 + 17   -1020  -  10    + 130
1020 - 1020 +17    -  10   + 130
1020 + 17  +130  - 1020  - 10

akkor a következő eredményeket kapjuk sorrendben: 0, 120, 137 és –10. Itt 1020+17 a számítógép számára 1020, mivel  az első példában  a 17 túl kicsi, hogy a 1020 mellett még szerepet játsszon. A 1020 levonása után aztán a hátsó jegyekből lesznek az elsők.
Ez a levágás semmi esetre se kell hogy olyan nyilvánvaló legyen, mint a példánkban. A következő 

a x + b y = 1         vagy konkrét példán:   64919121   x – 159018721 y = 1   
c x + d y = 0                                             41869520,5 x – 102558961 y = 0

egyszerű egyenletrendszer megoldása csak négy lebegőpontos műveletet igényel. Hagyományos aritmetikával az

x = 102558961, y = 41869520,5

megoldás adódik. De ebből egy számjegy sem helyes, mert a valódi megoldás:

x = 205117992, y = 83739041.

Ugyanis (a d – b c) -vel kell osztani. Ez a szám ebben az esetben nagyon kicsi, a levágás miatt nagyon pontatlan, és ezt a pontatlanságot a vele való osztás rendkívül felnagyítja.

•	Gyakran nincs elég betekintésünk a feladat struktúrájába, hogy a megfelelő algoritmust megtaláljuk.
•	Az azonos nagyságú számok összeadása és kivonása általában nem egymás után következik, mint a példánkban, hanem ezeket számos más művelet választja el, úgy hogy az egymásnak megfelelő műveletpárt nem lehet azonosítani.
•	Az, hogy a levágás bekövetkezik-e, az olykor a beolvasott adatoktól függ, amiket a programozó még nem ismer, amikor a szoftvert írja.
•	Elméletileg lehetséges minden egyes aritmetikai művelet hibáját megbecsülni, és annak továbbterjedését követni a számítás hátralévő részére. Több milliárd művelet esetén  ami modern számítógépeken másodpercek alatt végbemegy  egy ilyen hibaanalízis azonban gyakorlatilag kivitelezhetetlen. 
•	Bizonyos helyzetekben, mindenek előtt egyenletek iteratív megoldása során, közel azonos mennyiségek kivonása csaknem rendszeresen előfordul, éspedig pont akkor, amikor az eljárás már nagyon közel jutott a megoldáshoz (lásd keretes írásunkat).

Egy szokásos segítség az úgy nevezett négyszeres pontosságú aritmetika. Elővigyázatosságból nem egy, hanem két lebegőpontos számnyi helyet használunk egy szám tárolására. (Egy 16 jegyű decimális szám hagyományosan dupla pontosságú, vagy dupla hosszú, mert egyszeres pontosságú számábrázolás is van fele akkora tárigénnyel.) Az így ábrázolható 32 számjegyből a baloldali 16 jegy az első tárhelyre kerül, a jobboldaliak a másodikba. A helyiértékeknek megfelelő összeadást két négyszeres pontosságú számra egy megfelelően programozott utasítássor oldja meg. Ha 16 értékes jegy elveszne a levágás miatt, a négyszeres pontossággal még mindig marad 16 használható hely (a szokásos eljárásnál ilyen nincs).
De miért kell a jegyek számánál a duplázásnál megállnunk? Miért nem használunk nyolcszoros, tizenhatszoros, százszoros pontosságot, vagy reálisabban, mindig a maximális, még ésszerű pontosságot? Ennek az ötletnek a következetes továbbfejlesztése a lebegőpontos számoktól való eltéréshez vezet. Végül is van még fixpontos (fixed point) adatábrázolás, amit normális esetben egészszámokra használunk. A fixpontos aritmetikában nem veszíthetünk el értékes jegyet a normalizálás során, mert nincs normalizálás: a fixpontos összeadás mindig hibamentes. Csak a számítási tartomány nagyon korlátozott: a 232-nél nagyobb számok (ez kicsit több mint két milliárd) nem ábrázolhatók. 
Ezek a korlátok viszont nem elvi természetűek, hanem egy konvención alapulnak: évtizedek óta a fixpontos számokat nem szokták 4 byte-nál (32 bitnél) hosszabb helyen tárolni. Ennek a korábbi alapja, miszerint takarékoskodni kell a memóriával, a hardvertechnika fejlődésével idejét múlta. Egy modern lapkán minden további nélkül megvalósítható egy olyan regiszter, chipen belüli tárolóhely, amely egy fixpontos számot 2000 bináris hellyel, azaz körülbelül 600 decimális hellyel tud ábrázolni.
És elméletileg ez már az új koncepció: számoljunk extrém hosszú fixpontos számokkal, olyan hosszúakkal, hogy azok minden lebegőpontos számot meg tudjanak jeleníteni. Egy 600 jegyű decimális szám, amelynek 300 jegye van a tizedesvesszőtől balra és jobbra is, már elegendő: mivel a ma szokásos lebegőpontos számokban az exponens felülről és alulról korlátozott, és ezeket a tág határokat figyelembe venni nem nehéz. Olyan számokat, amelyek abszolút értéke nagyobb, mint 21024  10300, vagy kisebb, mint 2-1024  10-300, a lebegőpontos számok se tudnak ábrázolni, tehát a fixpontos regiszter sem kell, hogy helyet tartson fel nekik.
Amikor egy ilyen hosszú fixpontos regiszterhez egy szokásos lebegőpontos számot adunk, a különben szükséges normalizálás nem kell. Az exponens csak azt határozza meg, hogy a decimális jegyek a fixpontos szám mely jegyeihez igazodnak, és kezdődhet az összeadás (a sémát lásd alul).
 Érdekes módon ez az eljárás egy olyan technikát újít fel, amely a régi asztali kalkulátorokban szokásos volt (lásd a fényképet jobbra). Az eredményregiszter nagyon hosszú volt, és a beviteli regiszterekkel ellentétben fixpontos volt. Ebbe a regiszterbe tetszőlegesen sok különböző nagyságrendű összeadandót lehetett bevinni anélkül, hogy kerekítési vagy levágási hibák léptek volna fel. Az első elektronikus kalkulátorok aztán ezt az elvet költségcsökkentés végett a lebegőpontos számítás javára feladták. Ma ezek az okok már nem játszanak szerepet: már csak a kialakult technika tehetetlensége és az egyoldalúan képzett szakemberek restsége áll a korábbi elv alkalmazása útjában. 
A dupla cikk második része a skalárszorzással mint az ötödik alapművelettel, az intervallum-aritmetikával és hatékony iterációval foglalkozik.

Irodalomjegyzék

Ulrich Kulisch: Computer, Arithmetik und Numerik  ein Memorandum. Überblicke Mathematik 1998, Vieweg, Wiesbaden, 1998.

Ulrich Kulisch: Advanced Arithmetic for the Digital Computer  Design of Arithmetic Units. Electronic Notes of Theoretical Computer Science, 1999. http://www.elsevier.nl/locate/entcs/volume24.html

Christoph Pöppe: Rechnen mit garantierter Genauigkeit. Spektrum der Wissenschaft, 2000 szeptember, 54-61.

A Spektrum der Wissenschaften nyomán

Egyenletek iteratív megoldása

A Newton-eljárás

Hogyan oldunk meg egy egyenletet? Az iskolában azt tanultuk, hogy addig alakítjuk mindenféle algebrai műveletekkel, amíg az egyenlőségjel egyik oldalán kizárólag az ismeretlen változó lesz, a másik oldalán pedig ismert mennyiségek. A gyakorlatban az egyenletek ilyen eljárással ritkán kezelhetők: már olyan egyszerű egyenletek esetén, mint az x - 2 sin x = 0 egy formula alakú megoldás (az x = 0 kivételével) nem adható meg.
 A numerikus matematika működőképes megoldási módot ajánl ilyen feladatok megoldására tetszőleges pontossággal: a rendszeres próbálgatást. A megoldási eljárás prototípusa a Newton-eljárás. Az egyenlőséget az f(x) = 0 alakban kell felírni, ahol a függvény (a példában x -2 sin(x)) a bal oldalra kerül. A gyakorlati feladatok legtöbbjében f differenciálható, azaz a hozzá tartozó görbe nem tartalmaz sem lyukakat, sem töréseket, tehát minden ponthoz egyértelmű érintő húzható.
A függvény zérushelyét  ami az egyenlet megoldását jelenti  közvetlenül nem tudjuk meghatározni, nem így viszont egy érintő zérushelyét. Tehát választunk egy indulópontot, x0-at, kiszámítjuk a hozzá tartozó görbebeli ponthoz húzott érintő zérusmetszetét, és ezt a pontot x1-nek nevezzük. Mivel a görbe érintője legalábbis az érintőpont körül jól közelíti a görbét, az érintő zérushelye szintén közel kerül a görbe zérushelyéhez. Kedvező körülmények között tehát x1 közelebb van a keresett zérushelyhez, mint x0. Ezután megismételjük az eljárást, tehát x1-ben érintőt húzunk a görbéhez, meghatározzuk annak x2 zérushelyét, és így tovább. Képletekben: a k-adik közelítőértékből, xk-ból a következőt, xk+1-et a következő képlettel határozzuk meg: xk+1 = xk – f(xk)/f’(xk), ahol f’ az f deriváltja. Bizonyos (gyakran teljesülő) feltételek mellett az így definiált xk sorozat konvergál f egy zérushelyéhez. 

A Newton-eljárás kiterjesztései többváltozós egyenletrendszereket oldanak meg; számos változata olyan egyenletrendszerek megoldását is lehetővé teszi, amin a szokásos Newton-eljárás nem működik. Amikor az eljárás már előrehaladt, xk tehát már a zérushely közelében van, f(xk) már szükségképpen nagyon közel van a nullához. Általában ilyen érték két közepes mennyiség különbségeként jön létre, emiatt pontosságvesztés léphet fel, éspedig annál inkább, minél jobb az xk közelítés. Mivel az elkövetett hiba továbbterjed, az eljárás a rossz közelítésen nem tud javítani. 

2. Ábra. A Newton-eljárás az f(x) = x –2 sin x függvényre, az x0 = 4 kezdeti értékkel. Már három lépés után az iteráció olyan közel kerül a zérushelyhez, hogy a különbséget az ábrán nem lehet látni.