function [B,N,s,Cost] = coshher(A,kmax)
%   COSHHER Matrix hyperbolic cosine function based on Hermites series.
%
%   coshher(A, kmax) computes the matrix hyperbolic 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, J.Peinado, Solving engineering models  
%       using hyperbolic matrix functions, Applied Mathematical Modelling,
%       40(4), Pages, 2016, 2837–2844.

%   Authors:  Emilio Defez, Jorge Sastre, Javier Ibáñez, Jesús Peinado.
%   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=[ 1.398832216450461e-4 %N= 1
        0.004597769511080    %N= 2
        0.090555111535508    %N= 4
        0.365307732936407    %N= 6
        1.153979762538118    %N= 9;
        2.299607443312973    %N=12
        4.203595854404349    %N=16
        6.388962983958694    %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 cosh(A/2^s)
n=size(A,1);
I=eye(n);
switch N
    case 1
        p=[ 0.99999999999999999925418136139761, 0.50000000061066301628067307999673];
        B= p(2)*A2+p(1)*I;
    case 2
        p=[ 1.0000000000000000000337436127884, 0.49999999999991380105760854050549, 0.041666691133049394252225758798027];
        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=[ 1.0000000000000000000000308912629, 0.49999999999999999905816020124351, 0.041666666666670495412667972923443, 0.0013888888842195419696818549313115, 0.000024803621038301388246430685953622];       
        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=[ 1.00000000000000000000000001343, 0.49999999999999999999993093672715, 0.041666666666666666717402804840633, 0.0013888888888888764648274624127412, 0.000024801587302891152877875063089696, 0.00000027557312838565989861701307877965, 0.0000000020890970570218508317527013413452];
        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.99999999999999999999999999999995, 0.5000000000000000000000000000686, 0.041666666666666666666666651192182, 0.001388888888888888888890129989008, 0.000024801587301587301540642848291987, 0.00000027557319223985984203746864389814, 0.0000000020876756987761527936043866464993, 0.000000000011470745668166196572989599669951, 0.000000000000047794508577958753543024578306947, 0.00000000000000015671223237176281115049903631036];
        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.5, 0.041666666666666666666666666666674, 0.001388888888888888888888888888566, 0.000024801587301587301587301594690365, 0.00000027557319223985890652547851813979, 0.0000000020876756987868098986561869210711, 0.0000000000114707455977297210910213320611, 0.000000000000047794773323885456630877273617183, 0.00000000000000015619206966157285775754420565243, 0.00000000000000000041103179508110062440201842891519, 0.00000000000000000000088965190643596246743929856203558, 0.0000000000000000000000016243646951790764664796873383104];
        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.0013888888888888888888888888888905, 0.000024801587301587301587301587289951, 0.00000027557319223985890652557323471158, 0.000000002087675698786809897920924447862, 0.000000000011470745597729724713943387296944, 0.000000000000047794773323873852924056881270712, 0.00000000000000015619206968586227551501481979102, 0.00000000000000000041103176233121646520093651090009, 0.00000000000000000000088967913924504175317682318948197, 0.0000000000000000000000016117375711072668227857749554299, 0.0000000000000000000000000024795962575970955316957684553507, 0.0000000000000000000000000000032798911906084442704401902754839, 3.7695458991065818244622945975955e-33, 3.8587160332091442507915772993069e-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.041666666666666666666666666666665, 0.0013888888888888888888888888889089, 0.00002480158730158730158730158720858, 0.00000027557319223985890652557341674902, 0.0000000020876756987868098979207132834196, 0.000000000011470745597729724714063827108929, 0.000000000000047794773323873852897272947120695, 0.00000000000000015619206968586227572083422386473, 0.00000000000000000041103176233121644901697554111702, 0.00000000000000000000088967913924505738009405739216273, 0.000000000000000000000001611737571096118447409136929898, 0.0000000000000000000000000024795962632247972572893304286897, 0.0000000000000000000000000000032798892370698401372656035829844, 3.7699876288151072439314615948527e-33, 3.800390755064235524085323962831e-36, 3.3871574950285589150510366916995e-39, 2.6882258398126194237681827851627e-42, 1.9114479785482753309485230438393e-45, 1.2541803967754682372429750495242e-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

