Numeric problem with Ada

Numeric problem with Ada

[CNAM.CNAM.FR!bortz]

(It seems that my first attempt to post this was a failure. I repost.
Sorry if you receive it twice.)

Hello, numeric addicts,

I've a strange problem with numeric expressions in Ada. The test
program (the original
 one was a calendar) is:

procedure TEST is
	Mz : integer := 7;
	Z2 : float;
	package FLOAT_TEXT_IO is new TEXT_IO.FLOAT_IO (float);
begin
	Z2 := 2.6 * float (Mz) - 0.2;
	FLOAT_TEXT_IO.PUT (Z2, Aft => 6, Exp => 0);
end TEST;        

The program displays 17.999998 (the mathematical result is 18). I
understand that, since
 float'DIGITS is 6 on all the machines I used (Vax/VMS and Sun3), I
cannot expect the
 rightmost digit (the 8th) to be correct. If I declare my own
floating-point type with 15
 digits, I can have my result (18) fine. But my problem comes from the
fact that the C
 compilers don't behave this way. The program:

main ()
{
	int mz = 7;
	float z2;
	z2 = 2.6 * mz - 0.2;
	printf ("%2.6f\n", z2);
}

displays 18.000000.
This is true also on both the Vax and the Sun3.
Now the problem is: where does the difference comes from? I've checked
the assembly code
 produced and the two programs use the same types (I was thinking that
C was using the
 double-precision multiplication, but it's not the case). I'm not sure the code
s are
 equivalent (it's too difficult for me to read assembly language), but I think 
they
 should.
Is there something in the Ada floating-point model that prevents the
Ada compiler to give
 better results?



Stephane Bortzmeyer           Conservatoire National des Arts et Metiers
	
bortzmeyer@vcnam.cnam.fr      Laboratoire d'Informatique
bortz@cnam.cnam.fr            292, rue Saint-Martin			
                              75141 Paris Cedex 03
                              France	

"C'est la nuit qu'il est beau de croire a la lumiere." E. Rostand

					
	

Re: Numeric problem with Ada

[Jan Kok]

In article <9111131948.AA20227@cnam.cnam.fr> bortz@CNAM.CNAM.FR writes:
SB>  ...
SB>procedure TEST is
SB>	Mz : integer := 7;
SB>	Z2 : float;
SB>	package FLOAT_TEXT_IO is new TEXT_IO.FLOAT_IO (float);
SB>begin
SB>	Z2 := 2.6 * float (Mz) - 0.2;
SB>	FLOAT_TEXT_IO.PUT (Z2, Aft => 6, Exp => 0);
SB>end TEST;        
SB>
SB>The program displays 17.999998 (the mathematical result is 18).  ...
SB>
SB>main ()
SB>{
SB>	int mz = 7;
SB>	float z2;
SB>	z2 = 2.6 * mz - 0.2;
SB>	printf ("%2.6f\n", z2);
SB>}
SB>
SB>displays 18.000000.
SB>
SB>Stephane Bortzmeyer           Conservatoire National des Arts et Metiers
	

   Dear Stephane,
The answer with the Ada program is correct. You might accept
this, but I suppose you want to understand the phenomenon.
You did not say that the internal representations of the results
with the two programs are different (I suppose they are not).
Could you try something like:
   	Z2 := float (18);
   	FLOAT_TEXT_IO.PUT (Z2, Aft => 6, Exp => 0);
I am suggesting this, because the explanation might be in
the PUT procedure. Conversion from internal (binary)
floating-point value to a decimal character string of finite
length requires a sophisticated algorithm, and we certainly
know about Ada I/O implementations that do not recognise all
exact decimal values. Then, while computing successive decimal
digits for the output string, some accuracy can get lost.

I can also offer you an explanation in the case that the internal
values of the Ada and C results are different, but then I would
assume that the Ada PUT procedure would show a larger difference from 18.
--Jan Kok.
-- 
Mail: Jan Kok, CWI (afd. NW), Postbus 4079, NL-1009 AB Amsterdam, Nederland
Net : jankok@cwi.nl
---------------------------------------------------------------
  Why is a computer?

Re: Numeric problem with Ada

[Jim Giles]

At least on the Sun/Sparc, C is compiled so that it _always_ uses
double precision for intermediate calculations - no matter how
explicitly you ask for just (float).  

   main()
   { 
      float z;
      int mz;
      mz = 7;
      z = 2.6 * mz - 0.2;
      if ( z == 18.0 ) printf ("%2.6f\n", z);
   }

This program always prints `18.000000' even if mz is set in a separately
compiled procedure.  The reason is that the calculation in the assignment
statement to `z' is always done in double and then rounded to float when
the assignment is done.  The same continues to be true even if the
expression is altered to:
      z = (float) 2.6 * (float) mz - (float) 0.2;

So, that's why C always gives 18.000000 on Sun/Sparc.  I suspect the
same may be true of other C implementations.

J. Giles

Re: Numeric problem with Ada

[Jean-Pierre Rosen]

It is likely that the C IO library performs some roundings behind your back...
Printing a value is *never* a guarantee of the actual value, especiallly with
not-rigourously-defined languages...

Re: Numeric problem with Ada

[psinntp!vitro.com!v7.vitro.com!vaxs09]

> I've a strange problem with numeric expressions in Ada.
...
> procedure TEST is
>       Mz : integer := 7;
>       Z2 : float;
>       package FLOAT_TEXT_IO is new TEXT_IO.FLOAT_IO (float);
> begin
>       Z2 := 2.6 * float (Mz) - 0.2;
>       FLOAT_TEXT_IO.PUT (Z2, Aft => 6, Exp => 0);
> end TEST;        
>
(Program displays 17.999998)
...
> main ()
> {
>       int mz = 7;
>       float z2;
>       z2 = 2.6 * mz - 0.2;
>       printf ("%2.6f\n", z2);
> }
(Program displays 18.000000)

I can't speak for your results on the Sun.  On the VAX, however, I was able
to reproduce your results and found an underlying reason.

Using VAX/VMS and DEC's Ada compiler and DEC's C compiler, the following
machine code was generated:

        ADA                                  C

        cvtlf   -80(fp),r0                   cvtld      #7,r0
        mulf2   #1717977382,r0               muld2      $CODE,r0
        subf2   #-858964148,r0               subd2      $CODE+8,r0
        movf    r0,-84(fp)                   cvtdf      r0,ap

Result: FFFF428F (single)                    FFFF428F FFFFFFFF (double)
        = 18.0 minus 2**-19                  00004290 (single)
        ~ 18.0 minus .0000019                = 18.0 exactly
        ~ 17.999998

As you can see, the C program performed the calculations in double
precision and converted the result back to single.  The Ada program
performed the calculations in single precision.

Examining the value actually stored reveals that the result in the
two cases differs by only the least significant bit (the byte order
of VAX F-floating format makes the longword integer representations
of the result appear much different).  The C program computed the result
of the calculation as 18.0 exactly.  Although the double precision
intermediate result was inexact, rounding to single precision
produced a result that happens to exactly match the mathematically
correct answer.

Can a C guru confirm or deny that common practice in C compilers is to
perform floating point calculations in double precision?

VAX F-floating format has

	(Numbering bits in longword with low order = bit 0 and high order 31)
        Sign bit as the bit #15
        Binary exponent in bits 7..14, encoded excess 128
        High order bit of mantissa not encoded.  Because of normalization,
         it is always set.
        Next seven bits of mantissa in 0..7
        Low order sixteen bits of mantissa in 16..31

00004290 = sign bit clear (positive)
           exponent = 5 (multiply mantissa by 2**5)
           binary mantissa = .100100000000000000000000

FFFF428F = sign bit clear (positive)
           exponent = 5 (multiply mantissa by 2**5)
           binary mantissa = .100011111111111111111111

The D-floating format is identical with additional 32 bits of mantissa.

	John Briggs		vaxs09@v7.vitro.com

Re: Numeric problem with Ada

[Thomas Buenermann]

K&R C (pre-ANSI, pre-Std. C) specified that all floating point
expressions will be propagated to double. And parameter passing
is an expression. Thus, floating point values will always be put
on the stack as double, and printf expects them this way.

(If the formal parameter is a float, the function must read the 
double from the stack and store it back as a float before entering 
the body of the function.)

- Thomas

#include <std/disclaimer.h>

Thomas Buenermann	thomas@hpugrca.grc.hp.com
HP-UX Kernel & 		Thomas BUENERMANN / HP 8340/00
Language Support	Hewlett-Packard GmbH, Support-Zentrum
CRC Germany		Berliner Str. 111, D-4030 Ratingen, FRG

Re: Numeric problem with Ada

[CNAM.CNAM.FR!bortz]

Thanks to all for the replies to my question about why C and Ada give
two different results in "equivalent" numeric programs.
Two kinds of solution were presented: those who said it was a problem
in the I/O library and those who were looking toward the use of
double-precision arithmetic by the C compiler. From the following
modification of my test program (to suppress problems related with
I/O), I think the second hypothesis is the good one.
Here is the new program:
Ada:
	if Z2 /= 18.0 then 
              TEXT_IO.PUT_LINE ("Z2 /= 18.0"); 
        else 
              TEXT_IO.PUT_LINE ("Z2 = 18.0"); 
        end if;
C:
	if (z2 != 18.0) printf ("z2 != 18.0\n"); else printf ("z2 = 18.0\n"); 

Ada finds that Z2 /= 18.0 and C that z2 = 18.0. 
So the difference is really in internal computations. Two replies by
E-mail were specially interesting:

pragma BEGIN_QUOTATION (hilfingr@tully.cs.berkeley.edu "Paul N. Hilfinger");
C implementations generally DO use double precision (the new ISO stanard
calls for it in this case; both of the constants 2.6 and 0.2 are
double-precision, lacking any suffix to the contrary).  When I tried
this on a Sun 3, the compiler generated double-precision instructions
and got the result you indicated (18.000000).  When I changed the
calculation of z2 to read

       z2 = (float) 2.6 * mz - (float) 0.2;

the result was printed as 17.999999 (I used gcc version 1.39 on a
Sun3).  Perhaps you misread the assembly language code (those mnemonics
ARE rather cryptic). 
pragma END_QUOTATION (hilfingr@tully.cs.berkeley.edu "Paul N. Hilfinger");

I agree with the possibility that I missed something in the assembly
code (I prefer to read Ada!).

pragma BEGIN_QUOTATION (Per-Gunnar.Johansson@imsa.hv.se "Per-Gunnar Johansson")
;
The computer store Your object in a binary representation, and if You have
six digits in the decimal system You have about 20 digits in the binary
system.
If You wamt to store a float-object the coputer transform it on the form
    number = m * (2 ** e), where 0.5 <= m < 1 or m = 0.
It is no problem with the number 7 (7 = 0.875 * (2 ** 3), because 
0.875 in the decimal system is exact 0.111 in the binary system.
BUT, the number 2.6 = 0.65 * (2 ** 2), and 0.65 in the decimal system
is 0.1010 0110 0110 0110 0110 ..... (and so on), in the binary system.
The representation in Your computer, with six decimal digits, is about 20
binary digits. It means that Your computer represent the number 2.6 with
m = 0.1010 0110 0110 0110 in the binary systemand it is about 0.649999619
in the decimal system, so (2.6 =) 0.649999619 * (2 ** 2) = 2.599998476.
In the same way (0.2 =) 0.199999809
In that way, You calculating operation
    2.6 * 7 - 0.2 =
  = 20599998476 * 7 - 0.199999809 =
  = 17.99998952
It is not 17.999998 but I think it maybe depend on the accurancy on my
pocket calculator.
pragma END_QUOTATION (Per-Gunnar.Johansson@imsa.hv.se "Per-Gunnar Johansson");


Stephane Bortzmeyer           Conservatoire National des Arts et Metiers
	
bortzmeyer@vcnam.cnam.fr      Laboratoire d'Informatique
bortz@cnam.cnam.fr            292, rue Saint-Martin			
                              75141 Paris Cedex 03
                              France	

"C'est la nuit qu'il est beau de croire a la lumiere." E. Rostand