=================================================================================
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:   getf("rangj.sci");stacksize(2000000);
// 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=rand(1,N,"normal");b=rand(1,N,"normal");c=rand(1,N,"normal");
  JAC=[JAC;kron(I,kron(b,c));kron(a,kron(I,c));kron(a,kron(b,I))];
 temp=ranknew;ranknew=rank(JAC);
 printf('ranknew=%i\n',ranknew)
 rankold=temp;
 end;
 D=ranknew;
elseif K==4,
 rankold=0;ranknew=1;
 while ranknew>rankold,
  a=rand(1,N,"normal");b=rand(1,N,"normal");c=rand(1,N,"normal");d=rand(1,N,"normal");
  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);
 printf('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;
endfunction

=================================================================================
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: stacksize(2000000);
// getf("rangjs.sci");for i=2:6,for j=3:4,M(i-1,j-2)=rangjs(i,j);end;end;M
// NB: preferably limited to d=3 or d=4 and to N<7 for reasons of memory space
getf("/users/comon/Documents/Scilab/fonctions/gsm.sci")
I=eye(N,N); JAC=[];
if d==3,
 rankold=0;ranknew=1;
 while ranknew>rankold,
  a=rand(1,N,"normal");
  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);
 printf('ranknew=%i\n',ranknew)
 rankold=temp;
 end;
 D=ranknew;
elseif d==4,
 rankold=0;ranknew=1;
 while ranknew>rankold,
  a=rand(1,N,"normal");
  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);
 printf('ranknew=%i\n',ranknew)
 rankold=temp;
 end;
 D=ranknew;
else
 Rbar=0;D=0;
end;
[lign,col]=size(JAC);Rbar=lign/N-1;
endfunction

=================================================================================
function [Rbar,D,F]=rangj3(N1,N2,N3)
// [Rbar,D,F]=rangj3(N1,N2,N3)
// P.Comon, 22 May 2006.
// Rbar is the generic rank of tensors of order 3 and dimensions N1xN2xN3, 
// 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 loading vectors in each term, F=F=Rbar*(N1+N2+N3-2)-D
// Call:   getf("rangj3.sci");stacksize(2000000);
// [Rbar,D,F]=rangj3(N1,N2,N3)
// Example: K=4;for I=2:9,for J=3:5,M(I-1,J-2)=rangj3(I,J,K);end;end;M
I1=eye(N1,N1);I2=eye(N2,N2);I3=eye(N3,N3);JAC=[];
 rankold=0;ranknew=1;
 while ranknew>rankold,
  a=rand(1,N1,"normal");b=rand(1,N2,"normal");c=rand(1,N3,"normal");
  JAC=[JAC;kron(I1,kron(b,c));kron(a,kron(I2,c));kron(a,kron(b,I3))];
 temp=ranknew;ranknew=rank(JAC);
 printf('ranknew=%i\n',ranknew)
 rankold=temp;
 end;
 D=ranknew;
[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*(N2+N3-1)-D
// rmin is a lower bound on the generic rank, if known. The default is 1.
// Call:   getf("rgindscal3.sci");stacksize(2000000);
// 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=rand(1,N2);c=rand(1,N3);
  JAC=[JAC;kron(I2,kron(b,c))+kron(b,kron(I2,c));kron(b,kron(b,I3))];
 temp=ranknew;ranknew=rank(JAC);
 printf('ranknew=%i\n',ranknew)
 rankold=temp;
 end;
 D=ranknew;
[lign,col]=size(JAC);Rbar=lign/(N2+N3)-1;F=Rbar*(N2+N3-1)-D;
endfunction

=================================================================================
function [Rbar,D,F]=rgindscal2z(N,K)
// [Rbar,D,F]=rgindscal2z(N,K)
// P.Comon, 1st July 2006.
// Rbar is the generic rank of tensors of order 3 and dimensions (N,N,K) 
// with symmetric slices across the two first modes, and doubly centered
// 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 loading vectors in each term, F=Rbar*(N+K-2)-D
// Call:   getf("rgindscal2z.sci");stacksize(2000000);
// for K=2:15, for N=2:5, M(K-1,N-1)=rgindscal2z(N,K);end;end;M
I3=eye(K,K);JAC=[];IZ=[eye(N-1,N-1),-ones(N-1,1)];
 rankold=0;ranknew=1;
 while ranknew>rankold,
  b=rand(1,N-1);bb=[b,-sum(b)];c=rand(1,K);
  JAC=[JAC;kron(IZ,kron(bb,c))+kron(bb,kron(IZ,c));kron(bb,kron(bb,I3))];
 temp=ranknew;ranknew=rank(JAC);
 printf('ranknew=%i\n',ranknew)
 rankold=temp;
 end;
 D=ranknew;
[lign,col]=size(JAC);Rbar=lign/(N-1+K)-1;F=Rbar*(N+K-2)-D;
endfunction

=================================================================================
function [Rbar,D,F]=rangj3z(I,J,K)
// [Rbar,D,F]=rangj3z(I,J,K)
// P.Comon, 20 sept 2006.
// Rbar is the generic rank of tensors of order 3 and dimensions (I,J,K) 
// with free matrix slices across the two first modes, but doubly centered
// 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 loading vectors in each term, F=Rbar*(I+J+K-2)-D
// Call:   
// K=3;for I=2:12, for J=2:5, M(I-1,J-1)=rangj3z(I,J,K);end;end;M
IK=eye(K,K);JAC=[];
IZ=[eye(I-1,I-1),-ones(I-1,1)];
JZ=[eye(J-1,J-1),-ones(J-1,1)];
 rankold=0;ranknew=1;
 while ranknew>rankold,
  a=rand(1,I-1);aa=[a,-sum(a)];b=rand(1,J-1);bb=[b,-sum(b)];c=rand(1,K);
  JAC=[JAC;kron(IZ,kron(bb,c));kron(aa,kron(JZ,c));kron(aa,kron(bb,IK))];
 temp=ranknew;ranknew=rank(JAC);
 printf('ranknew=%i\n',ranknew)
 rankold=temp;
 end;
 D=ranknew;
[lign,col]=size(JAC);Rbar=lign/(I+J+K-2)-1;F=Rbar*(I+J+K-3)-D;
endfunction


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