Merge pull request #605 from zoziha/add_diff

[stdlib_math] Add function `diff`

Author: gareth-nx

CHANGELOG.md

@@ -36,6 +36,8 @@ Features available from the latest git source
     [#488](https://github.com/fortran-lang/stdlib/pull/488)
   - new procedures `arg`, `argd` and `argpi`
     [#498](https://github.com/fortran-lang/stdlib/pull/498)
+  - new procedure `diff`
+    [#605](https://github.com/fortran-lang/stdlib/pull/605)
 
 Changes to existing modules
 

doc/specs/stdlib_math.md

@@ -649,3 +649,88 @@ program demo_math_all_close
     
 end program demo_math_all_close
 ```
+
+### `diff`
+
+#### Description
+
+Computes differences between adjacent elements of an array.
+
+#### Syntax
+
+For a rank-1 array
+```fortran
+y = [[stdlib_math(module):diff(interface)]](x [, n, prepend, append])
+```
+and for a rank-2 array
+```fortran
+y = [[stdlib_math(module):diff(interface)]](x [, n, dim, prepend, append])
+```
+
+#### Status
+
+Experimental.
+
+#### Class
+
+Pure function.
+
+#### Arguments
+
+`x`: The array to take a difference of.
+Shall be a `real/integer` and `rank-1/rank-2` array.
+This argument is `intent(in)`.
+
+`n`: How many times to iteratively calculate the difference.
+Shall be an `integer` scalar.
+This argument is `intent(in)` and `optional`, and has value of `1` by default.
+
+`dim`: The dimension of the input array along which to calculate the difference.
+Its value must be between `1` and `rank(x)`.
+Shall be an `integer` scalar.
+This argument is `intent(in)` and `optional` and has a value of `1` by default.
+
+`prepend`, `append`: Arrays to prepend or append to a along axis prior to performing the difference.
+The dimension and shape must match a except along axis.
+Shall be a `real/integer` and `rank-1/rank-2` array.
+This argument is `intent(in)` and `optional`, which is no value by default.
+
+Note: 
+- The `x`, `prepend` and `append` arguments must have the same `type`, `kind` and `rank`.
+- If the value of `n` is less than or equal to `0` (which is not recommended), the return value of `diff` is `x`.
+- If the value of `dim` is not equal to `1` or `2` (which is not recommended),
+`1` will be used by the internal process of `diff`.
+
+
+#### Result value
+
+Returns the finite difference of the input array.
+Shall be a `real/integer` and `rank-1/rank-2` array.
+When both `prepend` and `append` are not present, the result `y` has one fewer element than `x` alongside the dimension `dim`.
+
+#### Example
+
+```fortran
+program demo_diff
+
+    use stdlib_math, only: diff
+    implicit none
+    
+    integer :: i(7) = [1, 1, 2, 3, 5, 8, 13]
+    real    :: x(6) = [0, 5, 15, 30, 50, 75]
+    integer :: A(3, 3) = reshape([1, 7, 17, 3, 11, 19, 5, 13, 23], [3, 3])
+    integer :: Y(3, 2)
+    
+    print *, diff(i)        ! [0, 1, 1, 2, 3, 5]
+    print *, diff(x, 2)     ! [5.0, 5.0, 5.0, 5.0]
+    
+    Y = diff(A, n=1, dim=2)
+    print *, Y(1, :)        ! [2, 2]
+    print *, Y(2, :)        ! [4, 2]
+    print *, Y(3, :)        ! [2, 4]
+
+    print *, diff(i, prepend=[0]) ! [1, 0, 1, 1, 2, 3, 5]
+    print *, diff(i, append=[21]) ! [0, 1, 1, 2, 3, 5, 8]
+    
+end program demo_diff
+```

src/CMakeLists.txt

@@ -51,6 +51,7 @@ set(fppFiles
     stdlib_math_arange.fypp
     stdlib_math_is_close.fypp
     stdlib_math_all_close.fypp
+    stdlib_math_diff.fypp
     stdlib_string_type.fypp
     stdlib_string_type_constructor.fypp
     stdlib_strings_to_string.fypp

src/Makefile.manual

@@ -48,6 +48,7 @@ SRCFYPP = \
         stdlib_math_logspace.fypp \
         stdlib_math_is_close.fypp \
         stdlib_math_all_close.fypp \
+        stdlib_math_diff.fypp \
         stdlib_string_type.fypp \
         stdlib_string_type_constructor.fypp \
         stdlib_strings.fypp \
@@ -238,6 +239,7 @@ stdlib_math_is_close.o: \
 stdlib_math_all_close.o: \
     stdlib_math.o \
     stdlib_math_is_close.o
+stdlib_math_diff.o: stdlib_math.o
 stdlib_stringlist_type.o: stdlib_string_type.o \
     stdlib_math.o \
     stdlib_optval.o

src/stdlib_math.fypp

@@ -14,7 +14,7 @@ module stdlib_math
     public :: EULERS_NUMBER_QP
 #:endif
     public :: DEFAULT_LINSPACE_LENGTH, DEFAULT_LOGSPACE_BASE, DEFAULT_LOGSPACE_LENGTH
-    public :: arange, arg, argd, argpi, is_close, all_close
+    public :: arange, arg, argd, argpi, is_close, all_close, diff
 
     integer, parameter :: DEFAULT_LINSPACE_LENGTH = 100
     integer, parameter :: DEFAULT_LOGSPACE_LENGTH = 50
@@ -359,6 +359,28 @@ module stdlib_math
         #:endfor
         #:endfor
     end interface all_close
+    
+    !> Version: experimental
+    !>
+    !> Computes differences between adjacent elements of an array.
+    !> ([Specification](../page/specs/stdlib_math.html#diff))
+    interface diff
+        #:set RI_KINDS_TYPES = REAL_KINDS_TYPES + INT_KINDS_TYPES
+        #:for k1, t1 in RI_KINDS_TYPES
+        pure module function diff_1_${k1}$(x, n, prepend, append) result(y)
+            ${t1}$, intent(in) :: x(:)
+            integer, intent(in), optional :: n
+            ${t1}$, intent(in), optional :: prepend(:), append(:)
+            ${t1}$, allocatable :: y(:)
+        end function diff_1_${k1}$
+        pure module function diff_2_${k1}$(X, n, dim, prepend, append) result(y)
+            ${t1}$, intent(in) :: x(:, :)
+            integer, intent(in), optional :: n, dim
+            ${t1}$, intent(in), optional :: prepend(:, :), append(:, :)
+            ${t1}$, allocatable :: y(:, :)
+        end function diff_2_${k1}$
+        #:endfor
+    end interface diff
 
 contains
 

src/stdlib_math_diff.fypp

@@ -0,0 +1,139 @@
+!> Inspired by original code (MIT license) written in 2016 by Keurfon Luu ([email protected])
+!> https://github.com/keurfonluu/Forlab
+
+#:include "common.fypp"
+#:set RI_KINDS_TYPES = REAL_KINDS_TYPES + INT_KINDS_TYPES
+submodule (stdlib_math) stdlib_math_diff
+
+    implicit none
+
+contains
+
+    !> `diff` computes differences of adjacent elements of an array.
+    
+    #:for k1, t1 in RI_KINDS_TYPES
+    pure module function diff_1_${k1}$(x, n, prepend, append) result(y)
+        ${t1}$, intent(in) :: x(:)
+        integer, intent(in), optional :: n
+        ${t1}$, intent(in), optional :: prepend(:), append(:)
+        ${t1}$, allocatable :: y(:)
+        integer :: size_prepend, size_append, size_x, size_work
+        integer :: n_, i
+
+        n_ = optval(n, 1)
+        if (n_ <= 0) then
+            y = x
+            return
+        end if
+        
+        size_prepend = 0
+        size_append = 0
+        if (present(prepend)) size_prepend = size(prepend) 
+        if (present(append)) size_append = size(append)
+        size_x = size(x)
+        size_work = size_x + size_prepend + size_append
+        
+        if (size_work <= n_) then
+            allocate(y(0))
+            return
+        end if
+
+        !> Use a quick exit for the common case, to avoid memory allocation.
+        if (size_prepend == 0 .and. size_append == 0 .and. n_ == 1) then
+            y = x(2:) - x(1:size_x-1)
+            return
+        end if
+
+        block
+        ${t1}$ :: work(size_work)
+        if (size_prepend > 0) work(:size_prepend) = prepend
+        work(size_prepend+1:size_prepend+size_x) = x
+        if (size_append > 0) work(size_prepend+size_x+1:) = append
+        
+        do i = 1, n_
+            work(1:size_work-i) = work(2:size_work-i+1) - work(1:size_work-i)
+        end do
+
+        y = work(1:size_work-n_)
+        end block
+
+    end function diff_1_${k1}$
+
+    pure module function diff_2_${k1}$(x, n, dim, prepend, append) result(y)
+        ${t1}$, intent(in) :: x(:, :)
+        integer, intent(in), optional :: n, dim
+        ${t1}$, intent(in), optional :: prepend(:, :), append(:, :)
+        ${t1}$, allocatable :: y(:, :)
+        integer :: size_prepend, size_append, size_x, size_work
+        integer :: n_, dim_, i
+
+        n_ = optval(n, 1)
+        if (n_ <= 0) then
+            y = x
+            return
+        end if
+        
+        size_prepend = 0
+        size_append = 0
+        if (present(dim)) then
+            if (dim == 1 .or. dim == 2) then
+                dim_ = dim
+            else
+                dim_ = 1
+            end if
+        else
+            dim_ = 1
+        end if
+        
+        if (present(prepend)) size_prepend = size(prepend, dim_)
+        if (present(append)) size_append = size(append, dim_)
+        size_x = size(x, dim_)
+        size_work = size_x + size_prepend + size_append
+        
+        if (size_work <= n_) then
+            allocate(y(0, 0))
+            return
+        end if
+        
+        !> Use a quick exit for the common case, to avoid memory allocation.
+        if (size_prepend == 0 .and. size_append == 0 .and. n_ == 1) then
+            if (dim_ == 1) then
+                y = x(2:, :) - x(1:size_x-1, :)
+            elseif (dim_ == 2) then
+                y = x(:, 2:) - x(:, 1:size_x-1)
+            end if
+            return
+        end if
+        
+        if (dim_ == 1) then
+            block
+            ${t1}$ :: work(size_work, size(x, 2))
+            if (size_prepend > 0) work(1:size_prepend, :) = prepend
+            work(size_prepend+1:size_x+size_prepend, :) = x
+            if (size_append > 0) work(size_x+size_prepend+1:, :) = append
+            do i = 1, n_
+                work(1:size_work-i, :) = work(2:size_work-i+1, :) - work(1:size_work-i, :)
+            end do
+
+            y = work(1:size_work-n_, :)
+            end block
+            
+        elseif (dim_ == 2) then
+            block
+            ${t1}$ :: work(size(x, 1), size_work)
+            if (size_prepend > 0) work(:, 1:size_prepend) = prepend
+            work(:, size_prepend+1:size_x+size_prepend) = x
+            if (size_append > 0) work(:, size_x+size_prepend+1:) = append
+            do i = 1, n_
+                work(:, 1:size_work-i) = work(:, 2:size_work-i+1) - work(:, 1:size_work-i)
+            end do
+            
+            y = work(:, 1:size_work-n_)
+            end block
+            
+        end if
+
+    end function diff_2_${k1}$
+    #:endfor
+
+end submodule stdlib_math_diff