Java och C

Man måste tänka sig lite för hur man skriver numeriska program i Java och C. Här följer ett enkelt summaexempel. Antag att vi skall beräkna summan (med flyttal)

-1/4 - 1/3 -1/2 - 1 + 1 + 1/2 + 1/3 + 1/4

i Java. Summan borde bli noll, men på grund av avrundningsfel kommer den sannolikt inte att bli det. Här följer några olika loopar som beräknar summan: 

    int    n = 4;
    float  s;
    double d;

    s = 0.0f;
    for(int k = -n; k <= n; k++)
      if (k != 0) s += 1.0 / k;

Jag börjar med att sätta s till noll, 0.0f är en enkelprecisionsnolla. Att i stället använda 0.0 (dubbel precision) går inte bra eftersom javac då ger felmeddelandet: "possible loss of precision" eftersom s är en float-variabel. n har värdet fyra och vi hoppar över k = 0 i loopen. Värdet på s blir 7.4505806E-8. Det är viktigt att notera att 1.0 är i dubbel precision, k konverteras då automatiskt till double, divisionen och additionen utförs i double, men partialsumman avrundas till float. Det kan man lätt kontrollera genom att studera bytekoden , eller lite mindre utförligt (koden produceras med javap -c Program). Den intressanta delen ser ut så här: 

  14 fload_2    hämta s
  15 f2d        s -> double  (ty s är ju float)
  16 dconst_1   1.0 (täljaren)
  17 iload_3    hämta k
  18 i2d        k -> double
  19 ddiv       double division
  20 dadd       double add
  21 d2f        s -> float
  22 fstore_2   lagra s

I följande exempel tvingar jag täljaren att vara en float-etta (1.0f eller (float) 1.0). 

    s = 0.0f;
    for(int k = -n; k <= n; k++)
      if (k != 0) s += 1.0f / k;

Så här blir bytekoden: 

  14 fload_2    hämta s
  15 fconst_1   1.0f
  16 iload_3    hämta k
  17 i2f        k -> float
  18 fdiv       floating divide
  19 fadd       floating add
  20 fstore_2   lagra s

I detta fall görs alla beräkningar således i enkel precision och resultatet blir ett annat, -1.4901161E-7.

Det skall också uppföra sig på detta vis (enligt "The Java language specification" avsnitt 4.2.4).

Följande loop använder en double-variabel d för summan och resultatet , ungefär 2.8E-16, blir ju mycket bättre eftersom hela beräkningen utförs i dubbel precision. 

    d = 0.0d;
    for(int k = -n; k <= n; k++)
      if (k != 0) d += 1.0 / (double) k;  // 1.0 / k går också bra

Observera att följande inte ger något bättre resultat: 

    s = 0.0f;
    for(int k = -n; k <= n; k++)
      if (k != 0) s += 1.0 / k;
    d = s;  // konvertera till double

Eftersom vi mellanavrundar summan till enkel precision kan vi ju inte förvänta oss att konverteringen från float till double skall ge något. Vi får samma resultat som i första exemplet på denna sida. Ett enkelprecisionsresultat mitt i en dubbelprecisionsberäkning kan förstöra resultatet.

C uppför sig på samma sätt förutsatt att man inte tvingar kompilatorn att följa K & R-versionen av C (definitionen av C i boken av Brian Kernighan och Dennis Ritchie, "The C Programming Language", Prentice-Hall 1978). Så här står det i manualbladet för cc (Suns C-kompilator)

-X[a|c|s|t]

Specifies the degree of conformance to the ANSI/ISO C standard. Specifies one of the following:

a (ANSI)
ANSI C plus K&R C compatibility extensions, with semantic changes required by ANSI C. Where K&R C and ANSI C specify different semantics for the same construct, the compiler will issue warnings about the conflict and use the ANSI C interpretation. This is the default compiler mode.

c (conformance)
Strictly conformant ANSI/ISO C, without K&R C compatibility extensions. The compiler will issue errors and warnings for programs that use non-ANSI/ISO C constructs.

s (K&R C)
The compiled language includes all features compatible with (pre-ANSI) K&R C. The compiler tries to warn about all language constructs that have differing behavior between Sun ANSI/ISO C and the K&R C. Invokes cpp for processing. __STDC__ is not defined in this mode. (See the C Transition Guide for differences between ANSI/ISO and K&R C.)

Vi kompilerar de två looparna: 

  int    n = 4;
  int    k;
  float  s;

  s = 0.0f;
  for(k = -n; k <= n; k++)
    if(k) s += 1.0 / k;
  printf("s = %13.6e\n", s);

  s = 0.0f;
  for(k = -n; k <= n; k++)
    if(k) s += 1.0f / k;    /* OBS 1.0f */
  printf("s = %13.6e\n", s);

% cc -Xc prog.c     cc prog.c ger samma resultat
% a.out             exekvera
s =  7.450581e-08   utskrifter
s = -1.490116e-07

Vi kompilerar sedan så att vi kan studera assemblerkoden (jag har bara skrivit ut flyttalsoperationerna). %f2 etc. är flyttalsregister i FPUn. 

% cc -S -Xc prog.c  generera assemblerkod på filen fil.s

Första loopen

        fitod   %f2,%f2       int -> double
        fdivd   %f4,%f2,%f2   divide double
        faddd   %f6,%f2,%f2   add double
        fdtos   %f2,%f2       double -> single

Andra loopen

        fitos   %f2,%f2       int -> double
        fdivs   %f3,%f2,%f2   divide single
        fadds   %f4,%f2,%f2   add single

Nu kompilerar vi enligt K & R. 

% cc -Xs prog.     kompilera enligt K & R
% a.out            exekvera
s =  7.450581e-08  samma
s =  7.450581e-08

Om man ser på assemblerkoden visar det sig att båda looparna ser ut som första loopen ovan. I K & R C gjordes alla beräkningar i double precision, men detta har ändrats i ANSI-standarden. Ett program som fungerade för några år sedan (eller med en gammal kompilator) behöver alltså inte fungera nu.

Back