#	        
# --> MapleMultivariatePolynomial(R, X:list(name));
#
# Creates a Maple multivariate polynomials domain for R[x] where R is
# the integers or rationals.  The Maple sum of products representation
# is used for efficiency. This works over all coefficient domains, but
# the efficiency is only reached over R, Z and Zmod.
#
# Author DG: 1992
#

macro( MMPGcd = `MMP/Gcd` );
macro( MMPPrem = `MMP/Prem` );

MapleMultivariatePolynomial := proc() 
    local X,P,R,env,n; option remember;

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

    P := DistributedMultivariatePolynomial(R,MapleExponentVector(X, 'plex'));

    P[DomainName] := MapleDistributedMultivariatePolynomial;
    env := ['C' = R, 'D' = P];

    # Rep := Maple Univariate Polynomial

    P[Type] := subs(env, proc(a) local x;
       if type(a,{`*`,`+`}) then
          for x in a do if not D[Type](x) then RETURN(false) fi od;
          true
       elif type(a,name) and member(a, D[Variables]) then true
       elif type(a,anything^posint) then D[Type](op(1,a))
       else C[Type](a)
       fi
    end);

    P[List2Poly] := subs(env, proc(a) local m, s, c;
       c := 0:
       for m in convert(map(t -> op(2,t), a),set) do
          s:=map(t -> op(1,t), select( (t,m) -> op(2,t)=m, a, m));
          if nops(s)=1 then c := c+s[1]*m else c := c+C[`+`](op(s))*m fi
       od;
       c
    end);
    P[Poly2List] := subs(env, proc(a) local i,n,c,t;
       c := coeffs(a, D[Variables], t); c := [c]; t := [t]; n := nops(c);
       [seq( [c[i],t[i]], i=1..n )]
    end);
    P[TotalDegree] := subs(env, proc(a) degree(a, {op(D[Variables])} ) end );
    P[Degree] := subs(env, proc(a,x) 
       if not member(x, D[Variables]) then ERROR(`bad variable`,x) fi;
       degree(a,x)
    end);
    P[Coeff] := subs(env, proc(a,e) local i, c, t; 
       c := coeffs(a, D[Variables], t);
       if member(e, [t], i) then c[i] else C[0] fi
    end );
    P[Lcoeff] := subs(env, proc(a, t) local c, w;
       c := lcoeff(a, D[Variables], w);
       if nargs=2 then t := w fi; c
    end);
    P[Lterm] := subs(env, proc(a) local c, w;
       c := lcoeff(a, D[Variables], w); c*w
    end);

    P[0] := R[0];
    P[1] := R[1];
    P[Coerce] := x -> x;

########
# Z[x] #
########
    if isDomain(R,Integer)  then
        P[List2Poly] := subs(env, proc(a)
           collect(convert(map(convert, a, `*`),`+`), D[Variables],
		'distributed');
        end);
        P[`+`] := proc() convert([args],`+`) end;
        P[`-`] := proc(a,b) if nargs = 1 then -a else a-b fi end;
        P[`*`] := proc() expand(`*`(args)) end;
        P[`.`] := proc(a,b) a*b end;
        P[Resultant] := subs(env, proc(a,b,x) resultant(a,b,x) end);
        P[Output] := x -> x;
        P[Div] := proc(a,b) local q;
           if irem( icontent(a), icontent(b) ) <> 0 then RETURN( FAIL ) fi;
           if divide(a,b,q) then q else FAIL fi
        end;
        P[Content] := proc(a,p) local c;
           c := icontent(a);
           if nargs = 2 then if a = 0 then p := 0 else p := a/c fi fi;
           c
        end;
        P[Primpart] := proc(a, cu) local c,p;
           if a = 0 then c := 0; p := 0
           else c := icontent(a)*sign(a); p := a/c
           fi;
           if nargs = 2 then cu := c fi;
           p
        end;
        P[Input] := subs(env, proc(a)
           if type(a,polynom(integer,D[Variables])) then expand(a) else FAIL fi
        end);
        P[Type] := subs(env, proc(a) type(a,polynom(integer,D[Variables])) end);
        P[Gcd] := proc() MMPGcd(args) end;
# 	P[PseudoRem] := subs(env, proc(a,b)
# 	   MMPPrem(a,b,D[Variables],args[3..nargs]) 
# 	end);
    fi;

######################
# Q[x], Q is a field #
######################
    if isDomain(R,Rational) then
        P[List2Poly] := subs(env, proc(a)
           collect(convert(map(convert, a, `*`),`+`), D[Variables],
		'distributed');
        end);
        P[`+`] := proc() convert([args],`+`) end;
        P[`-`] := proc(a,b) if nargs = 1 then -a else a-b fi end;
        P[`*`] := proc() expand(`*`(args)) end;
        P[`.`] := proc(a,b) a*b end;
        P[Resultant] := subs(env, proc(a,b,x) resultant(a,b,x) end);
        P[Output] := x -> x;
        P[Div] := proc(a,b) local q; if divide(a,b,q) then q else FAIL fi end; 
        P[Content] := proc(a,p)
           if nargs=2 then p := a fi;
           if a=0 then 0 else 1 fi
        end;
        P[Primpart] := subs(env, proc(a,r) local u;
           if a=0 then
              if nargs=2 then r := 0 fi; a
           else
              u := D[UnitNormal](a);
              if nargs=2 then r := u[1] fi; u[2]
           fi
        end);
        P[Input] := subs(env, proc(a)
           if type(a,polynom(rational,D[Variables])) then expand(a) else FAIL fi
        end);
    P[Type] := subs(env, proc(a) type(a,polynom(rational,D[Variables])) end);
        P[Gcd] := proc() local g;
           g := MMPGcd(args); if g = 0 then 0 else g/lcoeff(g) fi;
        end;
# 	    P[PseudoRem] := subs(env, proc(a,b)
# 	       MMPPrem(a,b,D[Variables],args[3..nargs]) 
# 	    end);
# 	P[Prime] := proc(a) end: # {true,false}
# 	P[Factor] := proc(a) end: # [ u, [[p,p], [u,u] , ... , [u,u]] ]
# 	P[Sqrfree] := proc(a) end: # [ u, [[p,p], [u,u] , ... , [u,u]] ]
    fi;

#########
# Zn[x] #
#########
    if isDomain(R,Zmod) then n := R[Size]:
        P[`+`] := subs('N'=n, proc() modp(convert([args], `+`), N) end);
        P[`-`] := subs('N'=n, proc(a,b) 
           if nargs = 1 then modp(-a,N) else modp(a-b,N) fi 
        end);
        P[`*`] := subs('N'=n, proc() modp( expand(`*`(args)) ,N) end);
        P[`.`] := subs('N'=n, proc(a,b) modp( a*b , N) end);
        P[Resultant] := subs(N=n, proc(a,b,x) modp(Resultant(a,b,x), N) end);
        P[Output] := x -> x;
        if hasCategory(R, Field) then # n is prime
           P[Div] := subs('N'=n, proc(a,b) local q; 
              if modp(Divide(a,b,q), N) then q else FAIL fi 
           end);
           P[Content] := proc(a,p) 
              if nargs=2 then p := a fi; if a=0 then 0 else 1 fi 
           end;
           P[Primpart] := subs(env, proc(a,r) local u;
              if a=0 then if nargs=2 then r := 0 fi; a
              else u := D[UnitNormal](a); if nargs=2 then r := u[1] fi; u[2]
              fi
           end);
           P[Gcd] := subs('N'=n, proc() local a, b, s;
              s := {args} minus {0};
              if s = {} then RETURN(0) else a := s[1]; s := s minus {a} fi;
              for b in s while a <> 1 do a := modp(Gcd(a,b), N) od;
              a
           end);
        fi;
    fi;
        
    op(P)
end:

`MMP/Gcd` := proc() local a,b,s;
    s := {args} minus {0};
    if s = {} then RETURN(0) else a := s[1]; s := s minus {a} fi;
    divide(a,icontent(a),'a');
    a := sign(a)*a;
    for b in s while a <> 1 do a := gcd(a,b) od; a
end:

# # Compute the Pseudo-Remainder of a divided b in Q[x] or Z[x]
# `MMP/Prem` := proc(a,b,x,m,q) local db, lb, dd, d, n, r, dr, lr;
#     if b = 0 then ERROR(`division by zero`) fi;
#     db := degree(b,x);
#     lb := coeff(b,x,db);
#     dr := degree(a,x);
#     dd := max(0,dr-db+1);
#     d := b - lb*x^db;
#     r := a;
#     for n from 0 while dr >= db and r <> 0 do
#     lr := coeff(r,x,dr);
#     r := lb * (r - lr*x^dr) - expand(lr*x^(dr-db)*d);
#     # r := lb * r - expand( coeff(r,x,dr) * x^(dr-db) * b );
#     dr := degree(r,x)
#     od;
#     r := lb^(dd-n)*r;
#     if nargs > 3 then m := lb^dd fi;
#     if nargs > 4 then divide(eval(m)*a-r,b,q) fi;
#     r
# end:

save `MMP.m`;
quit

