### ITEG / imult(f,g)

from PyM import *

'''
The function imult(f,g,O) computes the intersection
multiplicity of the plane curves f(x,y)=0 and g(x,y)=0
at the point O =[ a,b], which by default is the origin O = [0,0].
The algoritm is extracted from Fulton's book Algebraic Curves 
(see http://www.math.lsa.umich.edu/~wfulton/CurveBook.pdf,
Section 3.3, including the Example on page 40).
The auxiliary functions used in the body of imult, namely
cyclic_matrix(v,k), vector_resultant(v,w) and resultant(f,g),
are defined first.
'''

def cyclic_matrix(v,k):
    K = K_(v)
    if k < 0:
        return 'cyclic_matrix: number of rows must be non negative'
    if k == 0:
        return null_matrix(K)
    r = list(v)+(k-1)*[0>>K]
    C = [r]
    for _ in range(1,k):
        r = [0>>K] + r[:-1]
        C += [r]
    return matrix(C)
    
def vector_resultant(v,w):
    n = len(v)
    m = len(w)
    Cv = cyclic_matrix(v,m-1)
    Cw = cyclic_matrix(w,n-1)
    R = stack(Cv,Cw)
    return det(R)

def resultant(f,g, vars = None):
    if vars == None:
        V = variables(product(variables(f)+variables(g)))
        if len(V) == 2:
            x,y = V
        elif len(V)==1:
            v = coeffs_inc(f)
            w = coeffs_inc(g)
            return vector_resultant(v,w)
        elif len(V) < 1:
            return 0
        else:
            x = V[-2]
            y = V[-1]
    else:
        if len(vars) == 1:
            x = vars[0]
            v = coeffs_inc(f,x)
            w = coeffs_inc(g,x)
            return vector_resultant(v,w)
        else:
            x = vars[-2]
            y = vars[-1]
    m = degree(f,y); n=degree(g,y)
    v = coeffs_inc(f,y); w = coeffs_inc(g,y)
    return vector_resultant(v,w)
    
    

def imult(f,g,O=[0,0]):
    if f==0 or g==0: return 'Infinity'
    if len(O)==3:
        a,b,c = O
        if c!=0: O = [a/c,b/c]
        else: return imultinfinity(f,g,O)
    
    [_,x,y] = PR(K_(f),'x','y')
    if evaluate(f,[x,y],O)!=0 or evaluate(g,[x,y],O)!=0: return 0
    a,b = O
    if a!=0 or b!=0:
        f = evaluate(f,[x,y],[x+a,y+b])
        g = evaluate(g,[x,y],[x+a,y+b])
    f = f + x*y - x*y
    g = g + x*y - x*y
    f0 = constant_coeff(f,y); r = degree(f0)
    g0 = constant_coeff(g,y); s = degree(g0)
    if f0==0:
        if g0==0:
            return 'Infinity'
        else:
            return trailing_degree(g0) + imult(f/y,g)
    else: # f0!=0
        if g0==0:
            return trailing_degree(f0)+imult(f,g/y)
        else: # g0!=0
            c0 = leading_coeff(f0)
            d0 = leading_coeff(g0)
            if r<=s:
                return imult(f,c0*g-d0*x**(s-r)*f)
            else:
                return imult(d0*f-c0*x**(r-s)*g,g)
    return 'error'

# If f is a polynomial, the following function homogenizes it
# with respect to a specified new variable 
def homogenize(f,z=''):
    if z=='': _,z=polynomial_ring(Q_,'z')
    n = total_degree(f)
    F = 0
    for k in range(n+1):
        h = hpart(f,k)
        if h!=0: F += h*z**(n-k)
    return collect(F,z)
    
## Intersetion multiplicity at an improper point
def imultinfinity(f,g,O):
    VF = variables(f)
    VG = variables(g)
    if (VF == VG):
        x = VF[0]
        y = VF[1]
    else:
        if len(VF) == 2:
            x = VF[0]
            y = VF[1]  
        else:
            x = VG[0]
            y = VG[1]
    a,b,_ = O
    _,z = polynomial_ring(Q_,'z')
    F = homogenize(f,z); G= homogenize(g,z)
    
    if b!=0:
        a = a/b; b = 1
        f = evaluate(F,[x,y,z],[x,1,y])
        g = evaluate(G,[x,y,z],[x,1,y])
        return imult(f,g,[a,0])
    a=1; b=0
    f = evaluate(F,[x,y,z],[1,x,y])
    g = evaluate(G,[x,y,z],[1,x,y])
    return imult(f,g,[0,0])




## Examples

# Fulton's example over Q_, the field of rational numbers,
# and also over the finite field F5. For the latter, comment
# the first line below and uncomment the second.
show('First example')
F = Q_     # in this case it yields 14
#F = Zn(5) # in this case it yields 18

[Fxy,x,y] = polynomial_ring(F,'x','y')

f = (x**2+y**2)**3 - 4*x**2*y**2
g = (x**2+y**2)**2 + 3*x**2*y - y**3

fi = -x**3 + y
gi = -x**3 + (1 - x**2)*y

show(imult(f,g))
show(resultant(f,g))

show('Second example')
# Intersection multiplicites  of a cuspidal and a nodal cubic

# f1=collect(y**2-x**3,y)
# g1=collect(y**2-x**2-x**3,y)
f1=-x**3+y**2
g1=-x**3-x**2+y**2

show(imult(f1,g1))

# The explanation of values is checked by rhe resultant
# of f and g, as explained in ITEG.
# 
show("Resultant R(f1,g1) =",resultant(f1,g1))

show('intersecton multiplicity at infinity')
# 
show(imult(f1,g1,[0,1,0]))

## Exercises 3.3 in Walker

show('Third Example')
# First pair (a)
af = x*(y**2-x)**2 - y**5
ag = -x**2 + y**4 + y**5

show(imult(af,ag))

show(imultinfinity(af,ag,[1,0,0]))

R = resultant(af,ag)/x**10

R = R/(x**2-4*x-1)

show(R)

show('Fourth Example')
bf = x**3 - 2*x*y - y**3
bg = -2*x**2 + 2*x**3 - 4*x**2*y - 3*x*y**2 - y**3

show(imult(bf,bg))

show(-resultant(bf,bg)/x**5)

show('Fifth example')
cf = x**4 + y**4 - y**2
cg = x**4 + y**4 - 2*y**3-2*x**2*y-x*y**2 + y**2

show(imult(cf,cg))

show(resultant(cf,cg)/x**10)

show(56**2-64*49)

show((7*x+4)**2)

show('Sixth example')
f = x**3+y**3-2*x*y 
g = 2*x**3-4*x**2*y+3*x*y**2+y**3-2*y**2

m0 = imult(f,g)
show('m0 =',m0)

m1 = imult(f,g,[1,1])
show('m1 =',m1)

R = resultant(f,g)
show('R =',R)
R0 = R/x**5
R1 = R0/(x-1)**3
show('R1 =',R1)
x2=-trailing_coeff(R1)/leading_coeff(R1)
show('x2 =',x2)  # 4/7

f2 = evaluate(f,[x],[x2])
g2 = evaluate(g,[x],[x2])


h2 = g2-f2
h2 = h2/leading_coeff(h2)
a,b,c = C = coeffs(h2)

D = b**2-4*a*c
SD = sqrt(D) >> Q_
y21,y22 = (-b + SD)/2>>Q_, (-b - SD)/2>>Q_
show(y21,y22)

show(evaluate(h2,[y],[y21]))
show(evaluate(h2,[y],[y22]))
show(evaluate(f,[x,y],[x2,y21]))
show(evaluate(g,[x,y],[x2,y21]))
show(imult(f,g, [x2,y21]))
show(imult(f,g, [x2,y22]))