// Software    : Scilab >= 5.5
// Function    : pdiv_inc()
// Author      : Samuel GOUGEON
// Copyright(C): Samuel GOUGEON - 2015
// License     : CeCILL-C : http://www.cecill.info/licences
// Version     : 1.0
// Date        : 2015-08-20
// Release     : http://fileexchange.scilab.org/toolboxes/449000
// History     :
//   2015-08-20 : First release on FileExchange
//
// DISCLAIMER: This macro is provided as is, with NO WARRANTY of any kind.
//  Use it to your own risks, awareness, and responsability.
//
function [q, r] = pdiv_inc(p, d)
//
// Polynomial division by increasing terms powers
//
// SYNTAX
// ------
//          pdiv_inc   // Displays this help
// [q, r] = pdiv_inc(p, d)
// [q, r] = pdiv_inc(f)
//
// PARAMETERS
// ----------
// p,d,q,r: Single polynomials with real or complex coefficients.
//          p : polynomial to be divided
//          d : divisor. Its first non-null coefficient of lowest power must 
//              be of power equal or lower than p's one.
//          q : quotient, with deg(q) <= deg(p)
//          r : remainder, with deg(r) <= deg(p)+deg(d)
//          such that: p = q*d + r. p and d are not necessarily coprime.
//
// f : Single rational, with coprime polynomials p = numer(f) and d = denom(f)
//
// DESCRIPTION
// -----------
// pdiv_inc(p, d) performs the Euclidian division p/d with polynomial terms 
// of p and d sorted by increasing powers. It returns the quotient q and the 
// remainder r, such that p = q*d + r. 
// Detailled example:
//   2 - x + 3x^2               | 1 + 2x^2
//   ---------------------------|------------
//   2     + 4x^2               | 2 - x - x^2
//     - x -  x^2               |
//     - x        - 2x^3        |
//         -  x^2 + 2x^3        |
//         -  x^2        - 2x^4 |
//                  2x^3 + 2x^4 |
//  Hence, (2 - x + 3x^2) = (1 + 2x^2)*(2 - x - x^2) +  (2x^3 + 2x^4)
//
// NOTES:
//  * If d==0, q=0 and r=p are returned whatever is p.
//  * if both p and d are limited to a decimal, there are set to p+0*x and d+0*x
//
// WARNING: If f is specified explicitly through a p/d ratio with p and d not 
//  coprime (with some common factors), these factors are first simplified, and
//  the p = q*d + r relationshift holds only with simplified p=>sp and d=>sd. 
//  Example: p = poly([1 2 2 3 3],'x'), d = poly([2 3 4],'x'), p/d
//          [q,r] = pdiv_inc(p/d), p-q*d-r, f = p/d; numer(f)-q*denom(f)-r
//
// REFERENCE
// ---------
//  Comments, scoring and bug reports are welcome on
//  http://fileexchange.scilab.org/toolboxes/4490000#new_comment
// 
// SEE ALSO
// --------
//  pdiv : polynomial division (by decreasing terms powers)
//
// EXAMPLES
// --------
// x = poly(0,"x");
//// 1) Simple example
// p =  1 + 3*x + 2*x^2 - 7*x^3
// d = -1 +   x - 2*x^2
// [q,r] = pdiv_inc(p,d),  p-q*d-r   // => 0
//
//// 2) With some non-leading null coefficients in p and/or d
// p = 1 + x +      3*x^3
// d = 1 +    2*x^2
// [q,r] = pdiv_inc(p,d),  p-q*d-r   // => 0
//
//// 3) With p not starting with x^0
// p =    x -   x^2 + 2*x^3
// d = -1   + 2*x^2
// [q,r] = pdiv_inc(p,d),  p-q*d-r   // => 0
//
//// 4) With p and d not starting with x^0
// p =       x^2 + 2*x^3
// d =  -x + 2*x^2 - x^3
// [q,r] = pdiv_inc(p,d),  p-q*d-r   // => 0
//
//// 5) with complex coefficients:
// p = -1 + (1+%i)*x - %i*x^2 ; 
// d = %i - x + (1-%i)*x^2    ;
// [q,r] = pdiv_inc(p,d);  
// sci2exp(q), sci2exp(r), sci2exp(p-q*d-r)   // => 0
//
//// 6) d must start at a degree lower or equal to p's starting degree
// p =  2 -  x + 3*x^2
// d =     2*x -   x^2
// [q,r] = pdiv_inc(p,d)  //=> Error
//
//// 7) Random example
// p = poly(grand(1,5,"uin",-9,9), "x", "coeff")
// d = poly(grand(1,7,"uin",-9,9), "x", "coeff")
// [q,r] = pdiv_inc(p,d),  clean(p-q*d-r)   // => 0
//

// --------------------------------------------------------------------------

    fname = "pdiv_inc"
    
    // DISPLAYING the HELP
    // -------------------
    if argn(2)==0
        head_comments(fname)
        q = []
        r = []
        return
    end
    
    // CHECKING PARAMETERS
    // -------------------
    if argn(2)>2
        msg = _("%s: Wrong number of input arguments: %d or %d expected.\n")
        error(msprintf(msg, fname, 1, 2))
    end
    
    if argn(2)==1
        if typeof(p)~="rational"
            msg = _("%s: Wrong type for argument %d: Rational fraction expected.\n")
            error(msprintf(msg, fname, 1))
        else
            d = denom(p)
            p = numer(p)
        end
    else
        if ~or(type(p)==[1 2])
            msg = _("%s: Wrong type for argument %d: Polynomial expected.\n")
            error(msprintf(msg, fname, 1))
        end
        if ~or(type(d)==[1 2])
            msg = _("%s: Wrong type for argument %d: Polynomial expected.\n")
            error(msprintf(msg, fname, 2))
        end
    end
    // In case of vector or matrix, only the first component will be processed:
    p = p(1)    
    d = d(1)

    // In case of scalar(s) p and/or d, getting/homogenezing the unknown's name
    v = ""
    if type(d)==2
        v = varn(d)
    else
        d = poly(d,"x","coeff")
    end
    if type(p)==2
        v = varn(p)
    else
        p = poly(p,"x","coeff")
    end
    if v=="", v = "x", end
    p = varn(p, v)
    d = varn(d, v)
    
    // PROCESSING
    // ==========
    // d==0 case
    if d==0
        q = 0
        r = p
        return
    end
    
    cop = coeff(p)
    cod = coeff(d)
    
    // Are there any low order null coefficients?
    ip = find(cop~=0)
    id = find(cod~=0)
    if id(1)>ip(1)
        msg = _("%s: p/d: while p starts with x^%d, d must start with x^(n<=%d)\n")
        error(msprintf(msg, fname,ip(1)-1, ip(1)-1))
    end
    xfact = poly(0,v)^(id(1)-1)
    cop(1:(id(1)-1)) = []
    cod(1:(id(1)-1)) = []
    degp = degree(p)-id(1)+1
    degd = degree(d)-id(1)+1
    cop = [cop zeros(1,degd)]
    
    // Euclidian algorithm
    i = 1
    F = []
    while i<=degp+1
        f = 0
        if cop(i)~=0
            f = cop(i)/cod(1)
        end
        F(i) = f
        cop(i:i+degd) = cop(i:i+degd) - f*cod
        i = i+1
    end
    
    // Final formatting
    q = poly(F, v,"coeff")
    r = poly(cop,v,"coeff")*xfact
endfunction
