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!