rtfblog

Šis rašinys yra labai techniškas ir mokslo populiarinimo prasme moksliškas. Programuotojams rekomenduoju jį perskaityti dėl abiejų priežasčių, visiems kitiems rekomenduoju neišsigąsti pirmosios ir perskaityti bent jau dėl antrosios.

Kažkada „seniai seniai“, viename forume įvyko įdomi keisto reiškinio analizė. Reiškinys: neįtikėtinai didelė skaičiavimo paklaida, atliekant trivialius skaičiavimus su slankiojo kablelio skaičiais. Šiaip jau tiek reiškinys, tiek kapstymo procesas gana nuobodūs, o tiems, kam tenka darbo reikalais su juo susidurti, tai dar ir labai skaudi tema.

Ir visgi visus tuos metus ta analizė mane kažkuo žavėjo. Kelis kartus perskaičiau forumo giją iš naujo, ir tik neseniai supratau kas gi ten tokio svarbaus. Ogi svarbu tai, kad ta analizė yra puiki iliustracija iškart keliems dalykams:

1. Prigimtiniam programų sistemų kompleksiškumui
2. Scientific method‘ui
3. Klaidų paieškos programose metodikai

Paskutiniam punktui tai, ko gero, ne toks jau chrestomatinis pavyzdys, kaip pirmiems dviems, bet visgi.

Taigi, problematika. Vienas jaunuolis ėmė ir paklausė kodėl 1020 kartų prie nulio pridedant po 0.1 gaunama ne 102, kaip kad norisi tikėtis, o 101.99901580810547, jei dirbama su float ir 102.00000151991844, jei dirbama su double. Paklaida tikrai neįtikėtinai didelė.

Kitas kolega jam paantrino, papasakodamas baisesnių atvejų:

    ------------------------------------------------------
    double stat = count_wins () / count_loses ();

    if (stat != count_wins () / count_loses ())
        // sąlyga _niekada_ neturėtų būti tenkinama, bet...
    ------------------------------------------------------

beigi:

    ------------------------------------------------------
    if (a + b + c != a + c + b)
        // sąlyga _niekada_ neturėtų būti tenkinama, bet...
    ------------------------------------------------------

Daugiau nevarginsiu, štai kodas, kuris išgavo aukščiau minėtas paklaidas (kodas minimaliai paredaguotas, išlaikant esmę, bet patrumpinant):

    ------------------------------------------------------
    1: float  f = 0.0f;
    2: double d = 0.0f;
    3:
    4: for (int i = 1020; --i; )
    5: {
    6:     f += 0.1f;
    7:     d += 0.1f;
    8: }
    9: printf ("\n  %.20f \t %.20f\n", f, d);
    ------------------------------------------------------

Patyręs C programuotojas iškart pamatys, kad čia triskart pasikartoja viena ir ta pati klaida: neteisingai naudojamasi double tipu.

Antroje ir septintoje eilutėse nepridera naudoti literalo 0.1f. Mat skaičius su kableliu C kalboje automatiškai traktuojamas kaip double tipo. Tuo tarpu postfiksas ‘f’ reikalingas tam, kad jis būtų traktuojamas kaip float. Koks iš to skirtumas išsiaiškinsim mažumėlę vėliau, o kol kas trečia tos pačios klaidos apraiška.

Devintoje eilutėje taip pat neteisingai nurodomas antro kintamojo tipas formate: %.20f reiškia, kad parametras yra float tipo, ir spausdinti jį reikia iki 20-ies skaitmenų po kablelio tikslumu. Bet juk kintamasis yra double tipo! Taip ir reikia nurodyti: %.20lf.

Naudojant float tipą, klaidų neprivelta, bet kadangi float yra ženkliai mažesnio tikslumo, paklaida didesnė.

Dviem raidelėm per daug, ir viena per mažai. Ir dar tokios „neesminės“… Kokią tai gali turėti įtaką tikslumui? Štai, originalių skaičiavimų paklaida nuo norimo rezultato:

    102.00000151991844000000  /  102 =
    = 1.0000000149011611764705882352941       (7 nuliai po kablelio)

O štai rezultatas ir jo paklaida ištaisius klaidas:

    101.99999999999848000000  /  102 =
    = 0.99999999999998509803921568627451      (13 devynetų po kablelio)

Pataisykit mane, jei klystu, bet mano skaičiavimais, milijoną kartų tiksliau.

Labai įdomu pasiaiškinti koks gi tiksliai gaunasi skirtumas skaičiavimuose, naudojant šiek-tiek-ne-visiškai-bet-visgi-nevisai suderinamus tipus. Būtinai tai padarysim, kad pamatytume kiek daug visko sumuojasi į bet kurios programos galutinį sudėtingumą. Bet pradžioj atlikim mokslinį eksperimentą.

Kaip ir daugelį mokslinių eksperimentų, šį sufleruoja ne specialistui netikėtas, bet patirties diktuojamas smalsumas: kas, jeigu šią trumpą, paprastą programą, sukompiliuosime su optimizacijomis?

Kam? Tuoj paaiškinsiu, pradžiai pažiūrėkim, kas gavosi.

Originalių skaičiavimų su float paklaida (prisiminkite, jokios klaidos ten nebuvo, paklaidą sąlygojo iš prigimties mažas float tipo tikslumas):

    101.99901580810547000000  /  102 =
    = 0.99999035105985754901960784313725      (5 devynetai po kablelio)

Optimizuotų skaičiavimų su float paklaida:

    102.00000182539225000000  /  102 =
    = 1.0000000178960024509803921568627       (7 nuliai po kablelio)

Ir pagaliau optimizuotų skaičiavimų su double paklaida:

    102.00000000000000000000  /  102 =
    = skaičiuokit patys ;-)

Kas įvyko? Kodėl vien tik optimizavimo įjungimas 100x sumažino paklaidą, susikaupiančią float kintamajame? Kodėl double kintamojo tikslumas šoktelėjo iki palubių, kokių, regis, jau ir nebesitikėjome?

Žinantis žmogus iškart paaiškintų, bet aš apsimesiu, kad nežinau ir parodysiu kaip sužinojau. Kad suprasti kokį skirtumą sukėlė optimizacijos, reikia žiūrėti, ką konkrečiai kompiliatorius iškrėtė vienu ir kitu atveju iš identiško išeities kodo.

(Keletas pastraipų N-14, nagrinėsim asemblerį :-))

Neoptimizuotas kodas (smarkiai apkastruotas, kad liktų tik esmė):

    ------------------------------------------------------
    LC1:
            .long   1036831949        // (float) 0.1f
    LC2:
            .long   -1717986918       // (double) 0.1 (apatiniai 4 baitai iš 8)
            .long   1069128089        // (double) 0.1 (viršutiniai 4 baitai iš 8)

            // ...
            // Ciklas:
            movl    $1020, -4(%ebp)   // "for (i = 1020;"
    L6:
            leal    -4(%ebp), %eax    // eax = 1020
            decl    (%eax)            // sumažinam i vienetu
            cmpl    $-1, -4(%ebp)     // ar i >= 0?
            jne     L9                // jei taip, vykdom ciklą
            jmp     L7                // jei ne, ciklą baigiam ir spausdinam
    L9:
            flds    -8(%ebp)          // parsiunčiam iš atminties kintamąjį f
            flds    LC1               // parsiunčiam iš atminties 0.1f
            faddp   %st, %st(1)       // sudedam juos
            fstps   -8(%ebp)          // dedam rezultatą atgal į atmintį
            fldl    -16(%ebp)         // parsiunčiam iš atminties kintamąjį d
            fldl    LC2               // parsiunčiam iš atminties 0.1
            faddp   %st, %st(1)       // sudedam juos
            fstpl   -16(%ebp)         // dedam rezultatą atgal į atmintį
            jmp     L6                // grįžtam prie ciklo sąlygos
    L7:
            // čia, jau po ciklo, spausdinami rezultatai
    ------------------------------------------------------

O štai optimizuotas:

    ------------------------------------------------------
    LC3:
            .long   1036831949        // (float) 0.1f
    LC4:
            .long   -1717986918       // (double) 0.1 (apatiniai 4 baitai iš 8)
            .long   1069128089        // (double) 0.1 (viršutiniai 4 baitai iš 8)

            // ...
            // Ciklas:
            movl    $1019, %ebx       // for (i = 1020

            // ties čia, st(0) registre yra kintamasis d,
            // o st(1) registre -- kintamasis f. Ciklas dar neprasidėjo

    L11:
            fxch    %st(1)            // sukeičiam reikšmes tarp st(0) ir st(1)
            decl    %ebx              // sumažinam i vienetu
            fadds   LC3               // pridedam 1.0f prie st(0) (t.y. f)
            fxch    %st(1)            // sukeičiam reikšmes tarp st(0) ir st(1)
            faddl   LC4               // pridedam 1.0  prie st(0) (t.y. d)
            cmpl    $-1, %ebx         // ar i >= 0?
            jne     L11               // jei taip, tęsiam ciklą

            // čia, jau po ciklo, spausdinami rezultatai
    ------------------------------------------------------

Įgudę pastebės greitai, neįgudę tik labai pasigilinę, bet visgi matosi vienas esminis skirtumas: neoptimizuota versija cikle nuolat siunčia duomenis iš kintamojo atmintyje į FPU (slankiojo kablelio skaičiavimų koprocesorių), sumuoja ir siunčia atgal į atmintį. Tuo tarpu optimizuota versija vieną kartą prieš ciklą sukiša duomenis į FPU registrų steką, ir visą ciklą „neišleidžia“ skaičių iš FPU.

Taigi, padarėm eksperimentą: pažiūrėjom kas vyksta, jei suoptimizuojam. Rezultatas: tikslumas ženkliai didėja. Paanalizavom eksperimento duomenis, matom skirtumus. Iškelkime hipotezę:

O gal FPU viduje visi skaičiavimai atliekami didesniu tikslumu, ir tik tam, kad tilptų į reikiamo dydžio kintamuosius atmintyje, jų tikslumas mažinamas?

Tikrinam hipotezę, ir štai ką turim:

„However, x87 does not perform operations according to strict IEEE-754 formats, since it uses wide registers internally“.

(IEEE-754 yra standartas, apibrėžiantis slankiojo kablelio skaičius technikoje. Ir jeigu būti visai tiksliems, tai modernūs x87-ieji FPU pagal nutylėjimą viduje atlieka skaičiavimus 80 bitų tikslumu (long double), bet tą galima pakeisti per FPU kontrolinius registrus).

Hipotezė pasitvirtino. Atkreipkime dėmesį, kad hipotezė buvo itin sėkminga: ji ne tik paaiškino eksperimento rezultatus, bet ir paaiškino anksčiau empiriškai pataisytas klaidas originaliam kode; paaiškino kodėl naudojant „beveik suderinamus“ kintamojo ir literalo tipus, paklaida būdavo didesnė, bei kaip gali nutikti a + b + c != a + c + b ir panašios „neįmanomybės“.

O dabar pabandykime šią hipotezę išplėsti būsimiems rezultatams: kas, jeigu optimizuotą kodą dar šiek tiek pajudinsime rankomis ir padarysime, kad cikle prie float tipo kintamojo (kuris, kaip jau žinome, visai ne float, kol FPU viduje) sumuosime ne mažo tikslumo 0.1f, o didelio tikslumo 0.1? Štai šiek tiek pakeistas kodas, kur skiriasi vos viena eilutė -- naudojamas tikslesnis literalas, bei jo dydį atitinkanti instrukcija:

    ------------------------------------------------------
            // kaip ir anksčiau, st(0) registre jau yra kintamasis d,
            // o st(1) registre -- kintamasis f

    L11:
            fxch    %st(1)
            decl    %ebx

            // skirtumas tik šioje eilutėje: faddl vietoj fadds,
            // LC4 (kuris rodo į 0.1) vietoj LC3 (kuris rodo į 0.1f)
            faddl   LC4

            fxch    %st(1)
            faddl   LC4
            cmpl    $-1, %ebx
            jne     L11
    ------------------------------------------------------

Dabar jau nenuostabu, kad rezultatai štai tokie:

    float:                          double:
  102.00000000000000000000       102.00000000000000000000

Matome, kad visus skaičiavimus atlikus didesniu tikslumu, atsakymas suskaičiuojamas be paklaidos. Tai, be abejo, yra tam tikra apgaulė, nes procesoriaus viduje skaičiavimai buvo atliekami didesniu tikslumu, nei nurodyta C kalboje, bet siekiant rezultato, nulinė paklaida yra lemiamas kriterijus.

Na ir pabaigai, iš kur išvis, net teoriškai gali rastis netikslumai, atliekant tokius trivialius skaičiavimus su tokiais paprastais skaičiais. Čia juk ne kokia nors šaknis iš trijų?

Ogi netikslumai atsiranda iš to, kad 0.1 ir sqrt(3) dvejetainiam pavidale yra vienodai neparankūs skaičiai:

    ------------------------------------------------------
               IEEE-754 single precision floating-point:
    0.1f:         0-01111011-10011001100110011001101
    sqrtf (3):    0-01111111-10111011011001111010111
    ------------------------------------------------------

Šis, paskutinis, pavyzdys, viso rašinio kontekste, turėtų paaiškinti tiems, kas dar nesuprato, kam universitete reikėjo klausyti viso to „beprasmio ir nenaudingo“ Mitašiūno kurso ;-)

Programavimas iš prigimties sudėtingas užsiėmimas, ir štai jums buvo pavyzdys. Tokias klaidas galima išgauti programuojant bet kuria kalba. Ar tai būtų Java, ar C#, ar Lisp, ar bet kuri kita kalba. Ir kaip matėte, tam, kad atkapstyti visus niuansus, reikėjo ir asemblerio žinių, ir šiek tiek procesorių architektūros, ir transliatorių bent minimaliai pažinoti.

Ačiū už dėmesius. Ypač tiems, kas ištvėrė šitą technomaninę kankynę iki galo :-)

Originali gija forume: http://gamedev.lt/viewtopic.php?f=18&t=297