function [F,cost] = expmC(A)
% OUR ALGORITHM 
%expmC   Matrix exponential.
%   expmC(H) is the matrix exponential of skew-symmetric matrix H. expmC computes exponentials using
%   a scaling and squaring algorithm with a Chebyshev approximation evaluated 
%   using a new decomposition.
%   
%   The choice of order and scaling is based solely on norm(A,1)
% % Changelog:
%   V1.18 Nov 2020 (improved degree 12 decomposition)
A=1i*A;
[m_vals, theta,cost_table] = cheb_parameter; % 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)
            F = cheb_approximant_of_degree(m_vals(method_selector));
            break;
        end
    end
else
    [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);
    F = cheb_approximant_of_degree(m_vals(method_selector ));
    for i = 1:s
        F = F*F;  % Squaring
    end   
end
cost               = s + cost_table(method_selector);
% End of expmC

%%%%Nested Functions%%%%
    function [m_vals, theta, cost_table] = cheb_parameter
        %    initialize M_VALS and THETA accordingly.
                m_vals     = [ 2 4 8 12 18];
                %product of  matrices
                cost_table = [ 1 2 3 4 5 ]; 
                % theta_m .
                theta = [ 1.38e-05  % m_vals = 2
                          2.92e-03  % m_vals = 4
                          0.1295    % m_vals = 8
                          0.636     % m_vals =12 
                          2.212 ];  % m_vals = 18
                                     
    end

    function E = cheb_approximant_of_degree(order)
%% Improved Paterson Stockmeyer scheme, with Pade at end point
        I = eye(size(A,1));
        
        if (order>=2)
            A2 = A*A; 
        end
        if (order>8)
            A3 = A*A2;
        end
        switch order
            case 2  
                alp_0 = 0.999999999999999999999811107325001498863;
                alp_1 = -0.999999999976195000000188892674999250568*1i;
                alp_2 = -0.499999999992065000000047223168749850113;
                    E = alp_0*I + alp_1*A+alp_2*A2;
            case 4  
                
                x1 = 0.16666657785001893215842467829973122728*1i;
                x2 = 0.04166664890333648869312416038897725009;
                alp_0 = 0.99999999999999999997309614471715549922;
                alp_1 = -0.99999999999981067844714105106755532083*1i;
                alp_2 = -0.49999999999994320353145192968957453728;               

                   
                    A4 = A2*(x1*A + x2*A2);
                    E  = alp_0*I + alp_1*A + alp_2*A2 + A4;
            case 8
                % P_8(A) with theta=0.1295
                x1 = 431/4000;    
                x2 = -0.02693906873598870733*1i;
                x3 = 0.66321004441662438593*1i;
                x4 = 0.54960853911436015786*1i;
                x5 = 0.16200952846773660904;
                x6 = -0.01417981805211804396*1i;
                x7 = -0.03415953916892111403;
                alp_0 = 0.99999999999999999928;
                alp_1 = -0.99999999999999233987*1i;
                alp_2 = -0.13549409636220703066;
                % Matrix products
                A4 = A2*(x1*A + x2*A2);
                A8 = (x3*A2 + A4)*(x4*I + x5*A + x6*A2 + x7*A4) ;
                E = alp_0*I + alp_1*A + alp_2*A2 + A8;
            case 12
                % P_12(A) with theta=0.636
                a01 = -6.26756985350202252845;
                a11 = 2.52179694712098096140*1i;
                a21 = 0.05786296656487001838;
                a31 = -0.07766686408071870344*1i;
                a02 = 0;
                a12 = 1.41183797496250375498*1i;
                a22 = 0;
                a32 = -0.00866935318616372016*1i;
                a03 = 2.69584306915332564689; 
                a13 = -1.35910926168869260391*1i; 
                a23 = -0.09896214548845831754; 
                a33 = 0.01596479463299466666*1i; 
                a04 = 0;
                a14 = 0.13340427306445612526*1i; 
                a24 = 0.02022602029818310774; 
                a34 = -0.00674638241111650999*1i;
                
                B1 = a01*I+a11*A+a21*A2+a31*A3;
                B2 = a02*I+a12*A+a22*A2+a32*A3;
                B3 = a03*I+a13*A+a23*A2+a33*A3;
                B4 = a04*I+a14*A+a24*A2+a34*A3;
                % Matrix products
                A6 = B3 + B4*B4;
                E  = B1 + (B2 + A6)*A6;
            case 18
                % P_18(A) with theta=2.212
                a01 = 0; 
                a11 = 3/25;
                a21 = -0.00877476096879703859*1i;
                a31 = -0.00097848453523780954;
                b01 = 0;
                b11 = -0.66040840760771318751*1i;
                b21 = -1.09302278471564897987;
                b31 = 0.25377155817710873323*1i;
                b61 = 0.00054374267434731225;
                b02 = -2.58175430371188142440;
                b12 = -1.73033278310812419209*1i;
                b22 = -0.07673476833423340755;
                b32 = -0.00261502969893897079*1i;
                b62 = -0.00003400011993049304;
                b03 = 2.92377758396553673559;
                b13 = 1.44513300347488268510*1i;
                b23 = 0.12408183566550450221;
                b33 = -0.01957157093642723948*1i;
                b63 = 0.00002425253007433925;
                b04 = 0;
                b14 = 0;
                b24 = -0.123953695858283131480*1i;
                b34 = -0.011202694841085592373;
                b64 = -0.000012367240538259896*1i;

                % Matrix products
                A6 = A3*A3;
                B1 = (a01*I + a11*A + a21*A2 + a31*A3);
                B2 = (b01*I + b11*A + b21*A2 + b31*A3 + b61*A6); 
                B3 = (b02*I + b12*A + b22*A2 + b32*A3 + b62*A6);
                B4 = (b03*I + b13*A + b23*A2 + b33*A3 + b63*A6);
                B5 = (b04*I + b14*A + b24*A2 + b34*A3 + b64*A6);
                
                A9 = B1*B5 + B4;
                E  = B2 + (B3 + A9)*A9;
                
        end
    end

%%%%Nested Functions%%%%
end
