Bug in gsl 2.8 formula

[AlexBuccheri]

Description of the Problem

Homebrew GSL 2.8 arm64 macOS bottle can segfault when called from a GNU Fortran program compiled with -ffpe-trap=zero.

In GSL 2.8, complex/math.c:gsl_complex_sqrt contains:

if (x >= y)
  {
    double t = y / x;
    w = sqrt (x) * sqrt (0.5 * (1.0 + sqrt (1.0 + t * t)));
  }
else
  {
    double t = x / y;
    w = sqrt (y) * sqrt (0.5 * (t + sqrt (1.0 + t * t)));
  }

For input 4 + 0i, x = 4 and y = 0, the source path should take the x >= y branch and compute y / x, i.e. 0 / 4.
However Homebrew bottle's optimized arm64 assembly speculatively emits both divisions before selecting the result (assembly produced by Codex 5.5):

_gsl_complex_sqrt:
...
fdiv d4, d2, d3   ; x / y, i.e. 4 / 0
...
fdiv d3, d3, d2   ; y / x

When the GNU Fortran caller enables FP traps with -ffpe-trap=zero, the speculative 4 / 0 traps even though that source branch should not be taken.

Minimal Reproduction

I don't appear to be allowed to upload .F90 files so copy-pasting:

CMakeLists.txt
gsl_sqrt_probe.c
gsl_sqrt_driver.F90

cmake_minimum_required(VERSION 3.20)

project(gsl_trap LANGUAGES C Fortran)

find_package(GSL REQUIRED)

add_library(gsl_call STATIC gsl_sqrt_probe.c)
target_link_libraries(gsl_call PUBLIC GSL::gsl)

add_executable(gsl_trap gsl_sqrt_driver.F90)
target_link_libraries(gsl_trap PRIVATE gsl_call)

#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>

void gsl_sqrt_probe(double re_in, double im_in, double *re_out, double *im_out)
{
  gsl_complex z = gsl_complex_rect(re_in, im_in);
  gsl_complex root = gsl_complex_sqrt(z);

  *re_out = GSL_REAL(root);
  *im_out = GSL_IMAG(root);
}

program gsl_sqrt_driver
  use, intrinsic :: iso_c_binding, only: c_double
  implicit none

  interface
    subroutine gsl_sqrt_probe(re_in, im_in, re_out, im_out) bind(c, name='gsl_sqrt_probe')
      import :: c_double
      implicit none
      real(c_double), value :: re_in
      real(c_double), value :: im_in
      real(c_double), intent(out) :: re_out
      real(c_double), intent(out) :: im_out
    end subroutine gsl_sqrt_probe
  end interface

  real(c_double) :: re_out
  real(c_double) :: im_out

  write (*, '(a)') 'Calling gsl_complex_sqrt(4 + 0i)'
  call gsl_sqrt_probe(4.0_c_double, 0.0_c_double, re_out, im_out)
  write (*, '(f12.6,1x,sp,f12.6,a)') re_out, im_out, 'i'

end program gsl_sqrt_driver

README

Assumes Homebrew GSL is in the $PATH

Configure with trapping:

cmake -B build --fresh \
-DCMAKE_BUILD_TYPE=Debug \
-DCMAKE_Fortran_COMPILER=/opt/homebrew/bin/gfortran-15 \
-DCMAKE_Fortran_FLAGS="-ffpe-trap=zero"

Configure without trapping:

cmake -B build --fresh \
-DCMAKE_BUILD_TYPE=Debug \
-DCMAKE_Fortran_COMPILER=/opt/homebrew/bin/gfortran-15

Build:

cmake --build build

Run:

./build/gsl_trap

Solution

I think it's an Apple Clang problem. Compiling GSL with Apple Clang and the CFLAG -ffp-exception-behavior=strict resolves the problem. Compiling with GCC is also an option (I did not test this).

Note, it would be great to be able report these Issues on the core page without having to tick off the mandatory requirements. I appreciate they're deliberate barriers to entry, but I'm also not going to remove installations from requested directories in order to prove there's a problem.

[carlocab]

Is this Homebrew-specific, or can you reproduce this with GSL built without Homebrew?