Oppolzer - Informatik / Stanford Pascal Compiler


Home       Lebenslauf       Schwerpunkte       Kenntnisse       Seminare       Kunden       Projekte       Produkte       Blog       Stanford Pascal       Kontakt

The Stanford Pascal Compiler / Evolution Steps

Back to Compiler main page

Differences on floating point computations and rounding

When comparing the results of floating point operations on Hercules/VM and Windows, I discovered many differences. Most details of this difficult analysis and work are documented here:

Stanford Compiler Facebook page.

The biggest problem was that the output values were implicitly rounded on the PC (using the sprintf() function from ANSI C), but not always correctly. On Hercules/VM, there was no rounding at all.

For example: when writing the value 12.35 with only one decimal digit, Stanford Pascal on VM would always write 12.3, but Stanford Pascal on the PC would write 12.3 or 12.4, depending on the precision of the computed value (if there has been some arithmetic before the output operation, the value may be a little bit lower than "true" 12.35, so 12.3 is printed instead of 12.4, which would be correct).

To make the floating point output comparable, there is only one solution (which I know from my long experience with this topic, when working for a big European insurance company):

the output has to be rounded to the needed decimal precision,
but before the rounding a correction has to be made
that compensates the error introduced by the computations before.

This solution, BTW, will not work in every pathological case - a specialist on this topic once told me: "what you do is pretty good, but you can do what you want: I will always find a counter example where it will not work".

The "correction" in fact means that - depending on the platform - a value is added to the original value, that depends on the size on the original value; for example: if x is the original value, then the correction value is x * (16 ** (- 13)).

To do this, I implemented a new function in the Pascal runtime that does this kind of "corrective" rounding. This function is called ROUNDX, which means "round extended". It has a second parameter which tells the decimal position where the rounding has to take place; -2 for example means: at the 2nd decimal position. Values from 4 to -15 are accepted for the second parameter.

The ROUNDX function needs the FLOOR function as defined by the C runtime; so I added this function, too.

On Hercules/VM, ROUNDX is implemented in PASLIBX.PAS in Pascal (calling sequence generated by the compiler using CALLLIBRARYFUNC). FLOOR is a true CSP function which is implemented in PASMONN using ASSEMBLER (only a very small number of ASSEMBLER instructions). See the ROUNDX1 source code (from PASLIBX.PAS) below.

On Windows etc., roundx is a C function, and floor is already available, because it is a ANSI C function.

Because the "classical" ROUND function from Pascal worked differently on the both platforms, too, I changed its implementation, too. Now on both platforms, ROUND (x) is implemented as ROUNDX (x, 0), followed by TRUNC ... because ROUNDX gives a REAL result, but ROUND needs INTEGER. With this change, now ROUND works the same on all platforms (in most cases).

Note: I decided to do the corrective action only for ROUND and ROUNDX, but not for TRUNC and FLOOR.

The ROUNDX1 Source Code (from PASLIBX)


local function ROUNDX1 ( WERT : REAL ; BEREICH : INTEGER ) : REAL ; /**********************************************************/ /* */ /* roundx.c */ /* */ /* Rundungsfunktion neu mit geaenderter Logik; */ /* die Korrekturkonstante wird anhand der Groessen- */ /* ordnung des Ausgangswertes bestimmt (Ausgangs- */ /* wert durch (16 hoch 13); damit wird bei beiden */ /* Plattformen mindestens eine 1 an der letzten */ /* Ziffernposition dazuaddiert). */ /* */ /* Autor: Bernd Oppolzer */ /* April 1995 */ /* */ /**********************************************************/ var FAKTOR : REAL ; TEST : REAL ; RUNDKONST : REAL ; const FAKTTAB : array [ 0 .. 19 ] of REAL = ( 10000.0 , 1000.0 , 100.0 , 10.0 , 1.0 , 10.0 , 100.0 , 1000.0 , 10000.0 , 100000.0 , 1000000.0 , 10000000.0 , 100000000.0 , 1000000000.0 , 10000000000.0 , 100000000000.0 , 1000000000000.0 , 10000000000000.0 , 100000000000000.0 , 1000000000000000.0 ) ; begin (* ROUNDX1 *) FAKTOR := FAKTTAB [ 4 - BEREICH ] ; if WERT < 0.0 then TEST := - WERT else TEST := WERT ; if TEST < 1.0E-55 then begin ROUNDX1 := 0.0 ; return end (* then *) ; /************************************************/ /* */ /* 4 * (16 hoch 12) = 1125899906842624.0 */ /* 8 * (16 hoch 12) = 2251799813685248.0 */ /* 12 * (16 hoch 12) = 3377699720527872.0 */ /* 16 hoch 13 = 4503599627370496.0 */ /* */ /************************************************/ RUNDKONST := TEST / 1125899906842624.0 ; if BEREICH < 0 then begin TEST := ( TEST + RUNDKONST ) * FAKTOR + 0.5 ; TEST := FLOOR ( TEST ) ; if WERT < 0.0 then ROUNDX1 := - TEST / FAKTOR else ROUNDX1 := TEST / FAKTOR end (* then *) else if BEREICH > 0 then begin TEST := ( TEST + RUNDKONST ) / FAKTOR + 0.5 ; TEST := FLOOR ( TEST ) ; if WERT < 0.0 then ROUNDX1 := - TEST * FAKTOR else ROUNDX1 := TEST * FAKTOR ; end (* then *) else begin TEST := ( TEST + RUNDKONST ) + 0.5 ; TEST := FLOOR ( TEST ) ; if WERT < 0.0 then ROUNDX1 := - TEST else ROUNDX1 := TEST ; end (* else *) end (* ROUNDX1 *) ;

Example program: the program didn't show 13 and -13 as result of ROUND before my changes:


begin (* HAUPTPROGRAMM *) /******************/ /* positiver wert */ /******************/ R := 12.5 ; DUMP ( ADDR ( R ) , PTRADD ( ADDR ( R ) , 7 ) ) ; R := R * 4724731.123456 ; R := R / 4724731.123456 ; DUMP ( ADDR ( R ) , PTRADD ( ADDR ( R ) , 7 ) ) ; WRITELN ( 'r vor round = ' , R : 15 : 7 ) ; I := ROUND ( R ) ; WRITELN ( 'i nach round = ' , I : 13 ) ; I := TRUNC ( R ) ; WRITELN ( 'i nach trunc = ' , I : 13 ) ; R := FLOOR ( R ) ; WRITELN ( 'r nach floor = ' , R : 15 : 7 ) ; R := 12.357 ; WRITELN ( 'r vor roundx = ' , R : 15 : 7 ) ; R := ROUNDX ( ( R ) , - 2 ) ; WRITELN ( 'r nach roundx = ' , R : 15 : 7 ) ; /******************/ /* negativer wert */ /******************/ R := - 12.5 ; DUMP ( ADDR ( R ) , PTRADD ( ADDR ( R ) , 7 ) ) ; R := R * 4724731.123456 ; R := R / 4724731.123456 ; DUMP ( ADDR ( R ) , PTRADD ( ADDR ( R ) , 7 ) ) ; WRITELN ( 'r vor round = ' , R : 15 : 7 ) ; I := ROUND ( R ) ; WRITELN ( 'i nach round = ' , I : 13 ) ; I := TRUNC ( R ) ; WRITELN ( 'i nach trunc = ' , I : 13 ) ; R := FLOOR ( R ) ; WRITELN ( 'r nach floor = ' , R : 15 : 7 ) ; R := - 12.357 ; WRITELN ( 'r vor roundx = ' , R : 15 : 7 ) ; R := ROUNDX ( ( R ) , - 2 ) ; WRITELN ( 'r nach roundx = ' , R : 15 : 7 ) ; end (* HAUPTPROGRAMM *) .

Now the program produces the correct (and the same) result on both platforms:


EXEC PASRUN TESTRND Dump Speicherbereich von 000302A0 bis 000302A7 000302A0: 41c80000 00000000 ........ ........ *.H......aaaaaaaa* Dump Speicherbereich von 000302A0 bis 000302A7 000302A0: 41c7ffff fffffffd ........ ........ *.G......aaaaaaaa* r vor round = 12.5000000 i nach round = 13 i nach trunc = 12 r nach floor = 12.0000000 r vor roundx = 12.3570000 r nach roundx = 12.3600000 Dump Speicherbereich von 000302A0 bis 000302A7 000302A0: c1c80000 00000000 ........ ........ *AH..........aaaa* Dump Speicherbereich von 000302A0 bis 000302A7 000302A0: c1c7ffff fffffffd ........ ........ *AG..........aaaa* r vor round = -12.5000000 i nach round = -13 i nach trunc = -12 r nach floor = -13.0000000 r vor roundx = -12.3570000 r nach roundx = -12.3600000 Ready; T=0.05/0.19 23:44:15

Back to Compiler main page