Merge pull request #1129 from MadhurKumar004/fix-1055

Fix: logspace generate wrong result for integer arguments

Author: loiseaujc

doc/specs/stdlib_math.md

@@ -234,7 +234,7 @@ For function calls where the `base` is specified, the `type` and `kind` of the r
 | "             " | "       " | `Integer` | `complex(KIND)` |
 | `Integer` | "        " | `real(KIND)` | `real(KIND)` |
 | "              " | "        " | `complex(KIND)` | `complex(KIND)` |
-| "              " | "        " | `Integer` | `Integer` |
+| "              " | "        " | `Integer` | `real(dp)` |
 
 #### Examples
 

src/math/stdlib_math.fypp

@@ -301,8 +301,8 @@ module stdlib_math
       integer, intent(in) :: end
       integer, intent(in) :: n
       integer, intent(in) :: base
-      ! integer endpoints + integer base = integer result
-      integer :: res(max(n, 0))
+      ! integer endpoints + integer base = real(dp) result
+      real(dp) :: res(max(n, 0))
     end function ${RName}$
 
 

src/math/stdlib_math_logspace.fypp

@@ -69,24 +69,30 @@ contains
   #:for k1 in REAL_KINDS
     #:set RName = rname("logspace", 1, "integer(int32)", "int32", "n_r" + str(k1) + "base")
     module procedure ${RName}$
-      integer :: exponents(max(n, 0))
-      exponents = linspace(start, end, n)
+      real(${k1}$) :: exponents(max(n, 0))
+      exponents = linspace(real(start, kind=${k1}$), real(end, kind=${k1}$), n)
       res = base ** exponents
+      if (n > 1) res(1) = base ** start
+      if (n > 0) res(n) = base ** end
     end procedure
 
     #:set RName = rname("logspace", 1, "integer(int32)", "int32", "n_c" + str(k1) + "base")
     module procedure ${RName}$
-      integer :: exponents(max(n, 0))
-      exponents = linspace(start, end, n)
+      complex(${k1}$) :: exponents(max(n, 0))
+      exponents = cmplx(linspace(real(start, kind=${k1}$), real(end, kind=${k1}$), n), 0.0_${k1}$)
       res = base ** exponents
+      if (n > 1) res(1) = base ** start
+      if (n > 0) res(n) = base ** end
     end procedure
   #:endfor
 
     #:set RName = rname("logspace", 1, "integer(int32)", "int32", "n_ibase")
     module procedure ${RName}$
-      integer :: exponents(max(n, 0))
-      exponents = linspace(start, end, n)
+      real(dp) :: exponents(max(n, 0))
+      exponents = linspace(real(start, kind=dp), real(end, kind=dp), n)
       res = base ** exponents
+      if (n > 1) res(1) = base ** real(start, kind=dp)
+      if (n > 0) res(n) = base ** real(end, kind=dp)
     end procedure
 
 

test/math/test_logspace.f90

@@ -31,7 +31,8 @@ subroutine collect_logspace(testsuite)
             new_unittest("logspace_default", test_logspace_default), &
             new_unittest("logspace_base_2", test_logspace_base_2), &
             new_unittest("logspace_base_2_cmplx_start", test_logspace_base_2_cmplx_start), &
-            new_unittest("logspace_base_i_int_start", test_logspace_base_i_int_start) &
+            new_unittest("logspace_base_i_int_start", test_logspace_base_i_int_start), &
+            new_unittest("logspace_int_overflow", test_logspace_int_overflow) &
             ]
 
     end subroutine collect_logspace
@@ -242,6 +243,43 @@ subroutine test_logspace_base_i_int_start(error)
 
     end subroutine
 
+    ! Regression test for issue #1055: integer start/end caused overflow in
+    ! intermediate integer arithmetic, producing garbage results.
+    subroutine test_logspace_int_overflow(error)
+        !> Error handling
+        type(error_type), allocatable, intent(out) :: error
+
+        integer, parameter :: n = 3
+        integer, parameter :: start = 10
+        integer, parameter :: end = 23
+        real(dp) :: expected_proportion
+        integer :: i
+
+        real(dp), allocatable :: x(:)
+
+        x = logspace(start, end, n)
+
+        expected_proportion = 10.0_dp ** ( real(end - start, dp) / real(n - 1, dp) )
+
+        call check(error, size(x), n, "Array not allocated to appropriate size")
+        if (allocated(error)) return
+        call check(error, x(1), 10.0_dp ** start, &
+            & "Initial value not equal to 10^start (possible integer overflow)", &
+            & thr=abs(10.0_dp ** start) * TOLERANCE_DP)
+        if (allocated(error)) return
+        call check(error, x(n), 10.0_dp ** end, &
+            & "Final value not equal to 10^end (possible integer overflow)", &
+            & thr=abs(10.0_dp ** end) * TOLERANCE_DP)
+        if (allocated(error)) return
+
+        do i = 1, n-1
+            call check(error, x(i + 1) / x(i), expected_proportion, &
+                & thr=abs(expected_proportion) * TOLERANCE_DP)
+            if (allocated(error)) return
+        end do
+
+    end subroutine
+
 
 end module