# $Source: /u/maple/research/lib/LREtools/src/RCS/ratbeke,v $
# $Notify: mvanhoei@daisy.uwaterloo.ca

macro(	ratbeke=`LREtools/ratbeke`,
	doit=`LREtools/ratbeke/doit`,
	find_rf=`LREtools/ratbeke/find_rf`,
	rf_vals=`LREtools/ratbeke/rf_vals`,
	cnd_count=`LREtools/ratbeke/cnd_count`,
	expansion=`LREtools/ratbeke/expansion`,
	subfields=`LREtools/ratbeke/subfields`,
	subf=`LREtools/ratbeke/subf`,

	g_p=`LREtools/g_p`,
	normal_poly=`LREtools/g_p/normal_p`,
	singularities=`LREtools/g_p/sing`,
	symprod1=`LREtools/g_p/symprod1`,
	g_factors=`LREtools/g_p/factors`,
	symprod1=`LREtools/g_p/symprod1`,
	p_curv=`LREtools/p_curv`,

	LOW=-3,
	MED=-2,
	HIG=-1
):

ratbeke:=proc(L1,x::name,tau::name,ext)
	local L,i,n,cnd,sing,essn,semi,gp,cnds,F,j,B,n_found;
	options trace;

	if nargs<4 or not type(ext,set) then
		return procname(L1,x,tau,
		 indets([args],{'RootOf','radical','nonreal'}))
	fi;
	Normalizer:=`if`(ext={},normal,`evala/Normal`);
	_Env_LRE_x:=x;
	# evala@Primpart can be slow, so we take primpart
	L:=collect(primpart(L1,tau,'co'),tau,normal_poly);
	i, n:=ldegree(L,tau), degree(L,tau);
	if L=0 then
		error `input must be non-zero`
	elif not has(L,tau) then
		return []
	elif i<>0 then
		userinfo(1,'LREtools',L,`has factor`,tau^i);
		return procname(collect(subs(x=x+i,L/tau^i)
		 ,tau,normal_poly),args[2..-1])
	fi;
	# Could catch factors in C[tau] by: gcd(L,subs(x=0,L)).
	cnd:=g_p(L,x,tau,infinity,ext);
	if cnd=NULL then
	userinfo(2,'LREtools',
	 `Point at infinity shows non-existence of hypergeometric solutions`);
	  return []
	fi;
	cnd:=[seq([op(i),nops(i[4])],i=[cnd])];
	cnd:=cnd_count(L,x,tau,2,cnd);
	if cnd<>[] then cnd:=cnd_count(L,x,tau,3,cnd) fi;
	if max(seq(i[5],i=cnd))>1 then cnd:=cnd_count(L,x,tau,5,cnd) fi;
	if cnd=[] then	
	userinfo(2,'LREtools',
	 `p-curvature shows non-existence of hypergeometric solutions`);
	  return []
	fi;
	sing:=singularities(L1,x,tau,ext,L,co);
	essn:=NULL;
	semi:=NULL;
	for i in sing do
		gp[i]:=g_p(L,x,tau,RootOf(i,x));
		gp[i]:=[gp[i], i,
			# max number with low growth rate:
			gp[i][2],
			# max number with medium growth rate:
			`if`(nops(g[i][1])=2,0,n-max(gp[i][2],gp[i][3])),
			# max number with high growth rate:
			gp[i][3] ];
		if nops(gp[i][1])>1 then
			# i is an essential singularity
			essn:=essn,i
		else
			# i is a (semi)-regular singularity
			semi:=semi,i
		fi
	od;
	essn:=`polytools/sort_poly`([essn],x);
	semi:=[semi];
	# Construct base fields + degree possible extension and cnd's:
	cnds:=NULL;
	for i in cnd do
		F:=indets(i,'RootOf') union ext;
		cnds:=cnds,seq([op(i),F,`algcurves/degree_ext`(F,ext),j]
		 ,j=1..i[5])
	od;
	cnds:=sort([cnds], proc(a,b)
		local i,j;
		i,j:=a[-1],b[-1];
		if i=j then
			i,j:=a[-2],b[-2];
			if i=j then
				i,j:=seq(length(i),i=[a,b])
			fi
		fi;
		evalb(i<j) end);
	doit(gp,L,x,tau,essn,semi,ext,n,cnds,{seq(i,i=1..n)},n_found,false)
end:


doit:=proc(gp,L,x,tau,essn,semi,ext,n,cnds,d_set,n_found,adj)
	local nf,tried_adj,res,cnd,F,d,essn_F,B,Fs,es,i,v,j,
	 sing,Sde,Str,gmin,gmax,de,tr,g,curr,r,f,tff,k,kk;

	res:=NULL;
  for cnd in cnds do
	if assigned(n_found[cnd[1..3]]) then
		nf:=n_found[cnd[1..3]]
	else
		nf:=0
	fi;
	F:=cnd[-3]; # extension field, but need another extension of
		    # degree cnd[-1] on top of that.
	d:=cnd[-1]*cnd[-2]; # = total degree of extension
	if cnd[-1]>cnd[-4]-nf or not member(d,d_set) then
		next
	elif d=1 then
		essn_F:=[essn,F] # essn factored over F, and F
	elif d=n and member(1,{seq(igcd(degree(i,x),n),i=essn)}) then
		# All elements of essn should factor over the field,
		# and they don't, so:
		next
	elif cnd[-1]=1 then
		es:=NULL;
		for i in essn do
			# factor essn over F:
			v:=map(j->j[1],g_factors(i,x,F));
			if nops(v)=1 and d=n then
				v:=NULL;
				break
			fi;
			es:=es,op(v);
			for j in v do gp[j]:=gp[i] od
		od;
		if v=NULL then next fi;
		essn_F:=[[es], F]
	elif adj then
		return [res]
	# now we need to construct an extension of degree cnd[-1]
	elif igcd(mul(degree(i,x),i=essn),d)=1 then
		# either some of the essn factor over F or some of
		# them contribute to the degree cnd[-1] extension,
		# if not, then:
		next
	elif cnd[-1]=n and isprime(n) then
		v:={seq(F union {i},i=`evala/Subfields`(
		 {op(essn)},n,ext,x))};
		for F in v do
			i:=procname(op(subsop(9=
			 [subsop(-1=1,-2=n,-3=F,cnd)],[args])));
			if i<>[] then return i fi
		od;
		next
	else
		# Now construct subfields
		v:=essn;
		if cnd[-2]>1 then
			# these have already been factored before
			# but this shouldn't take extra time due to
			# options remember
			v:=`polytools/sort_poly`(map(i->i[1],
			 [seq(op(g_factors(i,x,F)),i=v)]),x);
		fi;
		if igcd(mul(degree(i,x),i=v),cnd[-1])=1 then
			# can't construct degree cnd[-1] extension
			next
		fi;
		r:=subfields(cnd[-1],essn,x,F);
# Somehow we have to throw away the ones that were actually defined
# over a smaller field (that were already found).
# Ook n_found lijkt mij twijfelachtig allemaal want je kunt ze dubbel
# tellen als je niet uitkijkt doordat sommige lichamen ook sommige
# geconjugeerden bevatten en die vind je dan vanzelf ook.
		for F in r do
			res:=res,op(procname(op(subsop(9=
			 [subsop(-1=1,-2=d,-3=F,cnd)],[args]))))
		od;
		next
	fi;

    for Fs in [essn_F] while cnd[-1] <= cnd[-4]-nf do
	sing, F:= [op(Fs[1]), op(semi)], Fs[2];
	# Now sing is the set of singularities and F is the field
	if F<>{} then Normalizer:=`evala/Normal` fi;
	Sde:=cnd[2];
	Str:=cnd[3];
	for i to nops(sing) do
		j:=gp[sing[i]];
		gmin[i]:=min(op(j[1]));
		gmax[i]:=max(op(j[1]));
		r:=`algcurves/degree_ext`(F, [ext,sing[i]]);
		if r>max(j[LOW],j[HIG],`if`(gmax[i]-gmin[i]=1,NULL,j[MED])) then
			# degree field extension is too high
			F:=0; break
		fi;
		if r>j[LOW] then
			# Valuation grow can't be minimal because
			# we need at least r hypergeom solutions
			# of this local type but we can have only
			# j[2]=rank(E_{p,r}) of minimal type
			gmin[i]:=`if`(r>j[MED],gmax[i],gmin[i]+1)
		fi;
		if r>j[HIG] then
			gmax[i]:=`if`(r>j[MED],gmin[i],gmax[i]-1)
		fi;
		de[i]:=degree(sing[i],x);
		tr[i]:=coeff(sing[i],x,de[i]-1); # i should be monic
		g[i]:=gmin[i];
		Sde:=Sde+g[i]*de[i];
		Str:=Str+g[i]*tr[i]
	od;
	if F=0 then next fi;
	while cnd[-1] <= cnd[-4]-nf do
		curr:=1;
		if Sde=0 then
			Str:=Normalizer(Str);
			if type(Str,integer) then
				r:=[seq(find_rf(gp[sing[i]],g[i],x,n,sing[i])
				 ,i=1..nops(sing))];
				# tff=subs(x=x+1,f)/f
				r,f,tff:=cnd[1]*mul(i[1],i=r),
					mul(i[2],i=r), mul(i[3],i=r);
				userinfo(6,'LREtools',`trying type`,r);
				j:=symprod1(L,r*tff,x,tau);
lprint(j);
				v:=LREtools[polysols](add(
				 coeff(j,tau,i)*A(x+i),i=0..n)
				,A(x),{},'output'='basis');
				if v<>0 and v<>NULL then
					userinfo(3,'LREtools',`found a solution`);
					j:=indets(v,'name') minus indets([j,A,N]);

					# now adjust the gp's
					for i to nops(Fs[1]) do
					  kk:=nops(j);       #    *`algcurves/degree_ext`(F,[ext,sing[i]])
					                     # *degree(sing[i],x)/degree(gp[sing[i]][-4],x);
					  for k in map(op,[indices(gp)]) do if gp[sing[i]][-4]=gp[k][-4] then
						if g[i]=min(op(gp[k][1])) then
							gp[k]:=subsop(LOW=gp[k][LOW] - kk ,gp[k])
						elif g[i]=max(op(gp[k][1])) then
							gp[k]:=subsop(HIG=gp[k][HIG] - kk ,gp[k])
						else
							gp[k]:=subsop(MED=gp[k][MED] - kk ,gp[k])
						fi
					fi od od;

					v:=[seq(Normalizer(f*coeff(v,i)),i=j)];
					nf:=nf+nops(v)*cnd[-1];
					res:=res,[r,v]
				fi
			fi
		fi;
		if nops(sing)=0 then
			break
		fi;
		while g[curr]=gmax[curr] do
			if curr=nops(sing) then
				curr:=0;
				break
			fi;
			g[curr]:=gmin[curr];
			Sde:=Sde+(gmin[curr]-gmax[curr])*de[curr];
			Str:=Str+(gmin[curr]-gmax[curr])*tr[curr];
			curr:=curr+1
		od;
		if curr=0 then
			break
		fi;
		g[curr]:=g[curr]+1;
		Sde:=Sde+de[curr];
		Str:=Str+tr[curr]		
	od
    od;
  n_found[cnd[1..3]]:=nf;
  od;
  [res]
end:

# Input: gp = g_p(L,x,tau,p) where p is a finite singularity
# 	 g is an element of gp[1]
#	 n is the order
# The hypergeometric solution we are looking for is
# a solution of: tau - r*tau(f)/f
# times some polynomial.
find_rf:=proc(gp::list,g::integer,x::name,n::integer,sing)
	local r,f,i,tff,so,gpn;
	gpn:=gp[6]-n-1;
	r,f:=op(rf_vals(gp[4],gp[5],g));
	r:=mul(normal_poly(subs(x=x-i-gpn,sing)),i=r);
	tff:=[seq(`if`(i=nops(f),0,f[i+1])-f[i],i=1..nops(f))];
	so:=proc(f,x,gpn,sing)
	  local i;
	  mul(`if`(f[i]<>0,normal_poly(subs(x=x-i-gpn,sing))^f[i],1),i=1..nops(f))
	end:
	f:=so(f,x,gpn,sing);
	tff:=so(tff,x,gpn,sing);
	[`if`(g<0,1/r,r),f,tff]
end:
#
# Here we will calculate the information necessary
# to construct r and f for find_rf.
rf_vals:=proc(gp4::list,gp5::list,g::integer)
	local v,i,ng,j;
	ng:=nops(gp4);
    if g=0 then
	# [r,f]
	[[],[seq(max(gp4[i],gp5[i]),i=1..ng)]]
    else
	if g<0 then
		j:=min(seq(`if`(gp4[i]<0,i,NULL),i=1..ng));
		v:=procname(
		  [seq(gp4[i]+`if`(i<j,0,1),i=1..ng)],
		  [seq(gp5[i]-`if`(i<j,1,0),i=1..ng)],
		  g+1)
	else
		j:=max(seq(`if`(gp5[i]<0,i,NULL),i=1..ng))+1;
		v:=procname(
		  [seq(gp4[i]-`if`(i<j,0,1),i=1..ng)],
		  [seq(gp5[i]+`if`(i<j,1,0),i=1..ng)],
		  g-1)
	fi;
	[[op(v[1]),j-1],v[2]]
    fi
end:


# F is the p_curvature polynomial in x and tau.
# R is the c*x^n*(1+d*(1/x))
# We want the number of roots with c*x^n+lower terms
expansion:=proc(F,x,tau,R,p,deg_R,co)
	local l,n;
	l:=lcoeff(F,tau);
	if l<>1 then
		n:=degree(F,tau);
		procname(Expand(tau^n+
		 add(coeff(F,tau,i)*l^(n-1-i)*tau^i,i=0..n-1)) mod p
		 ,x,tau,Expand(R*l) mod p,p)
	elif nargs=5 then
		procname(args,max(0,degree(R,x)))
	elif nargs=6 then
		add(procname(F,x,tau,R+l*x^(deg_R-1) mod p
		 ,p,deg_R,0),l=1..p)
	elif ldegree(R,x)<0 then
		0
	elif co=0 then
		procname(Expand(subs(tau=tau+R,F)) mod p
		 ,x,tau,0,p,deg_R-2,1)
	elif deg_R<0 then
		ldegree(F,tau)
	else
		l:=Expand(subs(tau=n*x^deg_R,F)) mod p;
		l:=subs(n=x,lcoeff(l,x));
		if not has(l,x) then
			return 0
		fi;
		n:=Factors(l,op(indets(F,RootOf))) mod p;
		add(degree(l[1],x)*procname(Expand(`mod/ReduceField`(
		 subs(tau=tau+RootOf(l[1],x)*x^deg_R,F),{x,tau},p)) mod p
		,x,tau,0,p,deg_R-1,1),l=n[2])
	fi
end:


# cnd[-1] = max number of hypergeomsols with that local type.
cnd_count:=proc(L,x,tau,p,cnds)
	local v,w,cnd,delta,pc,res,i,f;
	v:=[L,op(indets(L,{'RootOf','radical','nonreal','name'}))];
	# Could generate an error:
	w:=traperror(`mod/ReduceField`(v,{x,tau},p));
	if w=lasterror then
		return cnds
	fi;
	delta:=degree(L,tau)-degree(w[1],tau);
	if delta >= max(seq(cnd[-1],cnd=cnds)) then
		return cnds
	fi;
	pc:=`LREtools/p_curv`(w[1],x,tau,p,`dont factor`);
	res:=NULL;
	for cnd in cnds do
		if delta >= cnd[-1] then
			res:=res,cnd;
			next
		fi;
		f:=cnd[1]*expand(x^(-cnd[2])*(1+cnd[3]/x));
		f:=traperror(`mod/ReduceField`([subs(0=0,
		 seq(v[i]=w[i],i=2..nops(v)),f),pc],{x,tau},p));
		if f=lasterror then
			res:=res,cnd;
			next
		fi;
		i, f := op(f);
		i:=expansion(f,x,tau,i,p);
		if i=0 then next fi;
		res:=res, subsop(-1=min(cnd[-1],i),cnd)
	od;
	[res]
end:



# Give all subfields for f of degree n
subf:=proc(n,f,x,ext)
	local i;
	options remember;
	if n=1 then
		{ext}
	else
		{seq({i} union ext,i=`evala/Subfields`(f,n,ext,x))}
	fi
end:

# Input polynomials should be irreducible over Q(ext).
# Compute all subfields of degree n for a list f of polynomials,
# i.e. compute field extensions over which some of the f factor
# but over which no_factor do not factor.
subfields:=proc(n,f::list,x,ext,no_factor)
	local i,j,k,v;
	option remember;
	if n=1 then
		{ext}
	elif nargs=4 then
		procname(args,{})
	elif member(f[1],no_factor) then
		{}
	else
	  if nops(f)=1 then
		v:=subf(n,f[1],x,ext)
	  else
		v:=[op(f[2..-1]),f[1]];
		v:=`union`(seq(seq(procname(n/i,`if`(j=ext,v,
		 map(k -> k[1],map(op,map(g_factors,v,x,j))))
		,x,j,{op(no_factor),`if`(i=1,f[1],NULL)})
		 ,j=subf(i,f[1],x,ext)), i=numtheory['divisors'](n)))
	  fi;
		{seq(`if`({1,op(map(nops,map(g_factors,no_factor,x,i)))}
		 ={1},i,NULL),i=v)}
	fi
end:

#savelib('ratbeke','doit','find_rf','rf_vals', \
 'expansion','cnd_count','subf','subfields'):