// Software    : Scilab >= 5.5
// Function    : partfrac()
// License     : CeCILL-C : http://www.cecill.info/licences
// Author      : Samuel GOUGEON
// Copyright(C): Samuel GOUGEON - 2018
// Version     : 3.0
// Date        : 2015-08-27
// Release     : http://fileexchange.scilab.org/toolboxes/451000
//
// History     :
//  2018-11-01: 1.2: denom() and numer() replaced (deprecated in Scilab 6.0)
//                   Bug fixed when E is a scalar.
//                   Inconsistent notations fixed in the help.
//  2015-08-27: 1.1: 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 [I, R, F, T] = partfrac(r)
// 
// Partial fraction decomposition of a rational in C
//
// SYNTAX
// ------
//               partfrac    // Displays this help
//  I          = partfrac(r)
// [I, R]      = partfrac(r)
// [I, R, F]   = partfrac(r)
// [I, R, F, T]= partfrac(r)
//
// PARAMETERS
// ----------
// r : Single rational = polynomial fraction with real or complex coefficients.
// I : Integer part = non-fractional part of f (single polynomial).
// R : Remainder of r: rational with degree(numerator) < degree(denominator)
//
// F : (3xN) matrix describing the N terms of the r's decomposition:
//     F(1,:): real or complex coefficient of numerators.
//     F(2,:): real or complex value of the considered poles.
//     F(3,:): integer>0: powers m of the denominators (x-pole)^m
//     F is such that sum( F(1,:)./ ((x-F(2,:)).^F(3,:)) ) - r == 0.
//
// 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 r = 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 r, or remainder R of p/q : R = modulo(p,q)/q
//    such that degree(r.num))< degree(q) and  r = I + R
//  * the rational decomposition of R, aka partial fraction decomposition of r:
//    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 = F(1,:), pole = F(2,:), m = d(3,:)
//
// In addition, partfrac() may return as text the literal expression T of the 
// whole decomposition of r. 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 F or/and T is expected, polyroots() is required (=> See also)
//
// NOTE: pfss(r) 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");
// r = (3-x+x^2-4*x^3+2*x^4) / ((x-1)*(x-2)^2)
// [I, R, F, T] = partfrac(r); 
// I, R, F, T
// D = F(1,:)./((x-F(2,:)).^F(3,:))
// clean(R-sum(D))
// // Comparison with pfss():
// list2vec(pfss(r)).' 

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

    if argn(2)==0
        head_comments(fname)
        I = []
        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 = r.num
    D = r.den
    R = pmodulo(N,D)
    if N==R
        I = 0
    else
        I = (N-R)/D
        I = I.num / coeff(I.den)
    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(I)
    if type(I)==1 then
        tmp = ["" ; tmp]
    end
    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

