CODEWITHSUNDEEP
fortranfloating-point-conversion
user3337813
user3337813

Reputation: 476

Hexadecimal floating point in fortran

Is there an equivalent for the 'a' format specifier known from C in Fortran?

C Example: 

printf("%a\n",43.1e6); // 0x1.48d3bp+25

Exporting floating point numbers in hexadecimal format prevents rounding errors. While the rounding errors are usually negligible, it is still advantageous to be able to restore a saved value exactly. Note, that the hexadecimal representation produced by printf is portable and human readable.

How can I export and parse floating point numbers in Fortran like I do in C using the 'a' specifier?

Upvotes: 1

Views: 999

Answers (3)

Michal
Michal

Reputation: 91

Fortran 2018 standard introduced 'EX' descriptor which does exactly what '%a' in C. Recent Intel ifort compiler does support that, GNU Fortran - not yet as of version 12.2.0

Upvotes: 1

user4490638
user4490638

Reputation:

Another option would be to let the C library do your work for you and interface via C binding. This rather depends on a modern compiler (some F2003 features used).

module x
  use, intrinsic :: iso_c_binding
  private
  public :: a_fmt
  interface
     subroutine doit(a, dest, n) bind(C) 
       import
       real(kind=c_double), value :: a
       character(kind=c_char), intent(out) :: dest(*)
       integer, value :: n
     end subroutine doit
  end interface

  interface a_fmt
     module procedure a_fmt_float, a_fmt_double
  end interface a_fmt

contains
  function a_fmt_float(a) result(res)
    real(kind=c_float), intent(in) :: a
    character(len=:), allocatable :: res
    res = a_fmt_double (real(a, kind=c_double))
  end function a_fmt_float

  function a_fmt_double(a) result(res)
    real(kind=c_double), intent(in) :: a
    character(len=:), allocatable :: res
    character(len=30) :: dest
    integer :: n
    call doit (a, dest, len(dest))
    n = index(dest, achar(0))
    res = dest(1:n)
  end function a_fmt_double

end module x

program main
  use x
  implicit none
  double precision :: r
  integer :: i
  r = -1./3.d0
  do i=1,1030
     print *,a_fmt(r)
     r = - r * 2.0
  end do
end program main

#include <stdio.h>

void doit (double a, char *dest, int n)
{
  snprintf(dest, n-1, "%a", a);
}

Upvotes: 1

user4490638
user4490638

Reputation:

If you want to have full precision, the best way is to use unformatted files, such as this:

program main
  real :: r
  integer :: i
  r = -4*atan(1.)
  open(20,access="stream")
  write (20) r
  close(20)
end program main

(I used stream access, which is new to Fortran 2003, because it is usually less confusing than normal unformatted access). You can then use, for example, od -t x1 fort.20 to look at this as a hex dump.

You can also use TRANSFER to copy the bit pattern to an integer and then use the Z edit descriptor.

If you really want to mimic the %a specifier, you'll have to roll your own. Most machines now use IEEE format. Use TRANSFER for copying the pattern to an integer, then pick that apart using IAND (and multiplications or divisions by powers of two for shifting).

Upvotes: 1

Related Questions