=================================================================================function [Rbar,D,F]=rangj(N,K)
% [Rbar,D]=rangj(N,K)
% P.Comon, 21 May 2006.
% Rbar is the generic rank of square tensors of order K and dimension N, 
% and D is the dimension of the whole space of square tensors
% F is the dimension of the variety of solutions
% if a norm constraint is imposed on K-1 loading vectors in each term, F=Rbar*(K*N-K+1)-D
% Call:   
% for n=2:5,for k=3:4,[R,D]=rangj(n,k);G(n-1,k-2)=R;F(n-1,k-2)=R*(k*n-k+1)-D;printf('\n');end;end;G,F
% NB: this implementation limited to K=3 or K=4, and to N<6 for reasons of numerical complexity
% and memory requirements.

I=eye(N,N);JAC=[];
if K==3,
 rankold=0;ranknew=1;
 while ranknew>rankold,
  a=randn(1,N);b=randn(1,N);c=randn(1,N);
  JAC=[JAC;kron(I,kron(b,c));kron(a,kron(I,c));kron(a,kron(b,I))];
 temp=ranknew;ranknew=rank(JAC);
% temp=ranknew;[Q,U]=gsm(JAC);ranknew=rank(U);
 fprintf('ranknew=%i\n',ranknew)
 rankold=temp;
 end;
 D=ranknew;
elseif K==4,
 rankold=0;ranknew=1;
 while ranknew>rankold,
  a=randn(1,N);b=randn(1,N);c=randn(1,N);d=randn(1,N);
  JAC=[JAC;kron(I,kron(b,kron(c,d)));kron(a,kron(I,kron(c,d)));kron(a,kron(b,kron(I,d)));kron(a,kron(b,kron(c,I)))];
 temp=ranknew;ranknew=rank(JAC);
% temp=ranknew;[Q,U]=gsm(JAC);ranknew=rank(U);
 fprintf('ranknew=%i\n',ranknew)
 rankold=temp;
 end;
 D=ranknew;
else
  Rbar=0;D=0;
end;
[lign,col]=size(JAC);Rbar=lign/N/K-1;F=Rbar*(K*N-K+1)-D;

=================================================================================
function [Rbar,D]=rangjs(N,d)
% [Rbar,D]=rangjs(N,d)
% Rbar is the generic rank of symmetric tensors of order d and dimension N, 
% and D is the dimension of the whole space of symmetric tensors
% F=N*Rbar-D is the dimension of the variety of solutions
% P.Comon, 10 January 1996; modified for scilab 20 may 2006.
% Call:
% for i=2:5,for j=3:4,M(i-2,j-2)=rangjs(i,j);end;end;M
% NB: limited to d=3 or d=4 and to N<7 for reasons of numerical complexity

I=eye(N,N); JAC=[];
if d==3,
 rankold=0;ranknew=1;
 while ranknew>rankold,
  a=randn(1,N);
  JAC=[JAC;kron(I,kron(a,a))+kron(a,kron(I,a))+kron(a,kron(a,I))];
 temp=ranknew;[Q,U]=gsm(JAC);ranknew=rank(U);
 fprintf('ranknew=%i\n',ranknew)
 rankold=temp;
 end;
 D=ranknew;
elseif d==4,
 rankold=0;ranknew=1;
 while ranknew>rankold,
  a=randn(1,N);
  JAC=[JAC;kron(I,kron(a,kron(a,a)))+kron(a,kron(I,kron(a,a)))+kron(a,kron(a,kron(I,a)))+kron(a,kron(a,kron(a,I)))];
 temp=ranknew;[Q,U]=gsm(JAC);ranknew=rank(U);
 fprintf('ranknew=%i\n',ranknew)
 rankold=temp;
 end;
 D=ranknew;
else
 Rbar=0;D=0;
end;
[lign,col]=size(JAC);Rbar=lign/N-1;
F=Rbar*N-D;

=================================================================================
function [Rbar,D,F]=rangj3(N1,N2,N3,rmin)
% [Rbar,D,F]=rangj3(N1,N2,N3,rmin)
% P.Comon, 22 May 2006.
% Rbar is the generic rank of tensors of order 3 and dimensions (N1,N2,N3) 
% and D is the dimension of the whole space of square tensors
% F is the dimension of the variety of solutions
% if a norm constraint is imposed on K-1 loading vectors in each term, F=Rbar*(N1+N2+N3-2)-D
% rmin is a lower bound on the generic rank, if known. The default is 1.

if nargin==3, rmin=1;end; % yet unused
I1=eye(N1,N1);I2=eye(N2,N2);I3=eye(N3,N3);JAC=[];

 rankold=0;ranknew=1;
 while ranknew>rankold,
  a=randn(1,N1);b=randn(1,N2);c=randn(1,N3);
  JAC=[JAC;kron(I1,kron(b,c));kron(a,kron(I2,c));kron(a,kron(b,I3))];
 temp=ranknew;ranknew=rank(JAC);
% temp=ranknew;[Q,U]=gsm(JAC);ranknew=rank(U);
 fprintf('ranknew=%i\n',ranknew)
 rankold=temp;
 end;
 D=ranknew;
 Rbar=0;D=0;
[lign,col]=size(JAC);Rbar=lign/(N1+N2+N3)-1;F=Rbar*(N1+N2+N3-2)-D;

=================================================================================
function [Rbar,D,F]=rgindscal3(N2,N3)
% [Rbar,D,F]=rgindscal3(N2,N3)
% P.Comon, 29 June 2006.
% Rbar is the generic rank of tensors of order 3 and dimensions (N2,N2,N3) 
% with symmetric slices across the two first modes
% and D is the dimension of the whole tensor space
% F is the dimension of the variety of solutions
% if a norm constraint is imposed on K-1 loading vectors in each term, F=Rbar*(2*N2+N3-2)-D
% rmin is a lower bound on the generic rank, if known. The default is 1.
% Call:
% for I=2:10, for J=2:5, M(I-1,J-1)=rgindscal3(J,I);end;end;M

I2=eye(N2,N2);I3=eye(N3,N3);JAC=[];

 rankold=0;ranknew=1;
 while ranknew>rankold,
  b=randn(1,N2);c=randn(1,N3);
  JAC=[JAC;kron(I2,kron(b,c));kron(b,kron(I2,c));kron(b,kron(b,I3))];
 temp=ranknew;ranknew=rank(JAC);
% temp=ranknew;[Q,U]=gsm(JAC);ranknew=rank(U);
 fprintf('ranknew=%i\n',ranknew)
 rankold=temp;
 end;
 D=ranknew;
[lign,col]=size(JAC);Rbar=lign/(2*N2+N3)-1;F=Rbar*(N2+N3-1)-D;

=================================================================================
% function [Rbar,D,F]=rgindscal2z(N,K)
% P.Comon, 1st July 2006.
% see Scilab function

=================================================================================
% function [Rbar,D,F]=rangj3z(I,J,K)
% P.Comon, 20 sept 2006.
% see Scilab function

=================================================================================