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
    q_shift

    %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^2-eye(n);
    endfor
    
    %end cos computation

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

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


endfunction