A polynomial generalization of the Euler characteristic for algebraic sets.

135 days ago by mmarco

This worksheet contains the code for computing the Euler characteristic and its polynomial generalization as described in the article of the same title.

 

This first function computes the Euler characteristic. The entry is an ideal on a polynomial ring.
def euler_characteristic(I): R=I.ring() if I.is_one(): return 0 if R.ngens()==1: return sum([j[0].degree() for j in I.gen().factor()]) J1=I.radical() if J1.is_homogeneous(): return 1 primdec=J1.associated_primes() J1=primdec[0] m=len(primdec) if m>1: J2=R.ideal(1) for j in [1..m-1]: J2=J2.intersection(primdec[j]) return euler_characteristic(J1)+euler_characteristic(J2)-euler_characteristic(J1+J2) P=J1.homogenize().hilbert_polynomial() if P.is_zero(): deg=0 else: deg=P.leading_coefficient()*P.degree().factorial() if deg==1: return 1 dim=J1.dimension() n=R.ngens() vars1=R.gens()[0:n-dim] vars2=R.gens()[n-dim:n] varpiv=vars1[-1] IH=J1.homogenize() S=IH.ring() JH=IH+S.ideal(S.gens()[n-dim:]) if JH.dimension()>0: det=0 while det==0: MH=random_matrix(R.base_ring(),n) det=MH.determinant() L=list(MH*vector(list(R.gens()))) return euler_characteristic(R.hom(L)(J1)) if dim==0: return deg M=matrix([[f.derivative(v) for v in vars1] for f in J1.gens()]) J=R.ideal(M.minors(n-dim)) K=(J+J1).elimination_ideal(vars1) S=PolynomialRing(R.base_ring(),vars2) H=R.hom([S(0) for j in vars1]+[S(j) for j in vars2]) C=deg-deg*euler_characteristic(H(K))+euler_characteristic(K+J1) return C 
       

This second function computes the polynomial F(V).
@parallel(7) @cached_function def FV(I,var='L'): FS=PolynomialRing(ZZ,var) L=FS.gen() R=I.ring() if R.ngens()==0: return 0 if I.is_zero(): return L^R.ngens() if I.is_one(): return 0 if R.ngens()==1: return FS(sum([j[0].degree() for j in I.gen().factor()])) J1=I.radical() if J1.is_homogeneous(): S1=PolynomialRing(R.base_ring(),R.gens()[0:-1]) H1=R.hom(list(S1.gens())+[S1(1)]) H2=R.hom(list(S1.gens())+[S1(0)]) resulp=FV([H1(J1),H2(J1)]) d=dict([[a[0][0][0],a[1]] for a in resulp]) [i1,i2]=[d[H1(J1)],d[H2(J1)]] return (L-1)*(i1)+i2 primdec=J1.associated_primes() J1=primdec[0] m=len(primdec) if m>1: J2=R.ideal(1) for j in [1..m-1]: J2=J2.intersection(primdec[j]) resulp=FV([J1,J2,J1+J2]) d=dict([[a[0][0][0],a[1]] for a in resulp]) [i1,i2,i3]=[d[J1],d[J2],d[J1+J2]] return i1+i2-i3 P=J1.homogenize().hilbert_polynomial() if P.is_zero(): deg=0 else: deg=P.leading_coefficient()*P.degree().factorial() dim=J1.dimension() if deg==1: return FS(L^dim) n=R.ngens() vars1=R.gens()[0:n-dim] vars2=R.gens()[n-dim:n] varpiv=vars1[-1] IH=J1.homogenize() S=IH.ring() JH=IH+S.ideal(S.gens()[n-dim:]) if JH.dimension()>0: det=0 while det==0: MH=random_matrix(R.base_ring(),n) det=MH.determinant() L=list(MH*vector(list(R.gens()))) return FV(R.hom(L)(J1)) if dim==0: return FS(deg) M=matrix([[f.derivative(v) for v in vars1] for f in J1.gens()]) J=R.ideal(M.minors(n-dim)) K=(J+J1).elimination_ideal(vars1) S=PolynomialRing(R.base_ring(),vars2) H=R.hom([S(0) for j in vars1]+[S(j) for j in vars2]) d=dict([[a[0][0][0],a[1]] for a in FV([H(K),K+J1])]) [i1,i2]=[d[H(K)],d[K+J1]] C=deg*FS(L^dim-i1)+i2 return C 
       

Now we show some examples to ilustrate the functions. Note that the actual timings can depend on the configuration and the workload of the server at every given moment.

Three examples of plane curves:
R.<x,y>=QQ[] time euler_characteristic(R.ideal(x^5+1)) 
        
5
Time: CPU 0.11 s, Wall: 0.12 s
5
Time: CPU 0.11 s, Wall: 0.12 s
time euler_characteristic(R.ideal(y^4+x^3-1)) 
        
-5
Time: CPU 0.68 s, Wall: 0.68 s
-5
Time: CPU 0.68 s, Wall: 0.68 s
time euler_characteristic(R.ideal(x^2+y^2+5*x^2*y^4+x*y-1)) 
        
-8
Time: CPU 35.07 s, Wall: 35.08 s
-8
Time: CPU 35.07 s, Wall: 35.08 s

A curve and a surface in \mathbb{C}^3:
S.<x,y,z>=QQ[] time euler_characteristic(S.ideal(x^5+y^2+2*x*y+1,3*x-5*y*x+y^2+1)) 
        
10
Time: CPU 0.14 s, Wall: 0.14 s
10
Time: CPU 0.14 s, Wall: 0.14 s
time euler_characteristic(S.ideal(x^5+y^2+2*x*y+1)) 
        
-3
Time: CPU 0.37 s, Wall: 0.37 s
-3
Time: CPU 0.37 s, Wall: 0.37 s

Now we show some examples of the computation of F(V).

The complex 2-sphere:
R.<x,y,z>=QQ[] time FV(R.ideal(x^2+y^2+z^2-1)) 
        
2*L^2 - 2*L + 2
Time: CPU 0.06 s, Wall: 0.45 s
2*L^2 - 2*L + 2
Time: CPU 0.06 s, Wall: 0.45 s

Another surface:
time FV(R.ideal(x^3+y^3+z^3-1)) 
        
3*L^2 - 6*L + 12
Time: CPU 0.07 s, Wall: 0.38 s
3*L^2 - 6*L + 12
Time: CPU 0.07 s, Wall: 0.38 s

The intersection and the union of the two:
time FV(R.ideal(x^3+y^3+z^3-1,x^2+y^2+z^2-1)) 
        
6*L - 15
Time: CPU 0.08 s, Wall: 62.31 s
6*L - 15
Time: CPU 0.08 s, Wall: 62.31 s
FV(R.ideal((x^3+y^3+z^3-1)*(x^2+y^2+z^2-1))) 
        
5*L^2 - 14*L + 29
5*L^2 - 14*L + 29

The Whitney umbrella:
time FV(R.ideal(x*y^2-z^2)) 
        
3*L^2 - 4*L + 2
Time: CPU 0.19 s, Wall: 2.61 s
3*L^2 - 4*L + 2
Time: CPU 0.19 s, Wall: 2.61 s

The complex 3-sphere:
S.<x,y,z,t>=QQ[] time FV(S.ideal(x^2+y^2+z^2+t^2-1)) 
        
2*L^3 - 2*L^2 + 2*L - 2
Time: CPU 0.08 s, Wall: 0.65 s
2*L^3 - 2*L^2 + 2*L - 2
Time: CPU 0.08 s, Wall: 0.65 s