
#
# MapleUnivariatePolynomial(R,x);
#
# Creates a Maple univariate polynomials domain for R[x] where R is
# the integers or rationals or integers modulo n.  I.e. R is one of the Gauss
# domains Integer, Rational, Zmodn.  The Maple sum of products representation
# is used for efficiency.
#
# Author MBM: 1991, 1992
#

macro( MUPGcd = `MUP/Gcd` );
macro( MUPRem = `MUP/Rem` );
macro( MUPPrem = `MUP/Prem` );

MapleUnivariatePolynomial := proc() local x,P,R,env,n; option remember;

    R := args[1];
    if nargs = 2 then x := args[2] fi;
    if not type(x,name) then ERROR(`2nd argument must be a name`) fi;

    P := UnivariatePolynomial(R);
    if nargs = 2 then P[Variable] := x fi;

    P[DomainName] := MapleUnivariatePolynomial;
    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 a = D[Variable] then true
    	elif type(a,anything^integer) then D[Type](op(1,a))
    	else C[Type](a)
    	fi
        end);

    if hasProperties(R,{UniquelyRepresented,CannonicalForm}) then
    	P[`=`] := proc(x,y) evalb(x=y) end;
    	P[`<>`] := proc(x,y) evalb(x<>y) end;
    	addProperties(P,{CannonicalForm,UniquelyRepresented})
    fi;

    P[0] := R[0];
    P[1] := R[1];
    P[Degree] := subs(env, proc(a) degree(a,D[Variable]) end );
    P[Ldegree] := subs(env, proc(a) ldegree(a,D[Variable]) end );
    P[Lcoeff] := subs(env, proc(a) lcoeff(a,D[Variable]) end );
    P[Tcoeff] := subs(env, proc(a) tcoeff(a,D[Variable]) end );
    P[Coeff] := subs(env, proc(a,i) coeff(a,D[Variable],i) end );
    P[Polynom] := subs(env, proc(a) local i,v;
    	v := D[Variable];
    	convert( [seq(a[i]*v^(i-1), i=1..nops(a))], `+` )
        end);
    P[Coerce] := proc(x) x end;
    P[Constant] := proc(x) x end;
    P[SetCoeffs] := subs(env, proc(a) {coeffs(a,D[Variable])} end);
    P[ListCoeffs] := subs(env, proc(a) local i,v;
    	v := D[Variable];
    	[seq( coeff(a,v,i), i=0..degree(a,v) )]
        end);

    if isDomain(R,Integer) or isDomain(R,Rational) then # Q[x] or Z[x]
        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[Diff] := subs(env, proc(a) diff(a,D[Variable]) end);
        P[PseudoRem] := subs(env, proc(a,b)
    	MUPPrem(a,b,D[Variable],args[3..nargs]) end);
        P[Resultant] := subs(env, proc(a,b) resultant(a,b,D[Variable]) end);
        P[Output] := proc(a) a end;
    fi;
    if isDomain(R, Zmod) then # Zmod(p)[x]
	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[Diff] := subs(env, 'N'=n, proc(a) diff(a,D[Variable]) end);
    	P[Div] := subs('N'=n, proc(a,b) local q;
		if modp(Divide(a,b,q), N) then q else FAIL fi
	    end);
	P[Gcd] := subs(env, 'N'=n, proc(a,b)
		if nargs = 0 then 0
		elif nargs = 1 then D[Normal](a)
		elif nargs > 2 then D[Gcd](D[Gcd](a,b),args[3..nargs])
		else modp(Gcd(a,b),N)
		fi
	    end);
	P[Resultant] := subs(env, N=n, proc(a,b)
		modp(Resultant(a,b,D[Variable]), N) end);
	P[Output] := x -> x;
    fi;
    if isDomain(R,Rational) then # Q[x]
        P[Div] := proc(a,b) local q; if divide(a,b,q) then q else FAIL fi end;
        P[Rem] := subs(env, proc(a,b)
    	MUPRem(a,b,D[Variable],args[3..nargs]) end);
        P[Gcd] := proc() local g;
    	    g := MUPGcd(args); if g = 0 then 0 else g/lcoeff(g) fi;
    	end;
        P[Input] := subs(env, proc(a)
    	    if type(a,polynom(rational,D[Variable])) then expand(a)
    		else FAIL fi
    	end);
        P[Type] := subs(env, proc(a) type(a,polynom(rational,D[Variable])) end);
    fi;
    if isDomain(R,Integer) then # Z[x]
        P[Gcd] := proc() MUPGcd(args) end;
        P[Primpart] := proc(a) local c,p;
    	    if a = 0 then c := 0; p := 0
    	    else c := icontent(a)*sign(a); divide(a,c,p);
    	    fi;
    	    if nargs = 2 then assign(args[2],c) fi;
    	    p
    	end;
        P[Content] := proc(a,p) local c;
    	    c := icontent(a);
    	    if nargs = 2 then if a = 0 then p := 0 else divide(a,c,p) fi fi;
    	    c
    	end;
        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[Input] := subs(env, proc(a)
    	    if type(a,polynom(integer,D[Variable])) then
    		expand(a) else FAIL fi
    	end);
        P[Type] := subs(env, proc(a) type(a,polynom(integer,D[Variable])) end);
    fi;
    op(P)
end:

`MUP/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 quotient/remainder of a divided b in Q[x] as fast as possible
`MUP/Rem` := proc(a,b,x) local db, lb, r, dr, q, t;
    if b = 0 then ERROR(`division by zero`) fi;
    db := degree(b,x);
    dr := degree(a,x);
    lb := coeff(b,x,db);
    r := a;
    q := 0;
    while dr >= db and r <> 0 do
        t := coeff(r,x,dr)/lb * x^(dr-db);
    q := t + q;
    r := r - expand( t*b );
    dr := degree(r,x)
    od;
    if nargs = 4 then assign(args[4],q) fi;
    r
end:

# Compute the Pseudo-Remainder of a divided b in Q[x] or Z[x]
`MUP/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 `MUP.m`;
quit
