(* Transportation problem kk X (nn-kk),   1995 Nov 14 by N.Takayama*)
(* 30 lines program to generate the fiber 
   of the kk X (nn-kk) transportation problem *)
(* cf. SIAM J MATH ANAL.  1991, p821-p846, by T.Sasaki*)
(* a[[1]] + ... + a[[kk]] = a[[kk+1]] + ... + a[[nn]] has to be satisfied. *)
(* Usage: set the data kk,nn,a and top on the line 2 and type in trans[] *)
(* ppp is a prime that has to be larger than the largest number in a[[i]] *)

Clear[kk,nn,u,start,a,top,ppp];  
ppp = Prime[1000];
(* Set your kk,nn and a. *)
kk = 3; nn = 6; a = {1,1,2,1,1,2}; top=1;   
(*kk = 3; nn = 6; a = {11,3,4,5,6,7}; top=1;*)  (* 12.5 s *) 
Upq[f_,p_,q_]:=Sum[u[k,p]*D[f,u[k,q]],{k,1,kk}]
Xpa[f_,p_,a_]:=-u[a,p]*f-Sum[ u[a,q]*Upq[f,p,q], {q,kk+1,nn}]
minIndex[a_,top_]:=Block[{k,mindex},mindex=0; 
           For[k=top,k<=nn,k++,If[a[[k]] <= a[[top]] && a[[k]] > 0,mindex=k]];
           Return[mindex]]
nonZeroIndex[a_,f_,t_]:=Block[{k},
           For[k=f,k<=t,k++,If[a[[k]] != 0,Return[k]]]; 
           Return[0]]
trans[]:=Block[{m,s,i,j,plan,f}, 
  If[Sum[a[[i]],{i,1,kk}] != Sum[a[[i]],{i,kk+1,nn}],Return["Invalid Input"]];
  plan = { };
  While[ top >= 1,
    m = minIndex[a,top]; s = a[[m]];
    If[ m <= kk, 
        i=m; j = nonZeroIndex[a,kk+1,nn]; a[[m]]=0; a[[j]] -= s,
        j=m; i = nonZeroIndex[a,1,kk];    a[[m]]=0; a[[i]] -= s];
    plan = Append[plan,{{i,j},s}]; top = nonZeroIndex[a,1,nn]];
  plan = Reverse[plan];
  Print[plan];
  f = 1;
  For[i=1, i<=Length[plan], i++,
    For[j=1,j<= plan[[i]][[2]],j++,
      f = Expand[Xpa[f,plan[[i]][[1]][[2]],plan[[i]][[1]][[1]]]];
      Print[Length[f]];(*f=Mod[f,ppp];*) (*Modulo calc to avoid big numbers *)
  ]]; 
  Return[f];
]





