diff --git a/doc/specs/fftpack.md b/doc/specs/fftpack.md
index b55891b..edbbae2 100644
--- a/doc/specs/fftpack.md
+++ b/doc/specs/fftpack.md
@@ -857,7 +857,7 @@ end program demo_dcosqi
 
 ### `dcosqf`
 
-#### Decsription
+#### Description
 
 Computes the fast fourier transform of quarter wave data. 
 That is, `dcosqf` computes the coefficients in a cosine series representation with only odd wave numbers. 
@@ -927,7 +927,7 @@ end program demo_dcosqf
 
 ### `dcosqb`
 
-#### Decsription
+#### Description
 
 Computes the fast fourier transform of quarter wave data.
 That is, `dcosqb` computes a sequence from its representation in terms of a cosine series with odd wave numbers. 
@@ -989,7 +989,7 @@ This initialization does not have to be repeated so long as `n` remains unchange
 program demo_dcosqb
     use fftpack, only: dcosqi, dcosqf, dcosqb
     real(kind=8) :: w(3*4 + 15)
-    real(kind=8) :: x(4) = [4, 3, 5, 10]
+    real(kind=8) :: x(4) = [1, 2, 3, 4]
     call dcosqi(4, w) 
     call dcosqf(4, x, w) 
     call dcosqb(4, x, w)    !! `x`: [1.0, 2.0, 3.0, 4.0] * 4 * n, n = 4, which is unnormalized.
@@ -1000,7 +1000,7 @@ end program demo_dcosqb
 
 #### Description
 
-Forward transform of quarter wave data.
+Forward cosine-transform of quarter wave data.
 
 #### Status
 
@@ -1087,6 +1087,286 @@ program demo_iqct
 end program demo_iqct
 ```
 
+## Sine transform with odd wave numbers
+
+### `dsinqi`
+
+#### Description
+
+Initializes the array `wsave` which is used in both `dsinqf` and `dsinqb`. 
+The prime factorization of `n` together with a tabulation of the trigonometric functions are computed and stored in `wsave`.
+
+#### Status
+
+Experimental
+
+#### Class
+
+Pure subroutine.
+
+#### Syntax
+
+`call [[fftpack(module):dsinqi(interface)]](n, wsave)`
+
+#### Arguments
+
+`n`: Shall be an `integer` scalar.
+This argument is `intent(in)`.  
+The length of the array to be transformed. 
+The method is most efficient when `n` is a product of small primes.
+
+`wsave`: Shall be a `real` and rank-1 array.
+This argument is `intent(out)`.  
+A work array which must be dimensioned at least `3*n+15`.
+The same work array can be used for both `dsinqf` and `dsinqb`
+as long as `n` remains unchanged. 
+Different `wsave` arrays are required for different values of `n`.
+The contents of `wsave` must not be changed between calls of `dsinqf` or `dsinqb`.
+
+#### Example
+
+```fortran
+program demo_dsinqi
+    use fftpack, only: dsinqi
+    real(kind=8) :: w(3*4 + 15)
+    call dsinqi(4, w)   !! Initializes the array `w` which is used in both `dsinqf` and `dsinqb`. 
+end program demo_dsinqi
+```
+
+### `dsinqf`
+
+#### Description
+
+Computes the fast fourier transform of quarter wave data. 
+That is, `dsinqf` computes the coefficients in a sine series representation with only odd wave numbers. 
+The transform is defined below at output parameter `x`.
+
+`dsinqf` is the unnormalized inverse of `dsinqb` since a call of `dsinqf` followed by a call of `dsinqb` will multiply the input sequence `x` by `4*n`.
+
+The array `wsave` which is used by subroutine `dsinqf` must be initialized by calling subroutine `dsinqi(n,wsave)`.
+
+#### Status
+
+Experimental
+
+#### Class
+
+Pure subroutine.
+
+#### Syntax
+
+`call [[fftpack(module):dsinqf(interface)]](n, x, wsave)`
+
+#### Arguments
+
+`n`: Shall be an `integer` scalar.
+This argument is `intent(in)`.  
+The length of the array `x` to be transformed. 
+The method is most efficient when `n` is a product of small primes.
+
+`x`: Shall be a `real` and rank-1 array.
+This argument is `intent(inout)`.  
+An array which contains the sequence to be transformed.
+```
+for i=1,...,n
+
+    x(i) = (-1)**(i-1)*x(n)
+
+    + the sum from k=1 to k=n-1 of
+
+    2*x(k)*sin((2*i-1)*k*pi/(2*n))
+
+    a call of dsinqf followed by a call of
+    dsinqb will multiply the sequence x by 4*n.
+    therefore dsinqb is the unnormalized inverse
+    of dsinqf.
+```
+
+`wsave`: Shall be a `real` and rank-1 array.
+This argument is `intent(in)`.  
+A work array which must be dimensioned at least `3*n+15`
+in the program that calls `dsinqf`. 
+The `wsave` array must be initialized by calling subroutine `dsinqi(n,wsave)` and a different `wsave` array must be used for each different value of `n`. 
+This initialization does not have to be repeated so long as `n` remains unchanged thus subsequent transforms can be obtained faster than the first.
+
+##### Warning
+
+`wsave` contains initialization calculations which must not be destroyed between calls of `dsinqf` or `dsinqb`.
+
+#### Example
+
+```fortran
+program demo_dsinqf
+    use fftpack, only: dsinqi, dsinqf
+    real(kind=8) :: w(3*4 + 15)
+    real(kind=8) :: x(4) = [1, 2, 3, 4]
+    call dsinqi(4, w) 
+    call dsinqf(4, x, w)    !! `x`: [12.0, -9.10, 2.62, -1.51]
+end program demo_dsinqf
+```
+
+### `dsinqb`
+
+#### Description
+
+Computes the fast fourier transform of quarter wave data.
+That is, `dsinqb` computes a sequence from its representation in terms of a sine series with odd wave numbers. 
+The transform is defined below at output parameter `x`.
+
+`dsinqb` is the unnormalized inverse of `dsinqf` since a call of `dsinqb` followed by a call of `dsinqf` will multiply the input sequence `x` by `4*n`.
+
+The array `wsave` which is used by subroutine `dsinqb` must be initialized by calling subroutine `dsinqi(n,wsave)`.
+
+#### Status
+
+Experimental
+
+#### Class
+
+Pure subroutine.
+
+#### Syntax
+
+`call [[fftpack(module):dsinqf(interface)]](n, x, wsave)`
+
+#### Arguments
+
+`n`: Shall be an `integer` scalar.
+This argument is `intent(in)`.  
+The length of the array `x` to be transformed. 
+The method is most efficient when `n` is a product of small primes.
+
+`x`: Shall be a `real` and rank-1 array.
+This argument is `intent(inout)`.  
+An array which contains the sequence to be transformed.
+```
+for i=1,...,n
+
+    x(i)= the sum from k=1 to k=n of
+
+    4*x(k)*sin((2k-1)*i*pi/(2*n))
+
+    a call of dsinqb followed by a call of
+    dsinqf will multiply the sequence x by 4*n.
+    therefore dsinqf is the unnormalized inverse
+    of dsinqb.
+```
+
+`wsave`: Shall be a `real` and rank-1 array.
+This argument is `intent(in)`.  
+A work array which must be dimensioned at least `3*n+15`
+in the program that calls `dsinqb`. 
+The `wsave` array must be initialized by calling subroutine `dsinqi(n,wsave)` and a different `wsave` array must be used for each different value of `n`. 
+This initialization does not have to be repeated so long as `n` remains unchanged thus subsequent transforms can be obtained faster than the first.
+
+##### Warning
+
+`wsave` contains initialization calculations which must not be destroyed between calls of `dsinqf` or `dsinqb`.
+
+#### Example
+
+```fortran
+program demo_dsinqb
+    use fftpack, only: dsinqi, dsinqf, dsinqb
+    real(kind=8) :: w(3*4 + 15)
+    real(kind=8) :: x(4) = [1, 2, 3, 4]
+    call dsinqi(4, w) 
+    call dsinqf(4, x, w) 
+    call dsinqb(4, x, w)    !! `x`: [1.0, 2.0, 3.0, 4.0] * 4 * n, n = 4, which is unnormalized.
+end program demo_dsinqb
+```
+
+### `qst`
+
+#### Description
+
+Forward sine-transform of quarter wave data.
+
+#### Status
+
+Experimental.
+
+#### Class
+
+Pure function.
+
+#### Syntax
+
+`result = [[fftpack(module):qst(interface)]](x [, n])`
+
+#### Argument
+
+`x`: Shall be a `real` and rank-1 array.
+This argument is `intent(in)`.  
+The data to transform.
+
+`n`: Shall be an `integer` scalar.
+This argument is `intent(in)` and `optional`.  
+Defines the length of the Fourier transform. If `n` is not specified (the default) then `n = size(x)`. If `n <= size(x)`, `x` is truncated, if `n > size(x)`, `x` is zero-padded.
+
+#### Return value
+
+Returns a `real` and rank-1 array, the Quarter-Sine Transform (QST) of `x`.
+
+#### Notes
+
+Within numerical accuracy, `x == iqst(qst(x))/(4*size(x))`.
+
+#### Example
+
+```fortran
+program demo_qst
+    use fftpack, only: qst
+    real(kind=8) :: x(4) = [1, 2, 3, 4]
+    print *, qst(x,3)      !! [7.4, -1.0, 0.54].
+    print *, qst(x)        !! [13.1, -1.62, 0.72, -0.52].
+    print *, qst(x,5)      !! [15.4, 2.57, -4.0, 4.4, -4.49].
+end program demo_qst
+```
+
+### `iqst`
+
+#### Description
+
+Unnormalized inverse of `qst`.
+
+#### Status
+
+Experimental.
+
+#### Class
+
+Pure function.
+
+#### Syntax
+
+`result = [[fftpack(module):iqst(interface)]](x [, n])`
+
+#### Argument
+
+`x`: Shall be a `real` array.
+This argument is `intent(in)`. 
+Transformed data to invert.
+
+`n`: Shall be an `integer` scalar.
+This argument is `intent(in)` and `optional`.  
+Defines the length of the Fourier transform. If `n` is not specified (the default) then `n = size(x)`. If `n <= size(x)`, `x` is truncated, if `n > size(x)`, `x` is zero-padded.
+
+#### Return value
+
+Returns a `real` and rank-1 array, the unnormalized inverse Quarter-Sine Transform (iQST).
+
+#### Example
+
+```fortran
+program demo_iqst
+    use fftpack, only: qst, iqst
+    real(kind=8) :: x(4) = [1, 2, 3, 4]
+    print *, iqst(qst(x))/(4.0*4.0)         !! [1.0, 2.0, 3.0, 4.0]
+    print *, iqst(qst(x), 3)/(4.0*3.0)      !! [1.77, 3.58, 5.16]
+end program demo_iqst
+```
+
 ## Cosine transform of a real even sequence
 
 ### `dcosti`
@@ -1113,7 +1393,7 @@ Pure subroutine.
 `n`: Shall be a `integer` scalar.
 This argument is `intent(in)`.  
 The length of the sequence to be transformed.  
-The method is most efficient when n-1 is a product of small primes.
+The method is most efficient when `n-1` is a product of small primes.
 
 `wsave`: Shall be a `real` and rank-1 array.
 This argument is `intent(out)`.  
@@ -1139,7 +1419,6 @@ Computes the discrete fourier cosine transform of an even sequence `x(i)`.
 The transform is defined below at output parameter `x`.
 
 `dcost` is the unnormalized inverse of itself since a call of `dcost` followed by another call of `dcost` will multiply the input sequence `x` by `2*(n-1)`. 
-The transform is defined below at output parameter `x`.
 
 The array `wsave` which is used by subroutine `dcost` must be initialized by calling subroutine `dcosti(n,wsave)`.
 
@@ -1295,7 +1574,209 @@ program demo_idct
     print *, idct(dct(x), 3)          !! (unnormalized): [7.0, 15.0, 23.0]
 end program demo_idct
 ```
+## Sine transform of a real odd sequence
+
+### `dsinti`
+
+#### Description
+
+Initializes the array `wsave` which is used in subroutine `dsint`. 
+The prime factorization of `n` together with a tabulation of the trigonometric functions are computed and stored in `wsave`.
+
+#### Status
+
+Experimental
+
+#### Class
+
+Pure subroutine.
+
+#### Syntax
+
+`call [[fftpack(module):dsinti(interface)]](n , wsave)`
+
+#### Arguments
+
+`n`: Shall be a `integer` scalar.
+This argument is `intent(in)`.  
+The length of the sequence to be transformed.  
+The method is most efficient when `n+1` is a product of small primes.
+
+`wsave`: Shall be a `real` and rank-1 array.
+This argument is `intent(out)`.  
+A work array which must be dimensioned at least `int(2.5*n+15)`. 
+Different `wsave` arrays are required for different values of `n`. 
+The contents of `wsave` must not be changed between calls of `dsint`.
+
+#### Example
+
+```fortran
+program demo_dsinti
+    use fftpack, only: dsinti
+    real(kind=8) :: w(int(2.5*4 + 15))
+    call dsinti(4, w)   !! Initializes the array `w` which is used in subroutine `dsint`. 
+end program demo_dsinti
+```
+
+### `dsint`
+
+#### Description
+
+Computes the discrete fourier sinine transform of an odd sequence `x(i)`. 
+The transform is defined below at output parameter `x`.
+
+`dsint` is the unnormalized inverse of itself since a call of `dsint` followed by another call of `dsint` will multiply the input sequence `x` by `2*(n+1)`.
+
+The array `wsave` which is used by subroutine `dsint` must be initialized by calling subroutine `dsinti(n,wsave)`.
+
+#### Status
+
+Experimental
+
+#### Class
+
+Pure subroutine.
+
+#### Syntax
+
+`call [[fftpack(module):dsint(interface)]](n, x, wsave)`
+
+#### Arguments
+
+`n`: Shall be a `integer` scalar.
+This argument is `intent(in)`.  
+The length of the sequence `x`. 
+The method is most efficient when `n+1` is a product of small primes.
+
+`x`: Shall be a `real` and rank-1 array.
+This argument is `intent(inout)`.
+An array which contains the sequence to be transformed.
+```
+for i=1,...,n
+
+    x(i)= the sum from k=1 to k=n
+
+        2*x(k)*sin(k*i*pi/(n+1))
+
+    a call of dsint followed by another call of
+    dsint will multiply the sequence x by 2*(n+1).
+    hence dsint is the unnormalized inverse
+    of itself.
+```
+
+`wsave`: Shall be a `real` and rank-1 array.
+This argument is `intent(in)`.  
+A work array which must be dimensioned at least `int(2.5*n+15)` in the program that calls `dsint`. 
+The `wsave` array must be initialized by calling subroutine `dsinti(n,wsave)` and a different `wsave` array must be used for each different value of `n`. 
+This initialization does not have to be repeated so long as `n` remains unchanged thus subsequent transforms can be obtained faster than the first.
+Contains initialization calculations which must not be destroyed between calls of `dsint`.
+
+#### Example
+
+```fortran
+program demo_dsint
+    use fftpack, only: dsinti, dsint
+    real(kind=8) :: x(4) = [1, 2, 3, 4]
+    real(kind=8) :: w(int(2.5*4 + 15))
+    call dsinti(4, w)
+    call dsint(4, x, w)     !! Computes the discrete fourier sinine (forward) transform of an even sequence, `x`(unnormalized): [15.4, -6.9, 3.6, -1.6]
+    call dsint(4, x, w)     !! Computes the discrete fourier sinine (backward) transform of an even sequence, `x`(unnormalized): [10.0, 20.0, 30.0, 40.0]
+end program demo_dsint
+```
+
+### `dst`
+
+#### Description
+
+Discrete fourier sine (forward) transform of an odd sequence. 
+
+#### Status
+
+Experimental.
+
+#### Class
+
+Pure function.
+
+#### Syntax
+
+`result = [[fftpack(module):dst(interface)]](x [, n])`
+
+#### Argument
+
+`x`: Shall be a `real` and rank-1 array.
+This argument is `intent(in)`.  
+The data to transform.
+
+`n`: Shall be an `integer` scalar.
+This argument is `intent(in)` and `optional`.  
+Defines the length of the Fourier transform. If `n` is not specified (the default) then `n = size(x)`. If `n <= size(x)`, `x` is truncated, if `n > size(x)`, `x` is zero-padded.
+
+#### Return value
+
+Returns a `real` and rank-1 array, the Discrete-Sine Transform (DST) of `x`.
+
+#### Notes
+
+Within numerical accuracy, `y == dst(idst(y))/2*(size(y) + 1)`.
+
+#### Example
+
+```fortran
+program demo_dst
+    use fftpack, only: dst
+    real(kind=8) :: x(4) = [1, 2, 3, 4]
+    print *, dst(x,3)      !! [9.7, -4.0, 1.7].
+    print *, dst(x)        !! [15.4, -6.9, 3.6, -1.6].
+    print *, dst(x,5)      !! [17.4, -1.7, -4.0, 5.2, -3.4].
+    print *, dst(dst(x))/(2*(4 + 1))   !! (normalized): [1.0, 2.0, 3.0, 4.0]
+end program demo_dst
+```
 
+### `idst`
+
+#### Description
+
+Unnormalized inverse of `dst`.
+In fact, `idst` and `dst` have the same effect, `idst` = `dst`.
+
+#### Status
+
+Experimental.
+
+#### Class
+
+Pure function.
+
+#### Syntax
+
+`result = [[fftpack(module):idst(interface)]](x [, n])`
+
+#### Argument
+
+`x`: Shall be a `real` array.
+This argument is `intent(in)`. 
+Transformed data to invert.
+
+`n`: Shall be an `integer` scalar.
+This argument is `intent(in)` and `optional`.  
+Defines the length of the Fourier transform. If `n` is not specified (the default) then `n = size(x)`. If `n <= size(x)`, `x` is truncated, if `n > size(x)`, `x` is zero-padded.
+
+#### Return value
+
+Returns a `real` and rank-1 array, the inverse Discrete-Cosine Transform (iDST) of `x`.
+
+#### Example
+
+```fortran
+program demo_idst
+    use fftpack, only: dst, idst
+    real(kind=8) :: x(4) = [1, 2, 3, 4]
+    print *, idst(dst(x))/(2*(4+1))   !! (normalized):   [1.0, 2.0, 3.0, 4.0]
+    print *, idst(idst(x))/(2*(4+1))  !! (normalized):   [1.0, 2.0, 3.0, 4.0]
+    print *, idst(dst(x), 3)          !! (unnormalized): [13.1, 23.5, 40.7]
+end program demo_idst
+```
 
 ## Utility functions
 
diff --git a/fpm.toml b/fpm.toml
index 03e3b07..7d8560c 100644
--- a/fpm.toml
+++ b/fpm.toml
@@ -72,6 +72,21 @@ name = "fftpack_iqct"
 source-dir = "test"
 main = "test_fftpack_iqct.f90"
 
+[[test]]
+name = "fftpack_dsinq"
+source-dir = "test"
+main = "test_fftpack_dsinq.f90"
+
+[[test]]
+name = "fftpack_qst"
+source-dir = "test"
+main = "test_fftpack_qst.f90"
+
+[[test]]
+name = "fftpack_iqst"
+source-dir = "test"
+main = "test_fftpack_iqst.f90"
+
 [[test]]
 name = "fftpack_dcost"
 source-dir = "test"
@@ -82,6 +97,16 @@ name = "fftpack_dct"
 source-dir = "test"
 main = "test_fftpack_dct.f90"
 
+[[test]]
+name = "fftpack_dsint"
+source-dir = "test"
+main = "test_fftpack_dsint.f90"
+
+[[test]]
+name = "fftpack_dst"
+source-dir = "test"
+main = "test_fftpack_dst.f90"
+
 # `fftpack` utility routines
 [[test]]
 name = "fftpack_fftshift"
diff --git a/src/Makefile b/src/Makefile
index 2a9f61f..76a5a24 100644
--- a/src/Makefile
+++ b/src/Makefile
@@ -59,7 +59,10 @@ SRCF90 = \
 	fftpack_ifftshift.f90\
 	fftpack_qct.f90\
 	fftpack_iqct.f90\
-	fftpack_dct.f90
+	fftpack_qst.f90\
+	fftpack_iqst.f90\
+	fftpack_dct.f90\
+	fftpack_dst.f90
 
 OBJF := $(SRCF:.f=.o)
 OBJF90 := $(SRCF90:.f90=.o)
@@ -82,6 +85,9 @@ fftpack_rfft.o: fftpack.o
 fftpack_irfft.o: fftpack.o
 fftpack_qct.o: fftpack.o
 fftpack_iqct.o: fftpack.o
+fftpack_qst.o: fftpack.o
+fftpack_iqst.o: fftpack.o
 fftpack_dct.o: fftpack.o
+fftpack_dst.o: fftpack.o
 fftpack_fftshift.o: fftpack.o
 fftpack_ifftshift.o: fftpack.o
\ No newline at end of file
diff --git a/src/fftpack.f90 b/src/fftpack.f90
index ba3ebb8..34a6805 100644
--- a/src/fftpack.f90
+++ b/src/fftpack.f90
@@ -16,9 +16,15 @@ module fftpack
     public :: dcosqi, dcosqf, dcosqb
     public :: qct, iqct
 
+    public :: dsinqi, dsinqf, dsinqb
+    public :: qst, iqst
+
     public :: dcosti, dcost
     public :: dct, idct
 
+    public :: dsinti, dsint
+    public :: dst, idst
+
     interface
 
         !> Version: experimental
@@ -153,6 +159,38 @@ pure subroutine dcosqb(n, x, wsave)
             real(kind=dp), intent(in) :: wsave(*)
         end subroutine dcosqb
 
+        !> Version: experimental
+        !>
+        !> Initialize `dsinqf` and `dsinqb`.
+        !> ([Specification](../page/specs/fftpack.html#dsinqi))
+        pure subroutine dsinqi(n, wsave)
+            import dp
+            integer, intent(in) :: n
+            real(kind=dp), intent(out) :: wsave(*)
+        end subroutine dsinqi
+
+        !> Version: experimental
+        !>
+        !> Forward transform of quarter wave data.
+        !> ([Specification](../page/specs/fftpack.html#dsinqf))
+        pure subroutine dsinqf(n, x, wsave)
+            import dp
+            integer, intent(in) :: n
+            real(kind=dp), intent(inout) :: x(*)
+            real(kind=dp), intent(in) :: wsave(*)
+        end subroutine dsinqf
+
+        !> Version: experimental
+        !>
+        !> Unnormalized inverse of `dsinqf`.
+        !> ([Specification](../page/specs/fftpack.html#dsinqb))
+        pure subroutine dsinqb(n, x, wsave)
+            import dp
+            integer, intent(in) :: n
+            real(kind=dp), intent(inout) :: x(*)
+            real(kind=dp), intent(in) :: wsave(*)
+        end subroutine dsinqb
+
         !> Version: experimental
         !>
         !> Initialize `dcost`. ([Specification](../page/specs/fftpack.html#dcosti))
@@ -173,6 +211,26 @@ pure subroutine dcost(n, x, wsave)
             real(kind=dp), intent(in) :: wsave(*)
         end subroutine dcost
 
+        !> Version: experimental
+        !>
+        !> Initialize `dsint`. ([Specification](../page/specs/fftpack.html#dsinti))
+        pure subroutine dsinti(n, wsave)
+            import dp
+            integer, intent(in) :: n
+            real(kind=dp), intent(out) :: wsave(*)
+        end subroutine dsinti
+
+        !> Version: experimental
+        !>
+        !> Discrete fourier sine transform of an odd sequence.
+        !> ([Specification](../page/specs/fftpack.html#dsint))
+        pure subroutine dsint(n, x, wsave)
+            import dp
+            integer, intent(in) :: n
+            real(kind=dp), intent(inout) :: x(*)
+            real(kind=dp), intent(in) :: wsave(*)
+        end subroutine dsint
+
     end interface
 
     !> Version: experimental
@@ -225,7 +283,7 @@ end function irfft_dp
 
     !> Version: experimental
     !>
-    !> Forward transform of quarter wave data.
+    !> Forward cosine-transform of quarter wave data.
     !> ([Specifiction](../page/specs/fftpack.html#qct))
     interface qct
         pure module function qct_dp(x, n) result(result)
@@ -237,7 +295,7 @@ end function qct_dp
 
     !> Version: experimental
     !>
-    !> Backward transform of quarter wave data.
+    !> Backward cosine-transform of quarter wave data.
     !> ([Specifiction](../page/specs/fftpack.html#iqct))
     interface iqct
         pure module function iqct_dp(x, n) result(result)
@@ -247,6 +305,30 @@ pure module function iqct_dp(x, n) result(result)
         end function iqct_dp
     end interface iqct
 
+    !> Version: experimental
+    !>
+    !> Forward sine-transform of quarter wave data.
+    !> ([Specifiction](../page/specs/fftpack.html#qst))
+    interface qst
+        pure module function qst_dp(x, n) result(result)
+            real(kind=dp), intent(in) :: x(:)
+            integer, intent(in), optional :: n
+            real(kind=dp), allocatable :: result(:)
+        end function qst_dp
+    end interface qst
+
+    !> Version: experimental
+    !>
+    !> Backward sine-transform of quarter wave data.
+    !> ([Specifiction](../page/specs/fftpack.html#iqst))
+    interface iqst
+        pure module function iqst_dp(x, n) result(result)
+            real(kind=dp), intent(in) :: x(:)
+            integer, intent(in), optional :: n
+            real(kind=dp), allocatable :: result(:)
+        end function iqst_dp
+    end interface iqst
+
     !> Version: experimental
     !>
     !> Discrete fourier cosine (forward) transform of an even sequence.
@@ -267,6 +349,26 @@ end function dct_dp
         module procedure :: dct_dp
     end interface idct
 
+    !> Version: experimental
+    !>
+    !> Discrete fourier sine (forward) transform of an odd sequence.
+    !> ([Specification](../page/specs/fftpack.html#dst))
+    interface dst
+        pure module function dst_dp(x, n) result(result)
+            real(kind=dp), intent(in) :: x(:)
+            integer, intent(in), optional :: n
+            real(kind=dp), allocatable :: result(:)
+        end function dst_dp
+    end interface dst
+
+    !> Version: experimental
+    !>
+    !> Discrete fourier sine (backward) transform of an even sequence.
+    !> ([Specification](../page/specs/fftpack.html#idst))
+    interface idst
+        module procedure :: dst_dp
+    end interface idst
+
     !> Version: experimental
     !>
     !> Shifts zero-frequency component to center of spectrum.
diff --git a/src/fftpack_dst.f90 b/src/fftpack_dst.f90
new file mode 100644
index 0000000..89ac15c
--- /dev/null
+++ b/src/fftpack_dst.f90
@@ -0,0 +1,36 @@
+submodule(fftpack) fftpack_dst
+
+contains
+
+    !> Discrete fourier sine transform of an odd sequence.
+    pure module function dst_dp(x, n) result(result)
+        real(kind=dp), intent(in) :: x(:)
+        integer, intent(in), optional :: n
+        real(kind=dp), allocatable :: result(:)
+
+        integer :: lenseq, lensav, i
+        real(kind=dp), allocatable :: wsave(:)
+
+        if (present(n)) then
+            lenseq = n
+            if (lenseq <= size(x)) then
+                result = x(:lenseq)
+            else if (lenseq > size(x)) then
+                result = [x, (0.0_dp, i=1, lenseq - size(x))]
+            end if
+        else
+            lenseq = size(x)
+            result = x
+        end if
+
+        !> Initialize FFT
+        lensav = int(2.5_dp*lenseq + 15)
+        allocate (wsave(lensav))
+        call dsinti(lenseq, wsave)
+
+        !> Discrete fourier sine transformation
+        call dsint(lenseq, result, wsave)
+
+    end function dst_dp
+
+end submodule fftpack_dst
diff --git a/src/fftpack_iqct.f90 b/src/fftpack_iqct.f90
index f116fef..3736477 100644
--- a/src/fftpack_iqct.f90
+++ b/src/fftpack_iqct.f90
@@ -2,7 +2,7 @@
 
 contains
 
-    !> Backward transform of quarter wave data.
+    !> Backward cosine-transform of quarter wave data.
     pure module function iqct_dp(x, n) result(result)
         real(kind=dp), intent(in) :: x(:)
         integer, intent(in), optional :: n
diff --git a/src/fftpack_iqst.f90 b/src/fftpack_iqst.f90
new file mode 100644
index 0000000..1fe2d26
--- /dev/null
+++ b/src/fftpack_iqst.f90
@@ -0,0 +1,36 @@
+submodule(fftpack) fftpack_iqst
+
+contains
+
+    !> Backward sine-transform of quarter wave data.
+    pure module function iqst_dp(x, n) result(result)
+        real(kind=dp), intent(in) :: x(:)
+        integer, intent(in), optional :: n
+        real(kind=dp), allocatable :: result(:)
+
+        integer :: lenseq, lensav, i
+        real(kind=dp), allocatable :: wsave(:)
+
+        if (present(n)) then
+            lenseq = n
+            if (lenseq <= size(x)) then
+                result = x(:lenseq)
+            else if (lenseq > size(x)) then
+                result = [x, (0.0_dp, i=1, lenseq - size(x))]
+            end if
+        else
+            lenseq = size(x)
+            result = x
+        end if
+
+        !> Initialize FFT
+        lensav = 3*lenseq + 15
+        allocate (wsave(lensav))
+        call dsinqi(lenseq, wsave)
+
+        !> Backward transformation
+        call dsinqb(lenseq, result, wsave)
+
+    end function iqst_dp
+
+end submodule fftpack_iqst
diff --git a/src/fftpack_qct.f90 b/src/fftpack_qct.f90
index 582918b..da9c62f 100644
--- a/src/fftpack_qct.f90
+++ b/src/fftpack_qct.f90
@@ -2,7 +2,7 @@
 
 contains
 
-    !> Forward transform of quarter wave data.
+    !> Forward cosine-transform of quarter wave data.
     pure module function qct_dp(x, n) result(result)
         real(kind=dp), intent(in) :: x(:)
         integer, intent(in), optional :: n
diff --git a/src/fftpack_qst.f90 b/src/fftpack_qst.f90
new file mode 100644
index 0000000..30debf9
--- /dev/null
+++ b/src/fftpack_qst.f90
@@ -0,0 +1,36 @@
+submodule(fftpack) fftpack_qst
+
+contains
+
+    !> Forward sine-transform of quarter wave data.
+    pure module function qst_dp(x, n) result(result)
+        real(kind=dp), intent(in) :: x(:)
+        integer, intent(in), optional :: n
+        real(kind=dp), allocatable :: result(:)
+
+        integer :: lenseq, lensav, i
+        real(kind=dp), allocatable :: wsave(:)
+
+        if (present(n)) then
+            lenseq = n
+            if (lenseq <= size(x)) then
+                result = x(:lenseq)
+            else if (lenseq > size(x)) then
+                result = [x, (0.0_dp, i=1, lenseq - size(x))]
+            end if
+        else
+            lenseq = size(x)
+            result = x
+        end if
+
+        !> Initialize FFT
+        lensav = 3*lenseq + 15
+        allocate (wsave(lensav))
+        call dsinqi(lenseq, wsave)
+
+        !> Forward transformation
+        call dsinqf(lenseq, result, wsave)
+
+    end function qst_dp
+
+end submodule fftpack_qst
diff --git a/test/Makefile b/test/Makefile
index 57a618a..1f9d30b 100644
--- a/test/Makefile
+++ b/test/Makefile
@@ -9,8 +9,13 @@ all: tstfft \
 	fftpack_dcosq \
 	fftpack_qct \
 	fftpack_iqct \
+	fftpack_dsinq \
+	fftpack_qst \
+	fftpack_iqst \
 	fftpack_dcost \
-	fftpack_dct
+	fftpack_dct \
+	fftpack_dsint \
+	fftpack_dst
 
 # Orginal test
 tstfft: tstfft.o
@@ -50,6 +55,18 @@ fftpack_iqct: test_fftpack_iqct.f90
 	$(FC) $(FFLAGS) $< -L../src -l$(LIB) -I../src -o $@.x
 	./fftpack_iqct.x
 
+fftpack_dsinq: test_fftpack_dsinq.f90
+	$(FC) $(FFLAGS) $< -L../src -l$(LIB) -I../src -o $@.x
+	./fftpack_dsinq.x
+
+fftpack_qst: test_fftpack_qst.f90
+	$(FC) $(FFLAGS) $< -L../src -l$(LIB) -I../src -o $@.x
+	./fftpack_qst.x
+
+fftpack_iqst: test_fftpack_iqst.f90
+	$(FC) $(FFLAGS) $< -L../src -l$(LIB) -I../src -o $@.x
+	./fftpack_iqst.x
+
 fftpack_dcost: test_fftpack_dcost.f90
 	$(FC) $(FFLAGS) $< -L../src -l$(LIB) -I../src -o $@.x
 	./fftpack_dcost.x
@@ -58,6 +75,14 @@ fftpack_dct: test_fftpack_dct.f90
 	$(FC) $(FFLAGS) $< -L../src -l$(LIB) -I../src -o $@.x
 	./fftpack_dct.x
 
+fftpack_dsint: test_fftpack_dsint.f90
+	$(FC) $(FFLAGS) $< -L../src -l$(LIB) -I../src -o $@.x
+	./fftpack_dsint.x
+
+fftpack_dst: test_fftpack_dst.f90
+	$(FC) $(FFLAGS) $< -L../src -l$(LIB) -I../src -o $@.x
+	./fftpack_dst.x
+
 # `fftpack` utility routines
 fftpack_fftshift: test_fftpack_fftshift.f90
 	$(FC) $(FFLAGS) $< -L../src -l$(LIB) -I../src -o $@.x
diff --git a/test/test_fftpack_dct.f90 b/test/test_fftpack_dct.f90
index ef1c6f8..7b5c949 100644
--- a/test/test_fftpack_dct.f90
+++ b/test/test_fftpack_dct.f90
@@ -27,7 +27,6 @@ end subroutine test_fftpack_dct
     subroutine test_fftpack_idct
         use fftpack, only: dct, idct
         use iso_fortran_env, only: dp => real64
-        real(kind=dp) :: eps = 1.0e-10_dp
         real(kind=dp) :: x(4) = [1, 2, 3, 4]
 
         call check(all(idct(dct(x))/(2.0_dp*(4.0_dp - 1.0_dp)) == [real(kind=dp) :: 1, 2, 3, 4]), &
diff --git a/test/test_fftpack_dsinq.f90 b/test/test_fftpack_dsinq.f90
new file mode 100644
index 0000000..afc9558
--- /dev/null
+++ b/test/test_fftpack_dsinq.f90
@@ -0,0 +1,32 @@
+program tester
+
+    call test_fftpack_dsinq_real
+    print *, "All tests in `test_fftpack_dsinq` passed."
+
+contains
+
+    subroutine check(condition, msg)
+        logical, intent(in) :: condition
+        character(*), intent(in) :: msg
+        if (.not. condition) error stop msg
+    end subroutine check
+
+    subroutine test_fftpack_dsinq_real
+        use fftpack, only: dsinqi, dsinqf, dsinqb
+        use iso_fortran_env, only: dp => real64
+        real(kind=dp) :: w(3*4 + 15)
+        real(kind=dp) :: x(4) = [1, 2, 3, 4]
+        real(kind=dp) :: eps = 1.0e-10_dp
+
+        call dsinqi(4, w)
+        call dsinqf(4, x, w)
+        call check(sum(abs(x - [13.137071184544091_dp, -1.6199144044217753_dp, &
+                                0.72323134608584416_dp, -0.51978306494828974_dp])) < eps, msg="`dsinqf` failed.")
+
+        call dsinqb(4, x, w)
+        call check(sum(abs(x/(4.0_dp*4.0_dp) - [1.0000000000000002_dp, 2.0_dp, 3.0000000000000009_dp, 4.0_dp])) < eps, &
+                   msg="`dsinqb` failed.")
+
+    end subroutine test_fftpack_dsinq_real
+
+end program tester
diff --git a/test/test_fftpack_dsint.f90 b/test/test_fftpack_dsint.f90
new file mode 100644
index 0000000..fa762fe
--- /dev/null
+++ b/test/test_fftpack_dsint.f90
@@ -0,0 +1,33 @@
+program tester
+
+    call test_fftpack_dsint_real
+    print *, "All tests in `test_fftpack_dsint` passed."
+
+contains
+
+    subroutine check(condition, msg)
+        logical, intent(in) :: condition
+        character(*), intent(in) :: msg
+        if (.not. condition) error stop msg
+    end subroutine check
+
+    subroutine test_fftpack_dsint_real
+        use fftpack, only: dsinti, dsint
+        use iso_fortran_env, only: dp => real64
+        real(kind=dp) :: w(3*4 + 15)
+        real(kind=dp) :: x(4) = [1, 2, 3, 4]
+        real(kind=dp) :: eps = 1.0e-10_dp
+
+        call dsinti(4, w)
+        call dsint(4, x, w)
+        call check(sum(abs(x - [15.388417685876266_dp, -6.8819096023558668_dp, &
+                                3.6327126400268046_dp, -1.6245984811645318_dp])) < eps, &
+                   msg="`dsinti` failed.")
+
+        call dsint(4, x, w)
+        call check(sum(abs(x/(2.0_dp*(4.0_dp + 1.0_dp)) - [real(kind=dp) :: 1, 2, 3, 4])) < eps, &
+                   msg="`dsint` failed.")
+
+    end subroutine test_fftpack_dsint_real
+
+end program tester
diff --git a/test/test_fftpack_dst.f90 b/test/test_fftpack_dst.f90
new file mode 100644
index 0000000..48071a3
--- /dev/null
+++ b/test/test_fftpack_dst.f90
@@ -0,0 +1,46 @@
+program tester
+
+    call test_fftpack_dst()
+    call test_fftpack_idst()
+    print *, "All tests in `test_fftpack_dst` passed."
+
+contains
+
+    subroutine check(condition, msg)
+        logical, intent(in) :: condition
+        character(*), intent(in) :: msg
+        if (.not. condition) error stop msg
+    end subroutine check
+
+    subroutine test_fftpack_dst
+        use fftpack, only: dst
+        use iso_fortran_env, only: dp => real64
+        real(kind=dp) :: eps = 1.0e-10_dp
+        real(kind=dp) :: x(3) = [1, 2, 3]
+
+        call check(sum(abs(dst(x, 2) - [5.1961524227066320_dp, -1.7320508075688772_dp])) < eps, &
+                   msg="`dst(x, 2)` failed.")
+        call check(all(dst(x, 3) == dst(x)), msg="`dst(x, 3)` failed.")
+        call check(sum(abs(dst(x, 4) - [10.686135667536481_dp, 0.72654252800536079_dp, &
+                                        -3.9757394903344245_dp, 3.0776835371752531_dp])) < eps, &
+                   msg="`dst(x, 4)` failed.")
+
+    end subroutine test_fftpack_dst
+
+    subroutine test_fftpack_idst
+        use fftpack, only: dst, idst
+        use iso_fortran_env, only: dp => real64
+        real(kind=dp) :: eps = 1.0e-10_dp
+        real(kind=dp) :: x(4) = [1, 2, 3, 4]
+
+        call check(sum(abs(idst(dst(x))/(2.0_dp*(4.0_dp + 1.0_dp)) - [real(kind=dp) :: 1, 2, 3, 4])) < eps, &
+                   msg="`idst(dst(x))/(2.0_dp*(4.0_dp+1.0_dp))` failed.")
+        call check(sum(abs(idst(dst(x), 2)/(2.0_dp*(2.0_dp + 1.0_dp)) - [2.4556173659421150_dp, 6.4288897274009438_dp])) < eps, &
+                   msg="`idst(dst(x), 2)/(2.0_dp*(2.0_dp+1.0_dp))` failed.")
+        call check(all(idst(dst(x, 2), 4)/(2.0_dp*(4.0_dp + 1.0_dp)) == &
+                       [0.28138871112761987_dp, 0.78475214007354754_dp, 1.1919817084376492_dp, 0.94029999396468544_dp]), &
+                   msg="`idst(dst(x, 2), 4)/(2.0_dp*(4.0_dp+1.0_dp))` failed.")
+
+    end subroutine test_fftpack_idst
+
+end program tester
diff --git a/test/test_fftpack_iqst.f90 b/test/test_fftpack_iqst.f90
new file mode 100644
index 0000000..b2f73ea
--- /dev/null
+++ b/test/test_fftpack_iqst.f90
@@ -0,0 +1,32 @@
+program tester
+
+    call test_fftpack_iqst()
+    print *, "All tests in `test_fftpack_iqst` passed."
+
+contains
+
+    subroutine check(condition, msg)
+        logical, intent(in) :: condition
+        character(*), intent(in) :: msg
+        if (.not. condition) error stop msg
+    end subroutine check
+
+    subroutine test_fftpack_iqst
+        use fftpack, only: qst, iqst
+        use iso_fortran_env, only: dp => real64
+        real(kind=dp) :: eps = 1.0e-10_dp
+
+        real(kind=dp) :: x(4) = [1, 2, 3, 4]
+
+        call check(sum(abs(iqst(qst(x))/(4.0_dp*4.0_dp) - [1.0000000000000002_dp, 2.0_dp, &
+                                                           3.0000000000000009_dp, 4.0_dp])) < eps, &
+                   msg="`iqst(qst(x)/(4.0_dp*4.0_dp)` failed.")
+        call check(sum(abs(iqst(qst(x), 2)/(4.0_dp*2.0_dp) - [4.0719298296065567_dp, 7.3784927944829333_dp])) < eps, &
+                   msg="`iqst(qst(x), 2)/(4.0_dp*2.0_dp)` failed.")
+        call check(sum(abs(iqst(qst(x, 2), 4)/(4.0_dp*4.0_dp) - [0.19134171618254481_dp, 0.5_dp, &
+                                                                 0.84462319862073310_dp, 1.0_dp])) < eps, &
+                   msg="`iqst(qst(x, 2), 4)/(4.0_dp*4.0_dp)` failed.")
+
+    end subroutine test_fftpack_iqst
+
+end program tester
diff --git a/test/test_fftpack_qst.f90 b/test/test_fftpack_qst.f90
new file mode 100644
index 0000000..37a3ab2
--- /dev/null
+++ b/test/test_fftpack_qst.f90
@@ -0,0 +1,31 @@
+program tester
+
+    call test_fftpack_qst()
+    print *, "All tests in `test_fftpack_qst` passed."
+
+contains
+
+    subroutine check(condition, msg)
+        logical, intent(in) :: condition
+        character(*), intent(in) :: msg
+        if (.not. condition) error stop msg
+    end subroutine check
+
+    subroutine test_fftpack_qst
+        use fftpack, only: qst
+        use iso_fortran_env, only: dp => real64
+        real(kind=dp) :: eps = 1.0e-10_dp
+
+        real(kind=dp) :: x(3) = [9, -9, 3]
+
+        call check(sum(abs(qst(x, 2) - [3.7279220613578570_dp, 21.727922061357859_dp])) < eps, &
+                   msg="`qst(x, 2)` failed.")
+        call check(sum(abs(qst(x, 3) - qst(x))) < eps, &
+                   msg="`qst(x,3)` failed.")
+        call check(sum(abs(qst(x, 4) - [-0.29634308371852036_dp, 1.6058089296547653_dp, &
+                                        27.061653052370481_dp, 25.159501038997192_dp])) < eps, &
+                   msg="`qst(x, 4)` failed.")
+
+    end subroutine test_fftpack_qst
+
+end program tester