The following code shows how to tackle the problem, by manipulating individual terms.
While this is a very direct approach, it is quite slow.
The method reducibleBy gives a test for divisibility.
def subst(f,x,c):
if c==1:
return sum([t/x for t in f.terms() if t.reducibleBy(x)])+\
sum([t for t in f.terms() if not t.reducibleBy(x)])
else:
#c==0
return sum([t for t in f.terms() if not t.reducibleBy(x)])