function [F, cost] = cosmsinmC(A)
%%  OUR ALGORITHM 
%cosmsinmC matrix cosine and sine simultaneously.
%  cosmsinmC(X) is matrix cosine and sine of X simultaneously.  cosmsinmC computes 
%  the matrix cosine sine functions using a scaling and squaring algorithm with 
%   a Chebyshev approximation evaluatedusing a new decomposition.   
%   The choice of order and scaling is based solely on norm(A,1)
% % Changelog:
%   V1.18 Nov 2020 (improved degree 24 decomposition)
   
      
   
A=1i*A;
[m_vals, theta] = cosmsinmT_new; % Initialization
normA = norm(A,1);

if normA <= theta(end)
    % no scaling and squaring is required.
    s = 0;
    for method_selector = 1:length(m_vals)
        if normA <= theta(method_selector)
            [C, S] = taylor_approximant_of_degree(m_vals(method_selector));
            break;
        end
    end
else
% Else we need to scale
    [t, s] = log2(normA/theta(end));
    s = s - (t == 0.5); % adjust s if normA/theta(end) is a power of 2.
    A = A/2^s;    % Scaling
    method_selector = length(m_vals);
    [C, S] = taylor_approximant_of_degree(m_vals(method_selector ));
   n = size(A);
I = eye(n); 
for k = 1:s
    Cold = C;
    C = I - 2*S^2;
    S = 2*Cold*S;
end  
end

% if method_selector==1 || method_selector==2
     cost               = 2*s + method_selector+2;
%     cost               = 2*s + method_selector+2;
% else
% cost               = 2*s + method_selector +3;
% end

% End of cosmsinm
F=C-1i*S;
%%%%Nested Functions%%%%
    function [m_vals, theta] = cosmsinmT_new

        %    initialize M_VALS and THETA accordingly.
       
                    m_vals = [5 8 9 16 24 25 ];
                cost_table = [3 4 5  6  7  8 ];
%               cost_table = [3 4 6 7];
                theta = [ 1.17e-2  % m_vals = 5, cost = 3
%                           0.1295  % m_vals = 8                         
                          0.068  % m_vals = 8, cost = 4                         
                          0.214  % m_vals = 9, cost = 5                         
                          0.7563  % m_vals = 16, cost = 6
%                          1.5867  % m_vals = 16
                          2.1556 % m_vals = 24, cost = 7
                          4.5743 % m_vals = 25, cost = 8
                         ];
           
    end

    function [Tc, Ts] = taylor_approximant_of_degree(order)
%% Improved Paterson Stockmeyer scheme, New Chebyshev Polynomials
        I = eye(size(A,1));
%         
        if (order>=4)
            B  = A*A;
            B2 = B*B;
        end
        switch order
            case 5  

               alpha_0 = 0.99999999999999988866;
                alpha_1 = -0.49999999998536031183;
                alpha_2 = 0.04166638147997997916;
                z0 = 0.99999999999999994433;
                z1 =-0.16666666666341340086;
                z2 = 0.00833328580219952161;
                   
                   Tc = alpha_0*I + alpha_1*B + alpha_2*B2;
                    Ts = A*(z0*I + z1*B + z2*B2);
                   
%                 alpha_0 = 0.999999999999999999973096144717155499;
%                 alpha_1 = -0.49999999999994320353145192968957453;
%                 alpha_2 = 0.04166664890333648869312416038897725;
%                 z0 = 0.99999999999981067844714105106;
%                 z1 = -0.1666665778500189321584246782;
%                    
%                    Tc = alpha_0*I + alpha_1*B + alpha_2*B2;
%                     Ts = A*(z0*I + z1*B);
                   
            case 8
                    x1 = -0.001388888719399421184316693373311;
                    x2 = 0.00002479003614491668409382551153693;
                    alpha_0 = 0.999999999999999999286304506154937745;
                    alpha_1 = -0.499999999999997872108440071891158249;
                    alpha_2 = 0.041666666665651566154678053376007374654;
                    z0 = 0.857217689470640124669891442563445;
                    z1 = -0.095275511395900472561711394270672;
                    z2 = 0.002384069087305688504634128482644273;
                    z3 = 0.142782310529352215307883173172932979;


                    B4 = B2*(x1*B + x2*B2);
                    Tc = alpha_0*I + alpha_1*B + alpha_2*B2 + B4;
                    Ts = A*(z0*I + z1*B + z2*B2 + z3*Tc);

           case 9
                    alpha_0 = 0.99999999999999989168;
                    alpha_1 = -0.49999999999988173685;
                    alpha_2 = 0.04166666664600636231;
                    alpha_3 = -0.00138888762558264513;
                    alpha_4 = 0.00002477005498155486;
                    
                   
                    z0 = 0.999999999999999945837;
                    z1 =-0.166666666666643012068;
                    z2 = 0.008333333330440664914;
                    z3 =-0.000198412554024823435;
                    z4 = 2.75257852630876250884*10^(-6); 


        B3 = B2*B;  B4 = B3*B;
        Tc = alpha_0*I + alpha_1*B + alpha_2*B2 + alpha_3*B3 + alpha_4*B4;
        Ts = A*(z0*I + z1*B + z2*B2 + z3*B3 + z4*B4);

                   
            case 16
 % lower values of ||_||_1
             
            x1 = 1/100; 
            x2 = -0.00008035854055477845; 
            x3 = -0.10743065643419630630;
            x4 = -0.12491372919298427513;
            x5 = 0.00130085397953037838;
            x6 = -0.00001633763177694857;
            x7 = 7.13215089463286614820*10^(-6);
            alpha_0=0.99999999999999999530;
            alpha_1=-0.49999999999999969795;
            alpha_2=0.028247102741817734721;
            z0 = 33/50;
            z1 = 0.00333333333335438849;
            z2 = -0.00583333333345309522;
            z3 = 0.02773310749258735833;
            z4 = 0.33999999999999886261;
            z5 =-0.00034915267907803119;
            z6 = 4.19573036995827807213*10^(-6);
            z7 =-2.63931697420854364428*10^(-6);
            z8 = -3.00240279002259730782*10^(-6);
           
            B4=B2*(x1*B+x2*B2);
            B8=(x3*B2+B4)*(x4*I+x5*B+x6*B2+x7*B4); 
            Tc=alpha_0*I+alpha_1*B+alpha_2*B2+B8;
            c24=(z5*I+z5*B+z6*B2+z7*B4+z8*Tc)*B4;
            Ts=A*(z0*I+z1*B+z2*B2+z3*B4+z4*Tc+c24);
            
            
          
            case 24

%                 a01 = 0.38272620931352327385;
%                 a11 = -0.03337071139879721851;
%                 a21 = 0.01894164507143601070;
%                 a31 = -0.00013131304047131598;
%                 a02 = 0;
%                 a12 =-0.54235659842328961975;
%                 a22 = 679/100000;
%                 a32 = -0.00002902999756981724;
%                 a03 = 0.78566773555140770914;
%                 a13 = -0.02578520551577453856;
%                 a23 = 0.00019815665089300452;
%                 a33 = -1.10083330495602029332*10^(-6);
%                 a04 = 0;
%                 a14 = -0.03931944346958836562;
%                 a24 = 0.00017839382197658767;
%                 a34 = -1.06908694221941432625*10^(-6);
% 
%                 z0 = -0.21521153162788200934;
%                 z1 = -0.06028551052776780907;
%                 z2 = 0.00035899274257347367;
%                 z3 = -9.88861032740280754157*10^(-6);
%                 z4 = 1.36330829895850807447;
%                 z5 = 2.75831913749065306033*10^(-14);
%                 z6 = 0.14398058430437755744;
%                 z7 = 0.00082817888777282902;
%                 z8 = 3.93755124535163288400*10^(-6);
%                 z9 = 1.67495949351675674068*10^(-8);
%                 z10 = 0.00015732251269792353;
%                 z11 = 1.07352968693173803803*10^(-15);

                a01 = 0.39272620931352327385;
                a11 =-0.08760637124112618048;
                a21 = 0.01962064507143601071;
                a31 =-0.00013421604022829771;
                a02 = 1/5;
                a12 =-0.54235659842328961975;
                a22 = 679/100000;
                a32 =-0.00002902999756981724;
                a03 = 0.68566773555140770915;
                a13 =-0.02578520551577453856;
                a23 = 0.00019815665089300452;
                a33 =-1.10083330495602029332*10^(-6);
                a04 = 0;
                a14 =-0.03931944346958836562;
                a24 = 0.00017839382197658767;
                a34 =-1.06908694221941432625*10^(-6);
                
z0= -0.01238438326981811663;
z1= -0.06180067679127220638;
z2=  0.00046275599640408615;
z3= -9.92990416300441584763*10^(-6);
z4=  1.26307934615308708610;
z5=  9.10439014880980346565*10^(-15);
z6=  0.14610549096048524519; 
z7=  0.00087697762149660844;
z8 = 4.12092186281469998191*10^(-6); 
z9 = 2.23743615053828476204*10^(-8);
z10 = 0.00033015662857238333; 
z11 = -2.405371071766852323329*10^(-7);
%                 

               B3=B2*B;
                
                K1 = a01*I+a11*B+a21*B2+a31*B3;
                K2 = a02*I+a12*B+a22*B2+a32*B3;
                K3 = a03*I+a13*B+a23*B2+a33*B3;
                K4 = a04*I+a14*B+a24*B2+a34*B3;
                % Matrix products
                B6 = K3 + K4*K4;
                Tc  = K1 + (K2 + B6)*B6;
                c48=(z6*I+z7*B+z8*B2+z9*B3+z10*B6+z11*Tc)*Tc;
                Ts=A*(z0*I+z1*B+z2*B2+z3*B3+z4*B6+z5*Tc+c48);

                
            case 25

                a01 = 0.39272620931352327385;
                a11 =-0.08760637124112618048;
                a21 = 0.01962064507143601071;
                a31 =-0.00013421604022829771;
                a02 = 1/5;
                a12 =-0.54235659842328961975;
                a22 = 679/100000;
                a32 =-0.00002902999756981724;
                a03 = 0.68566773555140770915;
                a13 =-0.02578520551577453856;
                a23 = 0.00019815665089300452;
                a33 =-1.10083330495602029332*10^(-6);
                a04 = 0;
                a14 =-0.03931944346958836562;
                a24 = 0.00017839382197658767;
                a34 =-1.06908694221941432625*10^(-6);


z0 = 2.852476503968736096649749199992;
z1 =-0.238389229843545097977339008724;
z2 = 0.01254735251131974478084283639006;
z3 =-0.0000318498423383495493918242125653;
z4 =-7.914119343579328111106336825077;
z5 =-0.4558495682876669453813535993974;
z6 =-2.349447231105943100695796743436;
z7 =-0.3431565053409967548586045129680;
z8 = 0.003795294092950146105075020944040;
z9 =-0.000015093120022447184555877288510115;
z10 = -17/1000;
z11 = 7.681457951181004729457438110146*10^(-9);
z12 =-2.7189617581026327876471588546302*10^(-11);
z13 = 0.45584956828766694538135359939740;
% 
% 
% z0 = 2.875059309042944078861115903401;
% z1 =-0.246418557581221659697383264703;
% z2 = 0.02070444297769069821348864052689;
% z3 =-0.0000878275806725757863404120786778;
% z4 =-8.323853723290007114120294241698;
% z5 =-0.8456220325712473444984179214171;
% z6 =-2.081180360179208201149409926722;
% z7 =-0.5546910175966516775517948108090;
% z8 = 0.006441849125435657015079887690553;
% z9 =-0.00002640821371338120179438694620187;
% z10 = -9/1000;
% z11 = 7.681457951181004729457438110146*10^(-9);
% z12 =-2.7189617581026327876471588546302*10^(-11);
% z13 = 0.8456220325712473444984179214171;


               B3=B2*B;
                
                K1 = a01*I+a11*B+a21*B2+a31*B3;
                K2 = a02*I+a12*B+a22*B2+a32*B3;
                K3 = a03*I+a13*B+a23*B2+a33*B3;
                K4 = a04*I+a14*B+a24*B2+a34*B3;
                % Matrix products
                B6 = K3 + K4*K4;
                Tc  = K1 + (K2 + B6)*B6;
                D5=B2*(z11*B2+z12*B3);
                c24=(z6*I+z7*B+z8*B2+z9*B3+D5+z13*B6)*(B6+z10*B);
                Ts=A*(z0*I+z1*B+z2*B2+z3*B3+z4*D5+z5*Tc+c24);
                
        end
    end

%%%%Nested Functions%%%%
end
