function cosA = cosm (A)
    if !(ismatrix(A) & issquare(A))
        error("cosm: expecting square matrix argument");
    endif;
    
    #dimention of A
    n=rows(A);
    #normA - the norm_inf of A
    normA=norm(A,inf);
    
    #shift A <- A-pi*q*eye
    normzero = normA;
    normceil = norm(A-pi*ceil(trace(A)/(pi*n))*eye(n),inf);
    normfloor = norm(A-pi*floor(trace(A)/(pi*n))*eye(n),inf);

    switch min([normzero,normceil,normfloor])
        case normzero
            q_shift=0;
        case normceil
            q_shift=ceil(trace(A)/(pi*n)); 
            A=A-pi*q_shift*eye(n);
            normA=normceil;
        case normfloor
            q_shift=floor(trace(A)/(pi*n)); 
            A=A-pi*q_shift*eye(n);
            normA=normfloor;
    endswitch
    #
    #end shift A

    #balancing
    #B=D^-1*A*D
    
    balancing_performed=0;
    [D,B]=balance(A);
    if (norm(B,inf)<normA)
        A=B;
        normA=norm(B,inf);
        balancing_performed=1;
    endif

    #end balancing

    #argument reduction
    #by factor 2^-m, such that 2^-m*normAinf <=1
    m=ceil(log2(normA));
    A=2^(-m)*A;
    
    #end argument reduction


    #padé approximant r88 of cos A
    q=[1.0,
     0.0,
     0.0221377935149955,
     0.0,
     0.000244665145201773,
     0.0,
     1.66685394532938e-6,
     0.0,
     6.23754467948672e-9];
    p=[1.0,
     0.0,
     -0.477862206485004,
     0.0,
     0.0308424350543707,
     0.0,
     -0.000587146544419633,
     0.0,
     3.42184348618221e-6];

    A2=A*A;
    A4=A2*A2;
    A6=A2*A4;
    A8=A4*A4;
    P=p(1)*eye(n)+p(3)*A2+p(5)*A4+p(7)*A6+p(9)*A8;
    Q=q(1)*eye(n)+q(3)*A2+q(5)*A4+q(7)*A6+q(9)*A8;
    R=Q^-1*P;
    #
    #end of padé approximation 

    #cos computation by iterative doubling
    #as proposed by Serbin & Blalock

    cosA=R;
    for i = 0:(m-1)
        cosA=2*cosA-eye(n);
    endfor
    
    #end cos computation

    #undo shifting
    cosA=(-1)^q_shift*cosA;

    #undo balacing
    
    if (balancing_performed)
        cosA=D*cosA*D^-1;
    endif


endfunction
