[Caml-list] float rounding

[Sean McLaughlin]

Hello,

I'm using Ocaml for an interval arithmetic application. I"m
curious about
what the Ocaml parser/compiler does to float constants. May I assume
that for any constant I enter, eg. 3.1415... (for N digits of pi), that
the compiler will give me a closest machine representable number?
i.e., if I bound a constant by the previous and next floating point
value to
that given me by the compiler,
will it always be the case that my original (mathematical) constant
lies in that interval?

Thanks,

Sean

_______________________________________________
Caml-list mailing list. Subscription management:
http://yquem.inria.fr/cgi-bin/mailman/listinfo/caml-list
Archives: http://caml.inria.fr
Beginner's list: http://groups.yahoo.com/group/ocaml_beginners
Bug reports: http://caml.inria.fr/bin/caml-bugs

[Florent Monnier]

> I'm using Ocaml for an interval arithmetic application. I"m
> curious about
> what the Ocaml parser/compiler does to float constants. May I assume
> that for any constant I enter, eg. 3.1415... (for N digits of pi), that
> the compiler will give me a closest machine representable number?
> i.e., if I bound a constant by the previous and next floating point
> value to
> that given me by the compiler,
> will it always be the case that my original (mathematical) constant
> lies in that interval?

For the concrete exemple of PI, perhaps you could use a mathematical
calculation to get its best value:

# 4. *. atan 1. ;;
- : float = 3.14159265358979312
(but the last number '2' is wrong)

For other constants, maybe there are equivalent things?

[Gerd Stolpmann]

Am Samstag, den 30.09.2006, 14:25 -0400 schrieb Sean McLaughlin:
> Hello,
>
> I'm using Ocaml for an interval arithmetic application. I"m
> curious about
> what the Ocaml parser/compiler does to float constants. May I assume
> that for any constant I enter, eg. 3.1415... (for N digits of pi), that
> the compiler will give me a closest machine representable number?
> i.e., if I bound a constant by the previous and next floating point
> value to
> that given me by the compiler,
> will it always be the case that my original (mathematical) constant
> lies in that interval?

Don't think so. float_of_string is a wrapper around the strtod C
function, and the standard for this function does not require that. I
found this interesting note about how different OS implement string to
float conversion:

http://www.wrcad.com/linux_numerics.txt

To be sure your constants are the best possible you probably have to use
float_of_bits...

Gerd
--
------------------------------------------------------------
Gerd Stolpmann * Viktoriastr. 45 * 64293 Darmstadt * Germany
ge...@gerd-stolpmann.de http://www.gerd-stolpmann.de
Phone: +49-6151-153855 Fax: +49-6151-997714
------------------------------------------------------------

[Xavier Leroy]

> I'm using Ocaml for an interval arithmetic application. I"m curious
> about what the Ocaml parser/compiler does to float constants.

The lexer, parser and most of the compilers actually keep them as
strings, just like they were given in the source code. Conversion to
IEEE binary representation occurs:

- For the bytecode compiler ocamlc, when the bytecode is generated.
The conversion is performed by the float_of_string function, which is
a thin wrapper around the C library strtod() function, as Gerd
Stolpmann said.

- For the native-code compiler ocamlopt, some ports go the
float_of_string route, but most output the literal as a string in the
generated assembly code, letting the assembler do the conversion.
The assembler is likely to use strtod() or atof(), though.

> May I assume
> that for any constant I enter, eg. 3.1415... (for N digits of pi), that
> the compiler will give me a closest machine representable number?
> i.e., if I bound a constant by the previous and next floating point
> value to that given me by the compiler, will it always be the case
> that my original (mathematical) constant lies in that interval?

As Gerd said, the C standards do not guarantee this, so it really
depends on how well-implemented your C library is.

- Xavier Leroy

[Christophe Raffalli]

Sean McLaughlin a écrit :

> Hello,
>
> I'm using Ocaml for an interval arithmetic application. I"m curious
> about
> what the Ocaml parser/compiler does to float constants. May I assume
> that for any constant I enter, eg. 3.1415... (for N digits of pi), that
> the compiler will give me a closest machine representable number?
> i.e., if I bound a constant by the previous and next floating point
> value to
> that given me by the compiler,
> will it always be the case that my original (mathematical) constant
> lies in that interval?

This is not really an answer to your question but I wrote an OCaml binding to the function
to change the rounding mode of the proc. You can download it from (the fifth item)
http://www.lama.univ-savoie.fr/~raffalli/?page=soft&lang=en

(four rounding modes are availaible: closest, toward +inf, -inf and zero. The mode far from zero
is not available)

And to answer your question, the default rounding mode is "closest" in OCaml like in C (I think).

--
Christophe Raffalli
Université de Savoie
Batiment Le Chablais, bureau 21
73376 Le Bourget-du-Lac Cedex

tél: (33) 4 79 75 81 03
fax: (33) 4 79 75 87 42
mail: Christoph...@univ-savoie.fr
www: http://www.lama.univ-savoie.fr/~RAFFALLI
---------------------------------------------
IMPORTANT: this mail is signed using PGP/MIME
At least Enigmail/Mozilla, mutt or evolution
can check this signature. The public key is
stored on www.keyserver.net
---------------------------------------------

[Christophe Raffalli]

Sean McLaughlin a écrit :

> Hello,
>
> I'm using Ocaml for an interval arithmetic application. I"m curious
> about
> what the Ocaml parser/compiler does to float constants. May I assume
> that for any constant I enter, eg. 3.1415... (for N digits of pi), that
> the compiler will give me a closest machine representable number?
> i.e., if I bound a constant by the previous and next floating point
> value to
> that given me by the compiler,
> will it always be the case that my original (mathematical) constant
> lies in that interval?
>

By the way, float constants need to be written in hexadecimal and this
is missing to
the printf/scanf functions (it is what man printf says at least) ...
just compute how many decimals you need to write the exact value of
2^{-n} as a decimal float constant (0,5 0,25 0,125 0,625e-1 0,3125e-1 ...).

Then, I am not sure the algorithm used to parse the constant will give
you the correct value and you will never know how many decimals you need
to write your constant.

The only problem with float and hexadecimal is to learn the hexadecimals
of pi, e, ln 2, sqrt 2, ... ;-)

Christophe Raffalli

[Olivier Andrieu]

Christophe Raffalli [Tuesday 3 October 2006] :

>
> Sean McLaughlin a écrit :
> > Hello,
> >
> > I'm using Ocaml for an interval arithmetic application. I"m
> > curious about what the Ocaml parser/compiler does to float
> > constants. May I assume that for any constant I enter,
> > eg. 3.1415... (for N digits of pi), that the compiler will give
> > me a closest machine representable number? i.e., if I bound a
> > constant by the previous and next floating point value to that
> > given me by the compiler, will it always be the case that my
> > original (mathematical) constant lies in that interval?
> >
>
> By the way, float constants need to be written in hexadecimal and
> this is missing to the printf/scanf functions (it is what man
> printf says at least) ... just compute how many decimals you need
> to write the exact value of 2^{-n} as a decimal float constant (0,5
> 0,25 0,125 0,625e-1 0,3125e-1 ...).

float_of_string (which uses strtod) already knows how to read a
hexadecimal float and there's a hack to have the C printf do the writing:
,----
| # external format_float: string -> float -> string = "caml_format_float" ;;
| external format_float : string -> float -> string = "caml_format_float"
| # let hex_float = format_float "%a" ;;
| val hex_float : float -> string = <fun>
| # hex_float 1.25 ;;
| - : string = "0x1.4p+0"
| # float_of_string "0x1.4p+0" ;;
| - : float = 1.25
| # hex_float 0.1 ;;
| - : string = "0x1.999999999999ap-4"
`----

--
Olivier

[David Monniaux]

Xavier Leroy a écrit :

>> I'm using Ocaml for an interval arithmetic application. I"m curious
>> about what the Ocaml parser/compiler does to float constants.
>>

>> - For the native-code compiler ocamlopt, some ports go the
>> float_of_string route, but most output the literal as a string in the
>> generated assembly code, letting the assembler do the conversion.
>> The assembler is likely to use strtod() or atof(), though.

> As Gerd said, the C standards do not guarantee this, so it really
> depends on how well-implemented your C library is.

In C, if you are really concerned about getting your exact constants up
to the last bit, you should use the hexadecimal float notation
introduced in C99. You should also use "%a" or "%A" to print out your
floating-point values in hexadecimal if you intend to reload them on
another IEEE-754 system in a portable way.

(Unfortunately on OCaml, %a is already used for another use... Xavier,
what do you think of adding the hexadecimal float conversion to the
family of printf functions?)

Anyway, "IEEE-754 compliant", as used by common vendors, is a fairly
vague claim. For instance, contrary to a commonly held belief, the same
C or Caml program ran on two different major platforms (say, IA32 vs
AMD64 or IA32 vs PowerPC) can give different results, possibly even
depending on compiler options, and this has nothing to do with how they
parse constants... You might wish to read my preprint:
http://www.di.ens.fr/~monniaux/biblio/pitfalls_of_floating_point.pdf

The only widespread language that I know about that has fixed semantics
for floating-point is Java if the 'strictfp' modifier is used.
(Originally, Java was supposed to have a fully fixed semantics, but they
realized that implementing it over the 80387 compatible floating-point
stack was difficult and that it also precluded some optimizations, so
they relaxed it by default and introduced 'strictfp'. Goodbye
predictability.)

[David Monniaux]

In the Astrée project, we use some small C functions and OCaml wrappers
to change rounding modes. It's not easy (some systems have
fpgetround/fpsetround, others fegetround/fesetround, some have none and
you have to use assembly, and old versions of GNU libc on IA32 have
fegetround/fesetround that change rounding for the 387 FPU but not for
SSE, which may break some external C libraries).
[ Hey, we should publish our libraries. ]

In addition, some C libraries don't work properly when rounding is not
"round to nearest". I filed bugs against FreeBSD libc (printf doesn't
work properly with some arguments) and GNU libc (pow() does not work
properly).

Rounding modes other than the default "round to nearest" tend to be
largely untested on many systems and libraries, because few people use
them. #pragma FENV_ACCESS is not even implemented on gcc;
-frounding-math is supposed to signal rounding modes other than round to
nearest, but the documentation states that it's not even sure it works...

Not that I want to advertise myself too much, but I think you should read
http://www.di.ens.fr/~monniaux/biblio/pitfalls_of_floating_point.pdf
if you are concerned about interactions between compilers and
floating-point.