Open ivan-pi opened 2 years ago
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.
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.
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.
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.
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.
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 https://github.com/j3-fortran/fortran_proposals/issues/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
@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.
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:
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.
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.
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 parameter
s, 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
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)
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
! ...
A simple module parameterization facility would avoid the need for file inclusion and preprocessing tricks.
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).
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.
@ivan-pi writes Nov. 12, 2021 12:52 PM EST:
.. @FortranFan .. Related to the current issue, I can't see how would it allow practitioners to select the right kind:
The proposal I mentioned above tries to make it possible for practitioners and implementations to offer a "grouping" of related named constants (say math), hence enum
's, with consistent floating-point representations.
I envision a program to define a single enum of suitable floating-point representation for a given set of constants, say for math:
integer, parameter :: HP = selected_real_kind( p=xx ) !<-- or ieee_selected_real_kind(..), etc.
enum :: math(real(hp))
enumerator :: PI = 3.14159265358979323846264338327950288.._hp
enumerator :: e = ..
..
end enum
And for the code to consume the constants as scoped enumerations, say enum_nameXenumerator where X can be %
, ::
, etc.
The key is to have convenient grouping of constants that are all tied together with consistent floating-point representation.
A simple module parameterization facility would avoid the need for file inclusion and preprocessing tricks.
I personally think the boat for module parameterization has long sailed, it is needless complication considering what has been introduced into the language with Fortran 2003 thru' 2018 with PDTs and GENERIC
statements. A better approach now will be to simply build on this in a way that will greatly simplify the consumption of generic types on the client side. Say as follows
generic, kind :: RK => real32, real64, real128 !<-- define a GENERIC for KINDs a la GENERIC interfaces currently
type, public :: solver(k, n)
integer, kind :: k = <RK> !<-- use some symbols to designate a GENERIC set, shown here with angle brackets
integer, len :: n = 0
real(k) :: t
real(k), dimension(n) :: x
end type solver
contains
subroutine solve(me,t0,x0,dt)
type(solver(k=<RK>,n=*)), intent(inout) :: me
real(<RK>),intent(in) :: t0
real(<RK>),dimension(me%n),intent(in) :: x0
real(<RK>),intent(in) :: dt
..
end subroutine
The processor should then be able to do the semantics by substitution and build up the generic interface solve
.
The advantage with this can be that the code on the client side to use such as a PDT
with a generic interface for procedures will be the same as that currently. It's only on the "library" side that code duplication is reduced greatly via a GENERIC mechanism.
That only works for integer kind values. It doesn't work for more general types, such as a data structure module parameterized over arbitrary types.
Actually having modules like:
use constants_sp, only: pi
use constants_dp, only: pi
use constants_qp, only: pi
would work also I think. But it doesn't feel "right" to me either.
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
This has also been suggested and discussed at length in one of the stdlib issues: https://github.com/fortran-lang/stdlib/issues/99
Current status
In most Fortran codes today, users rely upon statements like:
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 ofM_E
.M_LOG10E
: The logarithm to base 10 ofM_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
, andmath.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.
which accept an integer kind argument. This is still kind of verbose, but users can always use a parameter:
or an associate construct:
if they want a shorter variable name.