• Lessons learned from the implementation of FMPFR

    From Thomas Koenig@21:1/5 to All on Sun Jul 31 14:48:34 2022
    Here's something I wrote up on implementing FMPFR. Enjoy!

    This documents some of the design choices for, and lessons from, the implementation of MPFR, reflecting both strenghs and weaknesses of
    Fortran and the other tools that are used.

    The aim of the project was to make MPFR easy to use for a Fortran
    programmer by providing a "thick binding", using overloading, elemental functions and subroutines and user-defined I/O to make the code look as
    close as possible to "natural" Fortran code, and to facilitate
    conversion of existing code to use the library.

    The big strength of Fortran is that this is indeed possible.

    Why generate the code with a script?

    Automated generation of interfaces to C

    MPFR is a large library, with more than 300 functions. Picking out those functions which are required and translating them for hand would have
    been time-consuming and error-prone. Fortunately, the interface to MPFR
    is very regular, and its interfaces from the Texinfo source to the documentation, for example

    @deftypefun void mpfr_init2 (mpfr_t @var{x}, mpfr_prec_t @var{prec})

    are easy to find with grep and easy to parse with a script language.
    Also, functions which are not suitable to use with Fortran (for example
    those containing unsigned integers) can easily be removed. In this
    project, perl is used, but other script languages would be equally
    suitable as long as they have associative arrays.

    Generating Fortran code for similar functions

    MPFR contains many functions which serve a similar purpose, with
    different data types. As an example, mpfr_add_si adds a variable of type
    mpfr_t and an unsigned long, and mpfr_add_d adds an mpfr_t and a double.
    Also, operator overloading should work for both a+b and b+a where a is a multi-precision variable and b is one of Fortran's types.

    This is also an advantage when a change is made to these functions. It
    is then enough to change one place, and not change the customary n-1
    places when n would have been needed.

    Lack of generics

    A use would probably like to write a+3 instead of a + 3_c_long,
    especially since ISO C binding is not otherwise needed on the user side.
    This was solved by generating appropriate functions from the script.

    Why the configure scripts and preprocessing?

    While a Fortran processor that implements C binding has to provide both
    c_int and c_long, these might be the same, and an interface like

    interface foo
    module procedure foo_int
    module procedure foo_long
    end interface

    subroutine foo_int (a)
    integer(c_int) :: a
    ...
    subroutine foo_long
    integer(c_long) :: a

    would be rejected. There is no way known to the author do solve this
    issue within the bounds of the Fortran standard.

    So, there is a need for a preprocessing step. This was done using the de-facto-standard of using the C preprocessor, but the preprocessore
    needed to get its information from somewhere.

    For this, autoconf and the rest of the GNU autotools suggested itself - checking for sizes of C types is implemented there.

    Autoconf is now also used for checking features (not all compilers
    support user-defined I/O) and for getting information on the internal
    types of MPFR. Also, the testsuite uses autotools.

    For the author, autotools are not easy to use because of documentation
    which is aimed at describing the features, and less towards how to reach specific goals.

    Other issues

    User-defined I/O

    I found user-defined I/O to be difficult. For list-directed output, it
    is usable (if not implemented on all compiler versions I tried this on).
    Input appears to be poorly specified, and different compiler writers
    appear to have different ideas what is specified. I therefore had to
    leave that out.

    Generics with two interger arguments

    One design issue was how to design generics where two integer arguments
    (the precision and the rounding mode) could both be optional arguments,
    and the compiler could differentiate between them. The solution
    implemented is somewhat of a hack, but if there is a cleaner solution, I
    would like to hear about it.

    The constants for mpfr_rnd_t are specified as integer(int8), and the
    precision can be specified as integer(c_short), integer(c_int) or integer(c_long).

    For conversion from a string to a type(fmpfr). the interface then
    contains the functions

    function fun_set_str (s, rnd) result(rop)
    type (fmpfr) :: rop
    integer (kind=int8), optional :: rnd
    character(kind=c_char,len=*), intent(in) :: s
    !...
    function fun_set_str_long (s, prec, rnd) result(rop)
    type (fmpfr) :: rop
    character(kind=c_char,len=*), intent(in) :: s
    integer (c_long), intent(in) :: prec
    integer (kind=int8), intent(in), optional :: rnd

    plus two versions for integer(c_int) and integer(c_short).

    Pure functions and return codes

    Many MPFR functions return a value indicating the direction of rounding.
    In order to be able to use them in an elemental procedure, they have to
    be declared pure. However, an optimizting compiler will remove a pure
    function call like

    rc = mpfr_add_si (rop&mp, op1%mp, op2, rnd_val)

    when no further use is made of rc. While it is normal style in C to
    discard the return value of a function, Fortran C binding does not allow
    this. Therefore, the line above was replaced with

    call fmpfr_add_si (rop%mp, op1%mp, op2, rnd_val)

    where fmpfr_add_si is a glue function containing only the single call to mpfr_add_si with the same arguments, which is translated into a single
    jump.

    ELEMENTAL and RECURSIVE

    In Fortran 2008, a procedure could not both be elemental and recursive.
    This restriction was lifted in Fortran 2018, but that is not implemented
    by many compilers up to now.

    Array intrinsics

    It would have been nice to be able to use array intrinsics like MAXLOC
    or CSHIFT. Unfortunately, it is not possible to write these intrinsics
    in Fortran unless a version is written for each rank.

    The extended C interoperability of Fortran 2018 could be used to write
    such an implementation (as libgfortran shows), but this would be a much
    larger job than the current implementation, and the code from
    libgfortran cannot be used due to the difference in license.

    --- SoupGate-Win32 v1.05
    * Origin: fsxNet Usenet Gateway (21:1/5)
  • From Steven G. Kargl@21:1/5 to Thomas Koenig on Tue Aug 2 15:33:35 2022
    On Sun, 31 Jul 2022 14:48:34 +0000, Thomas Koenig wrote:

    Here's something I wrote up on implementing FMPFR. Enjoy!

    This documents some of the design choices for, and lessons from, the implementation of MPFR, reflecting both strenghs and weaknesses of
    Fortran and the other tools that are used.


    Thomas, This looks really interesting. I'll have to give it
    a workout on the weekend. I'm subscribed to the MPFR mailing
    list and saw an email you sent there. Are you planning to
    integrate/contribute this to MPFR? If gfortran (or other
    modern compiler is available), it would be nice to have a
    Fortran module is built when building MPFR?

    --
    steve

    --- SoupGate-Win32 v1.05
    * Origin: fsxNet Usenet Gateway (21:1/5)
  • From Thomas Koenig@21:1/5 to Steven G. Kargl on Tue Aug 2 16:09:54 2022
    Steven G. Kargl <sgk@REMOVEtroutmask.apl.washington.edu> schrieb:
    On Sun, 31 Jul 2022 14:48:34 +0000, Thomas Koenig wrote:

    Here's something I wrote up on implementing FMPFR. Enjoy!

    This documents some of the design choices for, and lessons from, the
    implementation of MPFR, reflecting both strenghs and weaknesses of
    Fortran and the other tools that are used.


    Thomas, This looks really interesting. I'll have to give it
    a workout on the weekend.

    Great! Please drop me a line if you encounter bugs.

    I'm subscribed to the MPFR mailing
    list and saw an email you sent there. Are you planning to integrate/contribute this to MPFR? If gfortran (or other
    modern compiler is available), it would be nice to have a
    Fortran module is built when building MPFR?

    I hadn't considered that, but if there's interest, that would
    indeed be cool.

    Licensing should not be a problem, I can easily make this a double
    license, either MIT or LGPL.

    --- SoupGate-Win32 v1.05
    * Origin: fsxNet Usenet Gateway (21:1/5)
  • From Arjen Markus@21:1/5 to All on Wed Aug 3 00:02:32 2022
    I read this with interest, but I do not quite understand the issue with the generics with two integers. Could you explain this further via the C interface? (I did not find that documentation immediately so gave up after 60 seconds ;)). Also, with respect
    to array intrinsics, you could simply start with one-dimensional versions, as I imagine that will be the most common use case.

    Regards,

    Arjen

    --- SoupGate-Win32 v1.05
    * Origin: fsxNet Usenet Gateway (21:1/5)
  • From Thomas Koenig@21:1/5 to Arjen Markus on Wed Aug 3 17:06:40 2022
    Arjen Markus <arjen.markus895@gmail.com> schrieb:

    I read this with interest, but I do not quite understand the
    issue with the generics with two integers. Could you explain this
    further via the C interface? (I did not find that documentation
    immediately so gave up after 60 seconds ;)).

    The point was that I wanted the user to be able to write

    - fmpfr ("1.2") (default digits, default rounding)
    - fmpfr ("1.2",100) (100 binary digits, default rounding)
    - fmpfr ("1.2", mpfr_rndz) (default digits, round towards towards zero)
    - fmpfr ("1.2", 100, mpfr_rndz) (100 binary digits, round towards zero)

    and I was a little bit stuck in my head because I thought that
    mpfr_rnd_t had to be an integer (it is a C enum), so I chose an 8-bit
    int for it.

    However, a better solution is actually rather simple: Make a
    derived type and declare parameters with the values of the enum.
    Duh. I think I will implement that next.

    Source code compatibility yes, but not yet binary compatibility :-)

    Nice thing is that I the code generation is in a script, I will
    only have to change a few places.

    Also, with respect
    to array intrinsics, you could simply start with one-dimensional
    versions, as I imagine that will be the most common use case.

    That would be relatively easy, but I would prefer a more general
    solution. Maybe a combination of SELECT RANK and pointer rank
    remapping would work for this to be doable in Fortran. Doing
    this in C would be a) no fun and b) more or less a duplication
    of libgfortran, which I would hate to do.

    --- SoupGate-Win32 v1.05
    * Origin: fsxNet Usenet Gateway (21:1/5)
  • From Steve Lionel@21:1/5 to Thomas Koenig on Wed Aug 3 19:32:23 2022
    On 8/3/2022 1:06 PM, Thomas Koenig wrote:
    and I was a little bit stuck in my head because I thought that
    mpfr_rnd_t had to be an integer (it is a C enum), so I chose an 8-bit
    int for it.

    However, a better solution is actually rather simple: Make a
    derived type and declare parameters with the values of the enum.
    Duh. I think I will implement that next.

    In Fortran 2023, you will be able to declare "real enums" that are typed
    for this. See https://j3-fortran.org/doc/year/21/21-110r1.txt
    --
    Steve Lionel
    ISO/IEC JTC1/SC22/WG5 (Fortran) Convenor
    Retired Intel Fortran developer/support
    Email: firstname at firstnamelastname dot com
    Twitter: @DoctorFortran
    LinkedIn: https://www.linkedin.com/in/stevelionel
    Blog: https://stevelionel.com/drfortran
    WG5: https://wg5-fortran.org

    --- SoupGate-Win32 v1.05
    * Origin: fsxNet Usenet Gateway (21:1/5)
  • From Thomas Koenig@21:1/5 to Steve Lionel on Thu Aug 4 05:11:23 2022
    Steve Lionel <steve@seesignature.invalid> schrieb:
    On 8/3/2022 1:06 PM, Thomas Koenig wrote:
    and I was a little bit stuck in my head because I thought that
    mpfr_rnd_t had to be an integer (it is a C enum), so I chose an 8-bit
    int for it.

    However, a better solution is actually rather simple: Make a
    derived type and declare parameters with the values of the enum.
    Duh. I think I will implement that next.

    In Fortran 2023, you will be able to declare "real enums" that are typed
    for this. See https://j3-fortran.org/doc/year/21/21-110r1.txt

    And I am glad that this is in F2023. It's a well thought-out feature
    (and not too hard to implement :-) which is very useful.

    --- SoupGate-Win32 v1.05
    * Origin: fsxNet Usenet Gateway (21:1/5)