Datorfabrikanternas strävan efter snabbhet leder ibland till oönskade resultat. Betrakta följande C-program:

#include <math.h>
#include <stdio.h>
main()
{
        double          r, p, q;

        r = 268435451.5;
        p = r - 2.0;
        q = r - 1.0;

        /* Check before taking the root */
        if ( q * q >= p * r )
                printf("%f\n", sqrt(q * q - p * r));
}

När man exekverar (på en IBM RS/6000) får man:

syk-0606 > a.out
NaNQ

dvs en (Quit) Not-a-Number (man skiljer på quit och signaling NaN) och detta på grund av att argumentet till sqrt är negativt trots att vi först kontrollerar att argumentet är större än eller lika med noll! Anledningen till att detta inträffar är att IBMs Power2-CPU har en så kallad "fused multiply-add" vilket är en instruktion som kan beräkna a*b+c. Det sker ingen avrundning av a*b utan endast avrundning av slutresultatet (detta bryter mot IEEE 754). Så tittar man på en del av assemblerkoden framgår det att (# är mina kommentarer):

	
.main:                            
        mfspr   r0,LR
        stm     r30,-8(SP)
        st      r0,8(SP)
        stu     SP,-112(SP)
        l       r31,T.22.__STATIC(RTOC)
        l       r30,T.16.NO_SYMBOL(RTOC)
        lfd     fp1,0(r30)
        stfd    fp1,56(SP)
        lfs     fp2,8(r30)
        fa      fp1,fp1,fp2      # addition, p = r - 2
        stfd    fp1,64(SP)
        lfd     fp1,56(SP)
        lfs     fp2,12(r30)
        fa      fp2,fp1,fp2      # addition, q = r - 1
        stfd    fp2,72(SP)

        fm      fp3,fp2,fp2      # fp3 = fp2 * fp2 (= q * q)
        lfd     fp4,64(SP)
        lfd     fp1,56(SP)
        fm      fp1,fp1,fp4      # fp1 = fp1 * fp4 (= p * r)
        fcmpu   1,fp3,fp1        # compare, if-satsen
        cror    CR0_EQ,CR1_FEX,CR1_VX
        bc      BO_IF_NOT,CR0_EQ,__L7c

        fms     fp1,fp2,fp2,fp1  # fp1 = fp2 * fp2 - fp1 (= q * q - fp1)
        bl      .sqrt{PR}        # sqrt(fp1)

        cror    CR3_SO,CR3_SO,CR3_SO
        cal     r3,0(r31)
        stfd    fp1,80(SP)
        l       r5,84(SP)
        l       r4,80(SP)
        cror    CR3_SO,CR3_SO,CR3_SO
__L7c:                                  # 0x0000007c (H.10.NO_SYMBOL+0x7c)
        l       r0,120(SP)
        mtspr   LR,r0
        ai      SP,SP,112
        lm      r30,-8(SP)
        bcr     BO_ALWAYS,CR0_LT
        .long   0x00000000

fms är instruktionen för a*b+c och fp1, fp2, fp3, fp4 är flyttalsregister. I if-satsen beräknas q*q och p*r var för sig varefter en subtraktion utförs. När sqrt skall beräknas utnyttjas det redan beräknade p* r som ligger i fp1, varefter fms-instruktionen beräknar q * q - fp1 i ett svep.

Det visar sig att p * r beräknas som 72057591085137952 vilket överskrider det exakta värdet med 2.75. q*q blir likdant varför if-satsen jämför två lika stora tal.  När det kommer till sqrt så beräknas q*q - 72057591085137952 (detta delresultat beräknas exakt) som blir -1.75 vilket givetvis inte sqrt gillar.

Det är ett elände att likartade operationer utförs på olika sätt. En anledning till att standarden infördes var att man skulle slippa sådana konstigheter. Det hjälper inte om man iställer skriver

sqrt((q*q) - (p*r))

eftersom kompilatorn ignorerar dessa extra parenteser. Om man formulerar om if-satsen som

if ( q * q - p * r >= 0.0 )

fungerar det på grund av att fms-instruktionen nu även används i if-satsen.

fms-instruktionen är inte alltid av ondo. Innerproduktion mellan två vektorer

x(1)*y(1) + x(2)*y(2) + ... + x(n)*y(n)

kan beräknas på n+3 maskincykler och med ett litet fel om fms används.

Kör vi samma C-program på en Sun-dator får vi resultatet noll; programmet är inte portabelt med andra ord!

Back