##USAGE   C = powm ( B, A  {,CCHK=1} )
##
##	power function C = B^A
##	arbitrary (numerical) input types
##
##	square matrix A exponent possible if
##	* base B scalar  -or-
##	* base B square matrix commutating with exponent A  >>NEW<<
##CCHK	flag for commutator matrix check
##	ONLY relevant at equally sized  SQUARE MAT B,A  input.
##	     If you KNOW a priori, that 'commutate(A,B)==TRUE',
##	     you may set CCHK=0 to speed-up calculation.
##	     Especially useful if called within iterations.
##	WARN result is definitely INCORRECT in general, if
##	     CCHK is set to 0 , but 'commutate(A,B)==FALSE'
##	>>>  DON'T USE THIS OPTION  if you're in doubt <<<
##NOTE	syntax identical to C/C++ pow scalar-function
##SEE TOO commutate              (author -written)
##	expm, logm, operator ^ (octave built-in)
##AUTHOR  (c) Rolf Fabian 000120V <fabian@tu-cottbus.de>

#used-scripts
#commutate.m	author-written

#FIXME	There's still the problem, that 
#	for e.g.  B=ones(2)  i.e. non-diagonalizable
#	powm (B,A) hangs     e.g. for A=B
#	even octave:>  B^(non-integer scalar A)  FAILS

function C=powm(B,A,CCHK)

	if (nargin>3||nargin<2)
		usage ("C = powm ( B, A  {,CCHK=1} )"); C=[]; return;
	else
		if (nargin <3) CCHK=1; end 
		tag  ="powm: ";
	   	merr0="%s square matrix  %s required .. quit.\n";
		merr1=\
		"%s  B^A  undefined  for A not commutating with B  .. quit.\n";
	   	merr2="%s A,B square matrices,  but sizes mismatch .. quit.\n";

		[ar,ac]=size(A);
		[br,bc]=size(B);

		C=[]; 				#init output

		if (ar!=ac)
			fprintf(stderr,merr0,tag,'A'); return;
		end
		if (br!=bc)
			fprintf(stderr,merr0,tag,'B'); return;
		end

		if	( (BS=all([br,bc]==1)) )#scalar base
			if (B==0)
			    if (any(any(A)))
				 C=zeros(ar);	#overwrite octave NaN-output
			    else C=  eye(ar);	#shortcut
			    end
			else
			    C=B^A;
			end
		elseif	(!BS&& all([ar,ac]==1)) #MAT base + scalar exponent 
			if (any(any(B)))
				C=B^A;
			else			#here	scalar A
				if (A) C=zeros(br);
				else   C=eye(br);
				end
			end
		elseif	(ar==br)		#equally sized A,B
			if     ( !any(any(B)) )	#zero MAT base
			    if (  any(any(A)) )
				 C=zeros(ar);	#zeros(N)^MAT(N)  := zeros(N)
			    else C=  eye(ar);	#zeros(N)^zeros(N) =\
			    end			#  MAT(N)^zeros(N):=   eye(N)
			else
			    #call 'commutate' ONLY if REQUESTED by nonzero CCHK 
			    #REM  if called, 'commutate' checks for equal sizes too
			    if ( !CCHK || commutate(A,B) )
				 C = expm(A*logm(B));
			    else
			     	 fprintf(stderr,merr1,tag);
			    end
			end
		else
			fprintf(stderr,merr2,tag);
		end
	end
end
