subroutine bsplpp ( t, bcoef, n, k, scrtch, break, coef, l )
c  from  * a practical guide to splines *  by c. de Boor (7 may 92)
calls  bsplvb
c
converts the b-representation  t, bcoef, n, k  of some spline into its
c  pp-representation  break, coef, l, k .
c
c******  i n p u t  ******
c  t.....knot sequence, of length  n+k
c  bcoef.....b-spline coefficient sequence, of length  n
c  n.....length of  bcoef  and  dimension of spline space  spline(k,t)
c  k.....order of the spline
c
c  w a r n i n g  . . .  the restriction   k .le. kmax (= 20)   is impo-
c        sed by the arbitrary dimension statement for  biatx  below, but
c        is  n o w h e r e   c h e c k e d   for.
c
c******  w o r k   a r e a  ******
c  scrtch......of size  (k,k) , needed to contain bcoeffs of a piece of
c        the spline and its  k-1  derivatives
c
c******  o u t p u t  ******
c  break.....breakpoint sequence, of length  l+1, contains (in increas-
c        ing order) the distinct points in the sequence  t(k),...,t(n+1)
c  coef.....array of size (k,l), with  coef(i,j) = (i-1)st derivative of
c        spline at break(j) from the right
c  l.....number of polynomial pieces which make up the spline in the in-
c        terval  (t(k), t(n+1))
c
c******  m e t h o d  ******
c     for each breakpoint interval, the  k  relevant b-coeffs of the
c  spline are found and then differenced repeatedly to get the b-coeffs
c  of all the derivatives of the spline on that interval. the spline and
c  its first  k-1  derivatives are then evaluated at the left end point
c  of that interval, using  bsplvb  repeatedly to obtain the values of
c  all b-splines of the appropriate order at that point.
c
      integer k,l,n,   i,j,jp1,kmax,kmj,left,lsofar
      parameter (kmax = 20)
      real bcoef(n),break(l+1),coef(k,l),t(n+k),   scrtch(k,k)
     *                                      ,biatx(kmax),diff,factor,sum
c
      lsofar = 0
      break(1) = t(k)
      do 50 left=k,n
c                                find the next nontrivial knot interval.
         if (t(left+1) .eq. t(left))    go to 50
         lsofar = lsofar + 1
         break(lsofar+1) = t(left+1)
         if (k .gt. 1)                  go to 9
         coef(1,lsofar) = bcoef(left)
                                        go to 50
c        store the k b-spline coeff.s relevant to current knot interval
c                             in  scrtch(.,1) .
    9    do 10 i=1,k
   10       scrtch(i,1) = bcoef(left-k+i)
c
c        for j=1,...,k-1, compute the  k-j  b-spline coeff.s relevant to
c        current knot interval for the j-th derivative by differencing
c        those for the (j-1)st derivative, and store in scrtch(.,j+1) .
         do 20 jp1=2,k
            j = jp1 - 1
            kmj = k - j
            do 20 i=1,kmj
               diff = t(left+i) - t(left+i - kmj)
               if (diff .gt. 0.)  scrtch(i,jp1) =
     *                       (scrtch(i+1,j)-scrtch(i,j))/diff
   20          continue
c
c        for  j = 0, ..., k-1, find the values at  t(left)  of the  j+1
c        b-splines of order  j+1  whose support contains the current
c        knot interval from those of order  j  (in  biatx ), then comb-
c        ine with the b-spline coeff.s (in scrtch(.,k-j) ) found earlier
c        to compute the (k-j-1)st derivative at  t(left)  of the given
c        spline.
c           note. if the repeated calls to  bsplvb  are thought to gene-
c        rate too much overhead, then replace the first call by
c           biatx(1) = 1.
c        and the subsequent call by the statement
c           j = jp1 - 1
c        followed by a direct copy of the lines
c           deltar(j) = t(left+j) - x
c                  ......
c           biatx(j+1) = saved
c        from  bsplvb . deltal(kmax)  and  deltar(kmax)  would have to
c        appear in a dimension statement, of course.
c
         call bsplvb ( t, 1, 1, t(left), left, biatx )
         coef(k,lsofar) = scrtch(1,k)
         do 30 jp1=2,k
            call bsplvb ( t, jp1, 2, t(left), left, biatx )
            kmj = k+1 - jp1
            sum = 0.
            do 28 i=1,jp1
   28          sum = biatx(i)*scrtch(i,kmj) + sum
   30       coef(kmj,lsofar) = sum
   50    continue
      l = lsofar
	  if (k .eq. 1)                     return
	  factor = 1.
	  do 60 i=2,k
		 factor = factor*float(k+1-i)
		 do 60 j=1,lsofar
   60       coef(i,j) = coef(i,j)*factor
                                        return
      end