`UP/Subresultant` := proc(D,p,q) local c,d,u,v,r,g,h,du,dv,t,C;
	C := D[CoefficientRing];
	if p = D[0] or q = D[0] then RETURN(C[0]) fi; 
	du := D[Degree](p); dv := D[Degree](q);
	if du = 0 then RETURN(C[`^`](D[Coeff](p,0),dv)) fi;
	if du < dv then 
		r := `UP/Subresultant`(D,q,p);
		if irem(du*dv,2) = 1 then r := C[`-`](r) fi;
		RETURN( r )
	fi;
	# Assumes degree(p) >= degree(q);
	u := p; v := q; c := C[1]; g := C[1]; h := C[1];
	while dv > 0 do
	    d := du-dv;
	    c := C[`*`]((-1)^(du*dv),c);
	    r := D[PseudoRem](u,v);
	    u := v; v := r; du := dv; dv := D[Degree](r);
	    if d>0 then
		t := C[`*`](g,C[`^`](h,d));
		v := D[Div](v,D[Constant](t));
	    fi;
	    g := D[Coeff](u,du);
	    # if d=1 then h := g
	    if d>0 then h := C[Div](C[`^`](g,d),C[`^`](h,d-1)) fi
	od;
	C[`*`](c, C[Div](C[`^`](D[Coeff](v,0),du),C[`^`](h,du-1)))
end:
save `Subres.m`;
quit
