function [Z1,Z2,INZEE] = fft_cd( Z1 , Z2 , INZEE )

%MATLAB version of SUBROUTINE FFT on page 281 of conte/deboor
%**************   I N P U T   ************
%  Z1, Z2  complex N-vectors
%  INZEE  integer indicating whether  Z1  or  Z2  is to be transformed
%       = 1 , transform  Z1
%       = 2 , transform  Z2
%**************  O U T P U T  ************
%  Z1  or  Z2  contains desired transform (in the correct order)
%  INZEE  integer indicating whether  Z1  or  Z2  contains the transform
%       = 1 , transform is in Z1
%       = 2 , transform is in Z2
%**************  M E T H O D  ************
%  The integer  N := length(ZINZEE(:)) is divided into its prime factors
% (up to a point). For each such factor P , the P-transform of appropriate
% P-subvectors of  Z1  (or  Z2 )  is calculated in  F F T S T E P and stored
% in a suitable way in  Z2  (or  Z1 ). See text for details.
%
%  Reference:
%
%    Samuel Conte, Carl de Boor,
%    Elementary Numerical Analysis,
%    Third Edition,
%    SIAM, 2017,
%    ISBN: 978-1-611975-19-2.
%
PRIME = [2 3 5 7 11 13 17 19 23 29 31 37]; NEXTMX = length(PRIME);
AFTER = 1; BEFORE = length(Z1(:)); NEXT = 1;  

while NEXT<=NEXTMX
   if fix(BEFORE/PRIME(NEXT))*PRIME(NEXT) < BEFORE, NEXT = NEXT+1;
      if NEXT <= NEXTMX, continue
      else, NOW = BEFORE; BEFORE = 1; 
      end
   else, NOW = PRIME(NEXT); BEFORE = BEFORE/NOW; 
   end
   if INZEE==1, [Z2,AFTER,NOW,BEFORE] = fftstep(Z1,AFTER,NOW,BEFORE);
   else,        [Z1,AFTER,NOW,BEFORE] = fftstep(Z2,AFTER,NOW,BEFORE); 
   end
   INZEE = 3-INZEE; if BEFORE==1, return, end
   AFTER = AFTER*NOW;
end
return, end  % actually, this statement will never be executed


