find quotient of two multivariate polynomials (which are divisible)

edit

There really should be a version of reduce somewhere that gives you which combination of the generators allowed for the reduction, and if someone can find where then that really should be the preferred solution, but at the expense of some extra variables, you can fake it:

Given two polynomials f and g in variables [z1,...,zn], embed the ring into a larger one with variables [z1,...,zn,F,G], with an elimination order for z1,...zn (i.e., those variables should be larger than F and G). Then reduce f-F with respect to the principal ideal generated by g-G. The result should be a polynomial that is free of z1,...,zn, and expresses F as a polynomial in G.

Your data is fishy:

sage: denom.parent().term_order()
Negative weighted degree reverse lexicographic term order with weights (1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 4)

This 'local' term ordering causes reduce and (hence) divides to behave differently. In particular, it does not give the same result as with a well-ordering such as 'lex' in your ring dR:

sage: dR(denom).divides(dR(num))
False

So denom does not divide num in the ordinary sense.

Edit: In the 'local' sense, we have that denom is invertible, so the ideal generated by denom is the unit ideal, so of course num.reduce(Ideal([denom])) is zero.