diff --git a/sympy/polys/__init__.py b/sympy/polys/__init__.py
index 52977070a9..1367a39870 100644
--- a/sympy/polys/__init__.py
+++ b/sympy/polys/__init__.py
@@ -2,6 +2,7 @@
 
 from polytools import (
     Poly, poly,
+    poly_from_expr,
     parallel_poly_from_expr,
     degree, degree_list,
     LC, LM, LT,
diff --git a/sympy/polys/constructor.py b/sympy/polys/constructor.py
index de5788d995..f7ae63bc1d 100644
--- a/sympy/polys/constructor.py
+++ b/sympy/polys/constructor.py
@@ -180,7 +180,7 @@ def _construct_expression(coeffs, opt):
 
 def construct_domain(obj, **args):
     """Construct a minimal domain for the list of coefficients. """
-    opt = build_options(**args)
+    opt = build_options(args)
 
     if isinstance(obj, dict):
         monoms, coeffs = zip(*obj.items())
diff --git a/sympy/polys/domains/domain.py b/sympy/polys/domains/domain.py
index 5bcdb1dc82..8c338575e2 100644
--- a/sympy/polys/domains/domain.py
+++ b/sympy/polys/domains/domain.py
@@ -190,7 +190,8 @@ def unify(K0, K1, gens=None):
         if gens is not None:
             if (K0.is_Composite and (set(K0.gens) & set(gens))) or (K1.is_Composite and (set(K1.gens) & set(gens))):
                 raise UnificationFailed("can't unify %s with %s, given %s generators" % (K0, K1, tuple(gens)))
-        elif K0 == K1:
+
+        if K0 == K1:
             return K0
 
         if not K0.has_CharacteristicZero:
diff --git a/sympy/polys/numberfields.py b/sympy/polys/numberfields.py
index ba3f723ed5..0a18b9ebc5 100644
--- a/sympy/polys/numberfields.py
+++ b/sympy/polys/numberfields.py
@@ -74,7 +74,8 @@ def bottom_up_scan(ex):
 
                     alg = ex.base - coeff
 
-                    inverse = invert(elt.gen + coeff, elt)
+                    # XXX: turn this into eval()
+                    inverse = invert(elt.gen + coeff, elt).as_basic()
                     base = inverse.subs(elt.gen, alg).expand()
 
                     if ex.exp == -1:
diff --git a/sympy/polys/partfrac.py b/sympy/polys/partfrac.py
index 716ca11871..68ddc6c26e 100644
--- a/sympy/polys/partfrac.py
+++ b/sympy/polys/partfrac.py
@@ -1,11 +1,9 @@
 """Algorithms for partial fraction decomposition of rational functions. """
 
-from sympy.core import S, Symbol, Function, Lambda
-from sympy.polys import Poly, RootSum, cancel
+from sympy.core import S, sympify, Symbol, Function, Lambda
+from sympy.polys import Poly, RootSum, cancel, parallel_poly_from_expr
 from sympy.utilities import numbered_symbols, take, threaded
 
-from sympy.polys.polytools import _polify_basic
-
 @threaded
 def apart(f, x=None, full=False):
     """Compute partial fraction decomposition of a rational function.
@@ -25,16 +23,19 @@ def apart(f, x=None, full=False):
            y/(1 + x) - y/(2 + x)
 
     """
+    f = sympify(f)
+
     if f.is_Atom:
         return f
+    else:
+        P, Q = f.as_numer_denom()
 
     if x is None:
         gens = ()
     else:
         gens = (x,)
 
-    P, Q = f.as_numer_denom()
-    P, Q = _polify_basic(P, Q, *gens)
+    (P, Q), opt = parallel_poly_from_expr((P, Q), *gens)
 
     if P.is_multivariate:
         raise NotImplementedError("multivariate partial fraction decomposition")
diff --git a/sympy/polys/polyerrors.py b/sympy/polys/polyerrors.py
index 83e8b9f8c9..4dd805583d 100644
--- a/sympy/polys/polyerrors.py
+++ b/sympy/polys/polyerrors.py
@@ -54,6 +54,16 @@ class UnificationFailed(Exception):
 class GeneratorsNeeded(Exception):
     pass
 
+class ComputationFailed(Exception):
+
+    def __init__(self, func, nargs, exc):
+        self.func = func
+        self.nargs = nargs
+        self.exc = exc
+
+    def __str__(self):
+        return "%s(%s) failed without generators" % (self.func, ', '.join(map(str, self.exc.exprs[:self.nargs])))
+
 class GeneratorsError(Exception):
     pass
 
@@ -66,11 +76,25 @@ class MultivariatePolynomialError(PolynomialError):
 class OptionError(Exception):
     pass
 
+class FlagError(OptionError):
+    pass
+
 class PolificationFailed(PolynomialError):
 
-    def __init__(self, origs, exprs):
-        self.origs = tuple(origs)
-        self.exprs = tuple(exprs)
+    def __init__(self, origs, exprs, seq=False):
+        if not seq:
+            self.orig = origs
+            self.expr = exprs
+            self.origs = [origs]
+            self.exprs = [exprs]
+        else:
+            self.origs = origs
+            self.exprs = exprs
+
+        self.seq = seq
 
     def __str__(self): # pragma: no cover
-        return "can't construct polynomials from %s" % ', '.join(map(str, self.origs))
+        if not self.seq:
+            return "can't construct a polynomial from %s" % str(self.orig)
+        else:
+            return "can't construct polynomials from %s" % ', '.join(map(str, self.origs))
diff --git a/sympy/polys/polyoptions.py b/sympy/polys/polyoptions.py
index 3ce2eb66df..4e690212ce 100644
--- a/sympy/polys/polyoptions.py
+++ b/sympy/polys/polyoptions.py
@@ -6,6 +6,7 @@
     PolynomialError,
     GeneratorsError,
     OptionError,
+    FlagError,
 )
 
 from sympy.polys.monomialtools import monomial_key
@@ -21,6 +22,8 @@ class Option(object):
 
     option = None
 
+    is_Flag = False
+
     requires = []
     excludes = []
 
@@ -39,6 +42,11 @@ def postprocess(cls, options):
         """ """
         pass
 
+class Flag(Option):
+    """ """
+
+    is_Flag = True
+
 class BooleanOption(Option):
     """ """
 
@@ -64,9 +72,6 @@ def getter(self):
         setattr(Options, cls.option, getter)
         Options.__options__[cls.option] = cls
 
-class Flag(object):
-    """ """
-
 class Options(dict):
     """ """
 
@@ -153,6 +158,10 @@ class Expand(BooleanOption):
     requires = []
     excludes = []
 
+    @classmethod
+    def default(cls):
+        return True
+
 class Gens(Option):
     """ """
 
@@ -507,7 +516,18 @@ class Monic(BooleanOption, Flag):
     def default(cls):
         return True
 
-class Gen(Option, Flag):
+class All(BooleanOption, Flag):
+    """ """
+
+    __metaclass__ = OptionType
+
+    option = 'all'
+
+    @classmethod
+    def default(cls):
+        return False
+
+class Gen(Flag):
     """ """
 
     __metaclass__ = OptionType
@@ -517,16 +537,34 @@ class Gen(Option, Flag):
     requires = []
     excludes = []
 
+    @classmethod
+    def default(cls):
+        return 0
+
     @classmethod
     def preprocess(cls, gen):
-        if isinstance(gen, Basic):
+        if isinstance(gen, (Basic, int)):
             return gen
         else:
             raise OptionError("invalid argument for 'gen' option")
 
-def build_options(**args):
+def build_options(gens, args=None):
     """Construct options from keyword arguments or ... options. """
-    if len(args) != 1 or 'opt' not in args:
-        return Options((), args)
+    if args is None:
+        gens, args = (), gens
+
+    if len(args) != 1 or 'opt' not in args or gens:
+        return Options(gens, args)
     else:
         return args['opt']
+
+def allowed_flags(args, flags):
+    """ """
+    flags = set(flags)
+
+    for arg in args.iterkeys():
+        try:
+            if Options.__options__[arg].is_Flag and not arg in flags:
+                raise FlagError("'%s' flag is not allowed in this context" % arg)
+        except KeyError:
+            raise OptionError("'%s' is not a valid option" % arg)
diff --git a/sympy/polys/polytools.py b/sympy/polys/polytools.py
index ffcc69240f..eead382f2f 100644
--- a/sympy/polys/polytools.py
+++ b/sympy/polys/polytools.py
@@ -4,6 +4,10 @@
     S, Basic, I, Integer, Add, Mul, sympify, ask,
 )
 
+from sympy.core.sympify import (
+    SympifyError,
+)
+
 from sympy.core.decorators import (
     _sympifyit,
 )
@@ -39,7 +43,8 @@
     OperationNotSupported, DomainError,
     CoercionFailed, UnificationFailed,
     GeneratorsNeeded, PolynomialError,
-    PolificationFailed,
+    PolificationFailed, FlagError,
+    ComputationFailed,
 )
 
 from sympy.polys.polycontext import (
@@ -64,7 +69,7 @@
 from sympy.polys import polyoptions as options
 
 def init_poly_from_dict(dict_rep, **args):
-    opt = options.build_options(**args)
+    opt = options.build_options(args)
     return _init_poly_from_dict(dict_rep, opt)
 
 def _init_poly_from_dict(dict_rep, opt):
@@ -81,7 +86,7 @@ def _init_poly_from_dict(dict_rep, opt):
     return DMP(dict_rep, domain, len(gens)-1)
 
 def init_poly_from_list(list_rep, **args):
-    opt = options.build_options(**args)
+    opt = options.build_options(args)
     return _init_poly_from_list(list_rep, opt)
 
 def _init_poly_from_list(list_rep, opt):
@@ -100,7 +105,7 @@ def _init_poly_from_list(list_rep, opt):
     return DMP(list_rep, domain, len(gens)-1)
 
 def init_poly_from_poly(poly_rep, **args):
-    opt = options.build_options(**args)
+    opt = options.build_options(args)
     return _init_poly_from_poly(poly_rep, opt)
 
 def _init_poly_from_poly(poly_rep, opt):
@@ -134,7 +139,7 @@ def _init_poly_from_poly(poly_rep, opt):
     return (rep, gens or poly_rep.gens)
 
 def init_poly_from_basic(basic_rep, **args):
-    opt = options.build_options(**args)
+    opt = options.build_options(args)
     return _init_poly_from_basic(basic_rep, opt)
 
 def _init_poly_from_basic(basic_rep, opt):
@@ -244,6 +249,21 @@ def gen(self):
         """Return principal generator. """
         return self.gens[0]
 
+    @property
+    def domain(self):
+        """Get the ground domain of `self`. """
+        return self.get_domain()
+
+    @property
+    def zero(self):
+        """Return zero polynomial with `self`'s properties. """
+        return self.new(self.rep.zero(self.rep.lev, self.rep.dom), self.gens)
+
+    @property
+    def one(self):
+        """Return one polynomial with `self`'s properties. """
+        return self.new(self.rep.one(self.rep.lev, self.rep.dom), self.gens)
+
     def unify(f, g):
         """Make `f` and `g` belong to the same domain. """
         _, per, F, G = f._unify(g)
@@ -1489,26 +1509,110 @@ def __hash__(self):
     def __nonzero__(f):
         return not f.is_zero
 
+def poly_from_expr(expr, *gens, **args):
+    """Construct a polynomial from an expression. """
+    opt = options.build_options(gens, args)
+    orig, expr = expr, sympify(expr)
+
+    if not isinstance(expr, Basic):
+        raise PolificationFailed(orig, expr)
+    elif expr.is_Poly:
+        poly = Poly(expr, opt=opt)
+
+        opt['gens'] = poly.gens
+        opt['domain'] = poly.domain
+
+        if opt.polys is None:
+            opt['polys'] = True
+
+        if opt.frac:
+            return (poly, poly.one), opt
+        else:
+            return poly, opt
+    elif opt.frac:
+        numer, denom = expr.as_numer_denom()
+
+        if opt.expand:
+            numer = numer.expand()
+            denom = denom.expand()
+
+        try:
+            return parallel_poly_from_expr((numer, denom), opt=opt)
+        except PolificationFailed:
+            raise PolificationFailed(orig, numer/denom)
+    elif opt.expand:
+        expr = expr.expand()
+
+    try:
+        rep, gens = dict_from_basic(expr, opt=opt)
+    except GeneratorsNeeded:
+        raise PolificationFailed(orig, expr)
+
+    monoms, coeffs = zip(*rep.items())
+    domain = opt.domain
+
+    if domain is None:
+        domain, coeffs = construct_domain(coeffs, opt=opt)
+    else:
+        coeffs = map(domain.from_sympy, coeffs)
+
+    level = len(gens)-1
+
+    poly = Poly.new(DMP.from_monoms_coeffs(monoms, coeffs, level, domain), gens)
+
+    opt['gens'] = gens
+    opt['domain'] = domain
+
+    if opt.polys is None:
+        opt['polys'] = False
+
+    return poly, opt
+
 def parallel_poly_from_expr(exprs, *gens, **args):
     """Construct polynomials from expressions. """
-    opt = options.Options(gens, args)
+    opt = options.build_options(gens, args)
+
+    if len(exprs) == 2:
+        f, g = exprs
+
+        if isinstance(f, Poly) and isinstance(g, Poly):
+            f = Poly(f, opt=opt)
+            g = Poly(g, opt=opt)
+
+            f, g = f.unify(g)
+
+            opt['gens'] = f.gens
+            opt['domain'] = f.domain
+
+            if opt.polys is None:
+                opt['polys'] = True
+
+            return [f, g], opt
 
     origs, exprs = list(exprs), []
     _exprs, _polys = [], []
 
+    failed = False
+
     for i, expr in enumerate(origs):
         expr = sympify(expr)
 
-        if expr.is_Poly:
-            _polys.append(i)
-        else:
-            _exprs.append(i)
+        if isinstance(expr, Basic):
+            if expr.is_Poly:
+                _polys.append(i)
+            else:
+                _exprs.append(i)
 
-            if opt.expand:
-                expr = expr.expand()
+                if opt.expand:
+                    expr = expr.expand()
+        else:
+            failed = True
 
         exprs.append(expr)
 
+    if failed:
+        raise PolificationFailed(origs, exprs, True)
+
     if _polys:
         # XXX: this is a temporary solution
         for i in _polys:
@@ -1517,7 +1621,7 @@ def parallel_poly_from_expr(exprs, *gens, **args):
     try:
         reps, gens = parallel_dict_from_basic(exprs, opt=opt)
     except GeneratorsNeeded:
-        raise PolificationFailed(origs, exprs)
+        raise PolificationFailed(origs, exprs, True)
 
     coeffs_list, lengths = [], []
 
@@ -1557,37 +1661,6 @@ def parallel_poly_from_expr(exprs, *gens, **args):
 
     return polys, opt
 
-def NonStrictPoly(f, *gens, **args):
-    """Create a Poly instance with `strict` keyword set.  """
-    args = dict(args)
-    args['strict'] = False
-    return Poly(f, *gens, **args)
-
-def _polify_basic(f, g, *gens, **args):
-    """Cooperatively make polynomials out of `f` and `g`. """
-    if gens:
-        F = NonStrictPoly(f, *gens, **args)
-        G = NonStrictPoly(g, *gens, **args)
-
-        if not F.is_Poly or not G.is_Poly:
-            raise CoercionFailed(F, G) # pragma: no cover
-        else:
-            return F, G
-    else:
-        F = NonStrictPoly(f, **args)
-        G = NonStrictPoly(g, **args)
-
-        if F.is_Poly:
-            if G.is_Poly:
-                return F, G
-            else:
-                return F, Poly(g, *F.gens, **args)
-        else:
-            if G.is_Poly:
-                return Poly(f, *G.gens, **args), G
-            else:
-                raise CoercionFailed(F, G)
-
 def _update_args(args, key, value):
     """Add a new `(key, value)` pair to arguments dict. """
     args = dict(args)
@@ -1597,28 +1670,6 @@ def _update_args(args, key, value):
 
     return args
 
-def _filter_args(args, *keys):
-    """Filter the given keys from the args dict. """
-    if not keys:
-        return {}
-
-    keys, result = set(keys), {}
-
-    for key, value in args.items():
-        if key in keys:
-            result[key] = value
-
-    return result
-
-def _should_return_basic(*polys, **args):
-    """Figure out if results should be returned as basic. """
-    query = args.get('polys')
-
-    if query is not None:
-        return not query
-    else:
-        return not all(isinstance(poly, Poly) for poly in polys)
-
 def _keep_coeff(coeff, factors):
     """Return ``coeff*factors`` unevaluated if necessary. """
     if coeff == 1:
@@ -1631,623 +1682,681 @@ def _keep_coeff(coeff, factors):
         return Mul(coeff, factors, evaluate=False)
 
 def degree(f, *gens, **args):
-    """Returns degree of `f` in the given generator. """
-    try:
-        F = Poly(f, *_analyze_gens(gens), **args)
-    except GeneratorsNeeded:
-        raise GeneratorsNeeded("can't compute degree of %s without generators" % f)
+    """Return the degree of ``f`` in the given variable. """
+    options.allowed_flags(args, ['gen', 'polys'])
 
-    degree = F.degree(args.get('gen', 0))
+    try:
+        F, opt = poly_from_expr(f, *gens, **args)
+    except PolificationFailed, exc:
+        raise ComputationFailed('degree', 1, exc)
 
-    return Integer(degree)
+    return Integer(F.degree(opt.gen))
 
 def degree_list(f, *gens, **args):
-    """Returns a list of degrees of `f` in all generators. """
+    """Return a list of degrees of ``f`` in all variables. """
+    options.allowed_flags(args, ['polys'])
+
     try:
-        F = Poly(f, *_analyze_gens(gens), **args)
-    except GeneratorsNeeded:
-        raise GeneratorsNeeded("can't compute degrees list of %s without generators" % f)
+        F, opt = poly_from_expr(f, *gens, **args)
+    except PolificationFailed, exc:
+        raise ComputationFailed('degree_list', 1, exc)
 
     degrees = F.degree_list()
 
     return tuple(map(Integer, degrees))
 
 def LC(f, *gens, **args):
-    """Returns the leading coefficient of `f`. """
-    order = args.pop('order', None)
+    """Return the leading coefficient of ``f``. """
+    options.allowed_flags(args, ['polys'])
 
     try:
-        return Poly(f, *gens, **args).LC(order=order)
-    except GeneratorsNeeded:
-        raise GeneratorsNeeded("can't compute the leading coefficient of %s without generators" % f)
+        F, opt = poly_from_expr(f, *gens, **args)
+    except PolificationFailed, exc:
+        raise ComputationFailed('LC', 1, exc)
+
+    return F.LC(order=opt.order)
 
 def LM(f, *gens, **args):
-    """Returns the leading monomial of `f`. """
-    order = args.pop('order', None)
+    """Return the leading monomial of ``f``. """
+    options.allowed_flags(args, ['polys'])
 
     try:
-        f = Poly(f, *gens, **args)
-    except GeneratorsNeeded:
-        raise GeneratorsNeeded("can't compute the leading monomial of %s without generators" % f)
+        F, opt = poly_from_expr(f, *gens, **args)
+    except PolificationFailed, exc:
+        raise ComputationFailed('LM', 1, exc)
+
+    monom = Monomial(*F.LM(order=opt.order))
 
-    return Monomial(*f.LM(order=order)).as_basic(*f.gens)
+    return monom.as_basic(*opt.gens)
 
 def LT(f, *gens, **args):
-    """Returns the leading term of `f`. """
-    order = args.pop('order', None)
+    """Return the leading term of ``f``. """
+    options.allowed_flags(args, ['polys'])
 
     try:
-        f = Poly(f, *gens, **args)
-    except GeneratorsNeeded:
-        raise GeneratorsNeeded("can't compute the leading term of %s without generators" % f)
+        F, opt = poly_from_expr(f, *gens, **args)
+    except PolificationFailed, exc:
+        raise ComputationFailed('LT', 1, exc)
 
-    monom, coeff = f.LT(order=order)
+    monom, coeff = F.LT(order=opt.order)
 
-    return coeff*Monomial(*monom).as_basic(*f.gens)
+    return coeff*Monomial(*monom).as_basic(*opt.gens)
 
 def pdiv(f, g, *gens, **args):
-    """Polynomial pseudo-division of `f` and `g`. """
-    gens = _analyze_gens(gens)
+    """Compute polynomial pseudo--division of ``f`` and ``g``. """
+    options.allowed_flags(args, ['polys'])
 
     try:
-        F, G = _polify_basic(f, g, *gens, **args)
-    except CoercionFailed, (f, g):
-        raise GeneratorsNeeded("can't compute pseudo division of %s and %s without generators" % (f, g))
+        (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args)
+    except PolificationFailed, exc:
+        raise ComputationFailed('pdiv', 2, exc)
 
     q, r = F.pdiv(G)
 
-    if _should_return_basic(f, g, **args):
+    if not opt.polys:
         return q.as_basic(), r.as_basic()
     else:
         return q, r
 
 def prem(f, g, *gens, **args):
-    """Polynomial pseudo-remainder of `f` and `g`. """
-    gens = _analyze_gens(gens)
+    """Compute polynomial pseudo--remainder of ``f`` and ``g``. """
+    options.allowed_flags(args, ['polys'])
 
     try:
-        F, G = _polify_basic(f, g, *gens, **args)
-    except CoercionFailed, (f, g):
-        raise GeneratorsNeeded("can't compute pseudo remainder of %s and %s without generators" % (f, g))
+        (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args)
+    except PolificationFailed, exc:
+        raise ComputationFailed('prem', 2, exc)
 
     r = F.prem(G)
 
-    if _should_return_basic(f, g, **args):
+    if not opt.polys:
         return r.as_basic()
     else:
         return r
 
 def pquo(f, g, *gens, **args):
-    """Polynomial pseudo-quotient of `f` and `g`. """
-    gens = _analyze_gens(gens)
+    """Compute polynomial pseudo--quotient of ``f`` and ``g``. """
+    options.allowed_flags(args, ['polys'])
 
     try:
-        F, G = _polify_basic(f, g, *gens, **args)
-    except CoercionFailed, (f, g):
-        raise GeneratorsNeeded("can't compute pseudo quotient of %s and %s without generators" % (f, g))
+        (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args)
+    except PolificationFailed, exc:
+        raise ComputationFailed('pquo', 2, exc)
 
     q = F.pquo(G)
 
-    if _should_return_basic(f, g, **args):
+    if not opt.polys:
         return q.as_basic()
     else:
         return q
 
 def pexquo(f, g, *gens, **args):
-    """Polynomial exact pseudo-quotient of `f` and `g`. """
-    gens = _analyze_gens(gens)
+    """Compute polynomial exact pseudo--quotient of ``f`` and ``g``. """
+    options.allowed_flags(args, ['polys'])
 
     try:
-        F, G = _polify_basic(f, g, *gens, **args)
-    except CoercionFailed, (f, g):
-        raise GeneratorsNeeded("can't compute pseudo quotient of %s and %s without generators" % (f, g))
+        (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args)
+    except PolificationFailed, exc:
+        raise ComputationFailed('pexquo', 2, exc)
 
     q = F.pexquo(G)
 
-    if _should_return_basic(f, g, **args):
+    if not opt.polys:
         return q.as_basic()
     else:
         return q
 
 def div(f, g, *gens, **args):
-    """Polynomial division with remainder of `f` and `g`. """
-    gens = _analyze_gens(gens)
+    """Compute polynomial division of ``f`` and ``g``. """
+    options.allowed_flags(args, ['polys'])
 
     try:
-        F, G = _polify_basic(f, g, *gens, **args)
-    except CoercionFailed, (f, g):
-        raise GeneratorsNeeded("can't compute division of %s and %s without generators" % (f, g))
+        (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args)
+    except PolificationFailed, exc:
+        raise ComputationFailed('div', 2, exc)
 
     q, r = F.div(G)
 
-    if _should_return_basic(f, g, **args):
+    if not opt.polys:
         return q.as_basic(), r.as_basic()
     else:
         return q, r
 
 def rem(f, g, *gens, **args):
-    """Computes polynomial remainder of `f` and `g`. """
-    gens = _analyze_gens(gens)
+    """Compute polynomial remainder of ``f`` and ``g``. """
+    options.allowed_flags(args, ['polys'])
 
     try:
-        F, G = _polify_basic(f, g, *gens, **args)
-    except CoercionFailed, (f, g):
-        raise GeneratorsNeeded("can't compute remainder of %s and %s without generators" % (f, g))
+        (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args)
+    except PolificationFailed, exc:
+        raise ComputationFailed('rem', 2, exc)
 
     r = F.rem(G)
 
-    if _should_return_basic(f, g, **args):
+    if not opt.polys:
         return r.as_basic()
     else:
         return r
 
 def quo(f, g, *gens, **args):
-    """Computes polynomial quotient of `f` and `g`. """
-    gens = _analyze_gens(gens)
+    """Compute polynomial quotient of ``f`` and ``g``. """
+    options.allowed_flags(args, ['polys'])
 
     try:
-        F, G = _polify_basic(f, g, *gens, **args)
-    except CoercionFailed, (f, g):
-        raise GeneratorsNeeded("can't compute quotient of %s and %s without generators" % (f, g))
+        (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args)
+    except PolificationFailed, exc:
+        raise ComputationFailed('quo', 2, exc)
 
     q = F.quo(G)
 
-    if _should_return_basic(f, g, **args):
+    if not opt.polys:
         return q.as_basic()
     else:
         return q
 
 def exquo(f, g, *gens, **args):
-    """Computes polynomial exact quotient of `f` and `g`. """
-    gens = _analyze_gens(gens)
+    """Compute polynomial exact quotient of ``f`` and ``g``. """
+    options.allowed_flags(args, ['polys'])
 
     try:
-        F, G = _polify_basic(f, g, *gens, **args)
-    except CoercionFailed, (f, g):
-        raise GeneratorsNeeded("can't compute quotient of %s and %s without generators" % (f, g))
+        (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args)
+    except PolificationFailed, exc:
+        raise ComputationFailed('exquo', 2, exc)
 
     q = F.exquo(G)
 
-    if _should_return_basic(f, g, **args):
+    if not opt.polys:
         return q.as_basic()
     else:
         return q
 
 def half_gcdex(f, g, *gens, **args):
-    """Half extended Euclidean algorithm of `f` and `g`. """
-    gens = _analyze_gens(gens)
+    """Half extended Euclidean algorithm of ``f`` and ``g``. """
+    options.allowed_flags(args, ['auto', 'polys'])
 
     try:
-        F, G = _polify_basic(f, g, *gens, **args)
-    except CoercionFailed, (f, g):
-        try:
-            return f.half_gcdex(g)
-        except (AttributeError, TypeError): # pragma: no cover
-            raise GeneratorsNeeded("can't compute half extended GCD of %s and %s without generators" % (f, g))
+        (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args)
+    except PolificationFailed, exc:
+        f, g = exc.exprs
 
-    s, h = F.half_gcdex(G, **_filter_args(args, 'auto'))
+        if hasattr(f, 'half_gcdex'):
+            try:
+                return f.half_gcdex(g)
+            except (SympifyError, ValueError):
+                pass
+
+        raise ComputationFailed('half_gcdex', 2, exc)
+
+    s, h = F.half_gcdex(G, auto=opt.auto)
 
-    if _should_return_basic(f, g, **args):
+    if not opt.polys:
         return s.as_basic(), h.as_basic()
     else:
         return s, h
 
 def gcdex(f, g, *gens, **args):
-    """Extended Euclidean algorithm of `f` and `g`. """
-    gens = _analyze_gens(gens)
+    """Extended Euclidean algorithm of ``f`` and ``g``. """
+    options.allowed_flags(args, ['auto', 'polys'])
 
     try:
-        F, G = _polify_basic(f, g, *gens, **args)
-    except CoercionFailed, (f, g):
-        try:
-            return f.gcdex(g)
-        except (AttributeError, TypeError): # pragma: no cover
-            raise GeneratorsNeeded("can't compute extended GCD of %s and %s without generators" % (f, g))
+        (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args)
+    except PolificationFailed, exc:
+        f, g = exc.exprs
 
-    s, t, h = F.gcdex(G, **_filter_args(args, 'auto'))
+        if hasattr(f, 'gcdex'):
+            try:
+                return f.gcdex(g)
+            except (SympifyError, ValueError):
+                pass
+
+        raise ComputationFailed('gcdex', 2, exc)
+
+    s, t, h = F.gcdex(G, auto=opt.auto)
 
-    if _should_return_basic(f, g, **args):
+    if not opt.polys:
         return s.as_basic(), t.as_basic(), h.as_basic()
     else:
         return s, t, h
 
 def invert(f, g, *gens, **args):
-    """Invert `f` modulo `g`, if possible. """
-    gens = _analyze_gens(gens)
+    """Invert ``f`` modulo ``g`` when possible. """
+    options.allowed_flags(args, ['auto', 'polys'])
 
     try:
-        F, G = _polify_basic(f, g, *gens, **args)
-    except CoercionFailed, (f, g):
-        try:
-            return f.invert(g)
-        except (AttributeError, TypeError): # pragma: no cover
-            raise GeneratorsNeeded("can't compute inversion of %s modulo %s without generators" % (f, g))
+        (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args)
+    except PolificationFailed, exc:
+        f, g = exc.exprs
 
-    h = F.invert(G, **_filter_args(args, 'auto'))
+        if hasattr(f, 'invert'):
+            try:
+                return f.invert(g)
+            except (SympifyError, ValueError):
+                pass
+
+        raise ComputationFailed('invert', 2, exc)
+
+    h = F.invert(G, auto=opt.auto)
 
-    if _should_return_basic(f, g, **args):
+    if not opt.polys:
         return h.as_basic()
     else:
         return h
 
 def subresultants(f, g, *gens, **args):
-    """Computes subresultant PRS sequence of `f` and `g`. """
-    gens = _analyze_gens(gens)
+    """Compute subresultant PRS of ``f`` and ``g``. """
+    options.allowed_flags(args, ['polys'])
 
     try:
-        F, G = _polify_basic(f, g, *gens, **args)
-    except CoercionFailed, (f, g):
-        raise GeneratorsNeeded("can't compute subresultants of %s and %s without generators" % (f, g))
+        (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args)
+    except PolificationFailed, exc:
+        raise ComputationFailed('subresultants', 2, exc)
 
     result = F.subresultants(G)
 
-    if _should_return_basic(f, g, **args):
+    if not opt.polys:
         return [ r.as_basic() for r in result ]
     else:
         return result
 
 def resultant(f, g, *gens, **args):
-    """Computes resultant of `f` and `g` via PRS. """
-    gens = _analyze_gens(gens)
+    """Compute resultant of ``f`` and ``g``. """
+    options.allowed_flags(args, ['polys'])
 
     try:
-        F, G = _polify_basic(f, g, *gens, **args)
-    except CoercionFailed, (f, g):
-        raise GeneratorsNeeded("can't compute resultant of %s and %s without generators" % (f, g))
+        (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args)
+    except PolificationFailed, exc:
+        raise ComputationFailed('resultant', 2, exc)
 
     result = F.resultant(G)
 
-    if _should_return_basic(f, g, **args):
+    if not opt.polys:
         return result.as_basic()
     else:
         return result
 
 def discriminant(f, *gens, **args):
-    """Computes discriminant of `f`. """
-    F = NonStrictPoly(f, *_analyze_gens(gens), **args)
+    """Compute discriminant of ``f``. """
+    options.allowed_flags(args, ['polys'])
 
-    if not F.is_Poly:
-        raise GeneratorsNeeded("can't compute discriminant of %s without generators" % f)
+    try:
+        F, opt = poly_from_expr(f, *gens, **args)
+    except PolificationFailed, exc:
+        raise ComputationFailed('discriminant', 1, exc)
 
     result = F.discriminant()
 
-    if _should_return_basic(f, **args):
+    if not opt.polys:
         return result.as_basic()
     else:
         return result
 
 def cofactors(f, g, *gens, **args):
-    """Returns GCD of `f` and `g` and their cofactors. """
-    gens = _analyze_gens(gens)
+    """Compute GCD and cofactors of ``f`` and ``g``. """
+    options.allowed_flags(args, ['polys'])
 
     try:
-        F, G = _polify_basic(f, g, *gens, **args)
-    except CoercionFailed, (f, g):
-        try:
-            return f.cofactors(g)
-        except (AttributeError, TypeError): # pragma: no cover
-            raise GeneratorsNeeded("can't compute cofactors of %s and %s without generators" % (f, g))
+        (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args)
+    except PolificationFailed, exc:
+        f, g = exc.exprs
+
+        if hasattr(f, 'cofactors'):
+            try:
+                return f.cofactors(g)
+            except (SympifyError, ValueError):
+                pass
+
+        raise ComputationFailed('cofactors', 2, exc)
 
     h, cff, cfg = F.cofactors(G)
 
-    if _should_return_basic(f, g, **args):
+    if not opt.polys:
         return h.as_basic(), cff.as_basic(), cfg.as_basic()
     else:
         return h, cff, cfg
 
 def gcd(f, g, *gens, **args):
-    """Returns polynomial GCD of `f` and `g`. """
-    gens = _analyze_gens(gens)
+    """Compute GCD of ``f`` and ``g``. """
+    options.allowed_flags(args, ['polys'])
 
     try:
-        F, G = _polify_basic(f, g, *gens, **args)
-    except CoercionFailed, (f, g):
-        try:
-            return f.gcd(g)
-        except (AttributeError, TypeError): # pragma: no cover
-            raise GeneratorsNeeded("can't compute GCD of %s and %s without generators" % (f, g))
+        (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args)
+    except PolificationFailed, exc:
+        f, g = exc.exprs
+
+        if hasattr(f, 'gcd'):
+            try:
+                return f.gcd(g)
+            except (SympifyError, ValueError):
+                pass
+
+        raise ComputationFailed('gcd', 2, exc)
 
     result = F.gcd(G)
 
-    if _should_return_basic(f, g, **args):
+    if not opt.polys:
         return result.as_basic()
     else:
         return result
 
 def lcm(f, g, *gens, **args):
-    """Returns polynomial LCM of `f` and `g`. """
-    gens = _analyze_gens(gens)
+    """Compute LCM of ``f`` and ``g``. """
+    options.allowed_flags(args, ['polys'])
 
     try:
-        F, G = _polify_basic(f, g, *gens, **args)
-    except CoercionFailed, (f, g):
-        try:
-            return f.lcm(g)
-        except (AttributeError, TypeError): # pragma: no cover
-            raise GeneratorsNeeded("can't compute LCM of %s and %s without generators" % (f, g))
+        (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args)
+    except PolificationFailed, exc:
+        f, g = exc.exprs
+
+        if hasattr(f, 'lcm'):
+            try:
+                return f.lcm(g)
+            except (SympifyError, ValueError):
+                pass
+
+        raise ComputationFailed('lcm', 2, exc)
 
     result = F.lcm(G)
 
-    if _should_return_basic(f, g, **args):
+    if not opt.polys:
         return result.as_basic()
     else:
         return result
 
 def terms_gcd(f, *gens, **args):
-    """Remove GCD of terms from the polynomial `f`. """
+    """Remove GCD of terms from ``f``. """
+    options.allowed_flags(args, ['polys'])
+
     try:
-        f = Poly(f, *_analyze_gens(gens), **args)
-    except GeneratorsNeeded:
-        return f
+        F, opt = poly_from_expr(f, *gens, **args)
+    except PolificationFailed, exc:
+        return exc.expr
 
-    (J, f), dom = f.terms_gcd(), f.get_domain()
+    J, f = F.terms_gcd()
 
-    if dom.has_Ring:
-        if dom.has_Field:
+    if opt.domain.has_Ring:
+        if opt.domain.has_Field:
             denom, f = f.clear_denoms(convert=True)
 
         coeff, f = f.primitive()
 
-        if dom.has_Field:
+        if opt.domain.has_Field:
             coeff /= denom
     else:
-        coeff = 1
+        coeff = S.One
 
     term = Mul(*[ x**j for x, j in zip(f.gens, J) ])
 
     return _keep_coeff(coeff, term*f.as_basic())
 
 def trunc(f, p, *gens, **args):
-    """Reduce `f` modulo a constant `p`. """
-    F = NonStrictPoly(f, *_analyze_gens(gens), **args)
+    """Reduce ``f`` modulo a constant ``p``. """
+    options.allowed_flags(args, ['polys'])
 
-    if not F.is_Poly:
-        return F % p
+    try:
+        F, opt = poly_from_expr(f, *gens, **args)
+    except PolificationFailed, exc:
+        raise ComputationFailed('trunc', 1, exc)
+
+    result = F.trunc(sympify(p))
 
-    if _should_return_basic(f, **args):
-        return F.trunc(p).as_basic()
+    if not opt.polys:
+        return result.as_basic()
     else:
-        return F.trunc(p)
+        return result
 
 def monic(f, *gens, **args):
-    """Divides all coefficients by `LC(f)`. """
-    F = NonStrictPoly(f, *_analyze_gens(gens), **args)
+    """Divide all coefficients of ``f`` by ``LC(f)``. """
+    options.allowed_flags(args, ['auto', 'polys'])
 
-    if not F.is_Poly:
-        raise GeneratorsNeeded("can't compute monic polynomial of %s without generators" % f)
+    try:
+        F, opt = poly_from_expr(f, *gens, **args)
+    except PolificationFailed, exc:
+        raise ComputationFailed('monic', 1, exc)
 
-    G = F.monic(**_filter_args(args, 'auto'))
+    result = F.monic(auto=opt.auto)
 
-    if _should_return_basic(f, **args):
-        return G.as_basic()
+    if not opt.polys:
+        return result.as_basic()
     else:
-        return G
+        return result
 
 def content(f, *gens, **args):
-    """Returns GCD of polynomial coefficients. """
-    F = NonStrictPoly(f, *_analyze_gens(gens), **args)
+    """Compute GCD of coefficients of ``f``. """
+    options.allowed_flags(args, ['polys'])
 
-    if not F.is_Poly:
-        raise GeneratorsNeeded("can't compute content of %s without generators" % f)
-    else:
-        return F.content()
+    try:
+        F, opt = poly_from_expr(f, *gens, **args)
+    except PolificationFailed, exc:
+        raise ComputationFailed('content', 1, exc)
+
+    return F.content()
 
 def primitive(f, *gens, **args):
-    """Returns content and a primitive form of `f`. """
-    F = NonStrictPoly(f, *_analyze_gens(gens), **args)
+    """Compute content and the primitive form of ``f``. """
+    options.allowed_flags(args, ['polys'])
 
-    if not F.is_Poly:
-        raise GeneratorsNeeded("can't compute primitive part of %s without generators" % f)
+    try:
+        F, opt = poly_from_expr(f, *gens, **args)
+    except PolificationFailed, exc:
+        raise ComputationFailed('primitive', 1, exc)
 
     cont, result = F.primitive()
 
-    if _should_return_basic(f, **args):
+    if not opt.polys:
         return cont, result.as_basic()
     else:
         return cont, result
 
 def compose(f, g, *gens, **args):
-    """Returns functional composition `f(g)`. """
-    gens = _analyze_gens(gens)
+    """Compute functional composition ``f(g)``. """
+    options.allowed_flags(args, ['polys'])
 
     try:
-        F, G = _polify_basic(f, g, *gens, **args)
-    except CoercionFailed, (f, g):
-        raise GeneratorsNeeded("can't compute composition of %s and %s without generators" % (f, g))
+        (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args)
+    except PolificationFailed, exc:
+        raise ComputationFailed('compose', 2, exc)
 
     result = F.compose(G)
 
-    if _should_return_basic(f, g, **args):
+    if not opt.polys:
         return result.as_basic()
     else:
         return result
 
 def decompose(f, *gens, **args):
-    """Computes functional decomposition of `f`. """
-    F = NonStrictPoly(f, *_analyze_gens(gens), **args)
+    """Compute functional decomposition of ``f``. """
+    options.allowed_flags(args, ['polys'])
 
-    if not F.is_Poly:
-        raise GeneratorsNeeded("can't compute functional decomposition of %s without generators" % f)
+    try:
+        F, opt = poly_from_expr(f, *gens, **args)
+    except PolificationFailed, exc:
+        raise ComputationFailed('decompose', 1, exc)
 
     result = F.decompose()
 
-    if _should_return_basic(f, **args):
+    if not opt.polys:
         return [ r.as_basic() for r in result ]
     else:
         return result
 
 def sturm(f, *gens, **args):
-    """Computes the Sturm sequence of `f`. """
-    F = NonStrictPoly(f, *_analyze_gens(gens), **args)
+    """Compute Sturm sequence of ``f``. """
+    options.allowed_flags(args, ['auto', 'polys'])
 
-    if not F.is_Poly:
-        raise GeneratorsNeeded("can't compute Sturm sequence of %s without generators" % f)
+    try:
+        F, opt = poly_from_expr(f, *gens, **args)
+    except PolificationFailed, exc:
+        raise ComputationFailed('sturm', 1, exc)
 
-    result = F.sturm(**_filter_args(args, 'auto'))
+    result = F.sturm(auto=opt.auto)
 
-    if _should_return_basic(f, **args):
+    if not opt.polys:
         return [ r.as_basic() for r in result ]
     else:
         return result
 
 def gff_list(f, *gens, **args):
-    """Returns a list of greatest factorial factors of `f`. """
-    F = NonStrictPoly(f, *_analyze_gens(gens), **args)
+    """Compute a list of greatest factorial factors of ``f``. """
+    options.allowed_flags(args, ['polys'])
 
-    if not F.is_Poly:
-        return []
+    try:
+        F, opt = poly_from_expr(f, *gens, **args)
+    except PolificationFailed, exc:
+        raise ComputationFailed('gff_list', 1, exc)
 
     factors = F.gff_list()
 
-    if _should_return_basic(f, **args):
+    if not opt.polys:
         return [ (g.as_basic(), k) for g, k in factors ]
     else:
         return factors
 
 def gff(f, *gens, **args):
-    """Computes greatest factorial factorization of `f`. """
+    """Compute greatest factorial factorization of ``f``. """
     raise NotImplementedError('symbolic falling factorial')
 
 def sqf_norm(f, *gens, **args):
-    """Computes square-free norm of `f`. """
-    F = NonStrictPoly(f, *_analyze_gens(gens), **args)
+    """Compute square--free norm of ``f``. """
+    options.allowed_flags(args, ['polys'])
 
-    if not F.is_Poly:
-        raise GeneratorsNeeded("can't compute square-free norm of %s without generators" % f)
+    try:
+        F, opt = poly_from_expr(f, *gens, **args)
+    except PolificationFailed, exc:
+        raise ComputationFailed('sqf_norm', 1, exc)
 
     s, g, r = F.sqf_norm()
 
-    if _should_return_basic(f, **args):
+    if not opt.polys:
         return Integer(s), g.as_basic(), r.as_basic()
     else:
         return Integer(s), g, r
 
 def sqf_part(f, *gens, **args):
-    """Computes square-free part of `f`. """
-    F = NonStrictPoly(f, *_analyze_gens(gens), **args)
+    """Compute square--free part of ``f``. """
+    options.allowed_flags(args, ['polys'])
 
-    if not F.is_Poly:
-        raise GeneratorsNeeded("can't compute square-free part of %s without generators" % f)
+    try:
+        F, opt = poly_from_expr(f, *gens, **args)
+    except PolificationFailed, exc:
+        raise ComputationFailed('sqf_part', 1, exc)
 
     result = F.sqf_part()
 
-    if _should_return_basic(f, **args):
+    if not opt.polys:
         return result.as_basic()
     else:
         return result
 
 def sqf_list(f, *gens, **args):
-    """Returns a list of square-free factors of `f`. """
-    F = NonStrictPoly(f, *_analyze_gens(gens), **args)
-
-    if not F.is_Poly:
-        raise GeneratorsNeeded("can't compute square-free decomposition of %s without generators" % f)
+    """Compute a list of square--free factors of ``f``. """
+    options.allowed_flags(args, ['all', 'include', 'polys'])
 
-    all = args.get('all', False)
+    try:
+        F, opt = poly_from_expr(f, *gens, **args)
+    except PolificationFailed, exc:
+        raise ComputationFailed('sqf_list', 1, exc)
 
-    if not args.get('include', False):
-        coeff, factors = result = F.sqf_list(all)
+    if not opt.include:
+        coeff, factors = F.sqf_list(all=opt.all)
 
-        if _should_return_basic(f, **args):
+        if not opt.polys:
             return coeff, [ (g.as_basic(), k) for g, k in factors ]
         else:
             return coeff, factors
     else:
-        factors = F.sqf_list_include()
+        factors = F.sqf_list_include(all=opt.all)
 
-        if _should_return_basic(f, **args):
+        if not opt.polys:
             return [ (g.as_basic(), k) for g, k in factors ]
         else:
             return factors
 
-def sqf(f, *gens, **args):
-    """Returns square-free decomposition of `f`. """
-    frac = args.get('frac', False)
-    gens = _analyze_gens(gens)
-
-    def _sqf(f):
-        """Squaqre-free factor a true polynomial expression. """
-        F = NonStrictPoly(f, *gens, **args)
+def _inner_sqf(f):
+    """Helper function for :func:`sqf`. """
+    (coeff, factors), result = f.sqf_list(), S.One
 
-        if not F.is_Poly:
-            return (S.One, F)
+    for g, k in factors:
+        result *= g.as_basic()**k
 
-        (coeff, factors), result = F.sqf_list(), S.One
+    return coeff, result
 
-        for g, k in factors:
-            result *= g.as_basic()**k
+def sqf(f, *gens, **args):
+    """Compute square--free decomposition of ``f``. """
+    options.allowed_flags(args, ['frac', 'polys'])
 
-        return (coeff, result)
+    try:
+        F, opt = poly_from_expr(f, *gens, **args)
+    except PolificationFailed, exc:
+        return exc.expr
 
-    if not frac:
-        coeff, factors = _sqf(f)
+    if not opt.frac:
+        coeff, factors = _inner_sqf(F)
     else:
-        p, q = cancel(f).as_numer_denom()
+        p, q = F
 
-        coeff_p, factors_p = _sqf(p)
-        coeff_q, factors_q = _sqf(q)
+        cp, fp = _inner_sqf(p)
+        cq, fq = _inner_sqf(q)
 
-        coeff = coeff_p / coeff_q
-        factors = factors_p / factors_q
+        coeff, factors = cp/cq, fp/fq
 
     return _keep_coeff(coeff, factors)
 
 def factor_list(f, *gens, **args):
-    """Returns a list of irreducible factors of `f`. """
-    F = NonStrictPoly(f, *_analyze_gens(gens), **args)
+    """Compute a list of irreducible factors of ``f``. """
+    options.allowed_flags(args, ['include', 'polys'])
 
-    if not F.is_Poly:
-        raise GeneratorsNeeded("can't compute factorization of %s without generators" % f)
+    try:
+        F, opt = poly_from_expr(f, *gens, **args)
+    except PolificationFailed, exc:
+        raise ComputationFailed('factor_list', 1, exc)
 
-    if not args.get('include', False):
-        coeff, factors = result = F.factor_list()
+    if not opt.include:
+        coeff, factors = F.factor_list()
 
-        if _should_return_basic(f, **args):
+        if not opt.polys:
             return coeff, [ (g.as_basic(), k) for g, k in factors ]
         else:
             return coeff, factors
     else:
         factors = F.factor_list_include()
 
-        if _should_return_basic(f, **args):
+        if not opt.polys:
             return [ (g.as_basic(), k) for g, k in factors ]
         else:
             return factors
 
-@register_context
-def factor(f, *gens, **args):
-    """Returns factorization into irreducibles of `f`. """
-    frac = args.get('frac', False)
-
-    def _factor(f):
-        """Factor a true polynomial expression. """
-        F = NonStrictPoly(f, *gens, **args)
+def _inner_factor(f):
+    """Helper function for :func:`factor`. """
+    (coeff, factors), result = f.factor_list(), S.One
 
-        if not F.is_Poly:
-            return (S.One, F)
+    for g, k in factors:
+        result *= g.as_basic()**k
 
-        (coeff, factors), result = F.factor_list(), S.One
+    return coeff, result
 
-        for g, k in factors:
-            result *= g.as_basic()**k
+def factor(f, *gens, **args):
+    """Compute factorization into irreducibles of ``f``. """
+    options.allowed_flags(args, ['frac', 'polys'])
 
-        return (coeff, result)
+    try:
+        F, opt = poly_from_expr(f, *gens, **args)
+    except PolificationFailed, exc:
+        return exc.expr
 
-    if not frac:
-        coeff, factors = _factor(f)
+    if not opt.frac:
+        coeff, factors = _inner_factor(F)
     else:
-        p, q = cancel(f).as_numer_denom()
+        p, q = F
 
-        coeff_p, factors_p = _factor(p)
-        coeff_q, factors_q = _factor(q)
+        cp, fp = _inner_factor(p)
+        cq, fq = _inner_factor(q)
 
-        coeff = coeff_p / coeff_q
-        factors = factors_p / factors_q
+        coeff, factors = cp/cq, fp/fq
 
     return _keep_coeff(coeff, factors)
 
 def intervals(F, all=False, eps=None, inf=None, sup=None, strict=False, fast=False, sqf=False):
-    """Compute isolating intervals for roots of `f`. """
+    """Compute isolating intervals for roots of ``f``. """
     if not hasattr(F, '__iter__'):
         try:
             F = Poly(F)
@@ -2292,7 +2401,7 @@ def refine_root(f, s, t, eps=None, steps=None, fast=False, check_sqf=False):
     return F.refine_root(s, t, eps=eps, steps=steps, fast=fast, check_sqf=check_sqf)
 
 def count_roots(f, inf=None, sup=None):
-    """Return the number of roots of ``f`` in ``[inf, sup]`` interval. """
+    """Compute the number of roots of ``f`` in ``[inf, sup]`` interval. """
     try:
         F = Poly(f, greedy=False)
     except GeneratorsNeeded:
@@ -2301,7 +2410,7 @@ def count_roots(f, inf=None, sup=None):
     return F.count_roots(inf=inf, sup=sup)
 
 def real_roots(f, multiple=True):
-    """Return a list of real roots with multiplicities of ``f``. """
+    """Compute a list of real roots with multiplicities of ``f``. """
     try:
         F = Poly(f, greedy=False)
     except GeneratorsNeeded:
@@ -2310,7 +2419,7 @@ def real_roots(f, multiple=True):
     return F.real_roots(multiple=multiple)
 
 def nroots(f, maxsteps=50, cleanup=True, error=False):
-    """Compute numerical approximations of roots of `f`. """
+    """Compute numerical approximations of roots of ``f``. """
     try:
         F = Poly(f, greedy=False)
     except GeneratorsNeeded:
@@ -2319,7 +2428,9 @@ def nroots(f, maxsteps=50, cleanup=True, error=False):
     return F.nroots(maxsteps=maxsteps, cleanup=cleanup, error=error)
 
 def cancel(f, *gens, **args):
-    """Cancel common factors in a rational function `f`.  """
+    """Cancel common factors in a rational function ``f``.  """
+    options.allowed_flags(args, ['polys'])
+
     f = sympify(f)
 
     if type(f) is not tuple:
@@ -2330,11 +2441,9 @@ def cancel(f, *gens, **args):
     else:
         p, q = f
 
-    gens = _analyze_gens(gens)
-
     try:
-        F, G = _polify_basic(p, q, *gens, **args)
-    except CoercionFailed:
+        (F, G), opt = parallel_poly_from_expr((p, q), *gens, **args)
+    except PolificationFailed, exc:
         if type(f) is not tuple:
             return f
         else:
@@ -2345,13 +2454,13 @@ def cancel(f, *gens, **args):
     if type(f) is not tuple:
         return c*(P.as_basic()/Q.as_basic())
     else:
-        if _should_return_basic(p, q, **args):
+        if not opt.polys:
             return c, P.as_basic(), Q.as_basic()
         else:
             return c, P, Q
 
 def reduced(f, G, *gens, **args):
-    """Reduces a polynomial `f` modulo a set of polynomials `G`. """
+    """Reduce a polynomial ``f`` modulo a set of polynomials ``G``. """
     polys, opt = parallel_poly_from_expr([f] + list(G), *gens, **args)
 
     for i, poly in enumerate(polys):
@@ -2370,7 +2479,7 @@ def reduced(f, G, *gens, **args):
         return Q, r
 
 def groebner(F, *gens, **args):
-    """Computes reduced Groebner basis for a set of polynomials. """
+    """Compute a reduced Groebner basis for a set of polynomials. """
     args = _update_args(args, 'field', True)
 
     polys, opt = parallel_poly_from_expr(F, *gens, **args)
@@ -2390,16 +2499,14 @@ def groebner(F, *gens, **args):
 
 def symmetrize(f, *gens, **args):
     """Rewrite a polynomial in terms of elementary symmetric polynomials. """
+    options.allowed_flags(args, ['formal'])
+
     try:
-        f = Poly(f, *_analyze_gens(gens), **args)
-    except GeneratorsNeeded:
-        if args.get('formal', False):
-            return (f, S.Zero, {})
-        else:
-            return (f, S.Zero)
+        f, opt = poly_from_expr(f, *gens, **args)
+    except PolificationFailed, exc:
+        return ComputationFailed('symmetrize', 1, exc)
 
     from sympy.polys.specialpolys import symmetric_poly
-
     polys, symbols = [], numbered_symbols('s', start=1)
 
     gens, dom = f.gens, f.get_domain()
@@ -2451,17 +2558,19 @@ def symmetrize(f, *gens, **args):
 
     polys = [ (s, p.as_basic()) for s, p in polys ]
 
-    if args.get('formal', False):
+    if opt.formal:
         return (Add(*symmetric), f.as_basic(), dict(polys))
     else:
         return (Add(*symmetric).subs(polys), f.as_basic())
 
 def horner(f, *gens, **args):
-    """Apply Horner's rule to put a polynomial in Horner form. """
-    F = NonStrictPoly(f, *_analyze_gens(gens), **args)
+    """Rewrite a polynomial in Horner form. """
+    options.allowed_flags(args, [])
 
-    if not F.is_Poly:
-        return F
+    try:
+        F, opt = poly_from_expr(f, *gens, **args)
+    except PolificationFailed, exc:
+        return exc.expr
 
     form, gen = S.Zero, F.gen
 
diff --git a/sympy/polys/polyutils.py b/sympy/polys/polyutils.py
index 60c98eb003..ad46ce2327 100644
--- a/sympy/polys/polyutils.py
+++ b/sympy/polys/polyutils.py
@@ -21,7 +21,7 @@
 def _sort_gens(gens, **args):
     """Sort generators in a reasonably intelligent way. """
     from sympy.polys.polyoptions import build_options
-    opt = build_options(**args)
+    opt = build_options(args)
 
     gens_order, wrt = {}, None
 
@@ -272,7 +272,7 @@ def _dict_from_basic_no_gens(expr, opt):
 def parallel_dict_from_basic(exprs, **args):
     """Transform expressions into a multinomial form. """
     from sympy.polys.polyoptions import build_options
-    opt = build_options(**args)
+    opt = build_options(args)
 
     if opt.expand is not False:
         exprs = [ expr.expand() for expr in exprs ]
@@ -285,7 +285,7 @@ def parallel_dict_from_basic(exprs, **args):
 def dict_from_basic(expr, **args):
     """Transform an expression into a multinomial form. """
     from sympy.polys.polyoptions import build_options
-    opt = build_options(**args)
+    opt = build_options(args)
 
     if opt.expand is not False:
         expr = expr.expand()
diff --git a/sympy/polys/tests/test_polyoptions.py b/sympy/polys/tests/test_polyoptions.py
index 44f3ee6875..6c236a8852 100644
--- a/sympy/polys/tests/test_polyoptions.py
+++ b/sympy/polys/tests/test_polyoptions.py
@@ -3,7 +3,7 @@
 from sympy.polys.polyoptions import (
     Options, Expand, Gens, Wrt, Sort, Order, Field, Greedy, Domain,
     Split, Gaussian, Extension, Modulus, Symmetric, Strict, Auto,
-    Frac, Formal, Polys, Include, Monic, Gen,
+    Frac, Formal, Polys, Include, Monic, All, Gen,
 )
 
 from sympy.polys.monomialtools import monomial_lex_key
@@ -374,6 +374,21 @@ def test_Monic_postprocess():
 
     assert opt == {'monic': True}
 
+def test_All_preprocess():
+    assert All.preprocess(False) is False
+    assert All.preprocess(True) is True
+
+    assert All.preprocess(0) is False
+    assert All.preprocess(1) is True
+
+    raises(OptionError, "All.preprocess(x)")
+
+def test_Monic_postprocess():
+    opt = {'all': True}
+    All.postprocess(opt)
+
+    assert opt == {'all': True}
+
 def test_Gen_preprocess():
     pass
 
diff --git a/sympy/polys/tests/test_polytools.py b/sympy/polys/tests/test_polytools.py
index cabfc51d9b..ed26d497c1 100644
--- a/sympy/polys/tests/test_polytools.py
+++ b/sympy/polys/tests/test_polytools.py
@@ -1,7 +1,7 @@
 """Tests for user-friendly public interface to polynomial functions. """
 
 from sympy.polys.polytools import (
-    Poly, poly, _polify_basic,
+    Poly, poly,
     parallel_poly_from_expr,
     init_poly_from_dict,
     init_poly_from_list,
@@ -34,6 +34,7 @@
     OperationNotSupported,
     ExactQuotientFailed,
     PolificationFailed,
+    ComputationFailed,
     UnificationFailed,
     RefinementFailed,
     GeneratorsNeeded,
@@ -322,6 +323,14 @@ def test_Poly__gens():
     assert Poly((x-p)*(x-q), wrt='p', sort='q > x').gens == (p, q, x)
     assert Poly((x-p)*(x-q), wrt='q', sort='p > x').gens == (q, p, x)
 
+def test_Poly_zero():
+    assert Poly(x).zero == Poly(0, x, domain=ZZ)
+    assert Poly(x/2).zero == Poly(0, x, domain=QQ)
+
+def test_Poly_one():
+    assert Poly(x).one == Poly(1, x, domain=ZZ)
+    assert Poly(x/2).one == Poly(1, x, domain=QQ)
+
 def test_Poly__unify():
     raises(UnificationFailed, "Poly(x)._unify(y)")
 
@@ -829,7 +838,7 @@ def test_Poly_degree():
     assert degree(x*y**2, x, y) == 1
     assert degree(x*y**2, y, x) == 2
 
-    raises(GeneratorsNeeded, "degree(1)")
+    raises(ComputationFailed, "degree(1)")
 
 def test_Poly_degree_list():
     assert Poly(0, x).degree_list() == (-1,)
@@ -847,7 +856,7 @@ def test_Poly_degree_list():
 
     assert degree_list(x*y**2) == (1,2)
 
-    raises(GeneratorsNeeded, "degree_list(1)")
+    raises(ComputationFailed, "degree_list(1)")
 
 def test_Poly_total_degree():
     assert Poly(x**2*y+x**3*z**2+1).total_degree() == 6
@@ -1050,34 +1059,6 @@ def test_parallel_poly_from_expr():
 
     raises(PolificationFailed, "parallel_poly_from_expr([0, 1])")
 
-def test__polify_basic():
-    assert _polify_basic(x-1, x**2-1, x) == (Poly(x-1, x), Poly(x**2-1, x))
-    assert _polify_basic(Poly(x-1, x), x**2-1, x) == (Poly(x-1, x), Poly(x**2-1, x))
-    assert _polify_basic(x-1, Poly(x**2-1, x), x) == (Poly(x-1, x), Poly(x**2-1, x))
-    assert _polify_basic(Poly(x-1, x), Poly(x**2-1, x), x) == (Poly(x-1, x), Poly(x**2-1, x))
-
-    assert _polify_basic(x-1, x**2-1, x, y) == (Poly(x-1, x, y), Poly(x**2-1, x, y))
-    assert _polify_basic(Poly(x-1, x), x**2-1, x, y) == (Poly(x-1, x, y), Poly(x**2-1, x, y))
-    assert _polify_basic(x-1, Poly(x**2-1, x), x, y) == (Poly(x-1, x, y), Poly(x**2-1, x, y))
-    assert _polify_basic(Poly(x-1, x), Poly(x**2-1, x), x, y) == (Poly(x-1, x, y), Poly(x**2-1, x, y))
-
-    assert _polify_basic(x-1, x**2-1) == (Poly(x-1, x), Poly(x**2-1, x))
-    assert _polify_basic(Poly(x-1, x), x**2-1) == (Poly(x-1, x), Poly(x**2-1, x))
-    assert _polify_basic(x-1, Poly(x**2-1, x)) == (Poly(x-1, x), Poly(x**2-1, x))
-    assert _polify_basic(Poly(x-1, x), Poly(x**2-1, x)) == (Poly(x-1, x), Poly(x**2-1, x))
-
-    assert _polify_basic(1, x**2-1) == (Poly(1, x), Poly(x**2-1, x))
-    assert _polify_basic(1, x**2-1) == (Poly(1, x), Poly(x**2-1, x))
-    assert _polify_basic(1, Poly(x**2-1, x)) == (Poly(1, x), Poly(x**2-1, x))
-    assert _polify_basic(1, Poly(x**2-1, x)) == (Poly(1, x), Poly(x**2-1, x))
-
-    assert _polify_basic(x**2-1, 1) == (Poly(x**2-1, x), Poly(1, x))
-    assert _polify_basic(x**2-1, 1) == (Poly(x**2-1, x), Poly(1, x))
-    assert _polify_basic(Poly(x**2-1, x), 1) == (Poly(x**2-1, x), Poly(1, x))
-    assert _polify_basic(Poly(x**2-1, x), 1) == (Poly(x**2-1, x), Poly(1, x))
-
-    raises(CoercionFailed, "_polify_basic(1, 2)")
-
 def test_pdiv():
     f, g = x**2 - y**2, x - y
     q, r = x + y, 0
@@ -1119,10 +1100,10 @@ def test_pdiv():
     assert pquo(F, G, polys=False) == q
     assert prem(F, G, polys=False) == r
 
-    raises(GeneratorsNeeded, "pdiv(4, 2)")
-    raises(GeneratorsNeeded, "pexquo(4, 2)")
-    raises(GeneratorsNeeded, "pquo(4, 2)")
-    raises(GeneratorsNeeded, "prem(4, 2)")
+    raises(ComputationFailed, "pdiv(4, 2)")
+    raises(ComputationFailed, "pexquo(4, 2)")
+    raises(ComputationFailed, "pquo(4, 2)")
+    raises(ComputationFailed, "prem(4, 2)")
 
 def test_div():
     f, g = x**2 - y**2, x - y
@@ -1165,10 +1146,10 @@ def test_div():
     assert quo(F, G, polys=False) == q
     assert rem(F, G, polys=False) == r
 
-    raises(GeneratorsNeeded, "div(4, 2)")
-    raises(GeneratorsNeeded, "exquo(4, 2)")
-    raises(GeneratorsNeeded, "quo(4, 2)")
-    raises(GeneratorsNeeded, "rem(4, 2)")
+    raises(ComputationFailed, "div(4, 2)")
+    raises(ComputationFailed, "exquo(4, 2)")
+    raises(ComputationFailed, "quo(4, 2)")
+    raises(ComputationFailed, "rem(4, 2)")
 
 def test_gcdex():
     f, g = 2*x, x**2 - 16
@@ -1230,7 +1211,7 @@ def test_subresultants():
     assert subresultants(f, g, polys=True) == [F, G, H]
     assert subresultants(F, G, polys=False) == [f, g, h]
 
-    raises(GeneratorsNeeded, "subresultants(4, 2)")
+    raises(ComputationFailed, "subresultants(4, 2)")
 
 def test_resultant():
     f, g, h = x**2 - 2*x + 1, x**2 - 1, 0
@@ -1255,7 +1236,7 @@ def test_resultant():
     assert resultant(f, g, polys=True) == H
     assert resultant(F, G, polys=False) == h
 
-    raises(GeneratorsNeeded, "resultant(4, 2)")
+    raises(ComputationFailed, "resultant(4, 2)")
 
 def test_discriminant():
     f, g = x**3 + 3*x**2 + 9*x - 13, -11664
@@ -1280,7 +1261,7 @@ def test_discriminant():
     assert discriminant(f, polys=True) == G
     assert discriminant(F, polys=False) == g
 
-    raises(GeneratorsNeeded, "discriminant(4)")
+    raises(ComputationFailed, "discriminant(4)")
 
 def test_gcd():
     f, g = x**3 - 1, x**2 - 1
@@ -1410,7 +1391,7 @@ def test_monic():
     assert monic(f, polys=True) == G
     assert monic(F, polys=False) == g
 
-    raises(GeneratorsNeeded, "monic(4)")
+    raises(ComputationFailed, "monic(4)")
 
     assert monic(2*x**2 + 6*x + 4, auto=False) == x**2 + 3*x + 2
     raises(ExactQuotientFailed, "monic(2*x + 6*x + 1, auto=False)")
@@ -1424,7 +1405,7 @@ def test_content():
     F.content() == 2
     content(f) == 2
 
-    raises(GeneratorsNeeded, "content(4)")
+    raises(ComputationFailed, "content(4)")
 
     f = Poly(2*x, modulus=3)
 
@@ -1442,7 +1423,7 @@ def test_primitive():
     assert primitive(f, polys=True) == (2, G)
     assert primitive(F, polys=False) == (2, g)
 
-    raises(GeneratorsNeeded, "primitive(4)")
+    raises(ComputationFailed, "primitive(4)")
 
     f = Poly(2*x, modulus=3)
 
@@ -1471,8 +1452,8 @@ def test_compose():
     assert decompose(f, polys=True) == [G, H]
     assert decompose(F, polys=False) == [g, h]
 
-    raises(GeneratorsNeeded, "compose(4, 2)")
-    raises(GeneratorsNeeded, "decompose(4)")
+    raises(ComputationFailed, "compose(4, 2)")
+    raises(ComputationFailed, "decompose(4)")
 
     assert compose(x**2 - y**2, x - y, x, y) ==  x**2 - 2*x*y
     assert compose(x**2 - y**2, x - y, y, x) == -y**2 + 2*x*y
@@ -1489,7 +1470,7 @@ def test_sturm():
     assert sturm(f, polys=True) == [F, G]
     assert sturm(F, polys=False) == [f, g]
 
-    raises(GeneratorsNeeded, "sturm(4)")
+    raises(ComputationFailed, "sturm(4)")
     raises(DomainError, "sturm(f, auto=False)")
 
     f = Poly(S(1024)/(15625*pi**8)*x**5   \
@@ -1562,8 +1543,8 @@ def test_sqf():
     assert F.sqf_list_include() == [(G, 1), (H, 2)]
     assert sqf_list(f, include=True) == [(g, 1), (h, 2)]
 
-    raises(GeneratorsNeeded, "sqf_part(4)")
-    raises(GeneratorsNeeded, "sqf_list(4)")
+    raises(ComputationFailed, "sqf_part(4)")
+    raises(ComputationFailed, "sqf_list(4)")
 
     assert sqf(1) == 1
     assert sqf(1, frac=True) == 1
@@ -1585,6 +1566,12 @@ def test_sqf():
     assert sqf(6*x - 10) == Mul(2, 3*x - 5, evaluate=False)
 
     assert sqf((6*x - 10)/(3*x - 6), frac=True) == S(2)/3*((3*x - 5)/(x - 2))
+    assert sqf(Poly(x**2 - 2*x + 1), frac=True) == (x - 1)**2
+
+    f = 3 + x - x*(1 + x) + x**2
+
+    assert sqf(f) == 3
+    assert sqf(f, frac=True) == 3
 
 def test_factor():
     f = x**5 - x**3 - x**2 + 1
@@ -1606,7 +1593,7 @@ def test_factor():
     assert F.factor_list_include() == [(U, 1), (V, 2), (W, 1)]
     assert factor_list(f, include=True) == [(u, 1), (v, 2), (w, 1)]
 
-    raises(GeneratorsNeeded, "factor_list(4)")
+    raises(ComputationFailed, "factor_list(4)")
 
     assert factor(1) == 1
     assert factor(1, frac=True) == 1
@@ -1662,6 +1649,13 @@ def test_factor():
     assert factor(f) == y*(sin(x) + 1)/(pi**2 + 2*pi + 1)
     assert factor(f, frac=True) == y*(sin(x) + 1)/(pi + 1)**2
 
+    assert factor(Poly(x**2 - 2*x + 1), frac=True) == (x - 1)**2
+
+    f = 3 + x - x*(1 + x) + x**2
+
+    assert factor(f) == 3
+    assert factor(f, frac=True) == 3
+
 def test_intervals():
     assert intervals(0) == []
     assert intervals(1) == []
@@ -1992,15 +1986,9 @@ def test_groebner():
     assert sum([ q*g for q, g in zip(Q, G)]) + r == Poly(f, modulus=7)
 
 def test_symmetrize():
-    assert symmetrize(0) == (0, 0)
-    assert symmetrize(1) == (1, 0)
-
     assert symmetrize(0, x, y, z) == (0, 0)
     assert symmetrize(1, x, y, z) == (1, 0)
 
-    assert symmetrize(0, formal=True) == (0, 0, {})
-    assert symmetrize(1, formal=True) == (1, 0, {})
-
     s1 = x + y + z
     s2 = x*y + x*z + y*z
     s3 = x*y*z