function [F,m,s,nprods] = costaym(A)
%   COSTAYM  Matrix cosine function based on Taylor series.
%
%   costaym(A) computes the matrix cosine cos(A) using a scaling 
%   and squaring algorithm and Taylor approximation with maximum order equal to 32.
%   Function returns the matrix cosine, A; the appoximation
%   order used m; the scaling parameter s; and nprods in terms of 
%   number of matrix product evaluations.
%
%   References:
%   
%   [1] P. Alonso, J. Ibáñez, J. Sastre, J. Peinado, E. Defez, Efficient
%       and accurate algorithms for computing matrix trigonometric functions, 
%       J. Comput. Appl. Math., 309(2017), pp. 325-332.
%       http://dx.doi.org/10.1016/j.cam.2016.05.015.
%
%   Author: Javier Ibáñez.
%   Revised version: 2016/07/28.
%
%   Group of High Performance Scientific Computation (HiPerSC)
%   Universitat Politecnica de Valencia (Spain)
%   http://hipersc.blogs.upv.es
%  

kmax=7;%optimal order equal to m=16 (2m=32)
n=size(A,1);
I=eye(n);
Fa=[1 2 2 3 3  4  4];
M= [1 2 4 6 9 12 16];
%This theta values come from relativa backward errors 
Theta=[6.661338018806219e-16       %m=1 
     1.154075612730971e-07         %m=2
     2.491236564385514e-03         %m=4
     8.976968236812591e-02         %m=6
     1.189983654063290e+00         %m=9
     4.924177884630485e+00         %m=12
     1.606054585896760e+01];       %m=16
B=A*A;normB=norm(B,1);
%Computing the optimal values of m and s
for i=1:kmax
    if normB<=Theta(i)
        pB=compute_pot(B,Fa(i));
        F=tay_cos(pB,0,i);
        m=M(i);
        s=0;
        nprods=i;
        return
    end
end
m=M(kmax);
pB=compute_pot(B,Fa(kmax));
q=Fa(kmax);c=m/q;
pb2m1=norm1p(pB{q},c)^(1/m);
pb2m2=norm1pp( pB{q},c,B)^(1/(m+1));
alfa_min=max([pb2m1,pb2m2]);
s=max(0,ceil(log2(alfa_min/Theta(kmax))/2));
if isinf(s)
    s=30;
end
%Compute the Taylor approximation of A/4^s
F=tay_cos(pB,s,kmax);
%Recovering of cos(A) from cos(A/4^s)
for i=1:s
    F=2*F*F-I;
end
nprods=kmax+s;
end

function pB=compute_pot(B,q)
pB{1}=B;
for i=2:q
    pB{i}=pB{i-1}*B;
end
end

function [B]=tay_cos(A,s,k)
n = size(A{1},1); 
I = eye(n);
p=[ -2, 24, -720, 40320, -3628800, 479001600, -87178291200, 20922789888000, -6402373705728000, 2432902008176640000, -1124000727777607680000, 620448401733239439360000, -403291461126605635584000000, 304888344611713860501504000000, -265252859812191058636308480000000, 263130836933693530167218012160000000, -295232799039604140847618609643520000000, 371993326789901217467999448150835200000000, -523022617466601111760007224100074291200000000, 815915283247897734345611269596115894272000000000, -1405006117752879898543142606244511569936384000000000, 2658271574788448768043625811014615890319638528000000000, -5502622159812088949850305428800254892961651752960000000000, 12413915592536072670862289047373375038521486354677760000000000, -30414093201713378043612608166064768844377641568960512000000000000, 80658175170943878571660636856403766975289505440883277824000000000000, -230843697339241380472092742683027581083278564571807941132288000000000000, 710998587804863451854045647463724949736497978881168458687447040000000000000, -2350561331282878571829474910515074683828862318181142924420699914240000000000000, 8320987112741390144276341183223364380754172606361245952449277696409600000000000000, -31469973260387937525653122354950764088012280797258232192163168247821107200000000000000, 126886932185884164103433389335161480802865516174545192198801894375214704230400000000000000, -544344939077443064003729240247842752644293064388798874532860126869671081148416000000000000000, 2480035542436830599600990418569171581047399201355367672371710738018221445712183296000000000000000, -11978571669969891796072783721689098736458938142546425857555362864628009582789845319680000000000000000, 61234458376886086861524070385274672740778091784697328983823014963978384987221689274204160000000000000000, -330788544151938641225953028221253782145683251820934971170611926835411235700971565459250872320000000000000000, 1885494701666050254987932260861146558230394535379329335672487982961844043495537923117729972224000000000000000000, -11324281178206297831457521158732046228731749579488251990048962825668835325234200766245086213177344000000000000000000, 71569457046263802294811533723186532165584657342365752577109445058227039255480148842668944867280814080000000000000000000, -475364333701284174842138206989404946643813294067993328617160934076743994734899148613007131808479167119360000000000000000000, 3314240134565353266999387579130131288000666286242049487118846032383059131291716864129885722968716753156177920000000000000000000];
m=[1 2 4 6 9 12 16 20 25 30];

if s>0
    for i=1:length(A);
        A{i}=A{i}/4^(s*i);
    end
end
switch m(k)
    case 1
        B=A{1}/p(1) + I;
    case 2
        B=A{2}/p(2)+ A{1}/p(1) + I;
	case 4  
        B=A{2}/p(4) + A{1}/p(3) + I/p(2);
        B=B*A{2}    + A{1}/p(1) + I;
	case 6 
        B=A{2}/p(6) + A{1}/p(5) + I/p(4);
        B=B*A{2}    + A{1}/p(3) + I/p(2);
        B=B*A{2}    + A{1}/p(1) + I;
	case 9
        B=A{3}/p(9) + A{2}/p(8) + A{1}/p(7) + I/p(6);
        B=B*A{3}    + A{2}/p(5) + A{1}/p(4) + I/p(3);
        B=B*A{3}    + A{2}/p(2) + A{1}/p(1) + I;
	case 12
        B=A{3}/p(12) + A{2}/p(11) + A{1}/p(10) + I/p(9) ;
        B=B*A{3}     + A{2}/p(8)  + A{1}/p(7)  + I/p(6);
        B=B*A{3}     + A{2}/p(5)  + A{1}/p(4)  + I/p(3);
        B=B*A{3}     + A{2}/p(2)  + A{1}/p(1)  + I;
	case 16
        B=A{4}/p(16) + A{3}/p(15) + A{2}/p(14) + A{1}/p(13) + I/p(12);
        B=A{4}*B     + A{3}/p(11) + A{2}/p(10) + A{1}/p(9)  + I/p(8);
        B=A{4}*B     + A{3}/p(7)  + A{2}/p(6)  + A{1}/p(5)  + I/p(4);
        B=A{4}*B     + A{3}/p(3)  + A{2}/p(2)  + A{1}/p(1)  + I;
	case 20
        B=A{4}/p(20)  + A{3}/p(19) +A{2}/p(18)  + A{1}/p(17) + I/p(16);
        B=A{4}*B      + A{3}/p(15) + A{2}/p(14) + A{1}/p(13) + I/p(12);
        B=A{4}*B      + A{3}/p(11) + A{2}/p(10) + A{1}/p(9)  + I/p(8);
        B=A{4}*B      + A{3}/p(7)  + A{2}/p(6)  + A{1}/p(5)  + I/p(4);
        B=A{4}*B      + A{3}/p(3)  + A{2}/p(2)  + A{1}/p(1)  + I;
	case 25
        B=A{5}/p(25) + A{4}/p(24) + A{3}/p(23) + A{2}/p(22) + A{1}/p(21)  + I/p(20);
        B=A{5}*B     + A{4}/p(19) + A{3}/p(18) + A{2}/p(17) + A{1}/p(16)  + I/p(15);
        B=A{5}*B     + A{4}/p(14) + A{3}/p(13) + A{2}/p(12) + A{1}/p(11)  + I/p(10);
        B=A{5}*B     + A{4}/p(9)  + A{3}/p(8)  + A{2}/p(7)  + A{1}/p(6)   + I/p(5);
        B=A{5}*B     + A{4}/p(4)  + A{3}/p(3)  + A{2}/p(2)  + A{1}/p(1)   + I;
	case 30
        B=A{5}/p(30) + A{4}/p(29) + A{3}/p(28) + A{2}/p(27) + A{1}/p(26) + I/p(25);
        B=  A{5}*B   + A{4}/p(24) + A{3}/p(23) + A{2}/p(22) + A{1}/p(21) + I/p(20);
        B=  A{5}*B   + A{4}/p(19) + A{3}/p(18) + A{2}/p(17) + A{1}/p(16) + I/p(15);
        B=  A{5}*B   + A{4}/p(14) + A{3}/p(13) + A{2}/p(12) + A{1}/p(11) + I/p(10);
        B=  A{5}*B   + A{4}/p(9)  + A{3}/p(8)  + A{2}/p(7)  + A{1}/p(6)  + I/p(5);
        B=  A{5}*B   + A{4}/p(4)  + A{3}/p(3)  + A{2}/p(2)  + A{1}/p(1)  + I;    
end
end

function gama=norm1p(A,p)
%   gama=norm1p(A,p)
%
%   Estimates the 1-norm of A^p, where A is a square matrix.
%
%   Inputs:
%      A:  the input matrix
%      p:  power parameter
%
%   Outputs: 
%      gama: estimated 1-morm of A^p
%
%   References: 
%   [1] FORTRAN Codes for Estimating the One-Norm of a Real or Complex
%       Matrix, with Applications to Condition Estimation.
%       ACM Transactions on Mathematical Software, Vol. 14, No. 4, 
%       pp. 381-396, 1988 (Algorithm 4.1).
%
%
%   Revised version 1.0: 16/05/2013

n=size(A,1);
v=A*(ones(n,1)/n);
for i=1:p-1
    v=A*v;
end
gama=norm(v,1);
psi=sign(v);
x=A'*psi;
for i=1:p-1
    x=A'*x;
end
k=2;fin=0;
while fin==0
	[~,j]=max(abs(x));
    ej=zeros(n,1);
    ej(j)=1;
    v=A*ej;
    for i=1:p-1
        v=A*v;
    end
    gama1=gama;
    gama=norm(v,1);
    aux=max(abs(sign(v)-psi));
    if aux==0||(gama<=gama1)
        break
    end
    psi=sign(v);
    x=A'*psi;
	for i=1:p-1
        x=A'*x;
    end
    k=k+1;
    fin=norm(x,Inf)==x(j)||k>5;
end
for i=1:n
    x(i)=(1+(i-1)/(n-1))*(-1)^(i+1);
end
x=A*x;
for i=1:p-1
	x=A*x;
end
aux=2*norm(x,1)/(3*n);
if aux>gama
    gama=aux;
end

end

function gama=norm1pp(A,p,B)
%   gama=norm1pp(A,p,B)
%
%   Estimates the 1-norm of A^p*B, where A and B are square matrix.
%
%   Inputs:
%      A y B:  the input matrices
%      p:  power parameter
%
%   Outputs: 
%      gama: estimated 1-morm of A^p*B
%
%   References: 
%   [1] FORTRAN Codes for Estimating the One-Norm of a Real or Complex
%       Matrix, with Applications to Condition Estimation.
%       ACM Transactions on Mathematical Software, Vol. 14, No. 4, 
%       pp. 381-396, 1988 (Algorithm 4.1).
%
%
%   Revised version 1.0: 16/05/2013

n=size(A,1);
v=B*(ones(n,1)/n);
for i=1:p
    v=A*v;
end
gama=norm(v,1);
psi=sign(v);
x=A'*psi;
k=2;fin=0;
while fin==0
	[~,j]=max(abs(x));
    ej=zeros(n,1);
    ej(j)=1;
    v=B*ej;
    for i=1:p
        v=A*v;
    end
    gama1=gama;
    gama=norm(v,1);
    aux=max(max(abs(sign(v)-psi),abs(sign(v)+psi)));
    if aux==0||(gama<=gama1)
        break
    end
    psi=sign(v);
    x=B'*psi;
	for i=1:p
        x=A'*x;
    end
    k=k+1;
    fin=norm(x,Inf)==x(j)||k>5;
end
for i=1:n
    x(i)=(1+(i-1)/(n-1))*(-1)^(i+1);
end
x=B*x;
for i=1:p
	x=A*x;
end
aux=2*norm(x,1)/(3*n);
if aux>gama
    gama=aux;
end
end
