function [B,N,s,Cost] = cosher(A,kmax)
%   COSHER  Matrix cosine function based on Hermites series.
%
%   cosher(A, kmax) computes the matrix cosine of A using a scaling algorithm based
%   on the double angle formula and a Hermite approximation.
%   kmax is a vector index to determine the maximum Hermite approximation order 'N' to be used:
%     
%                                           kmax:  1 2 3 4 5  6  7  8
%      Max. Hermite approximation order (N(kmax)): 1 2 4 6 9 12 16 20
%
%      If kmax<2, kmax is set to 2. If kmax>8 or kmax is not provided, it is set to 8.
%
%   Function returns the cosine matrix B of the matrix A; the Hermite appoximation
%   order used N; the scaling parameter s; and the Cost in terms of 
%   number of matrix products.
%
%   References:
%
%   [1] E. Defez, J. Sastre, J. Ibáñez, P. Ruiz, Solving engineering models 
%       using matrix functions, in Proc. of the International Conference on 
%       Mathematical Modelling in Engineering & Human Behaviour 2011, Valencia, 
%       September 6-9th. 
%   
%   [2] E. Defez,J. Sastre, J. Ibáñez, P. Ruiz, Computing matrix functions 
%       arising in engineering models with orthogonal matrix polynomials, 
%       Mathematical and Computer Modelling 57(2013), pp. 1738–1743.
%
%  Authors: Emilio Defez, Jorge Sastre, Javier Ibáñez and Pedro Ruiz.
%  Revised version: 2015/11/25
%
%   Group of High Performance Scientific Computation (HiPerSC)
%   Universitat Politecnica de Valencia (Spain)
%   http://hipersc.blogs.upv.es


% Check input arguments
minkmax = 1;  % Minimum considered order (N = 1)
maxkmax = 8;  % Maximum considered order (N = 20)
if nargin<1 || nargin>2
    error('Incorrect number of input arguments.');
elseif nargin==1
    kmax = maxkmax;
elseif nargin==2
    if kmax<minkmax
        kmax = minkmax;
    elseif kmax>maxkmax
        kmax = maxkmax;
    end
end

saved_products=0;
u=eps/2; %the unit roundoff in the double precision floating-point
if max(max(abs(A)))<eps
    B=eye(size(A,1));
    N=0;s=0;Cost=0;
    return
end

vN=[1 2 4 6 9 12 16 20];%possible values of N

theta=[ 0.13988322173046763e-3 %N= 1
        0.45977704110066707e-2 %N= 2
        0.90556596644120163e-1 %N= 4
        0.36534325997941364    %N= 6
        1.15436374958047930    %N= 9;
        2.30098997117702760    %N=12
        4.20737031121960840    %N=16
        6.39599087275650820    %N=20
      ];

%Selecting scaling factor s and Hermite approximation order N
A2 = A*A;
nA2=norm(A2,1);
norma=sqrt(nA2);

for j=1:kmax
	if norma<=theta(j)
        break
    end
end
if (j==kmax)&&(norma>theta(kmax))
    s=max(0,ceil(log2(norma/theta(kmax))));
    A2=A2/4^s;
else
    s=0;
end

N=vN(j);

%Compute the Hermite approximation of cos(A/2^s)
n=size(A,1);
I=eye(n);
switch N
    case 1
        p=[ 0.99999999999999999925418135716616, -0.49999999938933698223561618618159];
        B= p(2)*A2+p(1)*I;
    case 2
        p=[ 1.0000005859725482336420265538639, -0.50000029298618791771002877253098, 0.041666666615778081852968015515132];
        B=p(3)*A2+p(2)*I;
        if norm(B,1)*nA2>u
            B= A2*B+p(1)*I;
        else
            B=p(1)*I;
            saved_products=saved_products+1;
        end
	case 4 
        A4=A2*A2;
        nA4=norm(A4,1);
        %p=exp(-1/L^2)*[(L^6 + L^4/2 + L^2/6 + 1/24)/L^8 + 1, - (L^4/2 + L^2/4 + 1/12)/L^6 - 1/2, (L^2/24 + 1/48)/L^4 + 1/24, - 1/(720*L^2) - 1/720, 1/40320];
        p=[ 0.99999999999999999999993819614656, -0.49999999999999999835974735024191, 0.041666666666660862426684695556712, -0.0013888888827271599215613898634118, 0.000024799251127513047521570072601744];
        B= p(5)*A4+p(4)*A2+p(3)*I;
        if norm(B,1)*nA4>u
            B = B*A4+p(2)*A2+p(1)*I;
        else
            B =      p(2)*A2+p(1)*I;
            saved_products=saved_products+1;
        end
	case 6  
        A4=A2*A2;
        nA4=norm(A4,1);
        %p=exp(-1/L^2)*[(L^10 + L^8/2 + L^6/6 + L^4/24 + L^2/120 + 1/720)/L^12 + 1, - (L^8/2 + L^6/4 + L^4/12 + L^2/48 + 1/240)/L^10 - 1/2, (L^6/24 + L^4/48 + L^2/144 + 1/576)/L^8 + 1/24, - (L^4/720 + L^2/1440 + 1/4320)/L^6 - 1/720, (L^2/40320 + 1/80640)/L^4 + 1/40320, - 1/(3628800*L^2) - 1/3628800, 1/479001600];
        p=[ 0.99999999999999999999999997306381, -0.4999999999999999999998746109367, 0.041666666666666666583281827020719, -0.0013888888888888704048704435039524, 0.000024801587299831261822688373209628, -0.00000027557311438423927587486667648646, 0.0000000020861066006143993630690678884352];
        B= p(7)*A4+p(6)*A2+p(5)*I;
        if norm(B,1)*nA4^2>u
            B = B*A4+p(4)*A2+p(3)*I;
        else
            B =      p(4)*A2+p(3)*I;
            saved_products=saved_products+1;
        end
        if norm(B,1)*nA4>u
            B=    B*A4+p(2)*A2+p(1)*I;
        else
            B=         p(2)*A2+p(1)*I;
            saved_products=saved_products+1;
        end
            
	case 9
        A4=A2*A2;
        A6=A4*A2;
        nA6=norm(A6,1);
        p=[ 0.99999999999999999999999999999991, -0.49999999999999999999999999987047, 0.041666666666666666666666639457973, -0.0013888888888888888888868566688076, 0.000024801587301587301516150551154419, -0.00000027557319223985755422594860420867, 0.0000000020876756987727132279078219166973, -0.000000000011470745510941823696972632902491, 0.000000000000047794469395500454186342367232531, -0.00000000000000015563538264585705314487460835329];
        B= p(10)*A6+p(9)*A4+p(8)*A2+p(7)*I;
        if norm(B,1)*nA6^2>u
            B=     B*A6+p(6)*A4+p(5)*A2+p(4)*I;
        else
            B=          p(6)*A4+p(5)*A2+p(4)*I;
            saved_products=saved_products+1;
        end
        if norm(B,1)*nA6>u
            B=     B*A6+p(3)*A4+p(2)*A2+p(1)*I;
        else
            B=          p(3)*A4+p(2)*A2+p(1)*I;
            saved_products=saved_products+1;
        end
	case 12
        A4=A2*A2;
        A6=A4*A2;
        nA6=norm(A6,1);
        p=[ 1.0, -0.49999999999999999999999999999972, 0.041666666666666666666666666649184, -0.001388888888888888888888888486931, 0.000024801587301587301587297086377657, -0.00000027557319223985890649734895097466, 0.0000000020876756987868097906515186838188, -0.000000000011470745597729465970296258131948, 0.000000000000047794773323468141859258551172228, -0.00000000000000015619206927001013731543890209762, 0.00000000000000000041103148766549802158860663684131, -0.00000000000000000000088956716860549123234855937451272, 0.0000000000000000000000015862344531539950028635673428217];
        B=p(13)*A6+p(12)*A4+p(11)*A2+p(10)*I;
        if norm(B,1)*nA6^3>u
            B=    B*A6+ p(9)*A4+ p(8)*A2+ p(7)*I;
        else
            B=          p(9)*A4+ p(8)*A2+ p(7)*I;
            saved_products=saved_products+1;
        end
        if norm(B,1)*nA6^2>u
            B=    B*A6+ p(6)*A4+ p(5)*A2+ p(4)*I;
        else
            B=          p(6)*A4+ p(5)*A2+ p(4)*I;
            saved_products=saved_products+1;
        end
        if norm(B,1)*nA6>u
            B=    B*A6+ p(3)*A4+ p(2)*A2+ p(1)*I;
        else
            B=          p(3)*A4+ p(2)*A2+ p(1)*I;
            saved_products=saved_products+1;
        end
	case 16
        A4=A2*A2;
        A6=A4*A2;
        A8=A4*A4;
        nA8=norm(A8,1);
        p=[ 1.0, -0.5, 0.041666666666666666666666666666667, -0.0013888888888888888888888888888924, 0.000024801587301587301587301587313235, -0.00000027557319223985890652557319204536, 0.0000000020876756987868098979211629888358, -0.000000000011470745597729724713851274231306, 0.000000000000047794773323873852974376947013898, -0.00000000000000015619206968586227476253490050706, 0.0000000000000000004110317623312164305266390046995, -0.00000000000000000000088967913924503731177904443252482, 0.0000000000000000000000016117375710824370281344578274338, -0.0000000000000000000000000024795962566230190309067916089332, 0.0000000000000000000000000000032798870453337818335080645666578, (-3.7695131567199169574552986591543)*10^(-33), 3.7402557704842414806645185093618*10^(-36)];
        B=p(17)*A8+p(16)*A6+p(15)*A4+p(14)*A2+p(13)*I;
        if norm(B,1)*nA8^3>u
            B=    B*A8+p(12)*A6+p(11)*A4+p(10)*A2+ p(9)*I;
        else
            B=         p(12)*A6+p(11)*A4+p(10)*A2+ p(9)*I;
            saved_products=saved_products+1;
        end
        if norm(B,1)*nA8^2>u
            B=    B*A8+ p(8)*A6+ p(7)*A4+ p(6)*A2+ p(5)*I;
        else
            B=          p(8)*A6+ p(7)*A4+ p(6)*A2+ p(5)*I;
            saved_products=saved_products+1;
        end
        if norm(B,1)*nA8>u
            B=    B*A8+ p(4)*A6+ p(3)*A4+ p(2)*A2+ p(1)*I;
        else
            B=          p(4)*A6+ p(3)*A4+ p(2)*A2+ p(1)*I;
            saved_products=saved_products+1;
        end
	case 20
        A4=A2*A2;
        A6=A4*A2;
        A8=A4*A4;
        nA8=norm(A8,1);
        p=[ 1.0, -0.5, 0.041666666666666666666666666666672, -0.0013888888888888888888888888889285, 0.000024801587301587301587301587390095, -0.00000027557319223985890652557319076233, 0.0000000020876756987868098979215267458516, -0.000000000011470745597729724713850288494422, 0.000000000000047794773323873852974394034180671, -0.00000000000000015619206968586227459362536429172, 0.00000000000000000041103176233121645226688094989606, -0.00000000000000000000088967913924505738009386714805963, 0.0000000000000000000000016117375710961184347757052441516, -0.0000000000000000000000000024795962632247972591963728241172, 0.0000000000000000000000000000032798892370698349556668126742649, (-3.7699876288149331978254983325809)*10^(-33), 3.8003907546092100129156600770326*10^(-36), (-3.3871574898324706154299582917937)*10^(-39), 2.6882142073996769410301338020973*10^(-42), (-1.9114227353554306784985342284553)*10^(-45), 1.1965913905684116817437346492255*10^(-48)];
        B= p(21)*A8+ p(20)*A6+ p(19)*A4+ p(18)*A2+ p(17)*I;
        if norm(B,1)*nA8^4>u
            B= B*A8     +p(16)*A6+ p(15)*A4+ p(14)*A2+ p(13)*I;
        else
            B=          p(16)*A6+ p(15)*A4+ p(14)*A2+ p(13)*I;
            saved_products=saved_products+1;
        end
        if norm(B,1)*nA8^3>u
            B= B*A8     +p(12)*A6+ p(11)*A4+ p(10)*A2+  p(9)*I;
        else
            B=           p(12)*A6+ p(11)*A4+ p(10)*A2+  p(9)*I;
            saved_products=saved_products+1;
        end
        if norm(B,1)*nA8^2>u
            B= B*A8     + p(8)*A6+  p(7)*A4+  p(6)*A2+  p(5)*I;
        else
            B=            p(8)*A6+  p(7)*A4+  p(6)*A2+  p(5)*I;
            saved_products=saved_products+1;
        end
        if norm(B,1)*nA8>u
            B= B*A8     + p(4)*A6+  p(3)*A4+  p(2)*A2+  p(1)*I;
        else
            B=            p(4)*A6+  p(3)*A4+  p(2)*A2+  p(1)*I;
            saved_products=saved_products+1;
        end
end

for i = 1:s % double angle phase
    B = 2*B*B-I;  
end
Cost=j+s-saved_products;
end

