median: some cleaning

Author: jvdp1

src/stdlib_stats_median.fypp

@@ -21,23 +21,27 @@ contains
         logical, intent(in), optional :: mask
         real(${o1}$) :: res
 
-        integer :: c, n
+        integer(kind = int64) :: c, n
         ${t1}$, allocatable :: x_tmp(:)
 
         if (.not.optval(mask, .true.) .or. size(x) == 0) then
           res = ieee_value(1._${o1}$, ieee_quiet_nan)
           return
         end if
+        
+        n = size(x, kind=int64)
+        c = floor( (n + 1) / 2._${o1}$, kind=int64 )
 
-        x_tmp = reshape(x, [size(x)])
+        x_tmp = reshape(x, [n])
 
         call sort(x_tmp)
 
-        n = size(x_tmp)
-        c = floor( (n + 1) / 2._${o1}$ )
-
-        if (mod(n, 2) == 0) then
+        if (mod(n, 2_int64) == 0) then
+          #:if t1[0] == 'r'
             res = sum(x_tmp(c:c+1)) / 2._${o1}$
+          #:else
+            res = sum( real(x_tmp(c:c+1), kind=${o1}$) ) / 2._${o1}$
+          #:endif
         else
             res = x_tmp(c)
         end if
@@ -56,9 +60,11 @@ contains
         real(${o1}$) :: res${reduced_shape('x', rank, 'dim')}$
 
         integer :: c, n
+        #:if rank > 1
         #:for fj in range(1, rank+1)
         integer :: j${"_" * fj}$
         #:endfor
+        #:endif
         ${t1}$, allocatable :: x_tmp(:)
 
         if (.not.optval(mask, .true.) .or. size(x) == 0) then
@@ -84,7 +90,12 @@ contains
                   call sort(x_tmp)
 
                   if (mod(n, 2) == 0) then
-                      res${reduce_subvector(rank, fi)}$ = sum(x_tmp(c:c+1)) / 2._${o1}$
+                      res${reduce_subvector(rank, fi)}$ = &
+                      #:if t1[0] == 'r'
+                          sum(x_tmp(c:c+1)) / 2._${o1}$
+                      #:else
+                          sum(real(x_tmp(c:c+1), kind=${o1}$) ) / 2._${o1}$
+                      #:endif
                   else
                       res${reduce_subvector(rank, fi)}$ = x_tmp(c)
                   end if
@@ -109,21 +120,25 @@ contains
         logical, intent(in) :: mask${ranksuffix(rank)}$
         real(${o1}$) :: res
 
-        integer :: c, n
-        ${t1}$, allocatable :: x_tmp(:)
+        integer(kind = int64) :: c, n
+        ${t1}$, allocatable   :: x_tmp(:)
 
         x_tmp = pack(x, mask)
 
         call sort(x_tmp)
 
-        n = size(x_tmp)
-        c = floor( (n + 1) / 2._${o1}$ )
+        n = size(x_tmp, kind=int64 )
+        c = floor( (n + 1) / 2._${o1}$, kind=int64 )
 
         if (n == 0) then
             res = ieee_value(1._${o1}$, ieee_quiet_nan)
-        else if (mod(n, 2) == 0) then
+        else if (mod(n, 2_int64) == 0) then
+          #:if t1[0] == 'r'
             res = sum(x_tmp(c:c+1)) / 2._${o1}$
-        else if (mod(n, 2) == 1) then
+          #:else
+            res = sum(real(x_tmp(c:c+1), kind=${o1}$)) / 2._${o1}$
+          #:endif
+        else if (mod(n, 2_int64) == 1) then
             res = x_tmp(c)
         end if
 
@@ -140,10 +155,12 @@ contains
         logical, intent(in) :: mask${ranksuffix(rank)}$
         real(${o1}$) :: res${reduced_shape('x', rank, 'dim')}$
 
-        integer :: c, n
+        integer(kind = int64) :: c, n
+        #:if rank > 1
         #:for fj in range(1, rank+1)
         integer :: j${"_" * fj}$
         #:endfor
+        #:endif
         ${t1}$, allocatable :: x_tmp(:)
 
         select case(dim)
@@ -159,15 +176,20 @@ contains
                                 mask${select_subvector(rank, fi)}$)
                   call sort(x_tmp)
 
-                  n = size(x_tmp)
-                  c = floor( (n + 1) / 2._${o1}$ )
+                  n = size(x_tmp, kind=int64)
+                  c = floor( (n + 1) / 2._${o1}$, kind=int64 )
 
                   if (n == 0) then
                       res${reduce_subvector(rank, fi)}$ = &
                           ieee_value(1._${o1}$, ieee_quiet_nan)
-                  else if (mod(n, 2) == 0) then
-                      res${reduce_subvector(rank, fi)}$ = sum(x_tmp(c:c+1)) / 2._${o1}$
-                  else if (mod(n, 2) == 1) then
+                  else if (mod(n, 2_int64) == 0) then
+                      res${reduce_subvector(rank, fi)}$ = &
+                      #:if t1[0] == 'r'
+                          sum(x_tmp(c:c+1)) / 2._${o1}$
+                      #:else
+                          sum(real(x_tmp(c:c+1), kind=${o1}$)) / 2._${o1}$
+                      #:endif
+                  else if (mod(n, 2_int64) == 1) then
                       res${reduce_subvector(rank, fi)}$ = x_tmp(c)
                   end if