subroutine fpintb(t,n,bint,nk1,x,y)
c  subroutine fpintb calculates integrals of the normalized b-splines
c  nj,k+1(x) of degree k, defined on the set of knots t(j),j=1,2,...n.
c  it makes use of the formulae of gaffney for the calculation of
c  indefinite integrals of b-splines.
c
c  calling sequence:
c     call fpintb(t,n,bint,nk1,x,y)
c
c  input parameters:
c    t    : real array,length n, containing the position of the knots.
c    n    : integer value, giving the number of knots.
c    nk1  : integer value, giving the number of b-splines of degree k,
c           defined on the set of knots ,i.e. nk1 = n-k-1.
c    x,y  : real values, containing the end points of the integration
c           interval.
c  output parameter:
c    bint : array,length nk1, containing the integrals of the b-splines.
c  ..
c  ..scalars arguments..
      integer n,nk1
      real*8 x,y
c  ..array arguments..
      real*8 t(n),bint(nk1)
c  ..local scalars..
      integer i,ia,ib,it,j,j1,k,k1,l,li,lj,lk,l0,min
      real*8 a,ak,arg,b,f,one
c  ..local arrays..
      real*8 aint(6),h(6),h1(6)
c  initialization.
      one = 0.1d+01
      k1 = n-nk1
      ak = k1
      k = k1-1
      do 10 i=1,nk1
        bint(i) = 0.0d0
  10  continue
c  the integration limits are arranged in increasing order.
      a = x
      b = y
      min = 0
      if (a.lt.b) go to 30
      if (a.eq.b) go to 160
      go to 20
  20  a = y
      b = x
      min = 1
  30  if(a.lt.t(k1)) a = t(k1)
      if(b.gt.t(nk1+1)) b = t(nk1+1)
c  using the expression of gaffney for the indefinite integral of a
c  b-spline we find that
c  bint(j) = (t(j+k+1)-t(j))*(res(j,b)-res(j,a))/(k+1)
c    where for t(l) <= x < t(l+1)
c    res(j,x) = 0, j=1,2,...,l-k-1
c             = 1, j=l+1,l+2,...,nk1
c             = aint(j+k-l+1), j=l-k,l-k+1,...,l
c               = sumi((x-t(j+i))*nj+i,k+1-i(x)/(t(j+k+1)-t(j+i)))
c                 i=0,1,...,k
      l = k1
      l0 = l+1
c  set arg = a.
      arg = a
      do 90 it=1,2
c  search for the knot interval t(l) <= arg < t(l+1).
  40    if(arg.lt.t(l0) .or. l.eq.nk1) go to 50
        l = l0
        l0 = l+1
        go to 40
c  calculation of aint(j), j=1,2,...,k+1.
c  initialization.
  50    do 55 j=1,k1
          aint(j) = 0.0d0
  55    continue
        aint(1) = (arg-t(l))/(t(l+1)-t(l))
        h1(1) = one
        do 70 j=1,k
c  evaluation of the non-zero b-splines of degree j at arg,i.e.
c    h(i+1) = nl-j+i,j(arg), i=0,1,...,j.
          h(1) = 0.0d0
          do 60 i=1,j
            li = l+i
            lj = li-j
            f = h1(i)/(t(li)-t(lj))
            h(i) = h(i)+f*(t(li)-arg)
            h(i+1) = f*(arg-t(lj))
  60      continue
c  updating of the integrals aint.
          j1 = j+1
          do 70 i=1,j1
            li = l+i
            lj = li-j1
            aint(i) = aint(i)+h(i)*(arg-t(lj))/(t(li)-t(lj))
            h1(i) = h(i)
  70    continue
        if(it.eq.2) go to 100
c  updating of the integrals bint
        lk = l-k
        ia = lk
        do 80 i=1,k1
          bint(lk) = -aint(i)
          lk = lk+1
  80    continue
c  set arg = b.
        arg = b
  90  continue
c  updating of the integrals bint.
 100  lk = l-k
      ib = lk-1
      do 110 i=1,k1
        bint(lk) = bint(lk)+aint(i)
        lk = lk+1
 110  continue
      if(ib.lt.ia) go to 130
      do 120 i=ia,ib
        bint(i) = bint(i)+one
 120  continue
c  the scaling factors are taken into account.
 130  f = one/ak
      do 140 i=1,nk1
        j = i+k1
        bint(i) = bint(i)*(t(j)-t(i))*f
 140  continue
c  the order of the integration limits is taken into account.
      if(min.eq.0) go to 160
      do 150 i=1,nk1
        bint(i) = -bint(i)
 150  continue
 160  return
      end