experimentalCAD
system:sage

{{{id=31|
reset()
///
}}}

{{{id=27|
def leadcoeff(pol,variable):
    r"""
    auxiliary function for computing the leading coefficient in the type of ring used
    """
    return pol.coefficient(variable^pol.degree(variable))
///
}}}

{{{id=45|
def sqrfr_factor(pol,x,basering):
    r"""
    returns the square free factors of the input polynomial pol
    """
    F=[]
    a1 = pol
    a2 = a1.gcd(a1.derivative(x))
    c1 = a1//a2
    a3 = a2.gcd(a2.derivative(x))
    c2 = a2//a3
    while c2!=1:
        a2=a3
        a3=a3.gcd(a3.derivative(x))
        c1=c2
        c2=a2//a3
        F.append(c1//c2)
    retF = []
    for i in range(len(F)):
        if basering(F[i]).degree()>0:
            retF.append(F[i])
    return retF
///
}}}

{{{id=79|
def getAklevel(polys,x):
    aklevel = []
    for i in range(len(polys)):
        if polys[i].variables()==[x]:
            aklevel.append(polys[i])
    return aklevel
///
}}}

{{{id=80|
def getroots(polys,basering):
    proots = []
    for i in range(len(polys)):
        proots.append(basering(polys[i]).roots())
    return proots
///
}}}

{{{id=0|
def pfs(polys,vars,basering):
    r"""
    input: list of polynomials and variables (needed for specifying the order) and basering (not sure whether needed)
    output: projection factor set w.r.t. the given order
    """
    inpolys = polys[:]
    outpolys = polys[:]
    z = vars[0]
    vars.remove(z)
    while len(vars)>0:    
# --------------------
#    k-level sets = set of polynomials with degree > 0 in the variable x[k]       
# --------------------
        aklevel = []
        help = polys[:]
        for i in range(len(polys)):
            if polys[i].degree(z)>0:
               aklevel.append(polys[i])
               help.remove(polys[i])
        polys = help[:]
# --------------------
#    compute projection wrt z using reduced brown-mccallum projection operator
#    projection-operator: proj(aklevel) = union(discr(p,xk),lc(p,xk),res(p,q,xk)) (all p and q in aklevel)        
# --------------------
        projak = []
        for i in range(len(aklevel)):
            m = aklevel[i].degree(z)
            discxk = aklevel[i].resultant(derivative(aklevel[i],z),z)*(-1)^(m*(m-1)/2)//leadcoeff(aklevel[i],z)
            projak.append(discxk)
            projak.append(leadcoeff(aklevel[i],z))
            for j in range(1,len(aklevel)):
                projak.append(aklevel[i].resultant(aklevel[j],z))
# --------------------
# update polynomials (k-level polynomials already removed - add projection factors)
# TODO: simplify polynomials (e.g.: [1-x^2] instead of [-4+4*x^2,x^2-1])
# --------------------             
        z = vars[0]
        vars.remove(z)
        pak = []
        for i in range(len(projak)):
            if basering(projak[i]).degree(z)>0:
                   sf = sqrfr_factor(projak[i],z,basering)
                   if sf == []:
                       pak.append(projak[i])
                   else:
                       pak = union(pak,sf)
        polys = union(polys,pak)
        outpolys = union(outpolys,pak)
    return union(outpolys)
///
}}}

{{{id=23|
R = PolynomialRing(QQ,(['x','y']))
R.inject_variables()
firstex = pfs([x^2+y^2-4,(x-1)*(y-1)-1],[y,x],R)
///
}}}

{{{id=69|
firstex
///
}}}

{{{id=72|
sum([implicit_plot(f,(x,-4,4),(y,-4,4)) for f in firstex]).show(aspect_ratio=1)
///
}}}

{{{id=14|
# three examples taken from Arnon, Collins + McCallum
///
}}}

{{{id=30|
R = PolynomialRing(QQ,(['x','y']))
R.inject_variables()
acc1 = [y^4-2*y^3+y^2-3*x^2*y+2*x^4]
acc1pfs = pfs(acc1,[x,y],R)
///
}}}

{{{id=53|
acc1pfs
///
}}}

{{{id=73|
sum([implicit_plot(f,(x,-2,2),(y,-0.1,2.4),plot_points=170) for f in acc1pfs]).show(aspect_ratio=1)
///
}}}

{{{id=29|
R = PolynomialRing(QQ,(['x','y','z']))
R.inject_variables()
ball = [x^2+y^2+z^2-1]
ballpfs = pfs(ball,[z,y,x],R)
///
}}}

{{{id=76|
ballpfs
///
}}}

{{{id=77|
implicit_plot3d(ballpfs[0],(x,-2,2),(y,-2,2),(z,-2,2))+implicit_plot3d(ballpfs[1],(x,-2,2),(y,-2,2),(z,-2,2))+implicit_plot3d(ballpfs[2],(x,-2,2),(y,-2,2),(z,-2,2))
///
}}}

{{{id=78|

///
}}}