s_bsplines_uniform__eval_basis_and_n_derivs Subroutine

  SLL_PURE subroutine s_bsplines_uniform__eval_basis_and_n_derivs( self, x, n, derivs, jmin )
    class(sll_t_bsplines_uniform), intent(in   ) :: self
    real(wp)                     , intent(in   ) :: x
    integer                      , intent(in   ) :: n
    real(wp)                     , intent(  out) :: derivs(0:,0:)
    integer                      , intent(  out) :: jmin

    integer  :: icell
    real(wp) :: x_offset

    integer  :: j, r
    real(wp) :: j_real
    real(wp) :: xx
    real(wp) :: temp
    real(wp) :: saved

    integer  :: k, s1, s2, rk, pk, j1, j2
    real(wp) :: d

    ! GFortran: to allocate on stack use -fstack-arrays
    real(wp) :: ndu (0:self%degree, 0:self%degree)
    real(wp) :: a   (0:1          , 0:self%degree)

    ! Check on inputs
    SLL_ASSERT( n >= 0 )
    SLL_ASSERT( size(derivs,1) == 1+n           )
    SLL_ASSERT( size(derivs,2) == 1+self%degree )

    ! 1. Compute cell index 'icell' and x_offset
    call s_bsplines_uniform__get_icell_and_offset( self, x, icell, x_offset )

    ! 2. Compute index range of B-splines with support over cell 'icell'
    jmin = icell - self%offset

    ! 3. Recursively evaluate B-splines (see "sll_s_uniform_bsplines_eval_basis")
    !    up to self%degree, and store them all in the upper-right triangle of ndu
    associate( spline_degree => self%degree, bspl => derivs )

      ndu(0,0) = 1.0_wp
      do j = 1, spline_degree
         j_real = real(j,wp)
         xx     = -x_offset
         saved  = 0.0_wp
         do r = 0, j-1
            xx       = xx + 1.0_wp
            temp     = ndu(r,j-1) * inv(j)
            ndu(r,j) = saved + xx * temp
            saved    = (j_real - xx) * temp
         end do
         ndu(j,j) = saved
      end do
      bspl(0,:) = ndu(:,spline_degree)

    end associate  ! spline_degree, bspl

    ! 4. Use equation 2.10 in "The NURBS Book" to compute n derivatives
    associate( deg => self%degree, bsdx => derivs )

      do r = 0, deg
         s1 = 0
         s2 = 1
         a(0,0) = 1.0_wp
         do k = 1, n
            d  = 0.0_wp
            rk = r-k
            pk = deg-k
            if (r >= k) then
               a(s2,0) = a(s1,0) * inv(pk+1)
               d = a(s2,0) * ndu(rk,pk)
            end if
            if (rk > -1) then
               j1 = 1
            else
               j1 = -rk
            end if
            if (r-1 <= pk) then
               j2 = k-1
            else
               j2 = deg-r
            end if
            do j = j1, j2
               a(s2,j) = (a(s1,j) - a(s1,j-1)) * inv(pk+1)
               d = d + a(s2,j) * ndu(rk+j,pk)
            end do
            if (r <= pk) then
               a(s2,k) = - a(s1,k-1) * inv(pk+1)
               d = d + a(s2,k) * ndu(r,pk)
            end if
            bsdx(k,r) = d
            j  = s1
            s1 = s2
            s2 = j
         end do
      end do

      ! Multiply result by correct factors:
      ! deg!/(deg-n)! = deg*(deg-1)*...*(deg-n+1)
      ! k-th derivatives are normalized, hence they should be divided by dx^k
      d = real(deg,wp) * self%inv_dx
      do k = 1, n
         bsdx(k,:) = bsdx(k,:) * d
         d = d * real(deg-k,wp) * self%inv_dx
      end do

    end associate  ! deg, bsdx

  end subroutine s_bsplines_uniform__eval_basis_and_n_derivs