function intgrf = fnint(f,ifa)
%FNINT Integrate a function.
%
%   FNINT(F)  returns the (representation of the) indefinite integral that 
%   at the left end point of the basic interval for (the function in) F.
%
%   FNINT(F,IFA)  provides the indefinite integral whose value at the 
%   left end point of the basic interval is IFA.
%
%   For example,
%
%      fnder(fnint(f));
%
%   is the same as F (except for round-off and, possibly, the multiplicity of
%   end knots).
%
%   See also FNDER.

%   Carl de Boor 25 feb 89
%   cb : Mar.21, 1991 (correct misprint in line 15)
%   cb : Sep.23, 1991 (add more knots at right end for B-form)
%   cb :  6 Feb  1992  (finish the last correction)
%   cb :  9 may '95 (use .' instead of ')
%   cb : 19 nov '95 (correct misprint; introduce second input argument)
%   cb : 29aug96 (include BB-form)
%   cb :  5 oct 97 (make use of ANS explicit, to help compilation)
%   cb : 27oct97 (switch to structures)
%   cb : 05may98 (standardize the help; accept older version of the forms)
%   cb : 11dec98 (use repmat to replicate)
%   Copyright 1987-2000 by C. de Boor and The MathWorks, Inc.
%   $Revision: 1.9 $

if ~isstruct(f), f = fn2fm(f); end

switch f.form(1)
case 'p'      % the function is in pp-representation:
   
   if length(f.order)>1
         error('Use FNDER with negative derivative order.'), end
      
   [breaks,coefs,l,k,d]=ppbrk(f);
   coefs=coefs./repmat([k:-1:1],d*l,1);
   if nargin==1, ifa = zeros(d,1); end
   if l<2
      intgrf=ppmak(breaks,[coefs ifa],d);
   else
      % evaluate each integrated polynomial at the right endpoint of its
      % interval (this is adapted from PPUAL)
      xs=diff(breaks(1:l));index=[1:l-1];
      if d>1
         xs = reshape(repmat(xs,d,1),1,(l-1)*d);
         index = reshape(1+repmat(d*index,d,1)+repmat([-d:-1].',1,l-1), ...
                        (l-1)*d,1);
      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([ifa junk].')).';
      else
         last=cumsum([ifa,vv]);
      end

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

case 'B'          % the function is in B-form or BB-form.
   
   if length(f.order)>1
      error('Use FNDER with negative derivative order.'), end
   
   % 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(d,needed)]; n = n+needed;
   end

   if nargin>1 % if a left-end value is specified, increase left-end knot
               % multiplicity to k+1, making the additional coefficients
               %  0 , then add IFA to all coefficients of the integral.
      needed = k - index(1);
      intgrf = spmak([t(ones(1,needed+1)) t t(n+k)], ...
    cumsum([ifa,zeros(d,needed),a.*(ones(d,1)*((t(k+[1:n])-t(1:n))/k))].').');
   else
      intgrf = spmak([t t(n+k)], ...
                     cumsum((a.*(ones(d,1)*((t(k+[1:n])-t(1:n))/k))).').');
   end
otherwise  % unknown representation
   error('Unknown function type encountered.')
end
