Intrinsics for mathematical constants

[ivan-pi]

This issue originates from one of the questions at the HPFPC (High-Performance computing with Fortran Promoting Consortium) symposium (https://site.hpfpc.org/home/events/parallel_fortran_sympo5)

In the requests & questions list we can read

Q27: The standard should include some physical constants such as pi. It is more helpful if they can used without using any modules, similarly to intrinsic procedures.

A27: This would not be a big addition to the language, but it has not been requested.

This has also been suggested and discussed at length in one of the stdlib issues: fortran-lang/stdlib#99

Current status

In most Fortran codes today, users rely upon statements like:

real(wp), parameter :: pi = 4*atan(1.0_wp) ! or
real(wp), parameter :: pi = 2*acos(0.0_wp)
real(wp), parameter :: e = exp(1.0_wp)

that are collected in a constants module, leading to duplication of effort and potential errors.

Personally, I often re-declare pi locally in each module (to lazy to use a module) and find myself wondering which trigonometric function am I supposed to use in my head or drawing a unit circle on a sheet of paper.

Other programming languages

Other languages are also known to provide such constants, e.g. in C the header file math.h provides the following:

• M_E: The base of natural logarithms.
• M_LOG2E: The logarithm to base 2 of M_E.
• M_LOG10E: The logarithm to base 10 of M_E.
• M_LN2: The natural logarithm of 2.
• M_LN10: The natural logarithm of 10.
• M_PI: Pi, the ratio of a circle’s circumference to its diameter.
• M_PI_2: Pi divided by two.
• M_PI_4: Pi divided by four.
• M_1_PI: The reciprocal of pi (1/pi).
• M_2_PI: Two times the reciprocal of pi.
• M_2_SQRTPI: Two times the reciprocal of the square root of pi.
• M_SQRT2: The square root of two.
• M_SQRT1_2: The reciprocal of the square root of two (also the square root of 1/2).

Julia also provides the constants for im (imaginary unit), pi, ℯ (the constant ℯ), Catalan's constant, Euler's constant, and the golden ratio, among others.

Python provides math.pi, math.e, and math.tau (2 times pi) to available precision (presumably C double in most implementations).

Solution

The solution could be provided in terms of intrinsic functions, e.g.

pure function math_pi(kind) real(pi)
   integer, kind :: kind  ! not valid Fortran
   real(kind=kind) :: pi

which accept an integer kind argument. This is still kind of verbose, but users can always use a parameter:

real(wp), parameter :: pi = math_pi(wp)

or an associate construct:

associate(pi => math_pi(wp), e => math_e(wp))
...
end associate

if they want a shorter variable name.

[ivan-pi]

Note in C, the constants are of type double, implemented as macro substitutions. A GNU C library extension also provides constants of type long double, by appending the a lowercase character l to their names.

Since the Fortran standard does not specify a preprocessor, intrinsic functions with a kind specifier seem like the cleanest solution. An optimizing compiler can easily unwind the function call and substitute the constant with the desired precision.

[klausler]

A more idiomatic interface would be to pass an unused real argument whose kind (not value) determines the kind of the result.. This interface could be implemented in current Fortran without adding intrinsics.

[certik]

As a user I think my preferred way is still to just import them from some module, such as math_constants or just constants:

use math_constants, only: pi

But I don't now how to do it to work with the different kinds properly. As a practical approach, just providing them in double precision would be enough for all my use cases. But I can see how this is not ideal, we really want this to work for all kinds well.

@klausler suggestion would work today:

real(wp), parameter :: x = 1
real(wp), parameter :: pi = math_pi(x)

But it's quite verbose.

[jacobwilliams]

This needs to be part of the generics feature in Fortran 202y. Not only do we need to be able to define routines and classes that can use any real (or integer, etc.) kind, but that needs to also trickle down to parameters. The kind in the parameter has to be able to inherit that somehow, without having to define multiple ones for each kind (which is basically impossible...which is the problem we have now), or give up on parameters altogether. Not having this is one of the reasons parameterized derived types are mostly useless.

[ivan-pi]

Since you mention generics, in C++ I sometimes see programmers use the following pattern:

template <typename T> T pi_const() {
  return static_cast<T>(3.14159265358979323846);
}

Usage looks like:

auto pi = pi_const<double>()

I guess the answer from @klausler would be one of the arguments against new intrinsic functions. It is similar to some of the intrinsic functions, like newline or kind (roughly). Instead of declaring an unused variable you can always just use a literal constant to achieve resolution of the correct interface, e.g:

module math
  implicit none
  private
  public :: math_pi
  interface math_pi
    module procedure math_pi_sp
    module procedure math_pi_dp
  end interface
contains
pure function math_pi_sp(a)
  integer, parameter :: sp = kind(1.0e0)
  real(sp), intent(in) :: a
  pi = 3.14159265358979323846264338327950288_sp
end function
pure function math_pi_dp(a) result(pi)
  integer, parameter :: dp = kind(1.0d0)
  real(dp), intent(in) :: a
  pi = 3.14159265358979323846264338327950288_dp
end function
end module

program circle
use math
implicit none
integer, parameter :: dp = kind(1.0d0)

real(dp), parameter :: pi = math_pi(1.0_dp)
real(dp) :: radius, area

write(*,*) "Enter circle radius"
read(*,*) radius
area = pi * radius**2
write(*,*) "Circle area = ", area
end program

This feels okay-ish. Perhaps this is the best available approach to use in stdlib.

[ivan-pi]

The kind in the parameter has to be able to inherit that somehow, without having to define multiple ones for each kind (which is basically impossible...which is the problem we have now), or give up on parameters altogether. Not having this is one of the reasons parameterized derived types are mostly useless.

@jacobwilliams, could you elaborate please? I fail to see how adding math constants relates to parameterized derived types. The first part of your answer relates to #91 if I understood correctly.

Apart from the boilerplate in the implementation, I don't think the Fortran approach:

real(dp), parameter :: pi = math_pi(1.0_dp)

is much worse than the C++ approach used in practice (C++ of course being a language famous for generic capabilities):

const double pi = pi_const<double>()

The C solution using macros is convenient, but limited to double precision. The same can be done easily with Fortran or Python:

from math import pi
print("pi = {}".format(pi))

include <math.h>
printf("pi = %f\n",M_PI);

use math, only: pi
write(*,*) "pi = ", pi
end

[FortranFan]

@ivan-pi writes Nov. 11, 2021, 6:11 AM EST:

.. The solution could be provided in terms of intrinsic functions ..

Please see the link in the original post of #46, an alternate solution I had proposed a couple of years ago toward such needs is a proper enumeration type facility in Fortran.

My proposal had "built-in" Generics i.e., the enumerator constants in an enumeration type can be of any intrinsic type as part of the definition of the enumeration type itself.

Here are some extracted bits from that proposal:

d) where the underlying type of enumerators in an integer with a
specified kind e.g.,

 enum :: MY_HASH_CODES(integer(kind=INT64))
    enumerator :: FOO = Z"d76aa478"
    enumerator :: BAR = Z"e8c7b756"
    ..
 end enum

e) an enumeration type where the underlying type of elements in the
enumerated list is that of the LOGICAL intrinsic type. An example
with the defaul logical kind is shown below though it is expected
other supported logical kinds can be used also.

 enum :: VALVE_STATE(logical)
    enumerator :: OPEN = .true.
    enumerator :: CLOSED = .false.
 end enum

f) an enumeration type where the underlying type of elements in the
enumerator list is that of the CHARACTER intrinsic type. An
example with the defaul characted kind is shown below though
it is expected other supported character kinds can be used also.

 enum :: LNG_CONSTITUENT(character(len=*))
    enumerator :: N2 = "Nitrogen"
    enumerator :: C1 = "Methane"
    enumerator :: C2 = "Ethane"
    enumerator :: C3 = "Propane"
    ..
 end enum

g) an enumeration type where the underlying type of elements in the
enumerator list is that of the REAL intrinsic type. An
example with a user-defined kind is shown below and it is
expected other supported real kinds can also be used similarly.

 integer, parameter :: R8 = selected_real_kind( p=12 )
 enum :: PHYS_CHEM_CONSTANTS(real(kind=R8)) ! SI units
    enumerator :: MU = 1.66053906660E-27_r8 ! Atomic mass constant
    enumerator :: NA = 6.02214076E23_r8     ! Avogadro number
    enumerator :: K = 1.380649E-23_r8       ! Boltzmann constant
    enumerator :: R = 8.314462618_r8        ! Molar gas constant
    ..
 end enum

..

a) "cast" the enumerator value to an object of an intrinsic type using
the intrinsic conversion functions e.g.,
..

 ii) REAL intrinsic
     real, parameter :: Rgas=real(PHYS_CHEM_CONSTANTS%R, kind(Rgas))

..

9. that it will be possible for an enumeration type to be the selector in
   an ASSOCIATE construct e.g., with reference to 1g above,

   type(PHYS_CHEM_CONSTANTS) :: k
   ..
   associate ( CONST => PHYS_CHEM_CONSTANTS )
   ..
   k = CONST%k
   ..
   end associate

10. that it is possible for an enumerator to be the selector in an
    ASSOCIATE construct e.g., with reference to 1e above,

    associate ( foo => MY_HASH_CODES%FOO )
    ..
    end associate

    associate ( R => real(PHYS_CHEM_CONSTANTS%R, kind=..) )
    ..
    Density = P/R/T
    ..
    end associate

Once such a comprehensive facility were to be added to the language, my vision was "intrinsic enumerations" a la "intrinsic named constants" in some "intrinsic modules" (along the lines of ISO_FORTRAN_ENV, IEEE_ARITHMETIC, etc.). Intrinsic enumerations for, say, MATH can then offer PI, E, and so forth for those users who would like to use them.

The proposal was seen as too comprehensive for the "minor revision" that is Fortran 202X and it failed to get any support.

But I maintain a proper enumeration type is the way to facilitate the use of many literal constants, be they of any intrinsic type, by the practitioners in their Fortran code.

[ivan-pi]

Thanks @FortranFan for the link. I like your proposed syntax for enumerators and think it would offer safer Fortran usage outside of the traditional scientific computing domain.

Related to the current issue, I can't see how would it allow practitioners to select the right kind:

enum :: math_sp(real(sp))
  enumerator :: pi = 4*atan(1.0_sp)
  enumerator :: e = exp(1.0_sp)
end enum

enum :: math_dp(real(dp))
  enumerator :: pi = 4*atan(1.0_dp)
  enumerator :: e = exp(1.0_dp)
end enum

In the stdlib thread on mathematical constants, some practitioners were against using a derived type as a "pseudo" enumerator, for the sole reason that it felt wrong to them to have to retrieve constants using the % operator. I think it's a case of de gustibus non disputandum est.

If I return to the original post, two questions that should be answered first are:

• Do we support adding mathematical constants (or functions/enums/derived types to retrieve them) and if yes, can they be in global namespace or should they be part of an intrinsic module?
• How do we allow practitioners to retrieve constants with the right kind, without resorting to expressions like 4*atan(1.0_wp) or exp(1.0_wp)?

If the answer from the majority is they are fine with a module, this can be done in stdlib with the approach @klausler suggested. If it is preferable to have pi built-in, we need to discuss further. By built-in I mean something along the following lines would be valid Fortran program:

real :: pi = 3.14 ! shadows intrinsic pi, be it function or literal constant

block
intrinsic :: pi 
print *, pi  ! or pi(1.0d0), prints 3.14... 
end block

I am afraid there would be opposition to having many short names like pi, e, and others. Other names like the Archimedes's constant or Euler's number are cumbersome. So prepending the names with M_ like C does, or math like Python, seems to be desirable.

A third alternative, what Julia does, is that the core mathematical constants are handled as special values of type "Irrational". These values get converted to floating point numbers without intermediate rounding when used in mathematical expressions.

This is an elegant solution which I guess could be done in Fortran too, but is cumbersome for implementation, and would need to be taken care of on the compiler side to guarantee the expressions get replaced with floating point numbers for minimum overhead.

[certik]

The most natural for me would be to redo how Fortran handles mixed precision. But that's a big task and probably not possible to do at this point.

But as a user I would like to just create a math_constants module with:

real(decimal), parameter :: pi = 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679 

which would be "exact" (in this case 100 digits) and then I would like to use it as:

real(sp) :: r4
real(dp) :: r8
r4 = sqrt(pi)*pi
r8 = sqrt(pi)*pi

and it would just work. It would do the operations with the single or double precision, say in here:

r4 = r4 * pi
r8 = r8 * pi

In the first line it would use a single precision multiplication, in the second line it would use a double precision multiplication. So there would be no runtime overhead.

In fact, as a user I would like exactly this behavior for any numbers such as:

r4 = r4 * 8.15
r8 = r8 * 8.15

As I said, this would require to rethink how mixed precision is done. As discussed elsewhere, in Fortran the precision is inferred from the right hand side, never from the left hand side. So 8.15 or "pi" would have to be a new type, "decimal" (with the exact decimal digits as you write them) and then it would be cast down to whatever accuracy is needed. However, I think this approach hits some obstacles how Fortran traditionally does things.

[jacobwilliams]

Consider this PDT example. Say you write to write a solver that works for any real kind. You get pretty far and then you realize that you have a bunch of parameters that need to be that same kind. It's basically impossible to do that. Functions are not really what you want (and you can't have functions for every real number). You just want to use parameters, which is the most natural way to do this with Fortran. PDT is basically a half-baked generic feature that never went anywhere and was never fixed. A real generic feature would account for this somehow.

I don't know what the syntax should look like, but I don't think we should have to use functions or derived types just to specify numeric parameters.

module test

use iso_fortran_env 

implicit none 

private

type,public  :: solver(rk, n)
    integer,kind :: rk = real64
    integer,len  :: n  = 0
    real(rk) :: t
    real(rk),dimension(n) :: x
   contains
   procedure :: solve => rk8
end type solver

contains 
    subroutine rk8(me,t0,x0,dt)
    implicit none
    class(solver(rk=real64,n=*)),intent(inout) :: me
    real(me%rk),intent(in) :: t0 
    real(me%rk),dimension(me%n),intent(in) :: x0 
    real(me%rk),intent(in) :: dt 

    real(me%rk),parameter :: coefficient_1 = 1.0_me%rk / 23.0_me%rk   ! not possible

   ...
    end subroutine rk8

end module test

[ivan-pi]

Here's an example of a PDT for performing LU factorization:

module lu_pdt

  implicit none 
  private

  public :: lu_workspace, factorize

  type :: lu_workspace(wp,n)
    integer, kind :: wp
    integer, len :: n
    real(wp) :: a(n,n)
    real(wp) :: b(n)
    integer :: ipiv(n)
    logical :: factorized = .false.
  end type

  interface factorize
    module procedure factorize_sp
    module procedure factorize_dp
  end interface

  integer, parameter :: sp = kind(1.0e0)
  integer, parameter :: dp = kind(1.0d0)

contains

  subroutine factorize_sp(this,info)
    type(lu_workspace(sp,*)), intent(inout) :: this
    integer, intent(out), optional :: info  

    integer :: info_
    external :: sgetrf

    if (.not. this%factorized) then
      call sgetrf(this%n,this%n,this%a,this%n,this%ipiv,info_)
      if (info_ == 0) then
        this%factorized = .true.
      end if
      if (present(info)) info = info_
    else
      return
    end if

  end subroutine

  subroutine factorize_dp(this,info)
    type(lu_workspace(dp,*)), intent(inout) :: this
    integer, intent(out), optional :: info  

    integer :: info_
    external :: dgetrf

    if (.not. this%factorized) then
      call dgetrf(this%n,this%n,this%a,this%n,this%ipiv,info_)
      if (info_ == 0) then
        this%factorized = .true.
      end if
      if (present(info)) info = info_
    else
      return
    end if

  end subroutine

end module

program main

  use lu_pdt
  implicit none

  integer, parameter :: sp = kind(1.0e0)
  integer, parameter :: dp = kind(1.0d0)

  type(lu_workspace(dp,:)), allocatable :: work
  integer :: info

  allocate(lu_workspace(dp,n=3) :: work)

  work%a = reshape(&
    [real(dp) :: 4, 2, -1, -3, 1, 2, 1, 3, -5], &
    [3,3])

  work%b = [real(dp) :: -10, 0, 17]

  call factorize(work,info)
  print *, "info = ", info

end program

I couldn't test it fully because I am missing LAPACK on this PC. But gfortran is able to compile it. If I make factorize a type-bound method then it breaks. You are right it's a half-baked feature, but I think we haven't even begun to explore what it offers in terms of simplified API'S.

With some include statements (read poor man's templates), you could make it support all allowable real kinds and even mixed-precision (apart from the calls to the LAPACK routines). Granted, it´s verbose and the methods wouldn't be type-bound (no chance for polymorphism), but it can still work to some degree.

If you look at some of the Julia packages, e.g. Pardiso.jl, they also use the style of passing a solver object to a function:

ps = PardisoSolver()

A = sparse(rand(10, 10))
B = rand(10, 2)
X = zeros(10, 2)
solve!(ps, X, A, B)

[ivan-pi]

In the context of your Runge Kutta PDT, the solve function would like like this:

type, public  :: solver(rk, n)
    integer, kind :: rk = real64
    integer, len  :: n  = 0
    real(rk) :: t
    real(rk), dimension(n) :: x
end type solver

generic :: solve => rk4_sp, rk4_dp

contains
    subroutine rk4_sp(me,t0,x0,dt)
    use prec, only: wp => sp
    include 'rk4.inc'
    end subroutine rk4_sp
    subroutine rk4_dp(me,t0,x0,dt)
    use prec, only: wp => dp
    include 'rk4.inc'
    end subroutine rk4_dp

and the include file would be

! rk4.inc
    type(solver(rk=wp,n=*)), intent(inout) :: me
    real(wp),intent(in) :: t0 
    real(wp),dimension(me%n),intent(in) :: x0 
    real(wp),intent(in) :: dt 

    real(wp),parameter :: coefficient_1 = 1.0_wp / 23.0_wp
    ! ...

[klausler]

A simple module parameterization facility would avoid the need for file inclusion and preprocessing tricks.

[jacobwilliams]

The problem with include, aside from the interface duplication you have to do, is that there's no way to make it work with any real kinds the compiler supports, without either knowing what they are in advance, or going outside of Fortran (some sort of preprocessor/introspection tricks).

[klausler]

The problem with include, aside from the interface duplication you have to do, is that there's no way to make it work with any real kinds the compiler supports, without either knowing what they are in advance, or going outside of Fortran (some sort of preprocessor/introspection tricks).

That's right, and that's why a module instantiation facility would need to be able to instantiate a module over a set of values. Easy to define and easy to implement.