#
# UnivariatePolynomial(R)
#
# The univariate polynoimal category over the ring R
#
# Author MBM: 1988 - 1992
#

macro(UFD=UniqueFactorizationDomain);

macro(UPInput=`UP/Input`, UPOutput=`UP/Output`, UPEqual=`UP/Equal`);
macro(UPAdd=`UP/Add`, UPMul=`UP/Mul`, UPInv=`UP/Inv`);
macro(UPShift=`UP/Shift`, UPEval=`UP/Eval`, UPPolyEval=`UP/PolyEval`);
macro(UPLess=`UP/Less`);
macro(UPNormal=`UP/Normal`, UPDiv=`UP/Div`);

macro(UPRem=readlib(`UP/Rem`,``.gauss.`/UP/Rem.m`));
macro(UPEuclideanGcd=readlib(`UP/EucGcd`,``.gauss.`/UP/EucGcd.m`));
macro(UPEuclideanRes=readlib(`UP/EucRes`,``.gauss.`/UP/EucRes.m`));
macro(UPPseudoRem=readlib(`UP/PRem`,``.gauss.`/UP/PRem.m`));
macro(UPReducedPRS=readlib(`UP/ReducedPRS`,``.gauss.`/UP/RedPRS.m`));
macro(UPSubresultant=readlib(`UP/Subresultant`,``.gauss.`/UP/Subres.m`));
macro(UPMGCD=readlib(`UP/ModGcd`,``.gauss.`/UP/ModGcd.m`));

UnivariatePolynomial := proc() local P,R,env,U,`?`;
	
	R := args[1];
	if nargs = 1 then P := newCategory(); else P := args[2] fi;

	addCategory(P,UnivariatePolynomial);
	if hasCategory(R,Field) then P := EuclideanDomain(P)
	elif hasCategory(R,UFD) then P := UFD(P)
	elif hasCategory(R,GcdDomain) then P := GcdDomain(P)
	elif hasCategory(R,IntegralDomain) then P := IntegralDomain(P)
	elif hasCategory(R,CommutativeRing) then P := CommutativeRing(P)
	elif hasCategory(R,Ring) then P := Ring(P)
	else ERROR(`1st argument must be a Ring`)
	fi;

	if hasCategory(R,OrderedSet) then
	    P := OrderedSet(P);
	fi;

	defOperation( Variable, Name, P );
	defOperation( Polynom, List(R) &-> P, P );
	defOperation( SetCoeffs, P &-> FiniteSet(R), P );
	defOperation( ListCoeffs, P &-> List(R), P );
	defOperation( ArrayCoeffs, P &-> Array(R), P );
	defOperations( {Degree,Ldegree}, P &-> Integer, P );
	defOperations( {Lcoeff,Tcoeff}, P &-> R, P );
	defOperation( Coeff, [P,Integer] &-> R, P );
	defOperation( Monomial, {[] &-> P, Integer &-> P}, P );
	defOperation( Constant, R &-> P, P );
	defOperation( Reductum, P &-> P, P );
	defOperation( Map, [R &-> R,P] &-> P, P );
	defOperation( Random, { [] &-> P, Integer &-> P,
				[Integer,[] &-> R] &-> P}, P );
	defOperation( Diff, P &-> P, P );
	defOperation( `.`, [R,P] &-> P, P ); # scalar multiplication
	defOperation( Eval, [P,R] &-> R, P ); # polynomial evaluation
	defOperation( PolyEval, [P,P] &-> P, P ); # polynomial evaluation
	defOperation( Shift, [P,Integer] &-> P, P );
	defOperation( CoefficientRing, Ring, P );

	env := ['D' = P, 'C' = R];
	P[CoefficientRing] := R;
	P[Characteristic] := R[Characteristic];
	P[Variable] := `?`;
	if hasCategory(R,OrderedSet) then
	    P[`<`] := subs(env, proc(a,b) UPLess(D,a,b) end);
	fi;


	P[SetCoeffs] := subs(env,
		proc(a) {op(D[ListCoeffs](a))}
	    end);
	P[ListCoeffs] := subs(env, proc(a) local i;
		[seq(D[Coeff](a,i), i = 0..D[Degree](a))]
	    end);
	P[ArrayCoeffs] := subs(env, proc(x) local a,l,n,i;
		l := D[ListCoeffs](x);
		n := nops(l);
		a := array(0..n-1);
		for i to n do a[i-1] := l[i] od;
		op(a);
	    end);

	P[EuclideanNorm] := subs(env, proc(a) D[Degree](a) end);
	P[Lcoeff] := subs(env, proc(a) D[Coeff](a,D[Degree](a)) end);
	P[Tcoeff] := subs(env, proc(a) D[Coeff](a,D[Ldegree](a)) end);
	P[Monomial] := subs(env, proc(n)
		if nargs=0 then RETURN( D[Monomial](1) ) fi;
		D[Polynom]([C[0]$n,C[1]])
	    end);
	P[Constant] := subs(env, proc(c) D[Polynom]([c]) end);
	P[Coerce] := subs(env, proc(n) D[Polynom]([C[Coerce](n)]) end);
	P[Reductum] := subs(env, proc(a) local c;
		c := D[ListCoeffs](a);
		D[Polynom]([c[1..nops(c)-1]])
	    end);
	P[Map] := subs(env, proc(f,a) local c,k;
		c := D[ListCoeffs](a);
		c := [seq(f(c[k],args[3..nargs]), k=1..nops(c))];
		D[Polynom](c)
	    end);
	P[RandomDegree] := rand(2..4);
	P[Random] := subs(env, proc() local n,f,l;
		if nargs = 0 then n := D[RandomDegree]() else n := args[1] fi;
		f := C[Random];
		l := f(); while l = C[0] do l := f() od;
		D[Polynom](['f()'$n,l])
	    end);
	P[Diff] := subs(env, proc(a) local d,k;
		d := D[Degree](a); if d = 0 then RETURN(D[0]) fi;
		D[Polynom]([seq(C[`*`](k,D[Coeff](a,k)), k=1..d)])
	    end);
		
	P[Input] := subs(env, proc(a) UPInput(D,a) end);
	P[Output] := subs(env, proc(a) UPOutput(D,a) end);
	P[`=`] := subs(env, proc(a,b) UPEqual(D,a,b) end);
	P[`+`] := subs(env, proc() UPAdd(D,args) end);
	P[`*`] := subs(env, proc() UPMul(D,args) end);
	P[Inv] := subs(env, proc(a) UPInv(D,a) end);
	P[`.`] := subs(env, proc(r,a) local c,k;
		if r = C[1] then RETURN(a) else c := D[ListCoeffs](a) fi;
		D[Polynom]( [seq(C[`*`](r,c[k]), k=1..nops(c))] )
	    end);

	P[Eval] := subs(env, proc(a,b) UPEval(D,a,b) end);
	P[PolyEval] := subs(env, proc(a,b) UPPolyEval(D,a,b) end);
	P[Shift] := subs(env, proc(a,n) UPShift(D,a,n) end);

	if hasCategory(R,IntegralDomain) then
	        defOperation( {PseudoRem,PseudoQuo}, [P,P,Name,Name] &-> P, P );
		defOperation( Resultant, [P,P] &-> R, P );
		if hasOperation(R,ModularHomomorphism) then
			defOperation( ModularHomomorphism,
				[] &-> [P &-> Union(P,FAIL),Integer], P );
			P[ModularHomomorphism] := subs(env,
			    proc() option remember;
				PolynomialHomomorphism(C,D)
			    end);
		fi;
		P[Unit] := subs(env, proc(a)
			D[Polynom]( [C[Unit](D[Lcoeff](a))] )
		    end);
		P[Normal] := subs(env, proc(a) UPNormal(D,a) end);
		P[Div] := subs(env, proc(a,b) UPDiv(D,a,b) end);
		P[PseudoRem] := subs(env, proc() UPPseudoRem(D,args) end);
		P[PseudoQuo] := subs(env, proc(a,b) local r,m,q;
			r := D[PseudoRem](a,b,m,q);
			if nargs > 2 then assign(args[3],m) fi;
			if nargs > 3 then assign(args[3],r) fi;
			q
		    end);
		P[Resultant] := subs(env, proc() UPResultant(D,args) end);
	fi;

	if hasCategory(R,GcdDomain) then
		defOperation( Content, [P,Name] &-> R, P );
		defOperation( Primpart, [P,Name] &-> P, P );
		P[Content] := subs(env, proc(a,p) local g;
			if a = D[0] then
			    if nargs = 2 then p := a fi; RETURN( C[0] )
			fi;
			g := C[Gcd]( op( D[SetCoeffs](a) ) );
			if nargs = 2 then p := D[Div](a,D[Constant](g)) fi;
			g
		    end);
		P[Primpart] := subs(env, proc(a,c) local p,t,u;
			t := D[Content](a,p);
			u := C[Unit](D[Lcoeff](p));
			p := D[`.`](C[Inv](u),p);
			if nargs = 2 then c := C[`*`](u,t) fi;
			p
		    end);
		P[Gcd] := subs(env, proc() UPReducedPRS(D,args) end)
	fi;

	if isDomain(R,Integer) then
		P[Gcd] := subs(env, proc() UPMGCD(D,args) end);
	fi;

	# if hasCategory(R,UFD) then
	#	P[Sqrfree] := ...
	# fi;

	if hasCategory(R,Field) then
		defOperation( Interpolate, [List(Integer),List(R)] &-> P, P );
		defOperation( Integrate, P &-> P, P );
		P[Rem] := subs(env, proc() UPRem(D,args) end);
		P[Gcd] := subs(env, proc() UPEuclideanGcd(D,args) end);
		P[Resultant] := subs(env, proc() UPEuclideanRes(D,args) end);
	fi;

	op(P)
end:


PolynomialHomomorphism := proc(C,D) local f,h,p;
	f := C[ModularHomomorphism](); p := f[2]; f := f[1];
	h := subs( {'I'=D,'g'=f,'m'=p}, proc(x) local a;
	        a := map(g,I[ListCoeffs](x));
	        if member(FAIL,a) then RETURN(FAIL) fi;
		modp1( Eval(ConvertIn(a),1231), m )
	    end);
	op(h),p
end:

RelativelyPrimeHeuristic := proc(C,D) local a,b,f,p;
	if not hasOperation(C,ModularHomomorphism) then RETURN(FAIL) fi;
	f := C[ModularHomomorphism](); p := f[2]; f := f[1];
	a := map(f,D[ListCoeffs](args[3]));
	b := map(f,D[ListCoeffs](args[4]));
	if member(FAIL,a) or member(FAIL,b) then RETURN(FAIL) fi;
	evalb( modp1( Gcd(ConvertIn(a),ConvertIn(b)), p ) = 1 )
end:

`UP/Input` := proc(D,a) local d,k,r,v,C;
	C := D[CoefficientRing];
	v := D[Variable];
	d := degree(a,v);
	if d = FAIL then RETURN(FAIL) fi;
	if ldegree(a,v) < 0 then RETURN(FAIL) fi;
	r := frontend(taylor,[a,v,d+1]);
	if r = lasterror then RETURN(FAIL) fi;
	r := [seq(coeff(r,v,k),k=0..d)];
	r := map(C[Input],r);
	if has(r,FAIL) then FAIL else D[Polynom](r) fi
end:

`UP/Output` := proc(D,a) local c,d,i,k,v,C;
	C := D[CoefficientRing];
	d := D[Degree](a);
	v := D[Variable];
	for k from 0 to d do
	    c[k] := C[Output](D[Coeff](a,k));
	    c[k] := c[k] * v^k
	od;
	convert([seq(c[d-i], i=0..d)],`+`);
end:


`UP/Shift` := proc(D,a,n) local C,t;
    C := D[CoefficientRing];
    if n < 0 then
        t := D[ListCoeffs](a);
        t := [t[1-n..nops(t)]];
        D[Polynom](t)
    elif n > 0 then
        t := C[0] $ n;
        D[Polynom]([t,op(D[ListCoeffs](a))])
    else a
    fi
end:

`UP/Equal` := proc(D,a,b) local C,d,i;
	C := D[CoefficientRing];
	d := D[Degree](a);
	if d <> D[Degree](b) then RETURN(false) fi;
	for i from 0 to d do
	    if C[`<>`](D[Coeff](a,i),D[Coeff](b,i)) then
		RETURN(false) fi
	od;
	true
end:

# polynomial evaluation : ($,$) -> $
`UP/PolyEval` := proc(D,a,b) local d,k,r;
	d := D[Degree](a);
	r := D[Constant](D[Coeff](a,d));
	for k from d-1 by -1 to 0 do
	    r := D[`+`](D[`*`](r,b),D[Constant](D[Coeff](a,k)))
	od;
	r
end:

# polynomial evaluation : ($,R) -> R
`UP/Eval` := proc(D,a,x) local d,k,r,C;
	C := D[CoefficientRing];
	d := D[Degree](a);
	r := D[Coeff](a,d);
	for k from d-1 by -1 to 0 do
	    r := C[`+`](C[`*`](r,x),D[Coeff](a,k))
	od;
	r
end:

# Add: UP* -> UP
`UP/Add` := proc(D) local b,c,i,j,m,n,C;
    C := D[CoefficientRing];
    m := 0;
    for b in [args[2..nargs]] do
	n := D[Degree](b);
	if n > m then m := n fi
    od;
    c := array(0..m);
    for i from 0 to m do c[i] := C[0] od;
    for b in [args[2..nargs]] do
	b := D[ListCoeffs](b);
	n := nops(b);
	for i from 0 to min(m,n-1) do c[i] := C[`+`](b[i+1],c[i]) od;
    od;
    D[Polynom]([seq(c[j], j=0..m)])
end:

`UP/Mul` := proc(D,a,b) local ca,cb,c,d,i,j,k,na,nb,s,C;
	C := D[CoefficientRing];
	if nargs = 1 then RETURN(D[1])
	elif nargs = 2 then RETURN(a)
	elif nargs > 3 then RETURN(D[`*`](D[`*`](a,b),args[4..nargs]))
	elif type(a,integer) then RETURN( D[`.`](C[Coerce](a),b) )
	elif a = D[0] or b = D[0] then RETURN( D[0] )
	elif a = D[1] then RETURN(b)
	elif b = D[1] then RETURN(a)
	fi;

	ca := D[ListCoeffs](a); na := nops(ca);
	if na = 1 then RETURN( D[`.`](ca[1],b) ) fi;
	if a = b then
	    d := 2*na; c := array(2..d); s := C[0];
	    for i from 2 to d do c[i] := s od;
	    for i to na do 
	        s := ca[i]; if s = C[0] then next fi;
	        for j to i-1 do
		    c[i+j] := C[`+`](c[i+j],C[`*`](s,ca[j]));
	        od
	    od;
	    for i from 3 to d-1 do c[i] := C[`*`](2,c[i]) od;
	    for i to d/2 do
		c[2*i] := C[`+`](c[2*i],C[`*`](ca[i],ca[i]))
	    od
	else
	    cb := D[ListCoeffs](b); nb := nops(cb);
	    d := na+nb; c := array(2..d); s := C[0];
	    for i from 2 to d do c[i] := s od;
	    for i to na do 
	        s := ca[i]; if s = C[0] then next fi;
	        for j to nb do
		    c[i+j] := C[`+`](c[i+j],C[`*`](s,cb[j]));
		od
	    od
	fi;
	D[Polynom]( [seq(c[k],k=2..d)] )
end:

`UP/Normal` := proc(D,a) local c,u,C;
	C := D[CoefficientRing];
	if a = D[0] then RETURN(a) fi;
	c := D[ListCoeffs](a);
	u := C[Unit](c[nops(c)]);
	c := map(C[`*`],c,C[Inv](u));
	D[Polynom](c)
end:

`UP/Inv` := proc(D,a) local b,C;
	C := D[CoefficientRing];
	# if type(a,integer) then RETURN( D[Inv](D[Coerce](a)) ) fi;
	if a = D[1] then RETURN(a) fi;
	if D[Degree](a) <> 0 then RETURN(FAIL) fi;
	b := C[Inv](D[Tcoeff](a));
	if b = FAIL then b else D[Polynom]([b]) fi
end:

# `<`: ($,$) -> Boolean
`UP/Less` := proc(D,a,b) local m,n,s,t,C;
	C := D[CoefficientRing];
	m := D[Degree](a);
	n := D[Degree](b);
	if m < n then RETURN(true) elif m > n then RETURN(false) fi;
	for i from 0 to min(m,n) do
	    s := D[Coeff](a,i); t := D[Coeff](b,i);
	    if C[`<>`](s,t) then RETURN( C[`<`](s,t) ) fi;
	od;
	false
end:

# Div : (UP,UP) -> Union(UP,FAIL)
`UP/Div` := proc(D,x,y) local a,b,da,db,dq,k,q,C;
	C := D[CoefficientRing];
	if y = D[0] then ERROR(`division by zero`) fi;
	if x = D[0] then RETURN(x) fi;
	if y = D[1] then RETURN(x) fi;
	if x = y then RETURN(D[1]) fi;
	da := D[Degree](x);
	db := D[Degree](y);
	dq := da-db; if dq < 0 then RETURN(FAIL) fi;
	a := D[ArrayCoeffs](x);
	b := D[ArrayCoeffs](y);
	q := array(0..dq);
	`UP/fill`(q,dq,C[0]);
	if not `UP/DivInPlace`(C,a,b,da,db,q) then RETURN(FAIL) fi;
	D[Polynom]([seq(q[k], k=0..dq)])
end:

# Computes q := Div(a,b) inplace
`UP/DivInPlace` := proc(C,a,b,da,db,q) local dr,e,i,j,lb;
	dr := da;
	e := da - db;
	lb := b[db];
	while e >= 0 do
	    q[e] := C[Div](a[dr],lb);
	    if q[e] = FAIL then RETURN(false) fi;
	    i := e;
	    for j from 0 to db-1 do
		a[i] := C[`-`](a[i],C[`*`](q[e],b[j]));
		i := i+1;
	    od;
	    dr := dr-1;
	    while dr >= 0 and a[dr] = C[0] do dr := dr-1 od;
	    e := dr - db
	od;
	evalb(dr = -1)
end:

`UP/fill` := proc(a,n,x) local i; for i from 0 to n do a[i] := x od end:

# Computes a := Rem(a,b) inplace
`UP/InPlaceRem` := proc(C,a,b,da,db,l) local dr,e,i,j,t;
	dr := da;
	e := da - db;
	while e >= 0 do
	    t := C[`*`](a[dr],l);
	    i := e;
	    for j from 0 to db-1 do
		a[i] := C[`-`](a[i],C[`*`](t,b[j]));
		i := i+1
	    od;
	    dr := dr-1;
	    while dr >= 0 and a[dr] = C[0] do dr := dr-1 od;
	    e := dr - db
	od;
	dr
end:

save `UP.m`;
quit
