Blogity blog blog


Skaičiavimų tikslumas: kvadratinės akys

2008-10-17, rtfb

Š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

 
20 komentarų:

avatar

O čia dar šiek tiek info:

http://www.python.org/doc/2.5.2/lib/module-decimal.html

Čia python’as, bet, manau, kitose kalbose irgi yra.


avatar

Ačiū už nuorodą, Daliau. Man ji priminė, kad neįdėjau nuorodos į kardinaliai priešingą šaltinį :-). Jei kas nori labai, Labai, LABAI pasigilinti į tai, ko čia tik paviršius paglostytas, skaitykit “What Every Computer Scientist Should Know About Floating-Point Arithmetic”:

http://web.cse.msu.edu/~cse320/Documents/FloatingPoint.pdf


Nuotaikingas bemaž 100 puslapių “skaitaliukas” :-)

-rtfb


avatar

perskaičiau bendram išsilavinimui. iki galo. patiko :) na tiesa, asemblerio dalį teko praleist


avatar

Ten specialiai tam reikalui kiekviena eilutė paaiškinta, galima skaityti vien tik komentarus :-)

Ir išvis, tu juk informatikos doktorantė :-P

-rtfb


avatar

kaip chemijos inžinerijos/biotechnologijos doktorantė turiu pasakyt, kad nieko nesupratau. Visiškai. Nors tas bendras tekstas tai įdomiai parašytas. Aš tiesiog matematikos nemoku;)

(paguoda ta, kad visą dieną skaičiausi apie afiniškumo žymes chimeriniams baltymams)


avatar

Eilinį kartą turėčiau tau padėkoti už ypač naudingą informaciją - ačiū!
O dabar lekiu apsitvarkyt savo kode… ;-)


avatar

Šiaip nieko nesupratau, nes deja Mitasiunas man nedeste, o ir neturejo , bet bėda iš tos lyg temos. Netiksliai Javoj apskaičiuoja trupmenas (float formatas). Tai ka man daryt, kad jas skaiciuotu tiksliai?


avatar

batty, apibūdink problemą tiksliau.

-rtfb


avatar

na kaip pvz:

public class Text2 { double h[]=new double[9]; double suma=0;

public Text2() {

    h[1]=0.6;
    h[2]=0.0005;
    h[3]=0.099;
    h[4]=0.015;
    h[5]=0.106;
    h[6]=0.096;
    h[7]=0.083;
    h[8]=0.0005;

    for (){
    suma=h[i]+suma;
}
    System.out.println(""+suma);    
}

public static void main(String args[]){ new Text2(); }
}

Turėtų ats būti 1, bet meta 0.9999999999999998 Negaliu apvalint, nes nezinau kokios trupmenos bus įvestos, kokio ilgio.


avatar

Oi, nemoku naudotis šitais daiktais :S žodžiu, sudedant double (vis gi ne float), atsakymas nera tikslus.
Kodel?? kaip tai pakeist?


avatar

Batty, naudojant FPU slankųjį kablelį, ko gero, tiksliau ir nepadarysi. Mirtinai prireikus didelio tikslumo, tektų imtis trupmenų arba programinio slankiojo kablelio bibliotekų. Bet patarimas nr. 1: pasitenkink tokiu tikslumu ;-)

-rtfb


avatar

Teks matyt kazka galvot kitaip, nes matai, man reikia susumuot tikimybes ir patikrint ar jos lygios 1. O jei apvalinsiu tą 0.9999999999999998 tarkim tūkstantųjų tikslumu, tai jei tikimybių suma bus 0,9999 iš tiesų, programa suapvalinus man mes, kad =1.

Anyway, ačiū už atsakymą ;)


avatar

Dirbant su float’ais yra labai bloga praktika tikrinti ar rezultatas lygus kokiam konkrečiam skaičiui. Paprastai nelygus ;)


avatar

Yep. Batty, lygink su kažkokia pasirinkta paklaida:

if (abs (result - 1.0) < 0.0001) // užskaitom, kad lygu

Kaip „protingai“ pasirinkti tą epsilon, yra ilga istorija, teks skaitytis čia: http://www.physics.ohio-state.edu/~dws/grouplinks/floating_point_math.pdf

Bet praktikoje dažnai pakanka „iš akies“ nustatyti kokią nors tikslinio skaičiaus aplinką ir tuo pasitenkinti.

-rtfb


avatar

Vis gi radau (kolkas atrodo veikiantį) sprendimą - sumuoti didėjančia tvarka skaičius :D


avatar

Hm, čia irgi teisingas pastebėjimas, kad veiksmų eiliškumas turi įtakos tikslumui. Ir visgi primygtinai rekomenduoju toleruoti tam tikrą paklaidą, kaip rodžiau anksčiau.

-rtfb


avatar


if (abs (result - 1.0) < 0.0001) // užskaitom, kad lygu

rtfb,
o kam tas abs čia? Kiek beskaičiuoju, tiek su neigiamais tiek su teigiamais skaliarais sąlyga tenkinama :-)


avatar

Bah, mistype:
sąlyga tenkinama arba ne, tačiau algoritmas puikiai veikia ir be abs (fabs)*


avatar

NightVisio? Nelabai supratau… Siūlai išmesti abs? Bet juk tada

if ((0.1 - 1.0) < 0.00001)

klaidingai įvertinama kaip true!

-rtfb


avatar

Ouch, my mistake. Sorry :-)


Mano nacis spam-filtras visus laiko botais.
Išspręskite šitą captcha jeigu jūs ne toks:

9 + 10 =