Open milancurcic opened 3 years ago
Thank you, Steve. While it may be trivial for you, I suspect it's not as easy for the average Fortran programmer, and especially not for a novice. We're targeting newcomers to the language as well, to make working with Fortran as easy as possible for them. It took you (an expert) 15 minutes, but it could take somebody else 30 minutes or even an hour. If stdlib helps even only 100 people not write their own log2
, that's already a huge win in my view. And some of these people could be introducing bugs to their code, or not be aware of the subtleties like the x->1 special case that you mention.
No matter how simple the implementation is, it's useful to have an off-the-shelf solution that you can just import and be done with it.
Would you like to contribute your implementation to stdlib?
I agree we should mine iso_fortran_env
for the available kinds, that's for another issue.
Steve, thanks, that's music to my ears. And if you agree, we can indicate in the commit message that it is co-authored by you, which GitHub will pick up and attribute proper credit. I will soon open an issue dedicated to mining iso_fortran_env
for supported kinds.
The "trivial" implementation also does not give equivalent results to log2
from libm for numbers other than those special cases. Consider
module m_libm
interface
function libm_log2(x) bind(C, name="log2")
use, intrinsic :: iso_c_binding
real(c_double) :: libm_log2
real(c_double), intent(in) :: x
end function libm_log2
function libm_log(x) bind(C, name="log")
use, intrinsic :: iso_c_binding
real(c_double) :: libm_log
real(c_double), intent(in), value :: x
end function libm_log
end interface
end module m_libm
module m_log2
contains
real(8) function my_log2(x)
real(8), intent(in) :: x
real(8), parameter :: invln2 = 1 / log(2.e0_8)
my_log2 = exponent(x) + log(fraction(x)) * invln2
end function my_log2
end module m_log2
program test
use m_libm
use m_log2
implicit none
integer, parameter :: n = 10000
real(8), parameter :: xmax = 1000.0_8
integer :: i
real(8) :: x,y1,y2,d
do
call random_number(x)
x = x * xmax
#ifdef LOG2
y1 = libm_log2(x)
y2 = my_log2(x)
#elif LOG
y1 = libm_log(x)
y2 = log(x)
#endif
d = y1 - y2
if (d /= 0.0_8) then
write(*,*) x, y1, y2
end if
end do
end program test
Very quickly produces differences:
$ gfortran -DLOG2 log.F90 -lm
$ ./a.out| head | sort -n
1.1585799099427252 0.21235755365654199 0.21235755365654196
38.211646926593843 5.2559405342531784 5.2559405342531793
81.461109688336862 6.3480395621989354 6.3480395621989363
207.84346787839991 7.6993535973111049 7.6993535973111058
224.22528693079769 7.8088051765375956 7.8088051765375965
281.33160252573532 8.1361278122684251 8.1361278122684233
595.08034940289781 9.2169406680421080 9.2169406680421098
609.06285212072703 9.2504473042184063 9.2504473042184046
634.89193775024694 9.3103672475828585 9.3103672475828603
758.29366600803792 9.5666128619771502 9.5666128619771520
The same test with -DLOG
seems to produce identical results, i.e. producing no output if I leave it running :
$ gfortran -DLOG log.F90 -lm
$ timeout 10m ./a.out
$
edit: add missing value
attribute
Thanks for the detailed reply (including pointing out some errors on my part). Hopefully this points to the utility of a high quality implementation of log2
for inclusion in stdlib. In terms of performance (and without doing any benchmarking...) I'll just note that this latest version does potentially introduce at least an additional branch over a direct implementation.
The previous comments show that the implementation is not trivial to get a high quality solution, but focusing on implementing log2
seems strange to me in the spirit of a consistent experience. Why would we not focus on a log
function where you can pass the base in as an optional argument as many other programming languages do?
By this logic there is no point in having any features in fortran that exist in any C standard library.
I don't know the full history of including a variety of Bessel functions, etc, but it does seem odd to have such an otherwise full set of special functions that isn't a full superset of what you get with math.h
@trippalamb there isn't much logic in specifying log2
in isolation, that is the just the one I stumbled across not existing which eventually lead to this issue.
edit2: to be clear, I don't actually feel very strongly one way or another about this feature/item in case the above seems too argumentative
Looks like we lost some replies to this issue. I will attempt to summarize some of the points made so they are not lost:
log(+-0) = -inf
with divide by zero exception, log(1)=+0
, log(x<0) = NaN + invalid exception
, log(+inf) = +inf
x
is close to 1math.h
is essentially everywhere, so just use iso_c_bind
log2(x) = log(x) / log(2.)
Implement base 2 logarithm function
log2
. Proposed by @brandongc in https://github.com/j3-fortran/fortran_proposals/issues/222. This looks like it belongs instdlib_math
.Prior art
numpy.log2
log2
Base.log2