# Author: Antonio Laface
# Date: May 5 2001
# E-mail: alaface@udec.cl
#-----------
#
# Here x is a vector with four entries, x=[d,m0,r,m], representing        
# polynomials of degree d trough a point of multiplicity m0 and r points
# of multiplicity r.
#
# virtual_dim(x) returns the virtual dimension of the vector x
#
# expected_dim(x) returns the expected dimension of the vector x
#
# system_product(x,y) returns the product of the strict transforms of x 
# and y over the blow up of P^2 along the r+1 points.
#
# system_difference(x,y,t) returns the difference of the strict transforms 
# of x+ty (t is a negative number).
#
# HC is the list of (-1)-curves, not compound, with multiplicity <= 2 on 
# P^2.
# HCN is the list of (-1)-curves, possibly compound, with multiplicity <= 
# 2 on P^2.
#
# effective_dim(x) return the effective dimension of the linear system z.
# This procedure works in a recursive way and is based on the conjecture 
# that a linear system is special if and only if its strict transform has
# product <= -2 with some -1 curve of the blow up.
#
# For a (k,b) degeneration we have the following:
#
# lp(x,k,b)  returns the dimension of the first system.
# lf(x,k,b)  returns the dimension of the second system.
# klp(x,k,b) returns the dimension of the first kernel-system.
# klf(x,k,b) returns the dimension of the second kernel-system.
# rp(x,k,b)  returns the dimension of the first restricted system.
# rf(x,k,b)  returns the dimension of the second restricted system.
# lo(x,k,b)  returns the dimension of the degenerated system.
#
# test(x) returns 1 if x is a (-1)-special system, return 0 if there exist
# a degeneration (k,b) such that lo(x)=effective_dim(x) (hence x is
# non-special) and returns -1 otherwise.
#
# analyze(deg,mult) analyze systems of degree <= deg and multiplicity = mult
# and prints out those systems ( x ) for which the degeneration does not work
# ( i.e. test(x) = -1 ).
#

virtual_dim := x -> binomial(x[1]+2,2)-binomial(x[2]+1,2)-x[3]*binomial(x[4]+1,2)-1:

expected_dim := x -> max(virtual_dim(x),-1):

system_product      := (x,y) -> x[1]*y[1]-x[2]*y[2]-min(x[3],y[3])*x[4]*y[4]:

system_difference   := (x,y,t) -> [x[1]+t*y[1],x[2]+t*y[2],x[3],x[4]+t*y[4]]:

HC  := [[2,0,5,1],[6,3,7,2],[e,e-1,2*e,1],[1,1,1,1],[1,0,2,1]]:
HCN := [[2,0,5,1],[6,3,7,2],[e,e-1,2*e,1],[t,t,t,1],[3,0,3,2]]:

effective_dim := proc(z) 

	local i, w, EC, ECN, ok; 
	global HC, HCN;  

	for i from 1 to 4 do 
		if z[i] < 0 then RETURN(-1):
		fi: 
	od;
	w := z:	ok := 0:
	EC := eval(HC,e=floor(z[3]/2)): 
	ECN := eval(HCN,[e=floor(z[3]/2),t=z[3]]):
  
	for i from 1 to 5 do 
		if system_product(z,EC[i]) <= -2 and z[3] = ECN[i,3] then
			w := system_difference(z,ECN[i],system_product(z,EC[i])): 
			i:=5: ok:=1:
		fi;
	od;

	if ok = 0 then 

		RETURN(expected_dim(w)):
	else 
		effective_dim(w):
	fi;
end: 

lp  := (x,k,b) -> effective_dim([x[1]-k,x[2],x[3]-b,x[4]]): 
lf  := (x,k,b) -> effective_dim([x[1],x[1]-k,b,x[4]]): 
klp := (x,k,b) -> effective_dim([x[1]-k-1,x[2],x[3]-b,x[4]]): 
klf := (x,k,b) -> effective_dim([x[1],x[1]-k+1,b,x[4]]): 
rp  := (x,k,b) -> lp(x,k,b) - klp(x,k,b) - 1:
rf  := (x,k,b) -> lf(x,k,b) - klf(x,k,b) - 1:
lo  := (x,k,b) ->
	if (rp(x,k,b) + rf(x,k,b) <= x[1] - k - 1) then
	
		klp(x,k,b) + klf(x,k,b) + 1: 
	else
		lp(x,k,b) + lf(x,k,b) - x[1] + k:
	fi:

test := proc(x) local b, k;  

	if effective_dim(x) > expected_dim(x) then
		RETURN(1):
	fi:
 
	for k from x[4]-1 to x[4] do
		for b from 1 to x[3]-1 do
			if lo(x,k,b) = effective_dim(x) then
				RETURN(0):
			fi:
		od:
	od:

	RETURN(-1):
end:


analyze := proc(deg,mult) local d,m,r_min,r_max; 

  for d from 8 to deg do

        for m from 0 to d-8 do

                r_max := ceil((binomial(d+2,2)-binomial(m+1,2))/binomial(mult+1,2)):
                r_min := r_max - 1: 
                
                while( test([d,m,r_max,mult]) = 1 ) do
                        r_max := r_max + 1:
                od:

                while( test([d,m,r_min,mult]) = 1 ) do
                        r_min := r_min - 1:
                od:

                if ( test([d,m,r_max,mult]) = -1 ) then
                        printf("L(%d,%d,%d^%d)\n",d,m,mult,r_max):

                fi:

                if ( test([d,m,r_min,mult]) = -1 ) then
                        printf("L(%d,%d,%d^%d)\n",d,m,mult,r_min):
                fi:
        od:
  od:

end: