// Software    : Scilab >= 5.5
// Function    : partfrac()
// License     : CeCILL-C : http://www.cecill.info/licences
// Author      : Samuel GOUGEON
// Copyright(C): Samuel GOUGEON - 2015
// Version     : 2.0
// Date        : 2015-08-27
// Release     : http://fileexchange.scilab.org/toolboxes/451000
//
// History     :
//  2015-08-27: 2.0: d's format changed. Reference to pfss() added.
//  2015-08-21: 1.0: First release on FileExchange
//  2008-10-31     : Request on http://bugzilla.scilab.org/6270
//
// DISCLAIMER: This macro is provided as is, with NO WARRANTY of any kind.
//  Use it to your own risks, awareness, and responsability.

function [E, R, F, T] = partfrac(r)
// 
// Partial fraction decomposition of a rational in C
//
// SYNTAX
// ------
//            partfrac    // Displays this help
//  i       = partfrac(f)
// [i,r]    = partfrac(f)
// [i,r,d]  = partfrac(f)
// [i,r,d,t]= partfrac(f)
//
// PARAMETERS
// ----------
// f : Single rational = polynomial fraction with real or complex coefficients.
// i : Integer part = non-fractional part of f (single polynomial).
// r : Remainder of f: rational with degree(numerator) < degree(denominator)
//
// d : (3xN) matrix describing the N terms of the r's decomposition:
//     d(1,:): real or complex coefficient of numerators.
//     d(2,:): real or complex value of the considered poles.
//     d(3,:): integer>0: powers m of the denominators (x-pole)^m
//     d is such that sum( d(1,:)./ ((x-d(2,:)).^real(d(3,:))) ) - r == 0.
//     NOTE: real() must be applied to d(3,:) only because d is complex-encoded.
//           See http://bugzilla.scilab.org/14116
//
// t : 3-rows column of Texts. write(%io(2),t) displays i+d in a comprehensive
//     way. The format of coefficients is set with format().
//
// DESCRIPTION
// -----------
// From a single polynomial fraction f = p/q with coprime polynomials p and q, 
// partfrac(f) extracts 
//  * the non fractional part i of f: polynomial such that 
//    0 <= degree(i) <= degree(p) - degree(q): i = p - modulo(p,q)
//  * the fractional part of f, or remainder r of p/q : r = modulo(p,q)/q
//    such that degree(numer(r))< degree(q) and  f = i + r 
//  * the rational decomposition of r, aka partial fraction decomposition of f:
//    The 'vector' of elementary rationals c/(x-pole)^m where c and poles are 
//    decimal or complex numbers and m are the multiplicities of poles:
//    c = d(1,:), pole = d(2,:), m = d(3,:)
//
// In addition, partfrac() may return as text the literal expression t of the 
// whole decomposition of f. This form shows the factorized forms of denominators
// where poles that are multiple appear with their powers. write(%io(2), t) may 
// be used to display it.
//
// DEPENDENCY: If d or/and t is expected, polyroots() is required (=> See also)
//
// NOTE: pfss(f) yields a list mixing the non-fractional part (as last element)
//  with the decomposition, and where terms of the decomposition are not always
//  elementary: their denominator may be of order 2 even for real poles. See the
//  example.
//
// REFERENCE
// ---------
//  Comments, scoring and bug reports are welcome on
//  http://fileexchange.scilab.org/toolboxes/451000#new_comment
//
// SEE ALSO
// --------
//  pfss      : Partial (uncomplete) fraction decomposition of a linear system
//  modulo    : Remainder after polynomial division
//  polyroots : Multiplicities and values of polynomial multiple roots:
//              http://fileexchange.scilab.org/toolboxes/362000
//  pdiv_inc  : Polynomial division with increasing terms powers:
//              http://fileexchange.scilab.org/toolboxes/449000
//
// EXAMPLE
// -------
// x = poly(0,"x");
// f = (3-x+x^2-4*x^3+2*x^4) / ((x-1)*(x-2)^2)
// [E, R, F, T] = partfrac(f); 
// E, R, F, T
// D = F(1,:)./((x-F(2,:)).^real(F(3,:)))
// clean(R-sum(D))
// // Comparison with pfss():
// list2vec(pfss(f)).' 

// -------------------------------------------------------------------------
    fname = "partfrac"

    if argn(2)==0
        head_comments(fname)
        E = []
        R = []
        F = []
        Dtxt = []
        return
    end
    // CHECKING PARAMETERS
    // ===================
    if argn(2)>1
        msg = _("%s: Wrong number of input arguments: %d expected.\n")
        error(msprintf(msg, fname, 1))
        
    elseif typeof(r)~="rational"
        msg = _("%s: Wrong type for argument %d: Rational expected.\n")
        error(msprintf(msg, fname, 1))
   elseif size(r,"*")>1
        msg = _("%s: Wrong size for argument %d: Single component expected.\n")
        error(msprintf(msg, fname, 1))
    end
    
    // PROCESSING
    // ==========
    // Non-fractional part:
    N = numer(r)
    D = denom(r)
    R = pmodulo(N,D)
    if N==R
        E = 0
    else
        E = (N-R)/D
        E = numer(E)/coeff(denom(E))
    end
    if argn(1)==1
        return
    elseif argn(1)==2
        R = R/D
        return
    else
        if ~isdef("polyroots","n")
            msg = _("%s: ""polyroots()"" must be installed and loaded:\n\t http://fileexchange.scilab.org/toolboxes/362000")
            error(msprintf(msg, fname))
        end
    end

    if degree(D)==1
        R = R/D
        F = R
        Dtxt = sci2exp(D)
        return
    end
        
    // Decomposition
    // -------------
    rootsmul = polyroots(D)
    // Expands the series of roots
    eroots = [] // expanded roots
    for i = 1:size(rootsmul,1)
        r = rootsmul(i,:)
        m = real(r(2))
        eroots = [eroots ones(1,m)*r(1)]
    end
    // Replicates the series
    Eroots = ones(length(eroots),1)*eroots
    
    x = poly(0,varn(D))
    k = 1
    for i = 1:size(rootsmul,1)
        v = rootsmul(i,1)
        m = real(rootsmul(i,2))
        for j=1:m
            n = find(eroots==v)
            Eroots(k,n(1:j)) = %nan // Cancels the term simplified with denominator
            k = k+1
        end
    end
    // Converts series of active roots in (x-root) polynomials 
    for i = 1:size(Eroots,1)
        j = find(~isnan(Eroots(i,:)))
        Eroots(i,j) = x-Eroots(i,j)
        j = find(isnan(Eroots(i,:)))
        Eroots(i,j) = 1
    end
    // Converts rows of (x-root) polynomials into full polynomials
    P = prod(Eroots,"c")
            
    // Computes coefficients of the decomposition
    d = degree(D);
    cR = coeff(R,0:d-1).'
    CP = coeff(P,0:d-1).'
    c = CP\cR   //

    // Builds F and T
    T = []
    F = []
    k = 1
    for i = 1:size(rootsmul,1)
        r = rootsmul(i,:)
        for m = 1:real(r(2))
            tmp = x-r(1)
            F = [F [c(k) r(1) m].']
            Ttmp = sci2exp(tmp)
            if m>1, Ttmp = msprintf("("+Ttmp+")^%d",m), end
            T = [ T [ string(c(k)) ; "" ; Ttmp ]]
            k = k+1
        end
    end
    
    // Output formatings
    // =================
    R = R/D
    // Formates the non-fractional part in text
    tmp = string(E)
    L = length(tmp)
    tmp(2) = tmp(2)+part(" ",1:(L(1)-L(2)))
    tmp(3) = part(" ",1:max(L))
    // Prepends it to T
    T = [ tmp justify(T,"c")]
    
    // Formates the remainder decomposition in text
    L = length(T(1,2:$))
    b = strsplit(part("-",ones(1:sum(L))), cumsum(L(1:$-1)))
    T(2,2:$) = b'
    // Concatenates all T components with "+"
    tmp = strcat(T(1,:),"   ")
    tmp(2) = strcat(T(2,:)," + ")
    tmp(3) = strcat(T(3,:),"   ")
    T = tmp
endfunction

