spline_roots Subroutine

    pure subroutine spline_roots(t,c,n,k,a,roots,iflag)

    implicit none

    integer(ip),intent(in)                        :: n      !! number of b-spline coefficients
    integer(ip),intent(in)                        :: k      !! order of the b-spline, `k >= 1`
    real(wp),dimension(n+k),intent(in)            :: t      !! knot array
    real(wp),dimension(n),intent(in)              :: c      !! b-spline coefficient array
    real(wp),intent(in)                           :: a      !! the value to solve for
    real(wp),dimension(:),allocatable,intent(out) :: roots  !! the roots, in increasing order
    integer(ip),intent(out)                       :: iflag  !! 0, or 4004 if \( s(x) = a \) on a whole interval

    real(wp),dimension(0:k-1) :: p   !! power basis coefficients of a piece, on \( u \in [0,1] \)
                                     !! (for the closed form roots of degree 1 and 2)
    real(wp),dimension(0:k-1) :: b   !! Bernstein coefficients of a piece minus `a`, on \( u \in [0,1] \)
    real(wp),dimension(k) :: cl      !! local b-spline coefficients minus `a`
    real(wp),dimension(3*k) :: work  !! work array for [[dbvalu]]
    real(wp),dimension(:),allocatable :: u  !! roots of a piece, on \( u \in [0,1] \)
    real(wp),dimension(:),allocatable :: xr !! roots of a piece, on \( x \in [t_l, t_{l+1}] \)
    logical,dimension(:),allocatable :: keep  !! roots to keep
    integer(ip) :: l     !! knot interval index: `t(l) <= x < t(l+1)`
    integer(ip) :: i     !! counter
    integer(ip) :: m     !! polynomial degree of the pieces: `k-1`
    integer(ip) :: inbv  !! initialization parameter for [[dbvalu]]
    integer(ip) :: ierr  !! status flag from [[dbvalu]]
    real(wp) :: h        !! width of the knot interval: `t(l+1) - t(l)`
    real(wp) :: val      !! spline value from [[dbvalu]]
    real(wp) :: ftol     !! absolute tolerance on \( |s(x)-a| \)

    real(wp),parameter :: utol = 16.0_wp*epsilon(1.0_wp)  !! accept roots just outside the interval
                                                           !! \( [0,1] \) (they are moved to the end)

    iflag = 0_ip
    allocate(roots(0))
    m = k - 1_ip  ! polynomial degree
    ftol = root_ftol * max(abs(a),maxval(abs(c)))
    inbv = 1_ip

    ! end points (a root here may not be a sign change after rounding)
    call dbvalu(t,c,n,k,0_ip,t(k),inbv,work,ierr,val)
    if (ierr==0_ip .and. abs(val-a)<=ftol) roots = [roots, t(k)]
    call dbvalu(t,c,n,k,0_ip,t(n+1_ip),inbv,work,ierr,val)
    if (ierr==0_ip .and. abs(val-a)<=ftol) roots = [roots, t(n+1_ip)]

    do l = k, n
        if (t(l)>=t(l+1_ip)) cycle  ! empty interval
        cl = c(l-k+1_ip:l) - a

        ! s(x) = a on the whole interval. Its left end is a root (it may
        ! not be found from the previous interval, where it may be a
        ! root of high multiplicity). The right end is found from the
        ! next interval (or is the end of the domain). The roots inside
        ! a flat region are removed below.
        if (maxval(abs(cl))<=ftol) then
            iflag = 4004_ip
            roots = [roots, t(l)]
            cycle
        end if

        ! convex hull: no root in this interval
        if (all(cl>ftol) .or. all(cl<-ftol)) cycle

        ! Bernstein form of the piece on u in [0,1], with x = t(l) + h*u
        h = t(l+1_ip) - t(l)
        b = bernstein_coefficients()
        ! root at the knot (the pieces on either side are computed
        ! separately, so it may not be a sign change after rounding)
        if (abs(b(0))<=ftol) roots = [roots, t(l)]

        select case (m)
        case(1_ip)
            p(0:1) = [b(0), b(1)-b(0)]
            u = linear_root()
        case(2_ip)
            p(0:2) = [b(0), 2.0_wp*(b(1)-b(0)), b(0)-2.0_wp*b(1)+b(2)]
            u = quadratic_roots()
        case(3_ip:)
            u = bernstein_roots()
        case default
            u = [real(wp) ::]  ! constant piece
        end select

        ! keep the roots in the interval (allowing for rounding at the ends)
        u = pack(u, u>=-utol .and. u<=1.0_wp+utol)
        xr = min(t(l+1_ip), t(l) + h*min(1.0_wp,max(0.0_wp,u)))
        if (t(l+1_ip)<t(n+1_ip)) then
            ! a root at the right end is the limit from the left. At an
            ! interior knot, s(x) is the value from the right, which is
            ! different if the spline is discontinuous there. A root at
            ! the knot is found from the next nonempty interval instead.
            xr = pack(xr, xr<t(l+1_ip))
        end if
        roots = [roots, xr]
    end do

    call sort_unique(roots, root_dup_tol(t(k),t(n+1_ip)))

    ! remove the roots inside regions where s(x) = a (only the ends are kept)
    if (iflag==4004_ip) then
        allocate(keep(size(roots)))
        do i = 1_ip, size(roots,kind=ip)
            keep(i) = .not. in_flat_region(t,c,n,k,a,ftol,roots(i))
        end do
        roots = pack(roots, keep)
    end if

    contains

        pure function bernstein_coefficients() result(bc)

        !! Bernstein coefficients of the piece on knot interval `l`, minus `a`.
        !!
        !! The `i`-th coefficient is the blossom of the piece evaluated at
        !! `m-i` copies of `t(l)` and `i` copies of `t(l+1)`, computed with
        !! de Boor's algorithm using these as the arguments at each level.
        !! Every step is a convex combination, so the coefficients are
        !! accurate to about `eps*max(|c|)` for any order (unlike the power
        !! basis computed from the derivatives, which loses accuracy at the
        !! right end of the interval for high orders).

        real(wp),dimension(0:m) :: bc  !! the Bernstein coefficients

        real(wp),dimension(0:m) :: w  !! de Boor triangle (`w(i)` is for coefficient `l-m+i`)
        real(wp) :: x       !! the argument at the current level
        real(wp) :: alpha   !! interpolation weight
        integer(ip) :: ib   !! index of the Bernstein coefficient
        integer(ip) :: r    !! level of the de Boor triangle
        integer(ip) :: jj   !! counter
        integer(ip) :: jc   !! b-spline coefficient index: `l-m+jj`

        do ib = 0_ip, m
            w = cl
            do r = 1_ip, m
                x = merge(t(l), t(l+1_ip), r<=m-ib)
                do jj = m, r, -1_ip
                    jc = l - m + jj
                    alpha = (x - t(jc)) / (t(jc+m+1_ip-r) - t(jc))
                    w(jj) = alpha*w(jj) + (1.0_wp-alpha)*w(jj-1_ip)
                end do
            end do
            bc(ib) = w(m)
        end do

        end function bernstein_coefficients

        pure function linear_root() result(ur)

        !! Root of a linear piece \( p_0 + p_1 u \).

        real(wp),dimension(:),allocatable :: ur  !! the root, on \( u \) (size 0 if \( p_1 = 0 \))

        if (p(1)/=0.0_wp) then
            ur = [-p(0)/p(1)]
        else
            allocate(ur(0))
        end if

        end function linear_root

        pure function quadratic_roots() result(ur)

        !! Roots of a quadratic piece \( p_0 + p_1 u + p_2 u^2 \)
        !! (numerically stable form of the quadratic formula).

        real(wp),dimension(:),allocatable :: ur  !! the real roots, on \( u \) (size 0, 1, or 2)

        real(wp),dimension(0:2) :: ps  !! the coefficients, scaled so the largest
                                       !! magnitude is 1 (so that the discriminant
                                       !! does not overflow or underflow)
        real(wp) :: disc  !! discriminant: \( p_1^2 - 4 p_2 p_0 \) (scaled)
        real(wp) :: q     !! \( -\frac{1}{2} \left( p_1 + \mathrm{sign}(p_1) \sqrt{disc} \right) \)
                          !! (the roots are \( q/p_2 \) and \( p_0/q \)) (scaled)

        if (p(2)==0.0_wp) then
            ur = linear_root()
            return
        end if
        allocate(ur(0))
        ps = p(0:2) / maxval(abs(p(0:2)))
        disc = ps(1)**2 - 4.0_wp*ps(2)*ps(0)
        if (disc<0.0_wp) return
        q = -0.5_wp * (ps(1) + sign(sqrt(disc),ps(1)))
        ur = [q/ps(2)]
        if (q/=0.0_wp) ur = [ur, ps(0)/q]

        end function quadratic_roots

        pure function bernstein_roots() result(ur)

        !! Isolate the roots of a piece by Bernstein subdivision,
        !! and refine each with a safeguarded Newton iteration.
        !!
        !! Subintervals are processed from a stack. The left half of a
        !! divided subinterval is processed first, so at most one right
        !! half per depth is waiting, and the stack never holds more
        !! than `max_depth+2` entries.

        real(wp),dimension(:),allocatable :: ur  !! the roots, on \( u \in [0,1] \)

        integer(ip),parameter :: max_depth = 52_ip  !! maximum number of subdivisions
                                                    !! (the interval width is then \( 2^{-52} \))

        real(wp),dimension(0:m,max_depth+2) :: bs   !! stack of Bernstein coefficients
        real(wp),dimension(max_depth+2) :: uas      !! stack of left ends of the subintervals
        real(wp),dimension(max_depth+2) :: ubs      !! stack of right ends of the subintervals
        integer(ip),dimension(max_depth+2) :: deps  !! stack of subdivision depths
        real(wp),dimension(0:m) :: bb     !! Bernstein coefficients on the current subinterval
        real(wp),dimension(0:m) :: left   !! Bernstein coefficients on the left half
        real(wp),dimension(0:m) :: right  !! Bernstein coefficients on the right half
        real(wp),dimension(0:m) :: tmp    !! work array for de Casteljau's algorithm
        integer(ip) :: top    !! top of the stack (0 if empty)
        integer(ip) :: nsign  !! number of sign changes in `bb` (an upper bound
                              !! on the number of roots in the subinterval)
        integer(ip) :: dep    !! subdivision depth of the current subinterval
        integer(ip) :: r      !! counter
        integer(ip) :: s      !! counter
        real(wp) :: ua        !! left end of the current subinterval
        real(wp) :: ub        !! right end of the current subinterval
        real(wp) :: um        !! midpoint of the current subinterval
        real(wp) :: sprev     !! sign of the previous nonzero coefficient (0 if none yet)

        allocate(ur(0))
        top = 1_ip
        bs(:,1) = b
        uas(1) = 0.0_wp
        ubs(1) = 1.0_wp
        deps(1) = 0_ip

        do while (top>0_ip)
            bb = bs(:,top)
            ua = uas(top)
            ub = ubs(top)
            dep = deps(top)
            top = top - 1_ip

            ! the piece is within the noise level of zero on this interval
            ! (a multiple root): keep it if the sign changes across the
            ! interval (odd multiplicity), otherwise it is found by db1root
            if (maxval(abs(bb))<=ftol) then
                if (bb(0)*bb(m)<0.0_wp) ur = [ur, 0.5_wp*(ua+ub)]
                cycle
            end if

            ! number of sign changes (Descartes' rule of signs),
            ! where coefficients within the noise level count as zero
            nsign = 0_ip
            sprev = 0.0_wp
            do r = 0_ip, m
                if (abs(bb(r))>ftol) then
                    if (sprev/=0.0_wp .and. sign(1.0_wp,bb(r))/=sprev) nsign = nsign + 1_ip
                    sprev = sign(1.0_wp,bb(r))
                end if
            end do

            if (nsign==0_ip) cycle
            if (nsign==1_ip .and. abs(bb(0))>ftol .and. abs(bb(m))>ftol) then
                ! exactly one simple root (and none at the ends)
                ur = [ur, newton(ua,ub)]
                cycle
            end if
            if (dep>=max_depth) then
                ur = [ur, 0.5_wp*(ua+ub)]  ! cluster of roots
                cycle
            end if

            ! subdivide at the midpoint (de Casteljau)
            tmp = bb
            do r = 1_ip, m
                left(r-1_ip) = tmp(0)
                right(m-r+1_ip) = tmp(m-r+1_ip)
                do s = 0_ip, m-r
                    tmp(s) = 0.5_wp*(tmp(s)+tmp(s+1_ip))
                end do
            end do
            left(m) = tmp(0)
            right(0) = tmp(0)
            um = 0.5_wp*(ua+ub)
            if (abs(tmp(0))<=ftol) ur = [ur, um]  ! root at the midpoint

            top = top + 1_ip
            bs(:,top) = right
            uas(top) = um
            ubs(top) = ub
            deps(top) = dep + 1_ip
            top = top + 1_ip
            bs(:,top) = left
            uas(top) = ua
            ubs(top) = um
            deps(top) = dep + 1_ip
        end do

        end function bernstein_roots

        pure function newton(ua,ub) result(u)

        !! Safeguarded Newton iteration for the single root of the
        !! piece in `[ua,ub]` (bisection when the Newton step leaves
        !! the bracket, or does not reduce it quickly enough).

        real(wp),intent(in) :: ua  !! left end of the bracket
        real(wp),intent(in) :: ub  !! right end of the bracket
        real(wp) :: u              !! the root, on \( u \in [u_a,u_b] \)

        integer(ip),parameter :: max_iter = 200_ip  !! maximum number of iterations

        integer(ip) :: iter  !! iteration counter
        real(wp) :: xl       !! end of the bracket where the piece is negative
        real(wp) :: xh       !! end of the bracket where the piece is positive
        real(wp) :: f        !! value of the piece at `u`
        real(wp) :: df       !! derivative of the piece at `u`
        real(wp) :: fl       !! value of the piece at `ua`
        real(wp) :: fh       !! value of the piece at `ub`
        real(wp) :: dx       !! the last step
        real(wp) :: dxold    !! the step before the last one
        real(wp) :: temp     !! previous value of `u`

        call bernstein_eval(ua,fl,df)
        call bernstein_eval(ub,fh,df)
        if (fl==0.0_wp) then
            u = ua; return
        else if (fh==0.0_wp) then
            u = ub; return
        else if (fl*fh>0.0_wp) then
            ! no sign change after rounding: take the closest end
            u = merge(ua,ub,abs(fl)<abs(fh)); return
        end if
        if (fl<0.0_wp) then
            xl = ua; xh = ub
        else
            xl = ub; xh = ua
        end if

        u = 0.5_wp*(ua+ub)
        dxold = abs(ub-ua)
        dx = dxold
        call bernstein_eval(u,f,df)
        do iter = 1_ip, max_iter
            if (((u-xh)*df-f)*((u-xl)*df-f)>0.0_wp .or. abs(2.0_wp*f)>abs(dxold*df)) then
                dxold = dx
                dx = 0.5_wp*(xh-xl)
                u = xl + dx
                if (xl==u) exit
            else
                dxold = dx
                dx = f/df
                temp = u
                u = u - dx
                if (temp==u) exit
            end if
            if (abs(dx)*h<=2.0_wp*epsilon(1.0_wp)*abs(t(l)+h*u)) exit
            call bernstein_eval(u,f,df)
            if (f==0.0_wp) exit
            if (f<0.0_wp) then
                xl = u
            else
                xh = u
            end if
        end do

        end function newton

        pure subroutine bernstein_eval(u,f,df)

        !! Value and derivative of the piece at `u`
        !! (from the Bernstein coefficients, by de Casteljau's algorithm).

        real(wp),intent(in)  :: u   !! the point, on \( u \in [0,1] \)
        real(wp),intent(out) :: f   !! value of the piece
        real(wp),intent(out) :: df  !! derivative of the piece with respect to `u`

        real(wp),dimension(0:m) :: w  !! de Casteljau triangle
        integer(ip) :: r   !! level of the triangle
        integer(ip) :: jj  !! counter

        w = b
        do r = 1_ip, m-1_ip
            do jj = 0_ip, m-r
                w(jj) = (1.0_wp-u)*w(jj) + u*w(jj+1_ip)
            end do
        end do
        ! the last level gives the value, and its two inputs the derivative
        df = real(m,wp)*(w(1)-w(0))
        f = (1.0_wp-u)*w(0) + u*w(1)

        end subroutine bernstein_eval

    end subroutine spline_roots