Matrix-exponentiation operator

Matrix-exponentiation operator
You are encouraged to solve this task according to the task description, using any language you may know.

Most all programming languages have a built-in implementation of exponentiation for integer and real only.

The following programs demonstrates how to implement a "complex number" matrix exponentiation (**) as an operator.

ALGOL 68

main:(

  INT default upb=3;
  MODE VEC = [default upb]COMPL;
  MODE MAT = [default upb,default upb]COMPL;
 
  OP * = (VEC a,b)COMPL: (
      COMPL result:=0;
      FOR i FROM LWB a TO UPB a DO result+:= a[i]*b[i] OD;
      result
    );
 
  OP * = (VEC a, MAT b)VEC: ( # overload vec times matrix #
      [2 LWB b:2 UPB b]COMPL result;
      FOR j FROM 2 LWB b TO 2 UPB b DO result[j]:=a*b[,j] OD;
      result
    );
 
  OP * = (MAT a, b)MAT: ( # overload matrix times matrix #
      [LWB a:UPB a, 2 LWB b:2 UPB b]COMPL result;
      FOR k FROM LWB result TO UPB result DO result[k,]:=a[k,]*b OD;
      result
    );

  OP IDENTITY = (INT upb)MAT:(
    [upb,upb] COMPL out;
    FOR i TO upb DO 
      FOR j TO upb DO
        out[i,j]:= ( i=j |1|0)
      OD
    OD;
    out
  );

  OP ** = (MAT base, INT exponent)MAT: (
    BITS binary exponent:=BIN exponent ;
    MAT out := IF bits width ELEM binary exponent THEN base ELSE IDENTITY UPB base FI;
    MAT sq:=base;

    WHILE 
      binary exponent := binary exponent SHR 1;
      binary exponent /= BIN 0 
    DO
      sq := sq * sq; 
      IF bits width ELEM binary exponent THEN out := out * sq FI
    OD;
    out
  );

  PROC compl matrix printf= (FORMAT compl fmt, MAT m)VOID:(
    FORMAT vec fmt = $"("n(2 UPB m-1)(f(compl fmt)",")f(compl fmt)")"$;
    FORMAT matrix fmt = $x"("n(UPB m-1)(f(vec fmt)","lxx)f(vec fmt)");"$;
    # finally print the result #
    printf((matrix fmt,m))
  );
 
  FORMAT compl fmt = $-z.z,+z.z"i"$; # width of 4, with no leading '+' sign, 1 decimals #
  MAT matrix=((sqrt(0.5)I0         , sqrt(0.5)I0        , 0I0),
              (        0I-sqrt(0.5),         0Isqrt(0.5), 0I0),
              (        0I0         ,         0I0        , 0I1));

  printf(($" matrix ** "g(0)":"l$,24));
  compl matrix printf(compl fmt, matrix**24); print(newline)
)

Output:

matrix ** 24:
(( 1.0+0.0i, 0.0+0.0i, 0.0+0.0i),
 ( 0.0+0.0i, 1.0+0.0i, 0.0+0.0i),
 ( 0.0+0.0i, 0.0+0.0i, 1.0+0.0i));