[Caml-list] Approximations when converting from string to float

[Caml-list] Approximations when converting from string to float

[Tung Vu Xuan]

[-- Attachment #1: Type: text/plain, Size: 264 bytes --]

Hi,

I am having a problem when converting from string to float as described in
this question:
http://stackoverflow.com/questions/40219852/string-to-float-in-ocaml-over-and-under-approximation
.
Your helps are mostly appreciated! Thanks a lot!

Best Regards,
Tung

[-- Attachment #2: Type: text/html, Size: 557 bytes --]

Re: [Caml-list] Approximations when converting from string to float

[Jacques-Henri Jourdan]

Hi,

I would say that the easiest is to use the ZArith library: you convert 
from tring to Z.t, then from Z.t to float.

By doing so,you get only an approximation, but you do not know whether 
this is a lower or upper rounding. You need then to convert back to Z.t 
and compare with the original number. Depending on whether this is a 
lower, exact or upper rounding, you can then get your interval by adding 
or removing an ulp....


-- 
JH

On 10/24/2016 04:01 PM, Tung Vu Xuan wrote:
> Hi,
>
> I am having a problem when converting from string to float as 
> described in this question:
> http://stackoverflow.com/questions/40219852/string-to-float-in-ocaml-over-and-under-approximation. 
>
> Your helps are mostly appreciated! Thanks a lot!
>
> Best Regards,
> Tung

RE: [Caml-list] Approximations when converting from string to float

[Soegtrop, Michael]

[-- Attachment #1: Type: text/plain, Size: 2936 bytes --]

Dear Tung,

converting strings to floats and keeping exact bounds is not that complicated. The steps are:


1.)    Convert the digits to an integer. You can stop when you have 1 or 2 digits more than the mantissa can hold. E.g. if you create a 32 or 64 bit float, a 32 bit or 64 bit integer is enough, because it is quite a bit longer than the float mantissa. The lower and upper bounds are the truncated value and the truncated value + 1. If the decimal number fits into the mantissa, the bounds are the same. For the truncated digits you just have to count them to correct the decimal exponent. Also you need to take care of the decimal point to get the proper decimal exponent.

2.)    Convert the integer to a float. If the integer is smaller than the mantissa length, this is trivial. If it is longer, you divide by 2 until it fits. The lower limit from step 1 is rounded down, the upper limit is rounded up while dividing. The number of divide steps is limited by the difference of the integer size you used and the floats mantissa size (that is the exponent size of the float +1). To make it faster, you can also do a binary intersection search of the right power of 2, but it might not be worth it, because dividing integers by 2 is cheap and the max loop count is small.

3.)    In order to handle the decimal point and decimal exponents, you need tables of 10^(2^n). Use the binary representation of the decimal exponent to find out which powers to multiply. For each power you need a lower and upper bound value (as a floating point number). You can compute these tables with libraries supporting higher precision numbers and store the numbers (best in hex float notation – which is standard C). In an IEEE floating point implementation, you can choose the rounding direction (up down middle). When you multiply the factors of the lower and upper bounds, you set the rounding direction accordingly. The C function for this is fesetround. The exponent tables are fairly small, because 10^(2^n) is larger than the maximal floating point number for rather modest n (for double precision, n=8 = 10^256 should be the largest value you need).

That’s it. Since you multiply several numbers (at most 10), each of which is rounded up /down, your range might be a bit larger than necessary. For getting the minimum bound you need multi precision arithmetic, but the method above works with just target precision arithmetic and gives reliable bounds.

Note that for integers and some other exactly representable numbers, this method gives identical lower and upper bounds.

Best regards,

Michael

Intel Deutschland GmbH
Registered Address: Am Campeon 10-12, 85579 Neubiberg, Germany
Tel: +49 89 99 8853-0, www.intel.de
Managing Directors: Christin Eisenschmid, Christian Lamprechter
Chairperson of the Supervisory Board: Nicole Lau
Registered Office: Munich
Commercial Register: Amtsgericht Muenchen HRB 186928

[-- Attachment #2: Type: text/html, Size: 9479 bytes --]

RE: [Caml-list] Approximations when converting from string to float

[Soegtrop, Michael]

[-- Attachment #1: Type: text/plain, Size: 551 bytes --]

P.S.:

a good intro to the not so well known hex floating point notation of C: http://www.exploringbinary.com/hexadecimal-floating-point-constants/. You should use this to get the 10^(2^n) constants right.

Best regards,

Michael
Intel Deutschland GmbH
Registered Address: Am Campeon 10-12, 85579 Neubiberg, Germany
Tel: +49 89 99 8853-0, www.intel.de
Managing Directors: Christin Eisenschmid, Christian Lamprechter
Chairperson of the Supervisory Board: Nicole Lau
Registered Office: Munich
Commercial Register: Amtsgericht Muenchen HRB 186928

[-- Attachment #2: Type: text/html, Size: 5232 bytes --]

Re: [Caml-list] Approximations when converting from string to float

[Jacques-Henri Jourdan]

[-- Attachment #1: Type: text/plain, Size: 3271 bytes --]

Hi all,

On 10/25/2016 01:20 PM, Soegtrop, Michael wrote:
>
> Dear Tung,
>
> converting strings to floats and keeping exact bounds is not that 
> complicated. The steps are:
>
> 1.)Convert the digits to an integer. You can stop when you have 1 or 2 
> digits more than the mantissa can hold. E.g. if you create a 32 or 64 
> bit float, a 32 bit or 64 bit integer is enough, because it is quite a 
> bit longer than the float mantissa. The lower and upper bounds are the 
> truncated value and the truncated value + 1. If the decimal number 
> fits into the mantissa, the bounds are the same. For the truncated 
> digits you just have to count them to correct the decimal exponent. 
> Also you need to take care of the decimal point to get the proper 
> decimal exponent.
>

The issue is that values of the remaining digits *do* matter for the 
right result to be computed, not just the number of them.

Imagine if the integer you get is 2^70 or 2^70-1. Then the lower bounds 
of the intervals have to be different (2^70 is exactly a double 
precision float, so 2^70-1 maps to the preceding float), but notice that 
these two integers differ only by their last digit.


> 2.)Convert the integer to a float. If the integer is smaller than the 
> mantissa length, this is trivial. If it is longer, you divide by 2 
> until it fits. The lower limit from step 1 is rounded down, the upper 
> limit is rounded up while dividing. The number of divide steps is 
> limited by the difference of the integer size you used and the floats 
> mantissa size (that is the exponent size of the float +1). To make it 
> faster, you can also do a binary intersection search of the right 
> power of 2, but it might not be worth it, because dividing integers by 
> 2 is cheap and the max loop count is small.
>
> 3.)In order to handle the decimal point and decimal exponents, you 
> need tables of 10^(2^n). Use the binary representation of the decimal 
> exponent to find out which powers to multiply. For each power you need 
> a lower and upper bound value (as a floating point number). You can 
> compute these tables with libraries supporting higher precision 
> numbers and store the numbers (best in hex float notation – which is 
> standard C). In an IEEE floating point implementation, you can choose 
> the rounding direction (up down middle). When you multiply the factors 
> of the lower and upper bounds, you set the rounding direction 
> accordingly. The C function for this is fesetround. The exponent 
> tables are fairly small, because 10^(2^n) is larger than the maximal 
> floating point number for rather modest n (for double precision, n=8 = 
> 10^256 should be the largest value you need).
>
> That’s it. Since you multiply several numbers (at most 10), each of 
> which is rounded up /down, your range might be a bit larger than 
> necessary. For getting the minimum bound you need multi precision 
> arithmetic, but the method above works with just target precision 
> arithmetic and gives reliable bounds.
>

That's the point : if you really want the tight interval, then you need 
arbitrary precision arithmetic.

Moreover, it is not possible to implement interval arithmetic with 
OCaml, since you cannot change the rounding mode without a bit of C...

-- 
JH

[-- Attachment #2: Type: text/html, Size: 9786 bytes --]

RE: [Caml-list] Approximations when converting from string to float

[Soegtrop, Michael]

[-- Attachment #1: Type: text/plain, Size: 2032 bytes --]

Dear Jacques-Henri,

if the number is cut off after n digits, the upper and lower bounds I suggested are the (truncated integer)*10^(#cut-off-digits) and (truncated integer+1)*10^(#cut-off-digits). The true number is obviously between these two. Since the integer has higher precision than the mantissa in the case of cut off, this is only a fraction of a mantissa bit. The errors you get by multiplying with the powers of 10 is likely larger in most cases.

In the case you mentioned, 2^70=1180591620717411303424 and 32 bit one would truncate after

1180591620 * 10^12

Maxima gives

is(1180591620 * 10^12 <= 2^70);
true

is(1180591621 * 10^12 >= 2^70);
true

As I said, this method doesn’t give the tightest bounds possible, but it gives you true upper and lower bounds without multi precision arithmetic. It gives 0 intervals e.g. for integers which fit in the mantissa, but not for large exact powers of 2.

> Moreover, it is not possible to implement interval arithmetic with OCaml, since you cannot change the rounding mode without a bit of C...

It should be possible to make the powers of 10 tables such that this works even with round to nearest, but it would likely be easier to make a C library for this.

Of cause it is a good question if I/O is the right place to save performance and waste precision. On the other hand you lose only 2 or 3 bits and you know what the precision is. I would think most applications are such that the required precision is not exactly an IEEE precision, so that these 2 or 3 bits are ok. I used this method once in an embedded platform, where multi precision arithmetic was not an option, and this was a good compromise.

Best regards,

Michael
Intel Deutschland GmbH
Registered Address: Am Campeon 10-12, 85579 Neubiberg, Germany
Tel: +49 89 99 8853-0, www.intel.de
Managing Directors: Christin Eisenschmid, Christian Lamprechter
Chairperson of the Supervisory Board: Nicole Lau
Registered Office: Munich
Commercial Register: Amtsgericht Muenchen HRB 186928

[-- Attachment #2: Type: text/html, Size: 9577 bytes --]

Re: [Caml-list] Approximations when converting from string to float

[Tung Vu Xuan]

[-- Attachment #1: Type: text/plain, Size: 3780 bytes --]

Hi,

Thanks for suggestions. Since I am using a library for interval arithmetic
[1], your idea becomes easier for me.

let intv_of_big_int bi =
    let n = Big_int.num_bits_big_int bi in
    if n <= 53 then
      let bi_float = Big_int.float_of_big_int bi in
      {low=bi_float;high=bi_float}
    else begin
      let n = n - 53 in
      (* Extract top 53 bits of x *)
      let top = Big_int.shift_right_big_int bi n in

      let top_float = Big_int.float_of_big_int top in
      (* interval [2, 2] *)
      let two_I = {low=2.;high=2.} in
      (* interval: [2, 2]**n *)
      let two_I_to_n = pow_I_i two_I n in

      (* compute upper bound of mantissa, represented as an intv *)
      (* interval arithmetic is needed because rounding error might occur *)
      let mantissa_upper_bound_intv = {low=top_float;high=top_float} +$
one_I in

      (* intv of mantissa *)
      let mantissa_intv = {low=top_float; high=mantissa_upper_bound_intv.high}
in

      (* compute final interval *)
      mantissa_intv *$ two_I_to_n
    end

[1] Alliot, J.M., Gotteland, J.B., Vanaret, C., Durand, N., Gianazza, D.:
Implementing an interval computation library for OCaml on x86/amd64
architectures. In: ICFP. ACM (2012)

Thanks again,
Tung.

On Tue, Oct 25, 2016 at 5:28 PM, Soegtrop, Michael <
michael.soegtrop@intel.com> wrote:

> Dear Jacques-Henri,
>
>
>
> if the number is cut off after n digits, the upper and lower bounds I
> suggested are the (truncated integer)*10^(#cut-off-digits) and (truncated
> integer+1)*10^(#cut-off-digits). The true number is obviously between
> these two. Since the integer has higher precision than the mantissa in the
> case of cut off, this is only a fraction of a mantissa bit. The errors you
> get by multiplying with the powers of 10 is likely larger in most cases.
>
>
>
> In the case you mentioned, 2^70=1180591620717411303424 and 32 bit one
> would truncate after
>
>
>
> 1180591620 * 10^12
>
>
>
> Maxima gives
>
>
>
> is(1180591620 * 10^12 <= 2^70);
>
> true
>
>
>
> is(1180591621 * 10^12 >= 2^70);
>
> true
>
>
>
> As I said, this method doesn’t give the tightest bounds possible, but it
> gives you true upper and lower bounds without multi precision arithmetic.
> It gives 0 intervals e.g. for integers which fit in the mantissa, but not
> for large exact powers of 2.
>
>
>
> > Moreover, it is not possible to implement interval arithmetic with
> OCaml, since you cannot change the rounding mode without a bit of C...
>
>
>
> It should be possible to make the powers of 10 tables such that this works
> even with round to nearest, but it would likely be easier to make a C
> library for this.
>
>
>
> Of cause it is a good question if I/O is the right place to save
> performance and waste precision. On the other hand you lose only 2 or 3
> bits and you know what the precision is. I would think most applications
> are such that the required precision is not exactly an IEEE precision, so
> that these 2 or 3 bits are ok. I used this method once in an embedded
> platform, where multi precision arithmetic was not an option, and this was
> a good compromise.
>
>
>
> Best regards,
>
>
>
> Michael
>
> Intel Deutschland GmbH
> Registered Address: Am Campeon 10-12, 85579 Neubiberg, Germany
> Tel: +49 89 99 8853-0, www.intel.de
> Managing Directors: Christin Eisenschmid, Christian Lamprechter
> Chairperson of the Supervisory Board: Nicole Lau
> Registered Office: Munich
> Commercial Register: Amtsgericht Muenchen HRB 186928
>



-- 
Tung Vu Xuan,
Japan Advanced Institute of Science and Technology,
Mobile: (+81) 080 4259 9135 -  (+84) 01689934381
Hometown: Bac Ninh
Email: toilatung90@gmail.com

[-- Attachment #2: Type: text/html, Size: 9138 bytes --]

Re: [Caml-list] Approximations when converting from string to float

[Jacques-Henri Jourdan]

[-- Attachment #1: Type: text/plain, Size: 4120 bytes --]

Ok.

But really, you should use ZArith, which is much more modern, faster and 
well-maintained than Big_int.

On 10/25/2016 06:58 PM, Tung Vu Xuan wrote:
> Hi,
>
> Thanks for suggestions. Since I am using a library for interval 
> arithmetic [1], your idea becomes easier for me.
>
> let intv_of_big_int bi =
>     let n = Big_int.num_bits_big_int bi in
>     if n <= 53 then
>       let bi_float = Big_int.float_of_big_int bi in
>       {low=bi_float;high=bi_float}
>     else begin
>       let n = n - 53 in
>       (* Extract top 53 bits of x *)
>       let top = Big_int.shift_right_big_int bi n in
>
>       let top_float = Big_int.float_of_big_int top in
>       (* interval [2, 2] *)
>       let two_I = {low=2.;high=2.} in
>       (* interval: [2, 2]**n *)
>       let two_I_to_n = pow_I_i two_I n in
>
>       (* compute upper bound of mantissa, represented as an intv *)
>       (* interval arithmetic is needed because rounding error might 
> occur *)
>       let mantissa_upper_bound_intv = {low=top_float;high=top_float} 
> +$ one_I in
>
>       (* intv of mantissa *)
>       let mantissa_intv = {low=top_float; 
> high=mantissa_upper_bound_intv.high} in
>       (* compute final interval *)
>       mantissa_intv *$ two_I_to_n
>     end
>
> [1] Alliot, J.M., Gotteland, J.B., Vanaret, C., Durand, N., Gianazza, 
> D.: Implementing an interval computation library for OCaml on 
> x86/amd64 architectures. In: ICFP. ACM (2012)
>
> Thanks again,
> Tung.
>
> On Tue, Oct 25, 2016 at 5:28 PM, Soegtrop, Michael 
> <michael.soegtrop@intel.com <mailto:michael.soegtrop@intel.com>> wrote:
>
>     Dear Jacques-Henri,
>
>     if the number is cut off after n digits, the upper and lower
>     bounds I suggested are the (truncated
>     integer)*10^(#cut-off-digits) and (truncated
>     integer+1)*10^(#cut-off-digits). The true number is obviously
>     between these two. Since the integer has higher precision than the
>     mantissa in the case of cut off, this is only a fraction of a
>     mantissa bit. The errors you get by multiplying with the powers of
>     10 is likely larger in most cases.
>
>     In the case you mentioned, 2^70=1180591620717411303424 and 32 bit
>     one would truncate after
>
>     1180591620 * 10^12
>
>     Maxima gives
>
>     is(1180591620 * 10^12 <= 2^70);
>
>     true
>
>     is(1180591621 * 10^12 >= 2^70);
>
>     true
>
>     As I said, this method doesn’t give the tightest bounds possible,
>     but it gives you true upper and lower bounds without multi
>     precision arithmetic. It gives 0 intervals e.g. for integers which
>     fit in the mantissa, but not for large exact powers of 2.
>
>     > Moreover, it is not possible to implement interval arithmetic
>     with OCaml, since you cannot change the rounding mode without a
>     bit of C...
>
>     It should be possible to make the powers of 10 tables such that
>     this works even with round to nearest, but it would likely be
>     easier to make a C library for this.
>
>     Of cause it is a good question if I/O is the right place to save
>     performance and waste precision. On the other hand you lose only 2
>     or 3 bits and you know what the precision is. I would think most
>     applications are such that the required precision is not exactly
>     an IEEE precision, so that these 2 or 3 bits are ok. I used this
>     method once in an embedded platform, where multi precision
>     arithmetic was not an option, and this was a good compromise.
>
>     Best regards,
>
>     Michael
>
>     Intel Deutschland GmbH
>     Registered Address: Am Campeon 10-12, 85579 Neubiberg, Germany
>     Tel: +49 89 99 8853-0, www.intel.de <http://www.intel.de>
>     Managing Directors: Christin Eisenschmid, Christian Lamprechter
>     Chairperson of the Supervisory Board: Nicole Lau
>     Registered Office: Munich
>     Commercial Register: Amtsgericht Muenchen HRB 186928
>
>
>
>
> -- 
> Tung Vu Xuan,
> Japan Advanced Institute of Science and Technology,
> Mobile: (+81) 080 4259 9135 -  (+84) 01689934381
> Hometown: Bac Ninh
> Email: toilatung90@gmail.com <mailto:toilatung90@gmail.com>


[-- Attachment #2: Type: text/html, Size: 11834 bytes --]

Re: [Caml-list] Approximations when converting from string to float

[Xavier Leroy]

On 10/25/2016 06:58 PM, Tung Vu Xuan wrote:

> Thanks for suggestions. Since I am using a library for interval arithmetic
> [1], your idea becomes easier for me.
> let intv_of_big_int bi = 
>     let n = Big_int.num_bits_big_int bi in
>     if n <= 53 then 
>       let bi_float = Big_int.float_of_big_int bi in
>       {low=bi_float;high=bi_float}
>     else begin
>       let n = n - 53 in
>       (* Extract top 53 bits of x *)
>       let top = Big_int.shift_right_big_int bi n in

Yes, that's how I would go about it, too.  But probably you don't need
the interval arithmetic in the end: both "top" and "top + 1" convert
to floats without rounding (*), and the scaling by 2^n can be achieved
with "ldexp".

(*) If "top" is 2^53 - 1, "top + 1" is 2^53, which converts exactly.
Otherwise, "top + 1" has at most 53 digits and converts exactly as well.

Jacques-Henri Jourdan adds:

> But really, you should use ZArith, which is much more modern, faster
> and well-maintained than Big_int.

Agreed.  And if you must use Big_int.float_of_big_int, make sure
you're using OCaml 4.03.0 or later, because older versions were flaky.

- Xavier Leroy

Re: [Caml-list] Approximations when converting from string to float

[Tung Vu Xuan]

[-- Attachment #1: Type: text/plain, Size: 436 bytes --]

Hi Xavier,


> But probably you don't need
> the interval arithmetic in the end: both "top" and "top + 1" convert
> to floats without rounding (*),


Yes, you're right :).


> and the scaling by 2^n can be achieved
> with "ldexp".
>

What is the semantics of "ldexp" in bounds cases? For example, I tried
"ldexp max_float max_int;;"  which returns "8.98846567431157854e+307".
Shouldn't we expect something like "infinity"?

Best,
Tung.

[-- Attachment #2: Type: text/html, Size: 1167 bytes --]

RE: [Caml-list] Approximations when converting from string to float

[Soegtrop, Michael]

[-- Attachment #1: Type: text/plain, Size: 968 bytes --]

Dear Tung,

if you are converting integers to floats (no decimal fractions) and if you can use an arbitrary precision integer library to convert from arbitrary size decimal to binary and do the truncation in the binary domain, there is no need for interval arithmetic – as Xavier also pointed out. The method I outlined is intended for limited precision integer and limited precision float arithmetic, in which case you have to throw away decimal digits and handle imprecise representations of powers of 10. In this case using a an interval arithmetic library e.g. to represent the powers of 10 would be a good solution.

Best regards,

Michael
Intel Deutschland GmbH
Registered Address: Am Campeon 10-12, 85579 Neubiberg, Germany
Tel: +49 89 99 8853-0, www.intel.de
Managing Directors: Christin Eisenschmid, Christian Lamprechter
Chairperson of the Supervisory Board: Nicole Lau
Registered Office: Munich
Commercial Register: Amtsgericht Muenchen HRB 186928

[-- Attachment #2: Type: text/html, Size: 3829 bytes --]