/*****************************************************************************
*                                                                            *
*   POWM.SPL     Copyright (C) 2019 DSP Development Corporation              *
*                               All Rights Reserved                          *
*                                                                            *
*   Author:      Randy Race                                                  *
*                                                                            *
*   Synopsis:    Computes the matrix power                                   *
*                                                                            *
*   Revisions:   10 May 2019  RRR  Creation                                  *
*                                                                            *
*****************************************************************************/


#if @HELP_POWM

    POWM

    Purpose: Raises a square matrix by a square matrix.

    Syntax:  POWM(a, b)

              a - Input scalar or square matrix, the base.

              b - Input scalar or square matrix, the exponent

    Returns: An N x N array.

    Example:
             W1: ravel(-1..2, 3)
             W2: powm(w1, w1);
             W3: powm(w2, inv(w1))

             W1 contains the array: 

              {{-1, 1}.
               { 0, 2}}

             W2 == 

              {{-1.000, 1,667}.
               { 0.000, 4.000}}

             W3 ==

              {{-1, 1}.
               { 0, 2}}

    Remarks:
             POWM(A) returns an array X such that X = A ^^ B by computing:

                X = expm(B *^ logm(A))

    See Also:
             ^^
             Eig
             Expm
             Logm
#endif


/* generalized matrix power */
powm(a, b, tol)
{
        local nra, nca, nrb, ncb;

        if (argc < 3)
        {
                if (argc < 2)
                {
                        if (argc < 1) error(sprintf("%s - matrix required", __FUNC__));

                        b = refseries(a);
                }

                tol = 10 * eps;
        }

        (nra, nca) = size(a);

        if (nra != nca) error(sprintf("%s - matrix must be square", __FUNC__));

        (nrb, ncb) = size(b);

        if (not(nrb == 1 && ncb == 1))
        {
                if (not(nra == 1 && nca == 1))
                {
                        if (not(nrb == nra && ncb == nca))
                        {
                                error(sprintf("%s - matrices must be the same size", __FUNC__));
                        }
                }
        }

        /* power */
        x = expm(b *^ logm(a));

        /* check if purely real */
        if (all(all(abs(imag(x)) < tol)))
        {
                x = real(x);
        }

        return(x);
}