N = {[32,32,32], [33,32,32], [32,33,32], [32,32,33], [10,10,10,10]};
len = length(N);
rep=2;
feps = 1e-13;

for testcmplx=0:1
  if (testcmplx)
    fprintf('\nComplex Fourier transforms\n');
  else
    fprintf('\nReal Fourier transforms\n');
  end


  for k=1:len
    nn = N{k};

    if (length(nn) < 3)
      fprintf('length([');
      for i=1:length(nn)
        if (i ~= 1)
          fprintf('x%2i', nn(i));
        else
          fprintf('%2i', nn(i));
        end
      end
      fprintf(']) < 3. Continuing\n');
      if (exist('OCTAVE_VERSION'))
        eval('fflush(stdout);');
      end
      continue;
    end

    fprintf('  Testing fft(');
    for i=1:length(nn)
      if (i ~= 1)
        fprintf('x%2i', nn(i));
      else
        fprintf('%2i', nn(i));
      end
    end
    fprintf(')       ');
    for j=1:(6-length(nn))
      fprintf('   ');
    end
    if (exist('OCTAVE_VERSION'))
      eval('fflush(stdout);');
    end

    clear a;
    for i=1:rep
      a = zeros(nn);
      a(:) = randn(prod(nn),1);
      if (testcmplx)
        a = a + 1i;
      end
    
      b = fft(a);
      for j=0:prod(nn(3:length(nn)))-1
	cmd1 = 'c = fft(a(:,:';
	cmd2 = 'd = b(:,:';
        div = 1;
	for k=3:length(nn)
	  kk = floor(mod(j / div, nn(k))) + 1;
	  div = div * nn(k);
          cmd1 = [cmd1 , ',' , num2str(kk)];
	  cmd2 = [cmd2 , ',' , num2str(kk)];
        end
	cmd1 = [cmd1, '));'];
	eval(cmd1);
        cmd2 = [cmd2, ');'];
	eval(cmd2);
        if (any(0.5*abs(d(:) - c(:))./abs(d(:)+c(:)) > feps))
          error('Bad fft(a)');
	end
      end
    end
    fprintf('Ok\n');
  end

  for k=1:len
    nn = N{k};

    if (length(nn) < 3)
      fprintf('length([');
      for i=1:length(nn)
        if (i ~= 1)
          fprintf('x%2i', nn(i));
        else
          fprintf('%2i', nn(i));
        end
      end
      fprintf(']) < 3. Continuing\n');
      if (exist('OCTAVE_VERSION'))
        eval('fflush(stdout);');
      end
      continue;
    end

    fprintf('  Testing fft(');
    for i=1:length(nn)
      if (i ~= 1)
        fprintf('x%2i', nn(i));
      else
        fprintf('%2i', nn(i));
      end
    end
    fprintf(',[],2)  ');
    for j=1:(6-length(nn))
      fprintf('   ');
    end
    if (exist('OCTAVE_VERSION'))
      eval('fflush(stdout);');
    end
  
    clear a;
    for i=1:rep
      a = zeros(nn);
      a(:) = randn(prod(nn),1);
      if (testcmplx)
        a = a + 1i;
      end
    
      b = fft(a,[],2);
      for j=0:prod(nn(3:length(nn)))-1
	cmd1 = 'c = transpose(fft(transpose(a(:,:';
	cmd2 = 'd = b(:,:';
        div = 1;
	for k=3:length(nn)
	  kk = floor(mod(j / div, nn(k))) + 1;
	  div = div * nn(k);
          cmd1 = [cmd1 , ',' , num2str(kk)];
	  cmd2 = [cmd2 , ',' , num2str(kk)];
        end
	cmd1 = [cmd1, '))));'];
	eval(cmd1);
        cmd2 = [cmd2, ');'];
	eval(cmd2);
        if (any(0.5*abs(d(:) - c(:))./abs(d(:)+c(:)) > feps))
          error('Bad fft(a,[],2)');
	end
      end
    end
    fprintf('Ok\n');
  end

  for k=1:len
    nn = N{k};

    if (length(nn) < 3)
      fprintf('length([');
      for i=1:length(nn)
        if (i ~= 1)
          fprintf('x%2i', nn(i));
        else
          fprintf('%2i', nn(i));
        end
      end
      fprintf(']) < 3. Continuing\n');
      if (exist('OCTAVE_VERSION'))
        eval('fflush(stdout);');
      end
      continue;
    end

    fprintf('  Testing fft(');
    for i=1:length(nn)
      if (i ~= 1)
        fprintf('x%2i', nn(i));
      else
        fprintf('%2i', nn(i));
      end
    end
    fprintf(',[],3)  ');
    for j=1:(6-length(nn))
      fprintf('   ');
    end
    if (exist('OCTAVE_VERSION'))
      eval('fflush(stdout);');
    end
  
    clear a;
    for i=1:rep
      a = zeros(nn);
      a(:) = randn(prod(nn),1);
      if (testcmplx)
        a = a + 1i;
      end
    
      b = fft(a,[],3);
      for j=0:prod(nn(3:length(nn)))-1
	cmd1 = ['c = transpose(fft(transpose(reshape(a(:,', ...
		num2str(mod(j,nn(2))+1), ',:'];
	cmd2 = ['d = b(:,', num2str(mod(j,nn(2))+1), ',:'];
        div = nn(2);
	for k=4:length(nn)
	  kk = floor(mod(j / div, nn(k))) + 1;
	  div = div * nn(k);
          cmd1 = [cmd1 , ',' , num2str(kk)];
	  cmd2 = [cmd2 , ',' , num2str(kk)];
        end
	cmd1 = [cmd1, '),',num2str(nn(1)),',',num2str(nn(3)),'))));'];
        cmd2 = [cmd2, ');'];
	eval(cmd1);
	eval(cmd2);
        if (any(0.5*abs(d(:) - c(:))./abs(d(:)+c(:)) > feps))
          error('Bad fft(a,[],3)');
	end
      end
    end
    fprintf('Ok\n');
  end

  for k=1:len
    nn = N{k};

    if (length(nn) < 3)
      fprintf('length([');
      for i=1:length(nn)
        if (i ~= 1)
          fprintf('x%2i', nn(i));
        else
          fprintf('%2i', nn(i));
        end
      end
      fprintf(']) < 3. Continuing\n');
      if (exist('OCTAVE_VERSION'))
        eval('fflush(stdout);');
      end
      continue;
    end

    fprintf('  Testing fft2(');
    for i=1:length(nn)
      if (i ~= 1)
        fprintf('x%2i', nn(i));
      else
        fprintf('%2i', nn(i));
      end
    end
    fprintf(')      ');
    for j=1:(6-length(nn))
      fprintf('   ');
    end
    if (exist('OCTAVE_VERSION'))
      eval('fflush(stdout);');
    end

    clear a;
    for i=1:rep
      a = zeros(nn);
      a(:) = randn(prod(nn),1);
      if (testcmplx)
        a = a + 1i;
      end
    
      b = fft2(a);
      for j=0:prod(nn(3:length(nn)))-1
	cmd1 = 'c = fft2(a(:,:';
	cmd2 = 'd = b(:,:';
        div = 1;
	for k=3:length(nn)
	  kk = floor(mod(j / div, nn(k))) + 1;
	  div = div * nn(k);
          cmd1 = [cmd1 , ',' , num2str(kk)];
	  cmd2 = [cmd2 , ',' , num2str(kk)];
        end
	cmd1 = [cmd1, '));'];
	eval(cmd1);
        cmd2 = [cmd2, ');'];
	eval(cmd2);
        if (any(0.5*abs(d(:) - c(:))./abs(d(:)+c(:)) > feps))
          error('Bad fft2(a)');
	end
      end
    end
    fprintf('Ok\n');
  end

  for k=1:len
    nn = N{k};

    if (length(nn) < 3)
      fprintf('length([');
      for i=1:length(nn)
        if (i ~= 1)
          fprintf('x%2i', nn(i));
        else
          fprintf('%2i', nn(i));
        end
      end
      fprintf(']) < 3. Continuing\n');
      if (exist('OCTAVE_VERSION'))
        eval('fflush(stdout);');
      end
      continue;
    end

    fprintf('  Testing fftn(');
    for i=1:length(nn)
      if (i ~= 1)
        fprintf('x%2i', nn(i));
      else
        fprintf('%2i', nn(i));
      end
    end
    fprintf(')      ');
    for j=1:(6-length(nn))
      fprintf('   ');
    end
    if (exist('OCTAVE_VERSION'))
      eval('fflush(stdout);');
    end

    clear a;
    for i=1:rep
      a = zeros(nn);
      a(:) = randn(prod(nn),1);
      if (testcmplx)
        a = a + 1i;
      end
    
      b = fftn(a);
      c = fft(a);
      for j=2:length(nn);
        c = fft(c,[],j);
      end

      if (any(0.5*abs(b(:) - c(:))./abs(b(:)+c(:)) > feps))
        max(0.5*abs(b(:) - c(:))./abs(b(:)+c(:)))
        error('Bad fftn(a)');
      end

    end
    fprintf('Ok\n');
  end












  for k=1:len
    nn = N{k};

    if (length(nn) < 3)
      fprintf('length([');
      for i=1:length(nn)
        if (i ~= 1)
          fprintf('x%2i', nn(i));
        else
          fprintf('%2i', nn(i));
        end
      end
      fprintf(']) < 3. Continuing\n');
      if (exist('OCTAVE_VERSION'))
        eval('fflush(stdout);');
      end
      continue;
    end

    fprintf('  Testing ifft(');
    for i=1:length(nn)
      if (i ~= 1)
        fprintf('x%2i', nn(i));
      else
        fprintf('%2i', nn(i));
      end
    end
    fprintf(')      ');
    for j=1:(6-length(nn))
      fprintf('   ');
    end
    if (exist('OCTAVE_VERSION'))
      eval('fflush(stdout);');
    end

    clear a;
    for i=1:rep
      a = zeros(nn);
      a(:) = randn(prod(nn),1);
      if (testcmplx)
        a = a + 1i;
      end
    
      b = ifft(a);
      for j=0:prod(nn(3:length(nn)))-1
	cmd1 = 'c = ifft(a(:,:';
	cmd2 = 'd = b(:,:';
        div = 1;
	for k=3:length(nn)
	  kk = floor(mod(j / div, nn(k))) + 1;
	  div = div * nn(k);
          cmd1 = [cmd1 , ',' , num2str(kk)];
	  cmd2 = [cmd2 , ',' , num2str(kk)];
        end
	cmd1 = [cmd1, '));'];
	eval(cmd1);
        cmd2 = [cmd2, ');'];
	eval(cmd2);
        if (any(0.5*abs(d(:) - c(:))./abs(d(:)+c(:)) > feps))
          error('Bad ifft(a)');
	end
      end
    end
    fprintf('Ok\n');
  end

  for k=1:len
    nn = N{k};

    if (length(nn) < 3)
      fprintf('length([');
      for i=1:length(nn)
        if (i ~= 1)
          fprintf('x%2i', nn(i));
        else
          fprintf('%2i', nn(i));
        end
      end
      fprintf(']) < 3. Continuing\n');
      if (exist('OCTAVE_VERSION'))
        eval('fflush(stdout);');
      end
      continue;
    end

    fprintf('  Testing ifft(');
    for i=1:length(nn)
      if (i ~= 1)
        fprintf('x%2i', nn(i));
      else
        fprintf('%2i', nn(i));
      end
    end
    fprintf(',[],2) ');
    for j=1:(6-length(nn))
      fprintf('   ');
    end
    if (exist('OCTAVE_VERSION'))
      eval('fflush(stdout);');
    end
  
    clear a;
    for i=1:rep
      a = zeros(nn);
      a(:) = randn(prod(nn),1);
      if (testcmplx)
        a = a + 1i;
      end
    
      b = ifft(a,[],2);
      for j=0:prod(nn(3:length(nn)))-1
	cmd1 = 'c = transpose(ifft(transpose(a(:,:';
	cmd2 = 'd = b(:,:';
        div = 1;
	for k=3:length(nn)
	  kk = floor(mod(j / div, nn(k))) + 1;
	  div = div * nn(k);
          cmd1 = [cmd1 , ',' , num2str(kk)];
	  cmd2 = [cmd2 , ',' , num2str(kk)];
        end
	cmd1 = [cmd1, '))));'];
	eval(cmd1);
        cmd2 = [cmd2, ');'];
	eval(cmd2);
        if (any(0.5*abs(d(:) - c(:))./abs(d(:)+c(:)) > feps))
          error('Bad ifft(a,[],2)');
	end
      end
    end
    fprintf('Ok\n');
  end

  for k=1:len
    nn = N{k};

    if (length(nn) < 3)
      fprintf('length([');
      for i=1:length(nn)
        if (i ~= 1)
          fprintf('x%2i', nn(i));
        else
          fprintf('%2i', nn(i));
        end
      end
      fprintf(']) < 3. Continuing\n');
      if (exist('OCTAVE_VERSION'))
        eval('fflush(stdout);');
      end
      continue;
    end

    fprintf('  Testing ifft(');
    for i=1:length(nn)
      if (i ~= 1)
        fprintf('x%2i', nn(i));
      else
        fprintf('%2i', nn(i));
      end
    end
    fprintf(',[],3) ');
    for j=1:(6-length(nn))
      fprintf('   ');
    end
    if (exist('OCTAVE_VERSION'))
      eval('fflush(stdout);');
    end
  
    clear a;
    for i=1:rep
      a = zeros(nn);
      a(:) = randn(prod(nn),1);
      if (testcmplx)
        a = a + 1i;
      end
    
      b = ifft(a,[],3);
      for j=0:prod(nn(3:length(nn)))-1
	cmd1 = ['c = transpose(ifft(transpose(reshape(a(:,', ...
		num2str(mod(j,nn(2))+1), ',:'];
	cmd2 = ['d = b(:,', num2str(mod(j,nn(2))+1), ',:'];
        div = nn(2);
	for k=4:length(nn)
	  kk = floor(mod(j / div, nn(k))) + 1;
	  div = div * nn(k);
          cmd1 = [cmd1 , ',' , num2str(kk)];
	  cmd2 = [cmd2 , ',' , num2str(kk)];
        end
	cmd1 = [cmd1, '),',num2str(nn(1)),',',num2str(nn(3)),'))));'];
        cmd2 = [cmd2, ');'];
	eval(cmd1);
	eval(cmd2);
        if (any(0.5*abs(d(:) - c(:))./abs(d(:)+c(:)) > feps))
          error('Bad ifft(a,[],3)');
	end
      end
    end
    fprintf('Ok\n');
  end

  for k=1:len
    nn = N{k};

    if (length(nn) < 3)
      fprintf('length([');
      for i=1:length(nn)
        if (i ~= 1)
          fprintf('x%2i', nn(i));
        else
          fprintf('%2i', nn(i));
        end
      end
      fprintf(']) < 3. Continuing\n');
      if (exist('OCTAVE_VERSION'))
        eval('fflush(stdout);');
      end
      continue;
    end

    fprintf('  Testing ifft2(');
    for i=1:length(nn)
      if (i ~= 1)
        fprintf('x%2i', nn(i));
      else
        fprintf('%2i', nn(i));
      end
    end
    fprintf(')     ');
    for j=1:(6-length(nn))
      fprintf('   ');
    end
    if (exist('OCTAVE_VERSION'))
      eval('fflush(stdout);');
    end

    clear a;
    for i=1:rep
      a = zeros(nn);
      a(:) = randn(prod(nn),1);
      if (testcmplx)
        a = a + 1i;
      end
    
      b = ifft2(a);
      for j=0:prod(nn(3:length(nn)))-1
	cmd1 = 'c = ifft2(a(:,:';
	cmd2 = 'd = b(:,:';
        div = 1;
	for k=3:length(nn)
	  kk = floor(mod(j / div, nn(k))) + 1;
	  div = div * nn(k);
          cmd1 = [cmd1 , ',' , num2str(kk)];
	  cmd2 = [cmd2 , ',' , num2str(kk)];
        end
	cmd1 = [cmd1, '));'];
	eval(cmd1);
        cmd2 = [cmd2, ');'];
	eval(cmd2);
        if (any(0.5*abs(d(:) - c(:))./abs(d(:)+c(:)) > feps))
          error('Bad ifft2(a)');
	end
      end
    end
    fprintf('Ok\n');
  end

  for k=1:len
    nn = N{k};

    if (length(nn) < 3)
      fprintf('length([');
      for i=1:length(nn)
        if (i ~= 1)
          fprintf('x%2i', nn(i));
        else
          fprintf('%2i', nn(i));
        end
      end
      fprintf(']) < 3. Continuing\n');
      if (exist('OCTAVE_VERSION'))
        eval('fflush(stdout);');
      end
      continue;
    end

    fprintf('  Testing ifftn(');
    for i=1:length(nn)
      if (i ~= 1)
        fprintf('x%2i', nn(i));
      else
        fprintf('%2i', nn(i));
      end
    end
    fprintf(')     ');
    for j=1:(6-length(nn))
      fprintf('   ');
    end
    if (exist('OCTAVE_VERSION'))
      eval('fflush(stdout);');
    end

    clear a;
    for i=1:rep
      a = zeros(nn);
      a(:) = randn(prod(nn),1);
      if (testcmplx)
        a = a + 1i;
      end
    
      b = ifftn(a);
      c = ifft(a);
      for j=2:length(nn);
        c = ifft(c,[],j);
      end

      if (any(0.5*abs(b(:) - c(:))./abs(b(:)+c(:)) > feps))
        max(0.5*abs(b(:) - c(:))./abs(b(:)+c(:)))
        error('Bad ifftn(a)');
      end

    end
    fprintf('Ok\n');
  end

end