#
# SparseDistributedMultivariatePolynomial(R,X)
#
# A sparse distributed multivariate polynomial domain R[X] where
# R is a ring and E is an exponent vector.
#
# The representation is a Maple list of elements from R and X, i.e.
#
#	[c[1], e[1], c[2], e[2], ..., c[n], e[n] ]
#
# Where the c[i] are from R and the e[i] are from X
# R must be a ring and X is an ordered abelian monoid hence defining an
# ordering on the terms: e[i] > e[j] for i > j
#
# Author MBM: 1991
#
SparseDistributedMultivariatePolynomial := proc()
local P,R,S,env; option remember;

    R := args[1];
    S := args[2];

    P := DistributedMultivariatePolynomial(R,S);
    P[DomainName] := SparseDistributedMultivariatePolynomial;

    env := ['C' = R, 'D' = P, 'E' = S];
    # Rep [ c[1], e[1], ..., c[n], e[n] ] where c[k] in R and e[k] in S

    P[0] := [R[0],S[0]];
    P[1] := [R[1],S[0]];

    P[`+`] := subs(env, proc(a) `SDMP/Add`(D,args) end);
    P[`*`] := subs(env, proc(a) `SDMP/Mul`(D,args) end);

    P[Lcoeff] := proc(a,t) if nargs=2 then t := a[2] fi; a[1] end;
    P[Lterm] := proc(a) [a[1], a[2]] end;

    P[Coeff] := subs(env, proc(a,e) local i;
        for i from 2 by 2 to nops(a) do
    	if a[i] = e then RETURN( a[i-1] ) fi
        od;
        C[0]
    end);

    P[Degree] := subs(env, proc(a,n) local d,i,j,m;
        if not member(n,D[Variables],i) then ERROR(`bad variable`,n) fi;
        m := 0;
        for j from 2 by 2 to nops(a) do
    	d := E[Degree](a[j],i);
    	if d > m then m := d fi;
        od;
        m
    end);

    P[Type] := subs(env, proc(a) local i, n;
        if not type(a,list) then RETURN(false) fi;
        n := nops(a); if irem(n,2)=1 then RETURN(false) fi;
        for i from 1 by 2 to n-1 do
            if not C[Type](a[i]) or not E[Type](a[i+1]) then RETURN(false) fi;
        od;
        true
    end);

    P[TermLess] := subs(env, proc(a,b) evalb(E[`<>=`](a[2],b[2]) = 1) end);
    P[List2Poly] := subs(env, proc(x) `SDMP/List2Poly`(D,x) end);
    P[Poly2List] := proc(x) local i;
    	[seq( [x[2*i-1],x[2*i]], i=1..nops(x)/2 )]
        end;

    if hasCategory(R,IntegralDomain) then
    	# implementations of PseudoRem, Res
    	# implementation of Unit, Div
    fi;

    if hasCategory(R,GcdDomain) then
    	# implementation of Gcd
    fi;

    op(P)
end:

macro( t = `SDMP/Temporary` );
`SDMP/List2Poly` := proc(P,x) local c,e,i,j,k,n,r,t,C,E;

    r := sort(x,P[TermLess]);
    n := nops(r);

    C := P[CoefficientRing];
    E := P[ExponentVector];

    i := 1;
    k := 1;
    while i <= n do
        c := r[i]; e := c[2]; c := c[1];
        for i from i+1 to n while E[`=`](e,r[i][2]) do 
        c := C[`+`](c,r[i][1])
        od;
    if C[`<>`](c,C[0]) then t[k] := c,e; k := k+1 fi;
    od;
    r := [ seq(t[j], j=1..k-1) ];
    if r = [] then P[0] else r fi

end:

`SDMP/Add` := proc(D,a,b) local i,j,k,l,m,n,r,comp,plus,equal,zero;

    if nargs = 2 then RETURN(a) fi;
    if nargs > 3 then RETURN( procname(D,procname(D,a,b),args[4..nargs]) ) fi;
    if a=D[0] then RETURN(b) fi:
    if b=D[0] then RETURN(a) fi:

    plus := D[CoefficientRing][`+`];
    comp := D[ExponentVector][`<>=`];
    zero := D[CoefficientRing][0]; 
    equal := D[CoefficientRing][`=`];

    m := nops(a); n := nops(b); i := 1; j := 1;
    for k by 2 while i < m and j < n do
        r := comp(a[i+1],b[j+1]);
        if r = 1 then
            t[k] := a[i]; t[k+1] := a[i+1];
            i := i + 2;
        elif r = -1 then
            t[k] := b[j]; t[k+1] := b[j+1];
            j := j + 2;
        else
            t[k] := plus(a[i],b[j]);
            if equal(t[k],zero) then k := k - 2
            else t[k+1] := a[i+1];
            fi;
            i := i + 2; j := j + 2;
        fi
    od;

    r := [ seq(t[l], l=1..k-1), a[i..m], b[j..n] ];
    if r = [] then D[0] else r fi;

end:

`SDMP/Mul` := proc(D,a,b) local i,j,k,m,n,c,e,r,mult,plus;

    if nargs = 2 then RETURN(a) fi;
    if nargs > 3 then RETURN(procname(D,procname(D,a,b),args[4..nargs])) fi;
    if type(a,integer) then RETURN(procname(D,D[Coerce](a),b)) fi;

    mult := D[CoefficientRing][`*`];
    plus := D[ExponentVector][`+`];

    m := nops(a);
    n := nops(b);
    
    r := [];
    for i by 2 to m do
        c := a[i];
        e := a[i+1];
        for j by 2 to n do t[j] := mult(c,b[j]); t[j+1] := plus(e,b[j+1]) od;
        c := [ seq(t[k], k=1..n) ];
        r := D[`+`](r,c);
    od;
    if r = [] then D[0] else r fi

end:

save `SDMP.m`;
quit
