function intgrf=fnint(f)
%FNINT	Indefinite integral of a function represented by splines.
%        intgrf=fnint(f)
%
%  returns the (representation of the) indefinite integral of the function 
%  contained in  f  (and in the same form). 

% C. de Boor /  Feb.25, 1989
% C. de Boor /  Mar.21, 1991 (correct misprint in line 15)
% C. de Boor /  Sep.23, 1991 (add more knots at right end for B-form)
% C. de Boor /   6 Feb  1992  (finish the last correction)
% Copyright (c) 1990-92 by Carl de Boor and The MathWorks, Inc.

if (f(1)==10), % the function is in pp-representation:
   [breaks,coefs,l,k,d]=ppbrk(f);
   coefs=coefs./(ones(d*l,1)*[k:-1:1]);
   if (l<2),
      intgrf=ppmak(breaks,[coefs zeros(d,1)],d);
   else
     % evaluate each integrated polynomial at the right endpoint of its 
	 % interval (adapted from  ppual )
      xs=diff(breaks(1:l));index=[1:l-1];
      if (d>1),
         ones(d,1)*xs; xs=ans(:)';
         1+d*ones(d,1)*index+[-d:-1]'*ones(1,l-1); index=ans(:);
      end
      vv=xs.*coefs(index,1)';
      for i=2:k,
         vv = xs.*(vv + coefs(index,i)');
      end
      if (d>1),
         junk=zeros(d,l-1);junk(:)=vv;last=(cumsum([zeros(1,d);junk']))';
      else,
         last=cumsum([0,vv]);
      end

      intgrf=ppmak(breaks,[coefs(:,1:k) last(:)],d);
   end

elseif (f(1)==11), % the function is in B-repr.
   % Set it up so that it would be correct on the interval [t(1),t(n+k)].
   % There is no way to make it correct everywhere since the integral of
   % a spline over the interval [t(1),t(n+k)] need not be zero.
   [t,a,n,k,d]=spbrk(f); 

   index = find(diff(t)>0); % increase multiplicity of last knot to  k
   needed = index(length(index)) - n; % =  k+1 - (n+k - index(length(index));
   if (needed > 0), 
	  t = [t t(n+k)*ones(1,needed)];a = [a zeros(1,needed)]; n = n+needed;
   end

   intgrf=spmak([t t(n+k)],cumsum((a.*(ones(d,1)*((t(k+[1:n])-t(1:n))/k)))')');
else,  % unknown representation
   error('unknown function type encountered')
end
