nonlin_polynomials.f90

! nonlin_polynomials.f90
 
 
! TO DO:
! - Multiplication (convolution)
! - Division (deconvolution)
! - C interface
 
module nonlin_polynomials
    use, intrinsic :: iso_fortran_env, only : int32, real64
    use linalg, only : eigen, solve_least_squares
    use ferror, only : errors
    use nonlin_constants, only : nl_invalid_input_error, nl_array_size_error, &
        nl_out_of_memory_error
    implicit none
 
private
public :: polynomial
public :: assignment(=)
public :: operator(+)
public :: operator(-)
public :: operator(*)
 
! ******************************************************************************
! INTERFACES
! ------------------------------------------------------------------------------
interface assignment(=)
    module procedure :: poly_equals
    module procedure :: poly_dbl_equals
    module procedure :: poly_equals_array
end interface
 
interface operator(+)
    module procedure :: poly_poly_add
end interface
 
interface operator(-)
    module procedure :: poly_poly_subtract
end interface
 
interface operator(*)
    module procedure :: poly_poly_mult
    module procedure :: poly_dbl_mult
    module procedure :: dbl_poly_mult
end interface
 
! ******************************************************************************
! TYPES
! ------------------------------------------------------------------------------
type polynomial
private
    real(real64), allocatable, dimension(:) :: m_coeffs
contains
    generic, public :: initialize => init_poly, init_poly_coeffs
    procedure, public :: order => get_poly_order
    procedure, public :: fit => poly_fit
    procedure, public :: fit_thru_zero => poly_fit_thru_zero
    generic, public :: evaluate => evaluate_real, evaluate_complex
    procedure, public :: companion_mtx => poly_companion_mtx
    procedure, public :: roots => poly_roots
    procedure, public :: get => get_poly_coefficient
    procedure, public :: get_all => get_poly_coefficients
    procedure, public :: set => set_poly_coefficient
 
    procedure :: evaluate_real => poly_eval_double
    procedure :: evaluate_complex => poly_eval_complex
    procedure :: init_poly
    procedure :: init_poly_coeffs
end type
 
contains
! ******************************************************************************
! POLYNOMIAL MEMBERS
! ------------------------------------------------------------------------------
    subroutine init_poly(this, order, err)
        ! Arguments
        class(polynomial), intent(inout) :: this
        integer(int32), intent(in) :: order
        class(errors), intent(inout), optional, target :: err
 
        ! Parameters
        real(real64), parameter :: zero = 0.0d0
 
        ! Local Variables
        integer(int32) :: n, istat
        class(errors), pointer :: errmgr
        type(errors), target :: deferr
 
        ! Initialization
        n = order + 1
        if (present(err)) then
            errmgr => err
        else
            errmgr => deferr
        end if
 
        ! Input Check
        if (order < 0) then
            ! ERROR: Negative order is not supported
            call errmgr%report_error("init_polynomial", &
                "A negative polynomial order is not supported.", &
                nl_invalid_input_error)
            return
        end if
 
        ! Process
        if (allocated(this%m_coeffs)) deallocate(this%m_coeffs)
        allocate(this%m_coeffs(n), stat = istat)
        if (istat /= 0) then
            ! ERROR: Out of memory
            call errmgr%report_error("init_polynomial", &
                "Insufficient memory available.", nl_out_of_memory_error)
            return
        end if
        this%m_coeffs = zero
    end subroutine
 
! ------------------------------------------------------------------------------
    subroutine init_poly_coeffs(this, c, err)
        ! Arguments
        class(polynomial), intent(inout) :: this
        real(real64), intent(in), dimension(:) :: c
        class(errors), intent(inout), optional, target :: err
 
        ! Local Variables
        integer(int32) :: i, n
        class(errors), pointer :: errmgr
        type(errors), target :: deferr
 
        ! Initialization
        n = size(c)
        if (present(err)) then
            errmgr => err
        else
            errmgr => deferr
        end if
 
        ! Initialize the polynomial
        call init_poly(this, n - 1, errmgr)
        if (errmgr%has_error_occurred()) return
 
        ! Populate the polynomial coefficients
        do i = 1, n
            call this%set(i, c(i))
        end do
    end subroutine
 
! ------------------------------------------------------------------------------
    pure function get_poly_order(this) result(n)
        class(polynomial), intent(in) :: this
        integer(int32) :: n
        if (.not.allocated(this%m_coeffs)) then
            n = -1
        else
            n = size(this%m_coeffs) - 1
        end if
    end function
 
! ------------------------------------------------------------------------------
    subroutine poly_fit(this, x, y, order, err)
        ! Arguments
        class(polynomial), intent(inout) :: this
        real(real64), intent(in), dimension(:) :: x
        real(real64), intent(inout), dimension(:) :: y
        integer(int32), intent(in) :: order
        class(errors), intent(inout), optional, target :: err
 
        ! Parameters
        real(real64), parameter :: one = 1.0d0
 
        ! Local Variables
        integer(int32):: j, n, ncols, flag
        real(real64), pointer, dimension(:,:) :: a
        class(errors), pointer :: errmgr
        type(errors), target :: deferr
 
        ! Initialization
        n = size(x)
        ncols = order + 1
        if (present(err)) then
            errmgr => err
        else
            errmgr => deferr
        end if
 
        ! Input Check
        if (size(y) /= n) then
            ! ERROR: Array size mismatch
            call errmgr%report_error("polynomial_fit", "Array size mismatch.", &
                nl_array_size_error)
            return
        else if (order >= n .or. order < 1) then
            ! ERROR: Requested order does not make sense
            call errmgr%report_error("polynomial_fit", "The requested " // &
                "polynomial order is not valid for this data set.", &
                nl_invalid_input_error)
            return
        end if
 
        ! Local Memory Allocation
        allocate(a(n,ncols), stat = flag)
        if (flag /= 0) then
            ! ERROR: Out of memory
            call errmgr%report_error("polynomial_fit", &
                "Insufficient memory available.", nl_out_of_memory_error)
            return
        end if
 
        ! Ensure the polynomial object is initialized and sized appropriately
        if (this%order() /= order) then
            call this%initialize(order, errmgr)
            if (errmgr%has_error_occurred()) return
        end if
 
        ! Populate A
        do j = 1, n
            a(j,1) = one
            a(j,2) = x(j)
        end do
        do j = 3, ncols
            a(:,j) = a(:,j-1) * x
        end do
 
        ! Solve: A * coeffs = y
        call solve_least_squares(a, y, err = errmgr)
        if (errmgr%has_error_occurred()) return
 
        ! Extract the coefficients from the first order+1 elements of Y
        this%m_coeffs = y(1:ncols)
    end subroutine
 
! ------------------------------------------------------------------------------
    subroutine poly_fit_thru_zero(this, x, y, order, err)
        ! Arguments
        class(polynomial), intent(inout) :: this
        real(real64), intent(in), dimension(:) :: x
        real(real64), intent(inout), dimension(:) :: y
        integer(int32), intent(in) :: order
        class(errors), intent(inout), optional, target :: err
 
        ! Parameters
        real(real64), parameter :: zero = 0.0d0
 
        ! Local Variables
        integer(int32):: j, n, ncols, flag
        real(real64), pointer, dimension(:,:) :: a
        class(errors), pointer :: errmgr
        type(errors), target :: deferr
 
        ! Initialization
        n = size(x)
        ncols = order
        if (present(err)) then
            errmgr => err
        else
            errmgr => deferr
        end if
 
        ! Input Check
        if (size(y) /= n) then
            ! ERROR: Array size mismatch
            call errmgr%report_error("polynomial_fit_thru_zero", &
                "Array size mismatch.", nl_array_size_error)
            return
        else if (order >= n .or. order < 1) then
            ! ERROR: Requested order does not make sense
            call errmgr%report_error("polynomial_fit_thru_zero", &
                "The requested polynomial order is not valid for this " // &
                "data set.", nl_invalid_input_error)
            return
        end if
 
        ! Local Memory Allocation
        allocate(a(n,ncols), stat = flag)
        if (flag /= 0) then
            ! ERROR: Out of memory
            call errmgr%report_error("polynomial_fit_thru_zero", &
                "Insufficient memory available.", nl_out_of_memory_error)
            return
        end if
 
        ! Ensure the polynomial object is initialized and sized appropriately
        if (this%order() /= order) then
            call this%initialize(order, errmgr)
            if (errmgr%has_error_occurred()) return
        end if
 
        ! Populate A
        a(:,1) = x
        do j = 2, ncols
            a(:,j) = a(:,j-1) * x
        end do
 
        ! Solve: A * coeffs = y
        call solve_least_squares(a, y, err = errmgr)
        if (errmgr%has_error_occurred()) return
 
        ! Extract the coefficients from the first order+1 elements of Y
        this%m_coeffs(1) = zero
        this%m_coeffs(2:ncols+1) = y(1:ncols)
    end subroutine
 
! ------------------------------------------------------------------------------
    elemental function poly_eval_double(this, x) result(y)
        ! Arguments
        class(polynomial), intent(in) :: this
        real(real64), intent(in) :: x
        real(real64) :: y
 
        ! Parameters
        real(real64), parameter :: zero = 0.0d0
 
        ! Local Variables
        integer(int32) :: j, order, n
 
        ! Initialization
        order = this%order()
        n = order + 1
        if (order == -1) then
            y = zero
            return
        else if (order == 0) then
            y = this%m_coeffs(1)
            return
        end if
 
        ! Process
        y =this%m_coeffs(n) * x + this%m_coeffs(order)
        do j = n - 2, 1, -1
            y = y * x + this%m_coeffs(j)
        end do
    end function
 
! ------------------------------------------------------------------------------
    elemental function poly_eval_complex(this, x) result(y)
        ! Arguments
        class(polynomial), intent(in) :: this
        complex(real64), intent(in) :: x
        complex(real64) :: y
 
        ! Parameters
        complex(real64), parameter :: zero = (0.0d0, 0.0d0)
 
        ! Local Variables
        integer(int32) :: j, order, n
 
        ! Initialization
        order = this%order()
        n = order + 1
        if (order == -1) then
            y = zero
            return
        else if (order == 0) then
            y = this%m_coeffs(1)
            return
        end if
 
        ! Process
        y =this%m_coeffs(n) * x + this%m_coeffs(order)
        do j = n - 2, 1, -1
            y = y * x + this%m_coeffs(j)
        end do
    end function
 
! ------------------------------------------------------------------------------
    pure function poly_companion_mtx(this) result(c)
        ! Arguments
        class(polynomial), intent(in) :: this
        real(real64), dimension(this%order(), this%order()) :: c
 
        ! Parameters
        real(real64), parameter :: zero = 0.0d0
        real(real64), parameter :: one = 1.0d0
 
        ! Local Variables
        integer(int32) :: i, n
 
        ! Process
        n = this%order()
        if (n == -1) return
        c = zero
        do i = 1, n
            c(i,n) = -this%m_coeffs(i) / this%m_coeffs(n + 1)
            if (i < n) c(i+1,i) = one
        end do
    end function
 
! ------------------------------------------------------------------------------
    function poly_roots(this, err) result(z)
        ! Arguments
        class(polynomial), intent(in) :: this
        complex(real64), dimension(this%order()) :: z
        class(errors), intent(inout), optional, target :: err
 
        ! Local Variables
        integer(int32) :: n
        real(real64), allocatable, dimension(:,:) :: c
 
        ! Initialization
        n = this%order()
 
        ! Quick Return
        if (n == 0) return
 
        ! Compute the companion matrix
        c = this%companion_mtx()
 
        ! Compute the eigenvalues of the companion matrix.  The eigenvalues are
        ! the roots of the polynomial.
        call eigen(c, z, err = err)
    end function
 
! ------------------------------------------------------------------------------
    function get_poly_coefficient(this, ind, err) result(c)
        ! Arguments
        class(polynomial), intent(in) :: this
        integer(int32), intent(in) :: ind
        class(errors), intent(inout), optional, target :: err
        real(real64) :: c
 
        ! Local Variables
        class(errors), pointer :: errmgr
        type(errors), target :: deferr
 
        ! Initialization
        c = 0.0d0
        if (present(err)) then
            errmgr => err
        else
            errmgr => deferr
        end if
 
        ! Quick Return
        if (this%order() == -1) return
 
        ! Input Check
        if (ind <= 0 .or. ind > this%order() + 1) then
            ! ERROR: Index out of range
            call errmgr%report_error("get_polynomial_coefficient", &
                "The specified index is outside the bounds of the " // &
                "coefficient array.", nl_invalid_input_error)
            return
        end if
 
        ! Get the coefficient
        c = this%m_coeffs(ind)
    end function
 
! ------------------------------------------------------------------------------
    pure function get_poly_coefficients(this) result(c)
        ! Arguments
        class(polynomial), intent(in) :: this
        real(real64), dimension(this%order() + 1) :: c
 
        ! Process
        if (this%order() == -1) return
        c = this%m_coeffs
    end function
 
! ------------------------------------------------------------------------------
    subroutine set_poly_coefficient(this, ind, c, err)
        ! Arguments
        class(polynomial), intent(inout) :: this
        integer(int32), intent(in) :: ind
        real(real64), intent(in) :: c
        class(errors), intent(inout), optional, target :: err
 
        ! Local Variables
        class(errors), pointer :: errmgr
        type(errors), target :: deferr
 
        ! Initialization
        if (present(err)) then
            errmgr => err
        else
            errmgr => deferr
        end if
 
        ! Quick Return
        if (this%order() == -1) return
 
        ! Input Check
        if (ind <= 0 .or. ind > this%order() + 1) then
            ! ERROR: Index out of range
            call errmgr%report_error("set_polynomial_coefficient", &
                "The specified index is outside the bounds of the " // &
                "coefficient array.", nl_invalid_input_error)
            return
        end if
 
        ! Process
        this%m_coeffs(ind) = c
    end subroutine
 
! ******************************************************************************
! OPERATORS
! ------------------------------------------------------------------------------
    subroutine poly_equals(x, y)
        ! Arguments
        class(polynomial), intent(inout) :: x
        class(polynomial), intent(in) :: y
 
        ! Local Variables
        integer(int32) :: i, ord
 
        ! Process
        ord = y%order()
        if (x%order() /= ord) call x%initialize(ord)
        do i = 1, ord + 1
            call x%set(i, y%get(i))
        end do
    end subroutine
 
! ------------------------------------------------------------------------------
    subroutine poly_dbl_equals(x, y)
        ! Arguments
        class(polynomial), intent(inout) :: x
        real(real64), intent(in) :: y
 
        ! Local Variables
        integer(int32) :: i, ord
 
        ! Process
        ord = x%order()
        do i = 1, ord + 1
            call x%set(i, y)
        end do
    end subroutine
 
! ------------------------------------------------------------------------------
    subroutine poly_equals_array(x, y)
        ! Arguments
        class(polynomial), intent(inout) :: x
        real(real64), intent(in), dimension(:) :: y
        call x%initialize(y)
    end subroutine
 
! ------------------------------------------------------------------------------
    function poly_poly_add(x, y) result(z)
        ! Arguments
        class(polynomial), intent(in) :: x, y
        type(polynomial) :: z
 
        ! Local Variables
        integer(int32) :: i, max_ord, x_ord, y_ord
 
        ! Initialization
        x_ord = x%order()
        y_ord = y%order()
        max_ord = max(x_ord, y_ord)
        call z%initialize(max_ord)
 
        ! Quick Return
        if (x_ord == -1 .and. y_ord == -1) return
        if (x_ord == -1 .and. y_ord /= -1) then
            do i = 1, max_ord + 1
                call z%set(i, y%get(i))
            end do
            return
        else if (x_ord /= -1 .and. y_ord == -1) then
            do i = 1, max_ord + 1
                call z%set(i, x%get(i))
            end do
            return
        end if
 
        ! Process
        if (x_ord > y_ord) then
            do i = 1, y_ord + 1
                call z%set(i, x%get(i) + y%get(i))
            end do
            do i = y_ord + 2, x_ord
                call z%set(i, x%get(i))
            end do
        else if (x_ord < y_ord) then
            do i = 1, x_ord + 1
                call z%set(i, x%get(i) + y%get(i))
            end do
            do i = x_ord + 2, y_ord + 1
                call z%set(i, y%get(i))
            end do
        else
            do i = 1, max_ord + 1
                call z%set(i, x%get(i) + y%get(i))
            end do
        end if
    end function
 
! ------------------------------------------------------------------------------
    function poly_poly_subtract(x, y) result(z)
        ! Arguments
        class(polynomial), intent(in) :: x, y
        type(polynomial) :: z
 
        ! Local Variables
        integer(int32) :: i, max_ord, x_ord, y_ord
 
        ! Initialization
        x_ord = x%order()
        y_ord = y%order()
        max_ord = max(x_ord, y_ord)
        call z%initialize(max_ord)
 
        ! Quick Return
        if (x_ord == -1 .and. y_ord == -1) return
        if (x_ord == -1 .and. y_ord /= -1) then
            do i = 1, max_ord + 1
                call z%set(i, y%get(i))
            end do
            return
        else if (x_ord /= -1 .and. y_ord == -1) then
            do i = 1, max_ord + 1
                call z%set(i, x%get(i))
            end do
            return
        end if
 
        ! Process
        if (x_ord > y_ord) then
            do i = 1, y_ord + 1
                call z%set(i, x%get(i) - y%get(i))
            end do
            do i = y_ord + 2, x_ord
                call z%set(i, x%get(i))
            end do
        else if (x_ord < y_ord) then
            do i = 1, x_ord + 1
                call z%set(i, x%get(i) - y%get(i))
            end do
            do i = x_ord + 2, y_ord + 1
                call z%set(i, -y%get(i))
            end do
        else
            do i = 1, max_ord + 1
                call z%set(i, x%get(i) - y%get(i))
            end do
        end if
    end function
 
! ------------------------------------------------------------------------------
    function poly_poly_mult(x, y) result(z)
        ! Arguments
        class(polynomial), intent(in) :: x, y
        type(polynomial) :: z
 
        ! Local Variables
        integer(int32) :: i, j, m, n
        real(real64) :: val
 
        ! Initialization
        n = x%order() + 1
        m = y%order() + 1
        call z%initialize(x%order() + y%order()) ! Sets z to all zeros
 
        ! Process
        do i = 1, n
            do j = 1, m
                val = z%get(i + j - 1) + x%get(i) * y%get(j)
                call z%set(i + j - 1, val)
            end do
        end do
    end function
 
! ------------------------------------------------------------------------------
    function poly_dbl_mult(x, y) result(z)
        ! Arguments
        class(polynomial), intent(in) :: x
        real(real64), intent(in) :: y
        type(polynomial) :: z
 
        ! Local Variables
        integer(int32) :: i, ord
 
        ! Process
        ord = x%order()
        call z%initialize(ord)
        do i = 1, ord + 1
            call z%set(i, x%get(i) * y)
        end do
    end function
 
! ------------------------------------------------------------------------------
    function dbl_poly_mult(x, y) result(z)
        ! Arguments
        real(real64), intent(in) :: x
        class(polynomial), intent(in) :: y
        type(polynomial) :: z
 
        ! Local Variables
        integer(int32) :: i, ord
 
        ! Process
        ord = y%order()
        call z%initialize(ord)
        do i = 1, ord + 1
            call z%set(i, y%get(i) * x)
        end do
    end function
 
! ------------------------------------------------------------------------------
 
! Example Polynomial Code (Coefficients go from lowest order to highest)
! src: https://github.com/JuliaMath/Polynomials.jl
!
! function *{T,S}(p1::Poly{T}, p2::Poly{S})
!     if p1.var != p2.var
!         error("Polynomials must have same variable")
!     end
!     R = promote_type(T,S)
!     n = length(p1)-1
!     m = length(p2)-1
!     a = zeros(R,m+n+1)
 
!     for i = 0:n
!         for j = 0:m
!             a[i+j+1] += p1[i] * p2[j]
!         end
!     end
!     Poly(a,p1.var)
! end
 
! ## older . operators, hack to avoid warning on v0.6
! dot_operators = quote
!     @compat Base.:.+{T<:Number}(c::T, p::Poly) = +(p, c)
!     @compat Base.:.+{T<:Number}(p::Poly, c::T) = +(p, c)
!     @compat Base.:.-{T<:Number}(p::Poly, c::T) = +(p, -c)
!     @compat Base.:.-{T<:Number}(c::T, p::Poly) = +(p, -c)
!     @compat Base.:.*{T<:Number,S}(c::T, p::Poly{S}) = Poly(c * p.a, p.var)
!     @compat Base.:.*{T<:Number,S}(p::Poly{S}, c::T) = Poly(p.a * c, p.var)
! end
! VERSION < v"0.6.0-dev" && eval(dot_operators)
 
 
! # are any values NaN
! hasnan(p::Poly) = reduce(|, (@compat isnan.(p.a)))
 
! function divrem{T, S}(num::Poly{T}, den::Poly{S})
!     if num.var != den.var
!         error("Polynomials must have same variable")
!     end
!     m = length(den)-1
!     if m == 0 && den[0] == 0
!         throw(DivideError())
!     end
!     R = typeof(one(T)/one(S))
!     n = length(num)-1
!     deg = n-m+1
!     if deg <= 0
!         return convert(Poly{R}, zero(num)), convert(Poly{R}, num)
!     end
 
!     aQ = zeros(R, deg)
!     # aR = deepcopy(num.a)
!     # @show num.a
!     aR = R[ num.a[i] for i = 1:n+1 ]
!     for i = n:-1:m
!         quot = aR[i+1] / den[m]
!         aQ[i-m+1] = quot
!         for j = 0:m
!             elem = den[j]*quot
!             aR[i-(m-j)+1] -= elem
!         end
!     end
!     pQ = Poly(aQ, num.var)
!     pR = Poly(aR, num.var)
 
!     return pQ, pR
! end
 
! div(num::Poly, den::Poly) = divrem(num, den)[1]
! rem(num::Poly, den::Poly) = divrem(num, den)[2]
 
! ==(p1::Poly, p2::Poly) = (p1.var == p2.var && p1.a == p2.a)
! ==(p1::Poly, n::Number) = (coeffs(p1) == [n])
! ==(n::Number, p1::Poly) = (p1 == n)
end module