/*									* 
 * diag([a1,a2,..])							*
 * ai is square matrix or one element					*
 * Jordan cell JF(labmda,n),labmda is eigenvalue,n is size of matrix 	*
 * jordan(mat) is jordan form but represented by list			*
 * we obtain matrix form ,dispJordan(this list)				*
 * minimalPoly(this list) display minimalpolynomial of that matrix	*
 * (M^^-1)AM=JF,I soon show maybe at this week end			*
 * I intended to write this program,after saw Barton's mat_exp.mc 	*
 * Partially because of Mathematica BUG ,				*
 * a:matix( [2,0,0,0,0,0,0,0],						*
 *        [1,2,0,0,0,0,0,0],						*
 *        [-4,1,2,0,0,0,0,0],						*
 *        [2,0,0,2,0,0,0,0],						*
 *        [-7,2,0,0,2,0,0,0],						*
 *        [9,0,-2,0,1,2,0,0],						*
 *        [-34,7,1,-2,-1,1,2,0],					*
 *        [145,-17,-16,3,9,-2,0,3])$	Mathematica returns wrong answer*
 * Jordan form,my program returns correct answer for this matrix        *
 * this program is GPL.							*
 */
 /* Make an n element Maxima list [a,a,...a].
  thank for Barton */
  load("eigen")$

  duplicate(a,n) := block([k],
  makelist(a,k,1,n))$
 

  diag(x):= block([ mmatrix,z1,ztemp,zero,n,ntemp,sumn,i,j,k,s,front,back],
    mode_declare([z1,ztemp,zero],list,[mmatrix],any,[n,ntemp,sumn,i,j,k,s,front,back],fixnum),
    n:length(x),sumn:0,front:0,sold:0,
    /*gard*/
    for i:1 thru n do (
      if (matrixp(x[i])) then 
      	x[i]:x[i] 
      else ( 
       	x[i]:matrix([x[i]]))
    ),
   /*sumn is matrix size */
    for i:1 thru n do sumn:length(x[i])+sumn,
    ztemp:[],z1:[],back:sumn,
    for i:1 thru n do (
      s:length(x[i]),
      back:back-s,front:front+sold,
      if (front=0) then
      	  zero:[]
      else (
       	  zero:duplicate(0,front)),
      for k:1 thru s do (
      	  ztemp:zero,
      	for j:1 thru s do (
      	    ztemp:endcons(x[i][k,j],ztemp)
      	),
      	ztemp:append(ztemp,duplicate(0,back)),
      	if (i=1 and k=1) then 
	    mmatrix:matrix(ztemp)
      	else (
	    mmatrix:addrow(mmatrix,ztemp))
       ),   
       sold:s 
     ),
  
  return(mmatrix))$
 
 /*Jordan cell JF(labmda,n),labmda is eigenvalue,n is size of matrix */
 JF(labmda,n):=block([mmatrix,ztemp,i],
    mode_declare([ztemp],list,[i,n],fixnum,[mmatrix],any),
    if (n <= 0) then (
       print("can not calc JF"),
       return(0)
    )else (
       if (n=1) then (
	 return(matrix([labmda]))
       ) else (
      	   front:-1,back:n-1,ztemp:[],
       	   for i:1 thru n do (
      	       front:front+1,back:back-1,
      	       if (front=0) then (
		  mmatrix:matrix(append([labmda,1],duplicate(0,back)))
      	       ) else ( 
	       	   if (front <= n-2) then (
      	      	       ztemp:append(duplicate(0,front),[labmda,1],duplicate(0,back)),
	      	       mmatrix:addrow(mmatrix,ztemp)
       	      	   ) else (
	      	       ztemp:append(duplicate(0,front),[labmda]),
	      	       mmatrix:addrow(mmatrix,ztemp)))
	    ),
 return(mmatrix))))$
	  

/*jordan(matrix),return diag(JF(labmda1,n1),..,JF(labmdak,nk))
or list [[labmda1,n1,..nl],[labmda2,m1,..mj],..] and M matrix  */

jordan(mat) := block( [JFLIST,eigenlist,multilist,zztemp,i,j,k,n,t,ii,neigen,tmpmulti,
     blocknum,pmulti,NJ,mattmp],
   mode_declare([i,j,k,n,t,ii,neigen,pmulti,tmpmulti,blocknum],fixnum,[JFLIST,eigenlist,
    multilist,NJ,zztemp],list,[mat,mattmp],any),
   n:length(mat),
   eigenlist:eigenvectors(mat),
   neigen:length(eigenlist[1][1]),
   multilist:eigenlist[1][2],
   JFLIST:[],
   for i:1 thru neigen do (
      pmulti:multilist[i],
      zztemp:[],NJ:duplicate(0,pmulti),tmpmulti:0,
     
      if (pmulti=1) then (
	 zztemp:cons(1,zztemp)
      ) else (
      	 mattmp:mat-(eigenlist[1][1][i])*ident(n),
	 blocknum:n-rank(mattmp),
	 if (blocknum=1) then (
	    zztemp:cons(pmulti,zztemp)
	 ) else (
	    if (blocknum=pmulti) then (
	       for ii:1 thru pmulti do zztemp:cons(1,zztemp)
	    ) else (
	       	for k:2  while (pmulti #tmpmulti) do (
       		   s:0,
		   t:n-rank(mattmp^^k),
		   if (k=2) then (
		      s:0
		   ) else (
		       for ii:1 thru k-2 do s:s+(k-ii)*NJ[ii]),
		   NJ[k-1]:k*blocknum-s-t,
		   for ii:1 thru NJ[k-1] do zztemp:cons(k-1,zztemp),
		   tmpmulti:tmpmulti+NJ[k-1]*(k-1)
		 )))),
		      
       JFLIST:endcons(cons(eigenlist[1][1][i],zztemp),JFLIST) 
     ),
  return(JFLIST))$
          
minimalPoly(a):=block([i,n,x],
  mode_declare([i,n],fixnum,[a],list,[x,T],any),
  n:length(a),T:1,
  for i:1 thru n do (
     T:T*(x - a[i][1])^(a[i][2])
  ),
  return(T))$

dispJordan(a):=block([i,n,p],
   mode_declare([i,n,j],fixnum,[a],list,[vtemp],any),
   n:length(a),vtemp:[],
   for i:1 thru n do (
      for j:1 thru length(a[i])-1 do (
       	 vtemp:endcons(JF(a[i][1],a[i][j+1]),vtemp)
      )
   ),
  
  return(diag(vtemp)))$