What is the best way to work with ratios of polynomials?

edit

What is the best way to work with ratios of polynomials with (real and complex) coefficients that can vary from very small to very big values? As I have to keep precision, I am working with symbolic expressions (type 'sage.symbolic.expression.Expression'), but as I operate with the expressions, they become bigger and bigger, and it is more and more difficult to plot them or find its roots, for example.

I perceived that making

# symbolic_expression is some symbolic polynomial with real or complex coefficients
symbolic_expression.polynomial(RR)
symbolic_expression.polynomial(CC)

simplifies the polynomials and eases computations but, as I understand, this converts the symbolic expressions to the real (RR) and complex (CC) rings, rounding the coefficients of the polynomials to machine precision and, in this manner, numerical errors would propagate in my computations. In this scenario, I would need to increase this conversion precision to an arbitrary value that guarantees correct results at the end of all computations.

I've read something about polynomial rings, symbolic ring, etc. I am studying this, but I still don't know what a ring is. Is symbolic ring the same thing as symbolic expressions? Does polynomial rings enable better performance and arbitrary precision?

Cordially

EDIT

I complement my question to answer to @B r u n o comment.

First, I would like to thank @B r u n o and @slelievre for the tips! I start to better understand things now.

Basically, I generate three polynomials. These polynomials are the basis of my work. I then operate on these polynomials generating other functions, e.g., as the ratio of these polynomials, and continue operating on the new functions, transforming them, plotting, finding the roots of the numerator and denominator polynomials, etc.

I'm slowly improving the code. I started with a plain script, then I created some functions to organize the code, in a procedural approach, and now I am moving to an object-oriented paradigm. This helps me better visualize and understand the problem. I even created a graphical interface to ease the interaction with the program, so as I can easily experiment with different input variables.

The code snippet below reflects the current state of the code it illustrates some inputs to test the code, how I generate my three basic polynomials and some derived functions. Sorry if the code seems silly :-/, I'm not used to develop software, I am still studying.

class FilterPrototype:

    def __init(self):
        pass

    def generate_polynomials(self):
        pass

    def generate_s21(self, e_poly, p_poly):
        return p_poly.polynomial(CC)/e_poly.polynomial(CC)

    def generate_s11(self, e_poly, f_poly):
        return f_poly.polynomial(CC)/e_poly.polynomial(CC)

    def generate_group_delay(self, s21):
        return -(1/s21(s=i*omega)*s21(s=i*omega).diff(omega)).imag()

    def plot(self, symb_express, bounds):
        pl = plot(symb_express, bounds[0], bounds[1])
        pl.show()
        return list(pl[0])

class ApproxArbitraryPhase(FilterPrototype):

    def __init__(self, spec_type, spec, e_order, p_order):
        self.spec_type = spec_type # type of specification
        self.spec = spec           # filter specification
        self.e_order = e_order     # order of the denominator polynomial
        self.p_order = p_order     # order of the numerator polynomial

    def generate_polynomials(self):
        listOmega = [] # List of frequencies (in radians per second)
        for a in range(self.e_order + 1):
            listOmega.append(a/self.e_order)

        if self.spec_type == "group delay":
            group_delay_spec(omega) = sage_eval(self.spec, locals={'omega':omega})
        elif self.spec_type == "phase":
            root.destroy()
            sys.exit("You cannot specify a phase. This is yet not implemented.")
        else:
            root.destroy()
            sys.exit("""The spec_type should be "group delay" or "phase".""")

        phase(omega) = group_delay_spec(omega).integral(omega)
        listPhases = [0]
        for a in range(1, self.e_order+1):
            listPhases.append(phase(listOmega[a]))

        alpha0 = listOmega[1]/tan(listPhases[1])
        listAlpha = [alpha0]
        for a in range(1, self.e_order):
            firstline = listOmega[a+1]^2 - listOmega[a]^2
            lastline = alpha0 - listOmega[a+1]/tan(listPhases[a+1])
            ic = 1
            while ic < a:
                lastline = listAlpha[ic] - ((listOmega[a+1]^2 - listOmega[ic]^2)/lastline)
                ic+=1
            listAlpha.append(firstline/lastline)

        ePolynomials = [1, s + listAlpha[0]]
        for a in range(2, self.e_order+1):
            ePolynomials.append(listAlpha[a-1]*ePolynomials[a-1] + (s^2 + listOmega[a-1]^2)*ePolynomials[a-2])

        ePoly = ePolynomials[self.e_order]
        pMod = sqrt(ePoly(s=i*omega)*ePoly(s=-i*omega))
        integrand = []
        for a in range(self.p_order+1):
            integrand.append(pMod*chebyshev_T(a, omega)/sqrt(1-omega^2))

        c = []
        for a in range(self.p_order+1):
            if a == 0:
                c.append((1/pi)*Rational(integrand[a].nintegral(omega, -1, 1)[0]))
            else:
                c.append((2/pi)*Rational(integrand[a].nintegral(omega, -1, 1)[0]))

        up = 0
        for a in range(self.p_order+1):
            up = up + c[a]*chebyshev_T(a, omega)

        eMod(omega) = (ePoly(s=i*omega)).polynomial(CC)
        den = lambda omega: up/abs(eMod(omega))
        pl = plot(den(omega), -1, 1) # to improve the precision of norm_factor, the plot_points=200 parameter of the plot function should be augmented
        plPoints = []
        for pair in list(pl[0]):
            plPoints.append(pair[1])
        normFactor = Rational(max(plPoints))

        pPoly(omega) = up/normFactor
        pPoly(s) = pPoly(omega=s/I)

        fMod2 = ePoly(s=s)*ePoly(s=-s) - pPoly(s=s)*pPoly(s=-s)

        solutions = fMod2.roots(ring=CC)
        fPoly = 0
        for solution in solutions:
            if solution[0].real() < 0 or (solution[0].real() == 0 and solution[0].imag() < 0):
                if fPoly == 0:
                    fPoly = (s - solution[0])
                else:
                    fPoly = fPoly*(s - solution[0])

        return (ePoly, fPoly, pPoly)


if __name__ == "__main__":
    # execute only if run as a script
    s = var("s")
    omega = var("omega", domain=RR)

    filtro = ApproxArbitraryPhase("group delay", "omega + 5/2", 4, 4)
#     filtro = ApproxArbitraryPhase("group delay", "omega + 9/2", 6, 6)
#     filtro = ApproxArbitraryPhase("group delay", "-omega + 27/5", 6, 6)
#     filtro = ApproxArbitraryPhase("group delay", "piecewise([([0,0.5], 2*omega + 143/20), ((0.5,1), 163/20), ([1,1], 163/20)])", 8, 8)
#     filtro = ApproxArbitraryPhase("group delay", "3/5*omega + 106/25", 6, 4)
    e,f,p = filtro.generate_polynomials()
    print(e)
    print(f)
    print(p(s))
    s21 = filtro.generate_s21(e, p(s))
    print(s21)
    s11 = filtro.generate_s11(e, f) # possui polinomio do numerador com coeficientes imaginarios
    print(s11)
    gd = filtro.generate_group_delay(s21)
    print(gd)
    points_s21 = filtro.plot(20*log(abs(s21(s=i*omega)),10), [-1.5, 1.5])
#     print(points_s21)
    points_s11 = filtro.plot(20*log(abs(s11(s=i*omega)),10), [-1.5, 1.5])
#     print(points_s11)
    points_gd = filtro.plot(gd, [-1.5, 1.5])
#     print(points_gd)