# $Source: /u/maple/research/lib/algcurves/src/RCS/monodromy,v $
# $Notify: deconinc@newton.colorado.edu $
# $Notify: mvanhoei@daisy.uwaterloo.ca $

# Maple V release 5 or 6?
R6:=evalb(traperror(type(0,nonreal))<>lasterror):

macro(
	display=plots['display'],
	puiseux=algcurves['puiseux'],
	convert_disjcyc = `group/listtocyc`,
	convert_permlist = `group/cyctolist`,
	mulperms=group['mulperms'],
	invperm=group['invperm'],
	EXT = {`if`(R6,nonreal,NULL), radical, algext},

	EXTRA1 = 3,
	EXTRA2 = 2,
	CLOSE = 4,
	CLOSE2 = 20,
	Increase_Digits = 1/16,
	Initial_tstep = 1/8,
	very_close=`algcurves/v_close`,

	Hurwitzsystem=`algcurves/monodromy`,
	Monodromy=`algcurves/Monodromy`,
	Linearmonodromy=`algcurves/Lmonodromy`,
	Generalmonodromy=`algcurves/Gmonodromy`,
	check_cstruct=`algcurves/check_cstruct`,
	Showpaths=`algcurves/Showpaths`,
	Permuted=`algcurves/Permuted`,
	Propagate=`algcurves/Propagate`,
	Analyticcontinuation=`algcurves/Acontinuation`,
	Continue=`algcurves/Continue`,
	Der=`algcurves/Der`,
	Distsort=`algcurves/Distsort`,
	Makepaths=`algcurves/Makepaths`,
	Sortnumbers=`algcurves/Sortnumbers`,
	Argsort=`algcurves/Argsort`
):


`algcurves/fsolve`:=proc(f,x)
local d,l,D,L,r,i;
	option `Copyright (c) 1999 Waterloo Maple Inc. All rights reserved.`;
	d,l:=degree(f,x),ldegree(f,x);
	D,L:=coeff(f,x,d),coeff(f,x,l);
	if d-l>2 and {seq(coeff(f,x,i),i=l+1..d-1)}={0} and (Im(D)<>0 or Im(L)<>0) then
		r:=evalf( (-L/D)^(1/(d-l)) );
		r,seq(evalf(r*exp(2*Pi*I*i/(d-l))),i=1..d-l-1),`if`(l>0,0,NULL)
	else
		fsolve(args)
	fi
end:

macro( fsolve=`algcurves/fsolve` ):

# This procedure computes the monodromy of a given plane Riemann surface. 
# The input of the procedure is the plane representation of the Riemann 
# surface and the names of the variables. Hence, "curve" is a polynomial in 
# two variables "x" and "y", such that the equation of the algebraic curve
# is curve=0.
# Considering the surface as a coveringsurface for y as a multivalued 
# function of x, the output of the procedure is a list containing the 
# following:
# 1) a non-special x value which does not correspond to a branch- or singular  
# point.
# 2) a list of pre-images of this x value, i.e. a list of y values,
# corresponding to the given x value. This effectively labels the sheets 
# of the Riemann surface for y(x).
# 3) A list of branchpoints with permutations giving how the sheets 
# interchange when one walks around the branchpoints. 
# The branchpoints are considered in order of
# increasing argument with respect to the non-special x value. 
# If the arguments coincide, branchpoints that are closer to the non-special
# x value are considered first.  
# Please note that the permutation of the sheets around infinity, if it is a
# branchpoint, is not given. It can easily be obtained as the inverse of the 
# composition of all the permutations around the finite branchpoints. Here 
# the composition must be done in the same order that they are given, i.e. in
# order of increasing argument.
# If as a fourth argument one gives `showpaths`, then maple graphs the paths
# that were used for the analytic continuation around the discriminant points.
# This option only results in a graph if the surface has more than two sheets, 
# because otherwise, no paths are constructed to find the monodromy. 

Monodromy:=proc(DI,F,x,y)
	options remember,
	    `Copyright (c) 1999 Waterloo Maple Inc. All rights reserved.`;
	if degree(F,y)=0 then
		ERROR(`curve is independent of the second variable`)
	elif degree(F,x)=0 and degree(F,y)<>1 then
		ERROR(`curve is independent of the first variable`)
	elif degree(F,y)=1 then
		Linearmonodromy(F,x,y)
	else if degree(F,x)<3 then
		userinfo(1,'algcurves',
		`Consider switching the roles of the variables. The curve has degree <=2 in the first variable.`);
	     fi;
		Generalmonodromy(args[2..nargs])
	fi
end:


# Given the monodromy of a Riemann surface, this procedure 
# computes a Hurwitz system for the surface. This is different 
# from the monodromy in that also the permutation around 
# infinity is included if infinity is a branchpoint. Also, 
# the permutations are rewritten in disjoint cycle form 
# (fixed points included).

Hurwitzsystem:=proc(curve,x::name,y::name)
	local F,sigma,C,i,b,mono;
	global `algcurves/A1`;
	option `Copyright (c) 1999 Waterloo Maple Inc. All rights reserved.`;
	description "Written by B. Deconinck with modifications by M. van Hoeij";
	if indets(curve,'name') minus {x,y}<>{} then
		ERROR(`Only 2 variables allowed`)
	elif indets(curve,float)<>{} then
		ERROR(`No floating point coefficients allowed`)
	elif member(group,[args[4..-1]]) then
		C:=procname(args[1..3],op(subs(group=NULL,[args[4..-1]])));
		'permgroup'(nops(C[2]),{seq(i[2],i=C[3])})
	elif args[nargs] <> `give paths` then
		C:=procname(args,`give paths`);
		[C[1],C[2],C[3]]
       else
	F:=collect(primpart(curve,{x,y}),{x,y},`distributed`,Normalizer);
	# at least 10 digits.
	Digits:=max(10,Digits);
	`algcurves/A1`:=[args[4..nargs],_Env_algcurves_opt];
	mono:=Monodromy(Digits,F,x,y);
	C:=[seq(mono[3][i][2],i=1..nops(mono[3]))];
	C:=map(convert_disjcyc,C);
	sigma:=[];
	for i from 1 to nops(C) do
		sigma:=mulperms(sigma,C[i]);
	od;
	check_cstruct(sigma,puiseux(F,x=infinity,y,`cycle structure`));
	b:=[seq(mono[3][i],i=1..nops(mono[3]))];
	if sigma<>[] then
		sigma:=convert_permlist(invperm(sigma),nops(mono[3][1][2]));
		b:=[op(b),[infinity,sigma]]
	fi;
	[mono[1],mono[2],[seq([i[1],convert_disjcyc(i[2])],i=b)],op(mono[4..-1])]
       fi
end:


# This procedure gives the monodromy information of the Riemann surface
# of a single-valued function.

Linearmonodromy:=proc(curve,x,y)
	local xvalue,preimages;
	option `Copyright (c) 1999 Waterloo Maple Inc. All rights reserved.`;
	xvalue:=0;
	while subs(x=xvalue,lcoeff(curve,y))=0 do
		xvalue:=xvalue-1
	od;
	preimages:=[fsolve(subs(x=xvalue,curve),y,complex)];
	[xvalue,preimages,[]]
end:

# This procedure gives the monodromy information in the general case

Generalmonodromy:=proc(curve,x,y)
	local discpoints,n,preimages,xvalue,paths,r,i,j,k,t,
	      orderedpt,permpoint,permutations,cyc_struct,
	      bpoints,dig,Pp,Q,rootof;
	global `algcurves/A1`;

	# covering number
	option `Copyright (c) 1999 Waterloo Maple Inc. All rights reserved.`;
	dig:=Digits;
    Digits:=dig+2;
	n:=degree(curve,y);
	# compute the discriminantpoints
	discpoints:=NULL;
	for i in evala(Factors(lcoeff(curve,y)^2*discrim(curve,y),indets(curve,EXT)))[2] do
		if i[2]>1 then
			r:=puiseux(curve,x=RootOf(i[1],x),y,`cycle structure`)
		else
			r:=[1$(n-2),2]
		fi;
		k:=fsolve(i[1],x,complex);
		# happens with fsolve(x^2+1,x,complex):
		if `algcurves/realcurve`(curve,x,y) and {k}<>map(conjugate,{k}) then
			k:=seq(op(`if`(Im(j)>0,[],[j,conjugate(j)])),j=k)
		fi;
		discpoints:=discpoints,k;
		for j in [k] do
			rootof[j]:=RootOf(i[1],x);
			cyc_struct[j]:=r
		od
	od;
	discpoints:=Sortnumbers([discpoints]);
	# radius of the circles used to walk around the discriminant points
	r:=table();
	for j in discpoints do
		r[j]:=min(seq(abs(i-j)/2.5,i={op(discpoints)} minus {j}));
	od;
	if nops(discpoints)=1 then r[discpoints[1]]:=1;
	fi;

if Im(discpoints[1])<>0 and `algcurves/realcurve`(curve,x,y) then
	# Make sure basepoint is real by adding a fake branchpoint.
	discpoints:=[evalf(Re(discpoints[1]-r[discpoints[1]])-1/10),op(discpoints)];
	r[discpoints[1]]:=-1/1000;
	cyc_struct[discpoints[1]]:=[1$n]
fi;

	# Choose the nonspecial x value to the left of all the branchpoints
	xvalue:=discpoints[1]-r[discpoints[1]];
	preimages:=[fsolve(subs(x=xvalue,curve),y,complex)];
	# make the paths along which the analytic continuation needs to be done
	paths:=Makepaths(t,discpoints,xvalue,r);
	if member('showpaths',`algcurves/A1`) then
		print(plot(Showpaths(paths,r,t),labels=[`Re`,`Im`]
		 ,scaling=CONSTRAINED,color=red,title=cat(
		 "Paths chosen for the analytic continuation for "
		 ,convert(y(x),string)),axes=NORMAL))
	fi;
	# Used for debugging.
	# if member('showindividualpaths',`algcurves/A1`) then
	#	print(plot(Showpaths(paths,r,t),labels=[`Re`,`Im`],scaling=CONSTRAINED,color=red,
	#	      title=cat(convert(curve,string),". Paths chosen for the analytic continuation")
	#	      ,axes=NORMAL));
	#	for j in discpoints do
	#	  if paths[j]<>[] then userinfo(1,'algcurves',`path for the analytic continuation to discriminant point`,j);
	#	  	pathplot:=plot({seq([Re(i),Im(i),t=0..1],i=paths[j])},scaling=CONSTRAINED,axes=NORMAL,color=blue):
	#	        datalist:=[seq([Re(i),Im(i)],i=discpoints)];
	#	        pointplot:=plot(datalist,style=POINT,color=red):
	#	        lineplot:=plot([[Re(xvalue),Im(xvalue)],[Re(j),Im(j)]],style=LINE,color=green):
	#	        print(display([pointplot,pathplot,lineplot]));         
	#	  fi;
	#	od;
	#fi;
	
	# do the analytic continuation over the semi-circles
	userinfo(1,'algcurves',`Starting the analytic continuation...`);
	# compute the permutations around the discriminant points,
	# relating the sheets to the indexing of the sheets above
	# xvalue.
	permutations:=table();
	bpoints:=Argsort([seq(`if`(cyc_struct[i]=[1$n],NULL,i)
	 ,i=discpoints)],xvalue);
	k:=0;
    Digits:=dig;
	for i in bpoints do
		k:=k+1;
	    if n=2 then
		permutations[i]:=[2,1]
	    elif Im(i)>0 and `algcurves/realcurve`(curve,x,y) then
		if not assigned(Pp) then
			Pp:=Permuted(preimages,map(conjugate,preimages))
		fi;
		Q:=Permuted(permutations[conjugate(i)],[seq(j,j=1..n)]);
		# Q is the inverse permutation for conjugate(i).
		permutations[i]:=[seq(Pp[Q[Pp[j]]],j=1..n)]
	    else
		orderedpt:=preimages;
		for j in paths[i] do
			orderedpt:=Propagate(Analyticcontinuation(curve,x,y,j,t,Digits),orderedpt);
		od;
		permpoint:=Propagate(Analyticcontinuation(curve,x,y,i-r[i]*exp(I*Pi*t),t,Digits),orderedpt);
		permpoint:=Propagate(Analyticcontinuation(curve,x,y,i+r[i]*exp(I*Pi*t),t,Digits),permpoint);
		permutations[i]:=Permuted(orderedpt,permpoint)
	    fi;
		check_cstruct(permutations[i],cyc_struct[i]);
		userinfo(2,'algcurves',`Computed the permutation around`,k,`of`,nops(bpoints),`branch points.`)
	od;
	[xvalue,preimages,[seq([i,permutations[i]],i=bpoints)],discpoints,paths,r,rootof]
end:

# Check if the cycle structure obtained by analytic continuation is the
# same as the one found by Puiseux expansions.
check_cstruct:=proc(perm,s)
	local a;
	option `Copyright (c) 1999 Waterloo Maple Inc. All rights reserved.`;
	if perm<>[] and not type(perm[1],list) then
		procname(convert_disjcyc(perm),s)
	else
		a:=map(nops,perm);
		if sort([1$(convert(s,`+`)-convert(a,`+`)),op(a)])<>sort(s) then
			ERROR( "Found wrong cycle structure" )
		fi
	fi;
	perm
end:

# This procedure is used to graph the paths.

Showpaths:=proc(ttable,r,t)
	local ind,i,j,pplot;
	option `Copyright (c) 1999 Waterloo Maple Inc. All rights reserved.`;
	ind:=[seq(op(i),i=indices(ttable))];
	ind:=[op({op(ind)} minus {'parameter'})];
	pplot:=[seq([Re(i)+r[i]*cos(2*Pi*t),Im(i)+r[i]*sin(2*Pi*t),t=0..1],i=ind)];
	for i in ind do
	    if i<>'parameter' then
		for j in ttable[i] do
			pplot:=[op(pplot),[Re(j),Im(j),t=0..1]];
		od;
	    fi;
	od;
	{op(pplot)}
end:

# This procedure gives the permutation that transforms 
# lijst1 to lijst2. Because of numerical roundoff, equalities
# between elements are never used.

Permuted:=proc(lijst1,lijst2)
	local k,outlijst,m,mm,j,r,S,d;
	option `Copyright (c) 1999 Waterloo Maple Inc. All rights reserved.`;
	outlijst:=NULL;
	S:={seq(k,k=1..nops(lijst1))};
	for k in lijst2 do
		m:=infinity;
		for j in S minus {outlijst} do
			mm:=abs(k-lijst1[j]);
			if mm<m then
				d:=m; # second smallest
				m:=mm;
				r:=j
			elif mm<d then
				d:=mm
			fi;
		od;
		outlijst:=outlijst,r;
		if m*10 > d then
		   ERROR("Computation inaccurate, use larger value for Digits")
		fi
	od;
	[outlijst]
end:

Propagate:=proc(R,lijst)
	local k,l;
	option `Copyright (c) 1999 Waterloo Maple Inc. All rights reserved.`;
	l:=R[-1][2];
	[seq(l[k],k=Permuted(R[1][2],lijst))]
end:

# This procedure computes the analytic continuation of 
# a set of function values on the Riemann surface, 
# along a given path. The output is a map from the initial
# points to the end points, given in the form of a table. 

Analyticcontinuation:=proc(curve,x,y,path,t,DI)
	local tstep,i,fibre,t1,dig,t_old,R,CL;
	global very_close;
	options remember, 
	    `Copyright (c) 1999 Waterloo Maple Inc. All rights reserved.`;
	# Not including the option remember results in going
	# through these slow calculations more than necessary.
	i:=exp(-I*Pi*t);
if has(path,i) then
	t1:=procname(curve,x,y,coeff(path,i,0)-coeff(path,i,1)*exp(I*Pi*t),t,DI);
	[seq([1-t1[-i][1],t1[-i][2]],i=1..nops(t1))]
else
	userinfo(6,'algcurves',`following path`,path);
	fibre:=[fsolve(subs(x= subs(t=0,path) ,curve),y,complex)];
	R:=[0,fibre];
	dig:=Digits;
	CL:=round(sqrt(max(1,dig-8)));
	if assigned(_Env_algcurves_CL) then CL:=_Env_algcurves_CL fi;
	tstep:=Initial_tstep / CL^2;
	t_old:=0;
	while t_old<1 do
		t1:=t_old + tstep;
		userinfo(9,'algcurves',`Digits, tstep, t1 `,Digits,tstep,t1);
		fibre:=Continue(CL,curve,x,y,path,t,t_old,t1,fibre);
		R:=R,fibre;
		fibre:=[fibre][-1][2];
		if very_close=1 and tstep<1/(4+CL) then
			tstep:=tstep*2;
			if Digits>dig then
				Digits:=Digits-3
			fi
		elif very_close=-1 then
			tstep:=tstep/2;
			if tstep < Increase_Digits / CL^2 then
				Digits:=Digits + 3
			fi
		fi;
		tstep:=min(tstep,1-t1);
		t_old:=t1
	od;
	[R]
fi
end:

# This procedure continues a given vector along a given
# path from t0 to t1. 

Continue:=proc(CL,curve,x,y,path,t,t0,t1,fibre,nested)
	local newy,neary,d,der,i,j,x0,x1,m,outy,
	      closeenough,root,mm,newfibre,Fx1;
	global very_close;
option `Copyright (c) 1999 Waterloo Maple Inc. All rights reserved.`;
if nargs=9 then
	RETURN(procname(args,0))
fi;
      Digits:=Digits+EXTRA1;
	x0:=evalf(subs(t=t0,path));
	x1:=evalf(subs(t=t1,path));
	Fx1:=subs(x=x1,curve);
      Digits:=Digits-EXTRA1;
if nargs>10 then
	newy:=args[11]
else
	newy:={fsolve(Fx1,y,complex)};
fi;
	newfibre:=newy;
      Digits:=Digits + EXTRA2;
	der:=(x1-x0)*subs(x=x0,Der(curve,x,y,Digits));
	neary:=[seq(i+subs(y=i,der),i=fibre)];
      Digits:=Digits - EXTRA2;
	# neary should be a nearby permutation of newy. It will
	# be used to reorder newy. 
	if nops(neary)<>nops(newy) then 
		ERROR(`Cannot determine all roots`)
	fi;
	outy:=[];
	closeenough:=true;
	very_close:=1;
	for i in neary do
		# d:=min(seq(evalf(abs(i-j)),j={op(neary)} minus {i}));
		m:=infinity;
		for j in newy do
			mm:=abs(i-j);
			if mm<m then
				d:=m;   # d = second closest
				m:=mm;  # m = closest
				root:=j;
			elif mm<d then
				d:=mm
			fi
		od;
		outy:=[op(outy),root];
		newy:=newy minus {root};
		if m>d/CLOSE2 / CL^2 then
			very_close:=0
		fi;
		if m>d/CLOSE / CL^2 then closeenough:=false;
		fi;
	od;
	newy:=fibre;
	if nested > 0 then
		very_close:=-1;
		Digits:=Digits + 3*nested
	fi;
	if closeenough then
		RETURN([t1,outy])
	fi;
	if nested > 1 then
		newfibre:=NULL;
		newy:=[fsolve(subs(x=x0,curve),y,complex)];
		newy:=[seq(newy[i],i=Permuted(newy,fibre))]
	else
		newy:=fibre
	fi;
	# The roots are not close enough
	# We'll repeat the continuation, but with
	# smaller stepsize
	userinfo(6,'algcurves',`nested, Digits, x0, x1`,nested,Digits,x0,x1);
	# i = intermediate
	i:=procname(args[1..7],(t0+t1)/2,newy,nested+1);
	i, procname(args[1..6],(t0+t1)/2,t1,[i][-1][2],nested+1,newfibre)
end:

# This procedure computes the first derivative y'(x) on the 
# Riemann surface, using implicit differentiation.

Der:=proc(curve,x,y,DI)
	options remember, 
	    `Copyright (c) 1999 Waterloo Maple Inc. All rights reserved.`;
	# unapply(-diff(curve,x)/diff(curve,y),x,y);
	evalf(collect(normal(-diff(curve,x)/diff(curve,y)),y),DI)
end:

# This procedure sorts a list of complex numbers 
# according to their distance from a given point x0.
# If two distances are equal, smallest argument comes
# first.

Distsort:=proc(points,x0)
local dec;
	option `Copyright (c) 1999 Waterloo Maple Inc. All rights reserved.`;
	dec:=(s,t)->evalb(
		abs(s-x0)<abs(t-x0) or
		(Im(s-x0)*Re(t-x0)<Im(t-x0)*Re(s-x0) and 
		 abs(s-x0)=abs(t-x0)
		)
		  	 );
	sort(points,dec)
end:

# This procedure gives as output a table of paths, one around each branchpoint, 
# given in the form of a table.

Makepaths:=proc(t,pnts,x0,r)
	local paths,i,j,used,connecting,p,q,
	      outpt,pts,inpt,ptj,ptjp1,ptjm1,linepaths,points,
	      pointseq,m,k,rampoints,n,tau,
	      circlepaths,connectionpt,newpath,rr,jj,
	      pj,minp,maxp,notseq,tau1,tau2,p1,p2,mm,nn,omega1,omega2;
	# make a graph with the discriminant points at the vertices
	option `Copyright (c) 1999 Waterloo Maple Inc. All rights reserved.`;
	points:=pnts;
	used:=[points[1]];
	pointseq:=table();
	pointseq[points[1]]:=[points[1]];
	connecting:=table();
	points:=Distsort(points,x0);
	for k from 2 to nops(points) do 
		q:=points[k];
		pointseq[q]:=[points[1]];
		# The path runs from the basepoint directly to
		# points[k], unless this straight line comes to close 
		# to one of the discriminant points. 
		for i in [seq(used[j],j=2..nops(used))] do
			# if the line from the last point in pointseq comes closer to 
			# point i than r[i], than this is too close, and point i is
			# added to pointseq. This way, we avoid getting too close to 
			# the discriminantpoints. 
			p:=op(nops(pointseq[q]),pointseq[q]);
			m:=p-r[p];
			if abs(p+r[p]-q)<abs(m-q) then 
				m:=p+r[p];
			fi;
			n:=q-r[q];
			if abs(q+r[q]-m)<abs(n-m) then 
				n:=q+r[q];
			fi;
			tau:=Re(evalc(conjugate(i-m))*(n-m))/abs(n-m)^2;
			# Closer than r[i]? Then reroute:
			if abs(Im(evalc(conjugate(i-m))*(n-m))/(n-m)) <abs(r[i]) and 
			   (tau>=0 and tau<=1)
			then 
				pointseq[q]:=pointseq[i];
			fi;
		od;
		# connect the last vertex using a point to the left
		# or a point to the right, whichever one is closest. 
		p:=op(nops(pointseq[q]),pointseq[q]);
		m:=p-r[p];
		if abs(p+r[p]-q)<abs(m-q) then 
			m:=p+r[p];
		fi;
		n:=q-r[q];
		if abs(q+r[q]-m)<abs(n-m) then 
			n:=q+r[q];
		fi;
		connecting[p,q]:=[m,n];
		pointseq[q]:=[op(pointseq[q]),q];
		used:=[op(used),q];
		# correcting for the points that are squeezed between
		# the path and the line from the starting point to the
		# discriminant point.
		
		notseq:=[op({op(points)} minus {op(pointseq[q])})];
		i:=1;
		while i<=nops(pointseq[q])-1 do
			m,n:=op(connecting[pointseq[q][i],pointseq[q][i+1]]);
			mm:=x0+(q-x0)*Re(evalc(conjugate(m-x0))*(q-x0))/abs(q-x0)^2;
			nn:=x0+(q-x0)*Re(evalc(conjugate(n-x0))*(q-x0))/abs(q-x0)^2;
			j:=1;
			while j>0 and j<=nops(notseq) do
			  pj:=notseq[j];
			  tau1:=Re(evalc(conjugate(pj-m))*(n-m))/abs(n-m)^2;
			  tau2:=Re(evalc(conjugate(pj-mm))*(nn-mm))/abs(nn-mm)^2;
			  if Re(n-m)<>0 then
			     omega1:=Im(n-m)/Re(n-m);
			     p1:=Im(pj)-omega1*Re(pj)-(Im(m)-omega1*Re(m));
			     omega2:=Im(q-x0)/Re(q-x0);
			     p2:=Im(pj)-omega2*Re(pj)-(Im(x0)-omega2*Re(x0));
			  else
			     omega1:=0;
			     p1:=Re(pj-m);
			     omega2:=Re(q-x0)/Im(q-x0);
			     p2:=Re(pj)-omega2*Im(pj)-(Re(x0)-omega2*Im(x0));
			  fi;
			  if (omega1*omega2>=0 and p1*p2<0 and tau1>=0 and tau1<=1 and tau2>=0 and tau2<=1) or 
			     (omega1*omega2<0 and p1*p2>0 and tau1>=0 and tau1<=1 and tau2>=0 and tau2<=1)
			  then
			  	pointseq[q]:=[seq(pointseq[q][jj],jj=1..i),pj,seq(pointseq[q][jj],jj=i+1..nops(pointseq[q]))];
			  	m:=pointseq[q][i]-r[pointseq[q][i]];
			  	if abs(pointseq[q][i]+r[pointseq[q][i]]-pj)<abs(m-pj) then
			  		m:=pointseq[q][i]+r[pointseq[q][i]];
			  	fi;
			  	n:=pj-r[pj];
			  	if abs(pj+r[pj]-m)<abs(n-m) then 
			  		n:=pj+r[pj];
			  	fi;
			  	connecting[pointseq[q][i],pj]:=[m,n];
			  	m:=pj-r[pj];
			  	if abs(pj+r[pj]-pointseq[q][i+2])<abs(m-pointseq[q][i+2]) then
			  		m:=pj+r[pj];
			  	fi;
			  	n:=pointseq[q][i+2]-r[pointseq[q][i+2]];
			  	if abs(pointseq[q][i+2]+r[pointseq[q][i+2]]-m)<abs(n-m) then 
			  		n:=pointseq[q][i+2]+r[pointseq[q][i+2]];
			  	fi;
			  	connecting[pj,pointseq[q][i+2]]:=[m,n];
			  	notseq:=[op({op(notseq)} minus {pj})];
			  	j:=0;
			  else j:=j+1;
			  fi;
			od;
			if j>nops(notseq) then i:=i+1 fi;
		od;
	od;
	
	# sort the discriminant points by argument
	# rearrange the radii accordingly
	rampoints:=Argsort(points,x0);
	
	# compose the paths from the non-special point around each
	# discriminant point.
	paths:=table();
	paths['parameter']:=t;
	for i in rampoints do
	if nops(pointseq[i])>1 then
		# the straight line pieces
		linepaths:=[];
		for j from 2 to nops(pointseq[i]) do
			pts:=connecting[pointseq[i][j-1],pointseq[i][j]];
			linepaths:=[op(linepaths),(1-t)*pts[1]+t*pts[2]];
		od;
		
		# the semi-circle pieces
		
		# the first semi-circular piece
		pts:=[subs(t=0,linepaths[1]),subs(t=1,linepaths[1])];
		rr:=r[pointseq[i][1]];
		if abs(pts[1]-pointseq[i][1]+rr)<10^(-Digits+2) 
		then
		  if Im(pts[2]-pts[1])<0 and Im(i-pointseq[i][1])>=0 then
		  	circlepaths:=[[pointseq[i][1]-rr*exp(-I*Pi*t),pointseq[i][1]+rr*exp(-I*Pi*t)]];
		  elif Im(pts[2]-pts[1])>0 and Im(i-pointseq[i][1])<0 then
		  	circlepaths:=[[pointseq[i][1]-rr*exp(I*Pi*t),pointseq[i][1]+rr*exp(I*Pi*t)]];
		  else circlepaths:=[[]];
		  fi;
		else
		  if Im(i-pointseq[i][1])>=0
		  then circlepaths:=[[pointseq[i][1]-rr*exp(-I*Pi*t)]];
		  else circlepaths:=[[pointseq[i][1]-rr*exp(I*Pi*t)]];
		  fi;
		fi;
		
		# all circular pieces but the first and last one
		for j from 2 to nops(pointseq[i])-1 do
			newpath:=[];
			ptj:=pointseq[i][j];
			ptjm1:=pointseq[i][j-1];
			ptjp1:=pointseq[i][j+1];
			inpt:=connecting[ptjm1,ptj][2];
			outpt:=connecting[ptj,ptjp1][1];
			# smaller argument
			if Im(i-x0)/Re(i-x0)<Im(ptj-x0)/Re(ptj-x0) then
			  if (inpt=ptj-r[ptj]) and (outpt=ptj+r[ptj]) then
			  	newpath:=[ptj-r[ptj]*exp(I*Pi*t)];
			  elif (inpt=ptj+r[ptj]) and (outpt=ptj-r[ptj]) then
			  	newpath:=[ptj+r[ptj]*exp(I*Pi*t)];
			  elif (inpt=ptj-r[ptj]) and (outpt=inpt) and (Im(ptjp1-ptjm1)>0) then
			  	newpath:=[ptj-r[ptj]*exp(I*Pi*t),ptj+r[ptj]*exp(I*Pi*t)];
			  elif (inpt=ptj+r[ptj]) and (outpt=inpt) and (Im(ptjp1-ptjm1)<0) then
			  	newpath:=[ptj+r[ptj]*exp(I*Pi*t),ptj-r[ptj]*exp(I*Pi*t)];
			  fi; 
			# bigger argument
			else
			  if (inpt=ptj-r[ptj]) and (outpt=ptj+r[ptj]) then
			  	newpath:=[ptj-r[ptj]*exp(-I*Pi*t)];
			  elif (inpt=ptj+r[ptj]) and (outpt=ptj-r[ptj]) then
			  	newpath:=[ptj+r[ptj]*exp(-I*Pi*t)];
			  elif (inpt=ptj-r[ptj]) and (outpt=inpt) and (Im(ptjp1-ptjm1)<0) then
			  	newpath:=[ptj-r[ptj]*exp(-I*Pi*t),ptj+r[ptj]*exp(-I*Pi*t)];
			  elif (inpt=ptj+r[ptj]) and (outpt=inpt) and (Im(ptjp1-ptjm1)>0) then
			  	newpath:=[ptj+r[ptj]*exp(-I*Pi*t),ptj-r[ptj]*exp(-I*Pi*t)];
			  fi;
			fi;
			circlepaths:=[op(circlepaths),newpath]
		od;
		# the last semi-circular piece
		connectionpt:=connecting[pointseq[i][-2],i][2];
		if abs(connectionpt-i+r[i])<10^(-Digits+2) then
			circlepaths:=[op(circlepaths),[]];
		elif Im(i-x0)>=0 then
			newpath:=[i+r[i]*exp(-I*Pi*t)];
			circlepaths:=[op(circlepaths),newpath]
		else 
			newpath:=[i+r[i]*exp(I*Pi*t)];
			circlepaths:=[op(circlepaths),newpath]
		fi;
		# put the two pieces together
		paths[i]:=op(circlepaths[1]);
		for j from 2 to nops(pointseq[i]) do
			paths[i]:=paths[i],linepaths[j-1];
			if circlepaths[j]<>[] then
				paths[i]:=paths[i],op(circlepaths[j])
			fi;
		od;
		paths[i]:=[paths[i]];
	else
		paths[i]:=[];
	fi;
	od;
	eval(paths)
end:



# This procedure sorts a list of complex numbers
# by real part first and then by imaginary part.
# When two numbers have the same real part, real
# numbers come before complex numbers.

Sortnumbers:=proc(ll)
	option `Copyright (c) 1999 Waterloo Maple Inc. All rights reserved.`;
	sort(ll,(s,t)->evalb(
		s=t or
		Re(s)<Re(t) or
		(Re(s)=Re(t) and type(s,realcons)) or 
		(Re(s)=Re(t) and Im(s)<Im(t) and not type(s,realcons) and not
			type(t,realcons))
	))
end:



# This procedure sorts a list of complex numbers ll
# according to smallest argument with respect to a
# reference point b to the left of all points in ll

Argsort:=proc(ll,b)
	local dec;
	option `Copyright (c) 1999 Waterloo Maple Inc. All rights reserved.`;
	dec:=(s,t)->evalb(s=b or
		 Im(s-b)*Re(t-b)<Im(t-b)*Re(s-b) or
		(Im(s-b)*Re(t-b)=Im(t-b)*Re(s-b) and 
		 Im(s-b)^2+Re(s-b)^2<Im(t-b)^2+Re(t-b)^2 ));
	sort(ll,dec)
end:

# Avoid integration around points with positive imaginary part. To turn this
# off, for debugging, set: _Env_algcurves_real:=false;
#
`algcurves/realcurve`:=proc(f,x,y)
options remember,`Copyright (c) 1999 Waterloo Maple Inc. All rights reserved.`;
	evalb(_Env_algcurves_real<>false and
	 #R5 code: map(Im,map(evalf,{coeffs(f,{x,y})}))={0})
	 not hastype(evalf(f),nonreal))
end:

#savelib('Monodromy','Linearmonodromy','`algcurves/realcurve`',\
 'Generalmonodromy','check_cstruct','Showpaths','Permuted','Propagate',\
 'Analyticcontinuation','Continue','Der','Distsort','`algcurves/fsolve`',\
 'Makepaths','Sortnumbers','Argsort','Hurwitzsystem' ):