logspace generates wrong results for integer arguments

[cneiss]

Description

The logspace routine seems to generate wrong results when called with integer arguments. Here is an example code (essentially taken from the stdlib specs page https://stdlib.fortran-lang.org/page/specs/stdlib_math.html#logspace-create-a-logarithmically-spaced-rank-one-array):

program example_logspace_int
  use stdlib_math, only: logspace
  use stdlib_kinds, only: dp
  implicit none

  integer, parameter :: start = 10
  integer, parameter :: end = 23
  integer, parameter :: n = 3

  real(dp) :: r(n) ! Integer values raised to real powers results in real values

  r = logspace(start, end, n)
  print *,  r
end program example_logspace_int

This produces 1410065408.0000000 1874919424.0000000 -159383552.00000000 as output, which is obviously wrong. I've tested with gfortran-14 and Intel oneapi ifx 2025; does not seem to depend on the compiler.

Expected Behaviour

The result should be 10000000000.000000 31622776601683792. 1.0000000000000001E+023

Version of stdlib

fffe0d7058180d22faba9a2d32377cd88f6aefb7

Platform and Architecture

Linux (OpenSUSE Leap 15.6)

Additional Information

Note: I can reproduce the numbers if I replace r = logspace(start, end, n) with

integer :: dum
dum = linspace(start, end, n)
r = 10**dum

where dum ist of type integer; with type real(dp) the correct results are obtained. Hence, the error seems to be related to a type mismatch in the logspace implementation.

[loiseaujc]

I've been able to reproduce the same problem with

program example_logspace_int
   use stdlib_math, only: logspace
   use stdlib_kinds, only: dp
   implicit none

   integer, parameter :: start = 9
   integer, parameter :: end = 10
   integer, parameter :: n = 2

   real(dp) :: r(n) ! Integer values raised to real powers results in real values

   r = logspace(start, end, n)
   print *, r
end program example_logspace_int

The issue appears as soon as 10 is in the range. It works as expected however if start and end are declared as real numbers.

Edit : Here is a bit more digging.

program example_logspace_int
   use stdlib_kinds, only: dp
   use stdlib_math, only: logspace, linspace
   implicit none
   integer, parameter :: start = 10
   integer, parameter :: end = 10
   integer, parameter :: n = 1
   integer :: exponents(n)
   real(dp) :: r(n)
   exponents = linspace(start, end, n)
   r = logspace(start, end, n)
   print *, 10**exponents
   print *, r
end program example_logspace_int

Both code snippet print the incorrect values 1410065408 and 1410065408.000000 while if I compute directly 10**10 (everything is integer), I get the following error

Error: Result of exponentiation at (1) exceeds the range of INTEGER(4)

I think this has to do with the default base begin declared as a regular integer. For some reason however, the compiler does not raise an error if the exponent is a size-1 array. gfortran bug?

[cneiss]

The problem in your second example is the line exponents = linspace(start, end, n) where exponents is of type integer, but linspace returns real (as it should). This is exactly what seems to happen in the logspace routine as well. The line print *, 10**exponents generates an integer overflow.

It's not a compiler bug IMO.

[loiseaujc]

Yep, that's correct. The possible bug I had in mind was the fact that gfortran correctly identifies that 10**10 exceeds the range, but that it raises no error for 10**exponents where exponents is a size-1 array filled with 10. The two computations should be identical.