diff --git a/doc/specs/fftpack.md b/doc/specs/fftpack.md
index be8b73f..3fe1062 100644
--- a/doc/specs/fftpack.md
+++ b/doc/specs/fftpack.md
@@ -809,13 +809,16 @@ end program demo_dzfftb
 
 ## DCT type-1 (DCT-1)
 
-### Initialize DCT-1: `dcosti`
+### Initialize DCT-1: `dcosti` or `dct_t1i`
 
 #### Description
 
 Initializes the array `wsave` which is used in subroutine `dcost`.
 The prime factorization of `n` together with a tabulation of the trigonometric functions are computed and stored in `wsave`.
 
+The two procedures are completely equivalent and expect the same arguments.
+It is a matter of personal preference which one you choose to use. 
+
 #### Status
 
 Experimental
@@ -851,18 +854,21 @@ program demo_dcosti
 end program demo_dcosti
 ```
 
-### Compute DCT-1: `dcost`
+### Compute DCT-1: `dcost` or `dct_t1`
 
 #### Description
 
 Computes the DCT-1 of the input real data.
 The transform is defined below at output parameter `x`.
 
+The two procedures are completely equivalent and expect the same arguments.
+It is a matter of personal preference which one you choose to use. 
+
 For real input data `x` of length `n`, the DCT-1 of `x` is equivalent, up to a
 scaling factor, to the DFT of the even extension of `x` with length `2*(n-1)`,
 where the first and last entries of the original data are not repeated in the
-extension. For example, the DCT-1 of input data *abcde* (size \[5\]) is
-equivalent to the DFT of data *abcdedcb* (size \[2*4=8\]).
+extension. For example, the DCT-1 of input data *abcde* (size \(5\)) is
+equivalent to the DFT of data *abcdedcb* (size \(2*4=8\)).
 
 Also, `dcost` is the unnormalized inverse of itself. This means that a call of
 `dcost` followed by another call of `dcost` will multiply the input sequence `x`
@@ -932,7 +938,7 @@ end program demo_dcost
 
 ## DCT of types 2, 3 (DCT-2, 3), a.k.a "Quarter" cosine transforms
 
-### Initialize DCT-2, 3: `dcosqi`
+### Initialize DCT-2, 3: `dcosqi` or `dct_t23i`
 
 #### Description
 
@@ -941,6 +947,9 @@ The prime factorization of `n` together with
 a tabulation of the trigonometric functions are computed and
 stored in `wsave`.
 
+The two procedures are completely equivalent and expect the same arguments.
+It is a matter of personal preference which one you choose to use. 
+
 #### Status
 
 Experimental
@@ -978,13 +987,16 @@ program demo_dcosqi
 end program demo_dcosqi
 ```
 
-### Compute DCT-3: `dcosqf`
+### Compute DCT-3: `dcosqf` or `dct_t3`
 
 #### Description
 
 Computes the DCT-3 of the input real data.
 The transform is defined below at output parameter `x`.
 
+The two procedures are completely equivalent and expect the same arguments.
+It is a matter of personal preference which one you choose to use. 
+
 Also, `dcosqf` (DCT-3) is the unnormalized inverse of `dcosqb` (DCT-2), since a
 call of `dcosqf` followed by a call of `dcosqb` will multiply the input sequence
 `x` by `4*n`.
@@ -1049,13 +1061,16 @@ program demo_dcosqf
 end program demo_dcosqf
 ```
 
-### Compute DCT-2: `dcosqb`
+### Compute DCT-2: `dcosqb` or `dct_t2`
 
 #### Description
 
 Computes the DCT-2 of the input real data.
 The transform is defined below at output parameter `x`.
 
+The two procedures are completely equivalent and expect the same arguments.
+It is a matter of personal preference which one you choose to use. 
+
 For real input data `x` of length `n`, the DCT-2 of `x` is equivalent, up to a
 scaling factor, to the DFT of the even extension of `x` with length `4*n`,
 where all the even-frequency entries are zero.
@@ -1076,7 +1091,7 @@ Pure subroutine.
 
 #### Syntax
 
-`call [[fftpack(module):dcosqf(interface)]](n, x, wsave)`
+`call [[fftpack(module):dcosqb(interface)]](n, x, wsave)`
 
 #### Arguments
 
diff --git a/src/fftpack.f90 b/src/fftpack.f90
index d02479a..cea79f9 100644
--- a/src/fftpack.f90
+++ b/src/fftpack.f90
@@ -17,6 +17,8 @@ module fftpack
     public :: dcosqi, dcosqf, dcosqb
     public :: dcosti, dcost
     public :: dct, idct
+    public :: dct_t1i, dct_t1
+    public :: dct_t23i, dct_t2, dct_t3
 
     public :: rk
 
@@ -125,7 +127,7 @@ end subroutine dzfftb
         !> Version: experimental
         !>
         !> Initialize `dcosqf` and `dcosqb`.
-        !> ([Specification](../page/specs/fftpack.html#dcosqi))
+        !> ([Specification](../page/specs/fftpack.html#initialize-dct-2-3-dcosqi-or-dct_t23i))
         pure subroutine dcosqi(n, wsave)
             import rk
             integer, intent(in) :: n
@@ -135,7 +137,7 @@ end subroutine dcosqi
         !> Version: experimental
         !>
         !> Forward transform of quarter wave data.
-        !> ([Specification](../page/specs/fftpack.html#dcosqf))
+        !> ([Specification](../page/specs/fftpack.html#compute-dct-3-dcosqf-or-dct_t3))
         pure subroutine dcosqf(n, x, wsave)
             import rk
             integer, intent(in) :: n
@@ -146,7 +148,7 @@ end subroutine dcosqf
         !> Version: experimental
         !>
         !> Unnormalized inverse of `dcosqf`.
-        !> ([Specification](../page/specs/fftpack.html#dcosqb))
+        !> ([Specification](../page/specs/fftpack.html#compute-dct-2-dcosqb-or-dct_t2))
         pure subroutine dcosqb(n, x, wsave)
             import rk
             integer, intent(in) :: n
@@ -156,7 +158,8 @@ end subroutine dcosqb
 
         !> Version: experimental
         !>
-        !> Initialize `dcost`. ([Specification](../page/specs/fftpack.html#dcosti))
+        !> Initialize `dcost`.
+        !> ([Specification](../page/specs/fftpack.html#initialize-dct-1-dcosti-or-dct_t1i))
         pure subroutine dcosti(n, wsave)
             import rk
             integer, intent(in) :: n
@@ -166,7 +169,7 @@ end subroutine dcosti
         !> Version: experimental
         !>
         !> Discrete fourier cosine transform of an even sequence.
-        !> ([Specification](../page/specs/fftpack.html#dcost))
+        !> ([Specification](../page/specs/fftpack.html#compute-dct-1-dcost-or-dct_t1))
         pure subroutine dcost(n, x, wsave)
             import rk
             integer, intent(in) :: n
@@ -245,7 +248,7 @@ end function irfft_rk
     !> Version: experimental
     !>
     !> Dsicrete cosine transforms.
-    !> ([Specification](../page/specs/fftpack.html#dct))
+    !> ([Specification](../page/specs/fftpack.html#simplified-dct-of-types-1-2-3-dct))
     interface dct
         pure module function dct_rk(x, n, type) result(result)
             real(kind=rk), intent(in) :: x(:)
@@ -258,7 +261,7 @@ end function dct_rk
     !> Version: experimental
     !>
     !> Inverse discrete cosine transforms.
-    !> ([Specification](../page/specs/fftpack.html#idct))
+    !> ([Specification](../page/specs/fftpack.html#simplified-inverse-dct-of-types-1-2-3-idct))
     interface idct
     pure module function idct_rk(x, n, type) result(result)
         real(kind=rk), intent(in) :: x(:)
@@ -268,6 +271,46 @@ pure module function idct_rk(x, n, type) result(result)
         end function idct_rk
     end interface idct
 
+    !> Version: experimental
+    !>
+    !> Initialize DCT type-1
+    !> ([Specification](../page/specs/fftpack.html#initialize-dct-1-dcosti-or-dct_t1i))
+    interface dct_t1i
+        procedure :: dcosti
+    end interface dct_t1i
+
+    !> Version: experimental
+    !>
+    !> Perform DCT type-1
+    !> ([Specification](../page/specs/fftpack.html#compute-dct-1-dcost-or-dct_t1))
+    interface dct_t1
+        procedure :: dcost
+    end interface dct_t1
+
+    !> Version: experimental
+    !>
+    !> Initialize DCT types 2, 3
+    !> ([Specification](../page/specs/fftpack.html#initialize-dct-2-3-dcosqi-or-dct_t23i))
+    interface dct_t23i
+        procedure :: dcosqi
+    end interface dct_t23i
+
+    !> Version: experimental
+    !>
+    !> Perform DCT type-2
+    !> ([Specification](../page/specs/fftpack.html#compute-dct-2-dcosqb-or-dct_t2))
+    interface dct_t2
+        procedure :: dcosqb
+    end interface dct_t2
+
+    !> Version: experimental
+    !>
+    !> Perform DCT type-3
+    !> ([Specification](../page/specs/fftpack.html#compute-dct-3-dcosqf-or-dct_t3))
+    interface dct_t3
+        procedure :: dcosqf
+    end interface dct_t3
+
     !> Version: experimental
     !>
     !> Shifts zero-frequency component to center of spectrum.
diff --git a/test/test_fftpack_dct.f90 b/test/test_fftpack_dct.f90
index f680058..cb4cf49 100644
--- a/test/test_fftpack_dct.f90
+++ b/test/test_fftpack_dct.f90
@@ -1,6 +1,6 @@
 module test_fftpack_dct
 
-    use fftpack, only: rk, dcosti, dcost, dct, idct, dcosqi, dcosqf, dcosqb
+    use fftpack
     use testdrive, only: new_unittest, unittest_type, error_type, check
     implicit none
     private
@@ -24,27 +24,53 @@ end subroutine collect_dct
 
     subroutine test_classic_dct(error)
         type(error_type), allocatable, intent(out) :: error
-        real(kind=rk) :: w(3*4 + 15)
+        real(kind=rk) :: w(3*4 + 15), w2(3*4 + 15)
         real(kind=rk) :: x(4) = [1, 2, 3, 4]
+        real(kind=rk) :: x2(4)
         real(kind=rk) :: eps = 1.0e-10_rk
-
+        
+        x2 = x
         call dcosti(4, w)
         call dcost(4, x, w)
-        call check(error, all(x == [real(kind=rk) :: 15, -4, 0, -1.0000000000000009_rk]), "`dcosti` failed.")
+        call check(error, sum(abs(x - [real(kind=rk) :: 15, -4, 0, -1.0000000000000009_rk])) < eps, &
+                   "`dcost` failed.")
+        if (allocated(error)) return
+        
+        call dct_t1i(4, w2)
+        call dct_t1(4, x2, w2)
+        call check(error, maxval(abs(x2-x)) < eps, "dct_t1 failed")
         if (allocated(error)) return
+
         call dcost(4, x, w)
-        call check(error, all(x/(2*3) == [real(kind=rk) :: 1, 2, 3, 4]), "`dcost` failed.")
+        call check(error, sum(abs(x/(2*3) - [real(kind=rk) :: 1, 2, 3, 4])) < eps, &
+                   "2nd `dcost` failed.")
+        if (allocated(error)) return
+
+        call dct_t1(4, x2, w2)
+        call check(error, maxval(abs(x2-x)) < eps, "2nd dct_t1 failed") 
+        if (allocated(error)) return
 
         x = [1, 2, 3, 4]
+        x2 = x
         call dcosqi(4, w)
         call dcosqf(4, x, w)
         call check(error, sum(abs(x - [11.999626276085150_rk, -9.1029432177492193_rk, &
                                        2.6176618435106480_rk, -1.5143449018465791_rk])) < eps, &
                    "`dcosqf` failed.")
         if (allocated(error)) return
+
+        call dct_t23i(4, w2)
+        call dct_t3(4, x2, w2)
+        call check(error, maxval(abs(x2-x)) < eps, "dct_t3 failed")
+        if (allocated(error)) return
+
         call dcosqb(4, x, w)
         call check(error, sum(abs(x/(4*4) - [real(kind=rk) :: 1, 2, 3, 4])) < eps, &
                    "`dcosqb` failed.")
+        if (allocated(error)) return
+
+        call dct_t2(4, x2, w2)
+        call check(error, maxval(abs(x2-x)) < eps, "dct_t2 failed")
 
     end subroutine test_classic_dct