
`UP/PseudoRemInPlace` := proc(C,a,b,da,db,m) local d,dr,i,j,k,t,l;
	dr := da;
	l := b[db];
	d := max(da-db+1,0);
	for k from 0 while dr-db >= 0 do
	    for j from 0 to dr-1 do a[j] := C[`*`](l,a[j]) od;
	    i := dr-db;
	    for j from 0 to db-1 do
		a[i] := C[`-`](a[i],C[`*`](a[dr],b[j]));
		i := i+1
	    od;
	    dr := dr-1;
	    while dr >= 0 and a[dr] = C[0] do dr := dr-1 od;
	od;
	if d-k > 0 then
	    t := C[`^`](l,d-k);
	    for j from 0 to dr do a[j] := C[`*`](t,a[j]) od;
	fi;
	for k from 0 to dr do a[k] := C[Div](a[k],eval(m)) od;
	m := C[`^`](l,d);
	dr
end:

`UP/ReducedPRS` := proc(D) local a,da,b,db,c,k,pm,r,dr,s,C;

	s := {args[2..nargs]} minus {D[0]};
	if s = {} then RETURN( D[0] ) fi;

	C := D[CoefficientRing];

	a := s[1];
	s := s minus {a};
	c := D[Content](a,'a');
	for b in s while a <> D[1] or c <> C[1] do


	    c := C[Gcd](c,D[Content](b,'b'));

	    # Heuristic test for Gcd(a,b) = 1
	    # if RelativelyPrimeHeuristic(C,D,a,b) then
	    #	a := D[1]; next fi;

	    da := D[Degree](a); a := D[ArrayCoeffs](a);
	    db := D[Degree](b); b := D[ArrayCoeffs](b);

	    pm := C[1];
	    while db > 0 do
		dr := `UP/PseudoRemInPlace`(C,a,b,da,db,'pm');
		r := op(a); a := op(b); b := op(r);
		da := db; db := dr;
	    od;

	    if db = 0 then a := D[1]
	    else
	        a := D[Polynom]([seq(a[k], k=0..da)]);
	        a := D[Primpart](a)
	    fi
	od;

	D[Normal](D[`.`](c,a))

end:

save `RedPRS.m`;
quit
