Implement constants module

[certik]

At the very least, it would contain the following constants:

• pi
• i_
• e_

An example implementation: https://github.com/certik/fortran-utils/blob/b43bd24cd421509a5bc6d3b9c3eeae8ce856ed88/src/constants.f90

Note about naming: The convention that we discussed in fortran-utils 10 years ago was that single letter constants contain underscore so that they do not clash with user variables ("e" and "i" are frequently used as loop variables). But we can definitely revisit this and choose a different convention.

[marshallward]

Here are the libm constants (and names) for reference:

https://www.gnu.org/software/libc/manual/html_node/Mathematical-Constants.html#Mathematical-Constants

[certik]

Julia seems to use Base.im, Base.MathConstants.pi, Base.MathConstants.ℯ, Base.MathConstants.eulergamma, ... (https://docs.julialang.org/en/v1/base/numbers/#Base.im).

[certik]

Matlab uses pi. Other constants: https://www.mathworks.com/help/matlab/constants-and-test-matrices.html.

[scivision]

to avoid namespace clashes, the constants could be in a derived type consisting mostly of parameter like

mc%pi  
mc%tau

[certik]

We had that discussion at #49 about using derived types to workaround the insufficiency of Fortran's namespaces. I personally think it's not a good idea, you can read the arguments at #49.

But if the majority in the community wants to use a derived type for this, we can, especially while it is still in experimental. We can always revisit later.

The proper solution, in my opinion, is j3-fortran/fortran_proposals#1 and j3-fortran/fortran_proposals#86. If both are implemented, then one could do something like this:

use, namespace :: stdlib, only: constants
...
constants%pi

Where constants is the module nested in stdlib. This would be equivalent to Python's from stdlib import constants.

[epagone]

I agree with @certik that the derived type solution is suboptimal, but I think it's the best that we can do with Fortran 2018. Once the standard is fixed (I wouldn't hold my breath meanwhile), we can revert to a better solution, similarly to what has been done with procedure optional argument default values and the optval function.

[nncarlson]

I think using a derived type for the constants is perfectly legitimate and if done "properly" acts exactly like a namespace and the user needn't even be aware that there is a derived type involved other than the "%" character which becomes part of the constant "name". For example,

module math_constants
  private
  type :: private_type
    real :: pi, e
  end type
  type(private_type), parameter, public :: mc = private_type(3.14, 2.71)
end module

use math_constants
print *, mc%pi, mc%e
end

Interestingly, for the NAG compiler this has 0 overhead compared to using individual parameters; the mc structure doesn't even appear in the C code for the print statement -- it's been completely unwound to the actual constants.

[scivision]

That's an interesting technique. Will have to check with other compiler disassemble to see if zero overhead holds there too.

[certik]

I propose we use this:

module math_constants
  private
  real, parameter, public :: pi = 3.14, e_ = 2.71
end module

use math_constants, only: pi, e_
print *, pi, e_
end

Both the module and user code is shorter and simpler.

Using the derived type approach, how do I import just pi? I don't think you can. So all your user code will have to always type mc%pi.

Here is an example from one of my codes

VeeG = 4*pi*neG / G2

Much more readable than:

VeeG = 4*mc%pi*neG / G2

[zjibben]

My personal preference is slightly toward @certik's bare parameters in a module, primarily for cleanliness in expressions. And very subjectively I realize, bundling constants into a single data object feels wrong. On the off chance there's a collision, one can always rename on use: use math_constants, only: mc_pi => pi. To me that's the superior compromise until true namespaces are available.

[nncarlson]

I want to clarify my earlier comment that I was not advocating for using the derived type approach in this case. I'm fairly ambivalent about it, and could go either way. I really just wanted to push back on the idea that this is misuse of derived types in general. I don't think it is. It provides a "namespace" like experience in lieu of having namespaces, and seems to do so without any overhead.

[nshaffer]

What precision should these constants be? Do we define, e.g., pi_sp, pi_dp, pi_qp and then expect the user to choose the one they want (giving them the chance to rename it)?

As for namespacing, I don't see the benefit here. Namespacing the constants seems like over-engineering, plus you lose the ability to just use pi or e_ as @certik mentioned.

[epagone]

What precision should these constants be? Do we define, e.g., pi_sp, pi_dp, pi_qp and then expect the user to choose the one they want (giving them the chance to rename it)?

How about we follow a pragmatic approach and implement only the highest possible precision type ("qp"?) with no suffix and then let Fortran implicitly take care of the conversions? Is it too much a "quick and dirty" solution?

[certik]

Regarding derived type versus just variables in a module: let's do both, so that we can all get what we want and move on. So let's do this:

module math_constants
private  
real, parameter, public :: pi = 3.14, e_ = 2.71  
 
type :: private_type  
    real :: pi = pi  
    real :: e = e_  
end type  
 
type(private_type), public, parameter :: mc = private_type()  
end module

Then this can be used both as:

use math_constants, only: pi, e_
print *, pi, e_
end

and as:

use math_constants, only: mc
print *, mc%pi, mc%e
end

I just tested it and it works.

There are essentially two camps here --- one side thinks it's an over engineering and an imperfect workaround for a fundamental deficiency of Fortran namespaces; the other side thinks it's worth using derived types as namespaces. The above approach gets both sides what they want, without forcing the other side to use an approach that feels wrong.

So let's try that and move on. We can do the same approach in #49.

---

Regarding the constant's precision, that's a very good point. If we set them to the highest precision available in the compiler, as in this code:

program test_pi
use stdlib_experimental_kinds, only: sp, dp, qp
implicit none
real(dp), parameter :: pi_dp    = 3.1415926535897932384626433832795_dp
real(qp), parameter :: pi_qp    = 3.1415926535897932384626433832795028841971_qp


real(dp) :: a
real(qp) :: b

a = pi_dp
print *, a
a = pi_qp
print *, a

b = pi_dp
print *, b
b = pi_qp
print *, b

end program

which prints:

   3.1415926535897931     
   3.1415926535897931     
   3.14159265358979311599796346854418516      
   3.14159265358979323846264338327950280      

Then it looks like things behave correctly. But unfortunately gfortran gives a warning:

test_pi.f90:13:4:

 a = pi_qp
    1
Warning: Change of value in conversion from ‘REAL(16)’ to ‘REAL(8)’ at (1) [-Wconversion]

Actually I do not even know how to get rid of this warning, as this also doesn't work:

test_pi.f90:13:9:

 a = real(pi_qp, dp)
         1
Warning: Change of value in conversion from ‘REAL(16)’ to ‘REAL(8)’ at (1) [-Wconversion]

So we need to figure this out.

[epagone]

Actually I do not even know how to get rid of this warning, as this also doesn't work:

Unfortunately, I had a quick look at the gcc manual and it doesn't seems possible :(

[certik]

A pragmatic solution: given that double precision is no doubt the most widely used precision, then we can have:

real(sp), parameter :: pi_sp = ...
real(dp), parameter :: pi_dp = ...
real(qp), parameter :: pi_qp = ...

real(dp), parameter :: pi = pi_dp

So one can use pi_sp, pi_dp, pi_qp directly, and one can also just use pi in the most common case of double precision.

[epagone]

Would it be complicated to just ignore the warnings with some scripting? I don't like the idea to bend the implementation to the quirks of individual compilers (especially for warnings).

[certik]

I checked Intel Fortran with -warn all and it does not warn about my code above. The PGI Fortran with -Minform=warn also does not warn (it does not support qp, so I changed qp -> dp and dp -> sp in the above example). Let's check other compilers too.

I agree if this is only gfortran's quirk, then we can simply disable this warning for GFortran as a workaround and report the bug, and simply declare constants as the highest precision supported in the compiler (typically qp).

[epagone]

Unfortunately, I can only test gfortran (...and lfortran, but I think there is someone more knowledgeable than me about it in this discussion :P), so I cannot tell what happens with NAG, IBM, PGI, etc...

[nncarlson]

NAG does not warn about these type conversions.

[marshallward]

I don't think gcc is wrong here. -Wconversion warns of implicit conversions, and this is an implicit conversion from real(16) to real(8).

However, it also seems to be balking on explicit conversions, e.g.

a = real(pi_qp, dp)

which does feel wrong to me.

Unfortunately I'm not sure the best way do this other than something like -Wall -Wno-conversion, but maybe there's a more explicit way to define conversions.

[epagone]

Unfortunately I'm not sure the best way do this other than something like -Wall -Wno-conversion, but maybe there's a more explicit way to define conversions.

Thanks @marshallward for finding the right gfortran flag (i.e. -Wno-conversion) and apologies for having missed that in the gcc manual (maybe I've looked in the wrong place).

I don't think gcc is wrong here. -Wconversion warns of implicit conversions, and this is an implicit conversion from real(16) to real(8).

I agree with you, but in this specific case, we don't need to worry about the warning because the implicit conversion it's intended and it is not caused by a distraction.

[nshaffer]

I wonder if there's a performance penalty for just providing the highest precision available on a compiler.

I tested the following:

program test
    use, intrinsic :: iso_fortran_env, only: dp => real32, qp => real128
    implicit none

    real(dp), parameter :: pi_dp = 3.1415926535897932384626433832795_dp
    real(qp), parameter :: pi_qp = 3.1415926535897932384626433832795028841971_qp

    real(dp) :: x

    x = 0.0_dp
    x = x + ___________
end program test

where I replaced the blank with

• pi_dp
• pi_qp
• real(pi_qp, dp)

After compiling and disassembly, e.g.,

$ gfortran test.f90
$ objdump --disassemble a.out > dp.s

Explicitly casting pi_qp down to dp produced identical assembly code as just using dp. Credit goes to gfortan's compile-time constant expression reduction, I expect.

Using pi_qp without the cast produced more complicated assembly code. The result is only a few instructions longer, but it introduces several more calls into what I think is the gfortran runtime library. That may or may not matter. Just some data.

Edit: If I compile with -O3, all three variants produce the same assembly code.

[nncarlson]

I wonder if there's a performance penalty for just providing the highest precision available on a compiler.

I like this idea, but there is (at least) one case where the user would have to be aware of the actual kind and explicitly down-cast if necessary: [ 0.0_dp, pi, 2*pi] would return an error if pi was quad kind

[jvdp1]

I like this idea, but there is (at least) one case where the user would have to be aware of the actual kind and explicitly down-cast if necessary: [ 0.0_dp, pi, 2*pi] would return an error if pi was quad kind

The user would be aware of the kind if stdlib provides pi_qp instead of pi.
Provding only pi_qp would be a solution half-way between providing only pi (or pi_) and providing all kinds (pi_sp,pi_dp, pi_qp)

[marshallward]

I think -O3 is just pre-computing all of the constant expressions because it knows that x is zero, producing identical assembly, including the real(pi_qp, dp) conversion.

If I take out the x = 0.0_dp line, so that x has random garbage, then it's a little more clear what is happening, including the conversion steps by GCC (via soft-fp).

Replacing pi_dp with pi_qp just adds this conversion code:

< 	movss	.LC1(%rip), %xmm0
< 	movq    %rbp, %rdi
< 	addss	12(%rsp), %xmm0
< 	movq    %r13, 24(%rsp)
< 	movl    $12, 32(%rsp)
< 	movss   %xmm0, 12(%rsp)
---
> 	movss	12(%rsp), %xmm0
> 	call	__extendsftf2@PLT
> 	movdqa	.LC1(%rip), %xmm1
> 	call	__addtf3@PLT
> 	call	__trunctfsf2@PLT
< 	movq    %rbp, %rdi
< 	movq    %r13, 24(%rsp)
< 	movss   %xmm0, 12(%rsp)
< 	movl    $12, 32(%rsp)

I don't thin this has anything to do with providing the quad precision numbers, it's just doing the necessary qp to dp conversion in x = x + pi_qp.

[urbanjost]

Curious what other compilers do when this is both compiled and run program testit use,intrinsic :: iso_fortran_env, only : int8, int16, int32, int64 implicit none write(*,*)int(huge(1_int64),kind=int32) end program testit Somtimes I want a warning even from int() and real() as gfortran does. I would prefer that the specific intrinsics for changing type not produce a warning but that I did get a warning over overflow/underflow, but the warning would make int() and real() IMPURE, I suppose. I had similar issues with defining a table of my own favorite constants in the past and so gave myself several choices (I really could not make up my mind and was new to modules and such anyway. So I did this module M_constants use, intrinsic :: iso_fortran_env, only : real32, real64, real128 implicit none real(kind=real128),parameter :: pi128=3.141592653589793238462643383279502884197169399375105820974944592307_real128 real(kind=real32),parameter :: pi32=real(pi128,kind=real32) real(kind=real64),parameter :: pi64=real(pi128,kind=real64) end module M_constants module M_constants_32 use M_constants, only : pi=> pi32 end module M_constants_32 module M_constants_64 use M_constants, only : pi=> pi64 end module M_constants_64 module M_constants_128 use M_constants, only : pi=> pi128 end module M_constants_128 so without TOO much maintenance I could use it several different ways: program testit call zero() call one() call two() call three() contains subroutine zero() use :: M_constants write(*,*)pi32/2 end subroutine zero subroutine one() use :: M_constants, only : pi=>pi32 write(*,*)pi/2 end subroutine one subroutine two() use :: M_constants_32, only : pi implicit none write(*,*)pi/2 end subroutine two subroutine three() use :: M_constants_128, only : pi implicit none write(*,*)pi/2 end subroutine three end program testit I suppose it is only my personal story as the discussion shows opinions vary, but when I look thru code wh ere I used it, I almost always ended up using case1, and otherwise case0 out of the four cases above. It ends up I like being able to tell the type by the name myself; although I initially thought I wanted to use 'PI' everywhere.

…

On January 13, 2020 at 1:13 PM Marshall Ward ***@***.***> wrote: I don't think gcc is wrong here. -Wconversion warns of implicit conversions, and this is an implicit conversion from real(16) to real(8). However, it also seems to be balking on explicit conversions, e.g. a = real(pi_qp, dp) which does feel wrong to me. Unfortunately I'm not sure the best way do this other than something like -Wall -Wno-conversion, but maybe there's a more explicit way to define conversions. — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub #99?email_source=notifications&email_token=AHDWN3NQOK4CVIVID772ICLQ5SVN3A5CNFSM4KENU6B2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEIZXDSQ#issuecomment-573796810 , or unsubscribe https://github.com/notifications/unsubscribe-auth/AHDWN3NAZLYVMNJUU4JFFCDQ5SVN3ANCNFSM4KENU6BQ .

[urbanjost]

propogating to the larger type in expressions is not the only issue, but using the constants in expressions or values passed as arguments or building arrays might be the bigger issue. I have sometimes wished there was a way to declare an argument to propogate to the type required by the called routine, like a built-in call to INT or REAL to the type required, sort of like K&R C propogated all the float values up before there were 10 000, types. Maybe even propogate down if the value fit. Metamorphic variables and CLASS(*) and generic routines give the functionality I wanted without the possible costs of a gather/scatter to the other type and so on nowadays, I suppose.

…

On January 13, 2020 at 1:32 PM nshaffer ***@***.***> wrote: I wonder if there's a performance penalty for just providing the highest precision available on a compiler. I tested the following: program test use, intrinsic :: iso_fortran_env, only: dp => real32, qp => real128 implicit none real(dp), parameter :: pi_dp = 3.1415926535897932384626433832795_dp real(qp), parameter :: pi_qp = 3.1415926535897932384626433832795028841971_qp real(dp) :: x x = 0.0_dp x = x + ___________ end program test where I replaced the blank with * pi_dp * pi_qp * real(pi_qp, dp) After compiling and disassembly, e.g., $ gfortran test.f90 $ objdump --disassemble a.out > dp.s Explicitly casting pi_qp down to dp produced identical assembly code as just using dp. Credit goes to gfortan's compile-time constant expression reduction, I expect. Using pi_qp without the cast produced more complicated assembly code. The result is only a few instructions longer, but it introduces several more calls into what I think is the gfortran runtime library. That may or may not matter. Just some data. — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub #99?email_source=notifications&email_token=AHDWN3KFGXNNY5X2MSP6573Q5SXUTA5CNFSM4KENU6B2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEIZZIGQ#issuecomment-573805594 , or unsubscribe https://github.com/notifications/unsubscribe-auth/AHDWN3N6UCX4QBYEBSJ4RWTQ5SXUTANCNFSM4KENU6BQ .

[certik]

If there is a performance penalty at runtime in Release mode (with all optimizations on) when using pi_qp in double precision expression, then that's not acceptable I think.

[jannisteunissen]

I am curious if any further progress has been made in the meantime. Personally, I think having suffixes that indicate the precision are worth it. The names become a bit longer, but when reading code it is certainly helpful to know what the precision of a number is. It would also prevent most name clashes.

[jvdp1]

I am curious if any further progress has been made in the meantime.

See this related comment

Personally, I think having suffixes that indicate the precision are worth it. The names become a bit longer, but when reading code it is certainly helpful to know what the precision of a number is. It would also prevent most name clashes.

I don't know how @MilanSkocic aims to to do that.