Naive approach

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)])