Base 2 logarithm

[milancurcic]

Implement base 2 logarithm function log2. Proposed by @brandongc in https://github.com/j3-fortran/fortran_proposals/issues/222. This looks like it belongs in stdlib_math.

Prior art

• numpy.log2
• MATLAB log2
• Julia Base.log2

[milancurcic]

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.

[milancurcic]

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.

[brandongc]

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

[brandongc]

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.

[trippalamb]

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?

[brandongc]

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

[brandongc]

Looks like we lost some replies to this issue. I will attempt to summarize some of the points made so they are not lost:

1. Naive comparison of results of floating point results should be avoided (Units in the Last Place ULP is preferred)
2. Special cases need care (log(+-0) = -inf with divide by zero exception, log(1)=+0, log(x<0) = NaN + invalid exception, log(+inf) = +inf
3. Need to be careful when x is close to 1
4. An argument against adding this kind of thing is that math.h is essentially everywhere, so just use iso_c_bind
5. Another argument against is that you already have a log function with different base so use log2(x) = log(x) / log(2.)