# Author: Mark van Hoeij, e-mail address:    hoeij@sci.kun.nl
# If you change/improve this program, or if you find a bug, please let me know

# Version: Dec 1995
# Contents:
#   Puiseux series computation
#   Normal forms for curves with genus 0 and 1 (in the case g=0 this is called
#     a parametrization and if g=1 then a Weierstrass normal form).
#   An algorithm by Harm Derksen and me to make nice tubeplots of
#     singularities knots
# References:
#   An algorithm for computing an integral basis in an algebraic
#     function field, J. Symbolic Computation, 18, 353-363 (1994)
#   Computing parametrizations of rational algebraic curves,
#     ISSAC '94 Proceedings, 187-190 (1994)
#   An algorithm for computing the Weierstrass normal form,
#     ISSAC '95 Proceedings, 90-95 (1995)

lprint(`Additional documentation and the most recent version of this program is`);
lprint(`available via WWW at  http://www.sci.kun.nl/math/compalg/IntBasis/`);
lprint();
lprint(`For help about Puiseux expansions type:  ?puiseux`);
lprint(`For help about the integral basis type:  ?integral_basis`);
lprint(`For help about genus and parametrization computation type: ?genus`);

macro(	bug=`genus/bug`,
	compute_x=`genus/compute_x`,
	cont_exp=`puiseux/cont_exp`,
	cont_exp_m1=`puiseux/cont_exp_m1`,
	degree_ext=`puiseux/degree_ext`,
	double_factors=`integral_basis/double_factors`,
	el_of=`puiseux/el_of`,
	ext_needed=`integral_basis/ext_needed`,
	ext_to_coeffs=`integral_basis/ext_to_coeffs`,
	find_point=`genus/find_point`,
	function_with_one_pole=`genus/function_with_one_pole`,
	g_conversion1=`puiseux/g_conversion1`,
	g_conversion2=`puiseux/g_conversion2`,
	g_denom=`puiseux/g_denom`,
	g_evala=`puiseux/g_evala`,
	g_evala_rem=`puiseux/g_evala_rem`,
	g_expand=`puiseux/g_expand`,
	g_ext_r=`puiseux/g_ext_r`,
	g_ext=`puiseux/g_ext`,
	g_factor=`puiseux/g_factor`,
	g_normal=`puiseux/g_normal`,
	g_numer=`puiseux/g_numer`,
	g_solve=`puiseux/g_solve`,
	g_zero_of=`puiseux/g_zero_of`,
	homogeneous=`genus/homogeneous`,
	monic=`puiseux/monic`,
	normal_tcoeff=`puiseux/normal_tcoeff`,
	singularities=`genus/singularities`,
	trancendental_ext=`puiseux/trancendental_ext`,
	truncate=`puiseux/truncate`,
	truncate_subs=`puiseux/truncate_subs`,
	v_ext_m=`puiseux/v_ext_m`
):

#############################
#  groundfield computation  #
#############################

g_conversion1:={}: # RootOf syntax -> my own syntax
g_conversion2:={}: # my syntax -> RootOf syntax
trancendental_ext:={}: # set of trancendental extensions

# Note: this code is written for constants field which are
# algebraic extensions of Q. For other constants fields it is not
# well tested, and certainly not optimized.

g_solve:=proc()
	if trancendental_ext={} then solve(evala(args))
	# evala needed because of a bug in solve
	else solve(evala(args),indets(args) minus trancendental_ext)
	fi
end:

bug:=proc()
	lprint(` Bug alert: please send this example to hoeij@sci.kun.nl`);
	ERROR(args)
end:

# Factorization over the groundfield
g_factor:=proc(ff,ext)
local f;
	f:=numer(ff);
	if ext=[] then factor(f)
	else
		subs(g_conversion1,evala(Factor(subs(g_conversion2,f)
		 ,op(subs(g_conversion2,ext)))))
	fi
end:

# g_ext: gives a list of the algebraic extensions.
g_ext_r:=proc(a)
local v,vv,i,tail;
options remember;
	v:=indets(a,RootOf);
	if nops(v)=0 then RETURN([]) fi;
	vv:={};
	for i in v do vv:=vv union indets(op(i),RootOf) od;
	tail:=g_ext_r(vv);
	v:=[op(v minus vv)];
	[op(v),op(tail)]
end:

######## NOTE ###########
# g_ext does not work properly if there are nested RootOfs that are treated
# the the alias procedure. The only way I know to "fix" this is by giving
# everyone the advice never (really never!) to use alias, because who knows
# what other bugs^H^H^H^Hfeatures the use of alias can introduce.

# Gives the algebraic extensions appearing in aa.
g_ext:=proc(a)
global g_conversion1,g_conversion2;
local v,i,result,ii,vv;
options remember;
	v:=g_ext_r(subs(g_conversion2,a));
	vv:=subs(g_conversion1,v);
	result:=NULL;
	for i from 0 to nops(g_conversion2) do if member(`puiseux/rootof`.i,vv) then
		result:=`puiseux/rootof`.i,result
	fi od;
	for i from nops(v) by -1 to 1 do if not member(subs(g_conversion1,v[i])
	 ,{seq(`puiseux/rootof`.ii,ii=0..nops(g_conversion2))}) then
		if g_conversion1={} then g_conversion1:=NULL fi;
		vv:=nops(g_conversion2);
		g_conversion1:=v[i]=`puiseux/rootof`.vv,g_conversion1;
		g_conversion2:={`puiseux/rootof`.vv=v[i],op(g_conversion2)};
		result:=subs(g_conversion1,v[i]),result
	fi od;
	[result]
end:

# Gives the algebraic degree of a over b
degree_ext:=proc(aa,bb)
local a,b,v,i,all,d,var;
options remember;
	a:=subs(g_conversion2,aa);
	b:=subs(g_conversion2,bb);
	v:=indets(a,RootOf) minus indets(b,RootOf);
	all:=[op(indets([a,b],RootOf))];
	all:={seq(all[i]=var[i],i=1..nops(all))};
	d:=1;
	for i in v do
		d:=d*degree(subs(all,op(i)),_Z)
	od;
	d
end:

# Some procedures are given twice. The first is meant to describe in Maple
# language what the second procedure does. The second procedures uses
# a different representation of algebraic numbers. I only wrote a replacement
# for evala(Expand()). For Normal and Factor I call the Maple procedures.
# If the evala(Expand()) ever gets faster then one can remove these second
# procedures and the code should still work.

# g_expand:=proc(a,ext) evala(Expand(a)) end:

# Give 3'rd argument x for normalizing the lowest coefficient
g_expand:=proc(a,ext)
	if type(a,polynom(anything,[op(ext),op(trancendental_ext)])) or
	(nargs<3 and type(a,polynom(anything,ext)))
		then g_evala(expand(a),ext)
	elif type(a,polynom(anything,ext)) # Note: nargs=3
		then normal_tcoeff(g_evala(expand(a),ext),args[3])
	elif nargs<3 then
		subs(g_conversion1,evala(Expand(subs(g_conversion2,a))))
	else
		normal_tcoeff(g_expand(a,ext),args[3])
	fi
end:

normal_tcoeff:=proc(a,x)
local c,i,r,d;
r:=a;
do
	d:=ldegree(r,x);
	c:=g_normal(coeff(r,x,d));
	if coeff(r,x,d)=c then RETURN(r) fi;
	r:=convert([c*x^d,seq(coeff(r,x,i)*x^i,i=d+1..degree(r,x))],`+`);
	if c<>0 then RETURN(r) fi
od
end:

# g_evala:=proc(a,ext) evala(Expand(a)) end:

g_evala:=proc(a,ext)
	local dummy,e;
	if nops(ext)=0 then RETURN(a)
	elif nops(ext)=1 then
		e:=ext[1];
		expand(convert([seq(coeff(a,e,dummy)*g_evala_rem(e^dummy)
		 ,dummy=0..degree(a,e))],`+`))
	else
	e:=g_evala(a,[ext[2..nops(ext)]]);
	g_evala(expand(convert([seq(coeff(e,ext[1],dummy)*g_evala_rem(ext[1]^dummy)
	 ,dummy=0..degree(e,ext[1]))],`+`)),[ext[2..nops(ext)]])
	fi
end:

g_evala_rem:=proc()
	options remember;
	expand(subs(g_conversion1,evala(Expand(subs(g_conversion2,args)))))
end:

# Bugfix: numer and denom sometimes factor the numerator and denominator.
# This causes: "ERROR could not compute coefficient" in some cases.
# Also expanding won't be useful anymore in g_normal.
# The following two procedures are meant to save some typing:
g_numer:=proc() expand(numer(args)) end:
g_denom:=proc() expand(denom(args)) end:

# works for both syntaxes
g_normal:=proc(aa)
	local a;
	if has(aa,RootOf) then 
		RETURN(evala(Normal(aa)))
	fi;
	a:=subs(g_conversion2,aa);
	if has(a,RootOf) then
		subs(g_conversion1,evala(Normal(a)))
	else
		normal(a)
	fi
end:

# Input : an irreducible polynomial kk in x, not necessarily monic
# Output: a zero of k
# If an algebraic extension is needed it will be placed in ext.
#g_zero_of:=proc(k,x,ext)
#	if degree(k,x)=1 then
#		ext:=NULL;
#		RETURN(evala(-coeff(k,x,0)/coeff(k,x,1)))
#	fi;
#	ext:=RootOf(k,x)
#end:

g_zero_of:=proc(k,x,ext)
	global  g_conversion1, g_conversion2;
	local a,vv;
	if degree(k,x)=1 then
		ext:=NULL;
		RETURN(g_normal(-coeff(k,x,0)/coeff(k,x,1)))
	fi;
	a:=RootOf(subs(g_conversion2,k),x);
	# Bugfix 7 Aug 1994: one must avoid having multiple
	# `puiseux/rootof`s for the same RootOf
	if not member(_Z,indets(op(subs(g_conversion1,a)))) then
		a:=subs(g_conversion1,a);
		ext:=a;
		RETURN(a)
	fi;
	if g_conversion1={} then g_conversion1:=NULL fi;
	vv:=nops(g_conversion2);
	g_conversion1:=a=`puiseux/rootof`.vv,g_conversion1;
	g_conversion2:={`puiseux/rootof`.vv=a,op(g_conversion2)};
	ext:=`puiseux/rootof`.vv
end:

# Gives the zeros of the factors, their multiplicities and algebraic extensions
v_ext_m:=proc(f,x)
local ext,nulp,i,result;
	if degree(f,x)=0 then RETURN({}) fi;
	if type(f,`^`) then
		nulp:=g_zero_of(op(f)[1],x,'ext');
		ext:=eval(ext);
		RETURN({[nulp,op(f)[2],[ext],degree(op(f)[1],x)]})
	fi;
	if type(f,`*`) then
		result:={};
		for i in {op(f)} do
			result:=result union v_ext_m(i,x)
		od;
		RETURN(result)
	fi;
	nulp:=g_zero_of(f,x,'ext');
	ext:=eval(ext);
	{[nulp,1,[ext],degree(f,x)]}
end:

# ext_to_coeffs does basically: coeffs(expression,RootOf( .. ))
# We need this procedure because coeffs doesn't always work this way.
# a must be evala'd
#ext_to_coeffs:=proc(a,ext)
#local r,i,j;
#	r:=indets(a,RootOf) minus indets(ext,RootOf);
#	if r={} then RETURN(a) fi;
#	for i in r do if not member(i,indets(r minus {i},RootOf)) then
#		r:=[coeffs(g_expand(subs(i=j,a)),j)];
#		RETURN(seq(ext_to_coeffs(j,ext),j=r))
#	fi od;
#	bug()
#end:

ext_to_coeffs:=proc(a,ext)
	local dummy,aa;
	aa:=(indets(a) minus indets(ext)) intersect
	{seq(`puiseux/rootof`.dummy,dummy=0..nops(g_conversion2))};
	coeffs(a,[op(aa)])
end:

########################
# Puiseux computation  #
########################

# The puiseux code is not optimally implemented yet, I should start
# using so-called "rational" (i.e. minimal alg extension) puiseux series,
# but it can take a while before this is implemented because it's not on
# top of my to-do list.

# Takes the lowest coefficients
truncate:=proc(aa,n,x,ext)
local dummy,a;
	a:=collect(aa,x);
	a:=g_expand(convert([seq(x^dummy*coeff(a,x,dummy),dummy=0..n-1)],`+`),ext)
end:

#truncate_subs:=proc(f,x,y,y_value,n,ext)
#	truncate(subs(y=y_value,f),n,x,ext)
#end:

truncate_subs:=proc(f,x,y,y_value,n,ext)
local ym,i,result;
	if degree(f,x)>=n and n>=1 then
		RETURN(truncate_subs(truncate(f,n,x,ext),x,y,y_value,n,ext))
	fi;
	result:=0;
	for i from 0 to degree(f,y) do
		if i=0 then ym:=1 else ym:=truncate(ym*y_value,n,x,ext) fi;
		result:=result+coeff(f,y,i)*ym
	od;
	result:=truncate(result,n,x,ext);
	# Make sure that ldegree(result,x) gives a right answer:
	normal_tcoeff(result,x)
end:

# Makes f monic, returns also q, the factor needed to multiply y by, to
# get it integral.
monic:=proc(f,y,ext,q)
	local dummy,ff,lc,qq;
	ff:=g_numer(g_normal(f));
	lc:=lcoeff(ff,y);
	if indets(lc,`name`) minus {op(ext),op(trancendental_ext)}={} then
		q:=1;
		RETURN(g_expand(g_normal(ff/lc),ext))
	fi;
	lc:=subs(g_conversion1,evala(Sqrfree(subs(g_conversion2,lc))));
	lc:=lc[2];
	lc:=g_expand(product(lc[dummy][1],dummy=1..nops(lc)),ext);
	ff:=monic(subs(y=y/lc,ff),y,ext,'qq');
	qq:=eval(qq);
	q:=g_expand(qq*lc,ext);
	ff
end:

# See ?puiseux
puiseux:=proc(aa:algfun(rational),eqn:name=algfun(rational),n:numeric)
global trancendental_ext;
local x,v,a,f,y,ext,q,verz,i,result,ma,dummy,dummy2,qx;
	if not type(aa,RootOf) then ERROR(`First argument must be a RootOf`) fi;
	x:=eval(op(1,eqn));
	v:=op(2,eqn);
	a:=subs(x=x+v,aa);
	ext:=g_ext(op(a));
	f:=subs(g_conversion1,op(a));
	f:=subs(_Z=y,f);
	f:=monic(f,y,ext,'q');
	q:=eval(q);
	qx:=ldegree(q,x);
	q:=g_expand(q/x^qx,ext);
	trancendental_ext:=indets(f) minus
		{x,y,seq(`puiseux/rootof`.i,i=0..nops(g_conversion2))};
	verz:=`puiseux/technical_answer`(f,x,y,n+qx,ext);
	ma:=ceil(max(seq(dummy[2]/dummy[3],dummy=verz)))+qx;
	q:=taylor(1/q,x=0,ma);
	q:=g_expand(sum('x^dummy2*coeff(q,x,dummy2)',dummy2=0..ma-1),ext);
	result:={};
	for i in verz do
		result:=result union {subs(x=x^(1/i[3]),
		 truncate(subs(x=x^i[3],q)*i[1],i[2],x,ext))/x^qx}
	od;
	subs(g_conversion2,subs(x=x-v,result))
end:

# This procedure computes Puiseux expansions.
# It's input must be in internal format.
# Output: a list, see cont_exp
`puiseux/technical_answer`:=proc(f,x,y,n,ext)
local result,i,verz;
options remember;
	verz:=cont_exp([0,0,1,ext,degree(f,y),1,0],f,x,y);
	result:={};
	for i in verz do
		result:=result union cont_exp_m1(i,f,x,y,n)
	od;
result
end:

# This procedure continues an expansion until it has multiplicity 1
# Input is a list v
# v[1] = The expansion so far
# v[2] = This expansion is determined modulo x^(v[2]/v[3])
# v[3] = Least common multiple of the denominators of the powers of x
# v[4] = The algebraic extensions in this expansion
# v[5] = The multiplicity of this factor
# v[6] = The degree of the algebraic extensions above the groundfield
# It will also return (not needed in the input)
# v[7] = sum of the valuations ( v(x^(1/d))=1/d ) of the differences
# of this expansion with the other expansions. This sum is needed for
# computing an integral basis in an algebraic function field.
# f    = minimal polynomial of y over L(x), must be monic
# Comment: this procedure is an awful hack, the reason is that at the time I
#  wrote this I didn't know what a Newton polygon was. Sorry about that.
cont_exp:=proc(v,f,x,y)
local t,n,a,i,ii,r,som,result,vv7,tabel,ld;
	if v[5]=1 then RETURN({v}) fi;
	result:={};
	r:=v[1]+a*x^v[2];
	# we try to find an equation for the unknown a
	i:=1;
	ii:=0;
	if v[1]=0 then ii:=g_expand(subs({x=x^v[3],y=r},f),v[4],x)
	else
		while ii=0 do
		   i:=i+v[3]+2;
		   ii:=truncate_subs(subs(x=x^v[3],f),x,y,r,v[2]+
		    v[3]*v[7]+i,v[4])
		od
	fi;
	vv7:=(ldegree(ii,x)-v[2])/v[3];
	r:=tcoeff(ii,x);
	ld:=ldegree(r,a);
	if degree(r,a)>0 then
		r:=g_expand(r*g_normal(1/lcoeff(r,a)),v[4],x);
		if v[5]=0 then r:=expand(r/a^ldegree(r,a)) fi;
		if degree(r,a)>1 then
			r:=g_factor(r,v[4])
		fi;
		r:=v_ext_m(r,a);
		for i in r do
			result:=result union cont_exp([v[1]+x^v[2]
			 *i[1],v[2]+1,v[3],[op(i[3]),op(v[4])],i[2],v[6]*i[4]
			 ,vv7],f,x,y)
		od
	fi;
	som:=0;
	for ii in result do
		som:=som+ii[5]*ii[6]/v[6]
	od;
	n:=2;
	if som<v[5] then tabel:={[0,1,v[5],som]} else tabel:={} fi;
	while tabel<>{} do for t from 1 to n-1 do
	if igcd(t,n)=1 and el_of(t/n,tabel) then
		r:=cont_exp([subs(x=x^n,v[1]),(v[2]-1)*n+t,v[3]*n,
		 v[4],0,v[6],vv7],f,x,y);
		result:=result union r[2];
		som:=0;
		for i in r[2] do
			som:=som+i[5]*i[6]/v[6]
		od;
		tabel:=el_of(t/n,tabel,r[1],som)
	fi od;
	n:=n+1;
	od;
if v[5]=0 then result:=[ld,result] fi;
result
end:

# used for deciding which t/n to use in cont_exp
el_of:=proc(q,tabel,rp,som)
local i,v,r;
	if nargs=2 then
		for i in tabel do if i[1]<q and i[2]>q then
			RETURN(true)
		fi od;
		RETURN(false)
	fi;
	for i in tabel do if i[1]<q and i[2]>q then v:=i fi od;
	r:=tabel minus {v};
	if v[3]>rp+som then r:=r union {[v[1],q,v[3],rp+som]} fi;
	if rp>v[4] then r:=r union {[q,v[2],rp,v[4]]} fi;
	r
end:

# This procedure continues (= computes more terms of) expansions that have
# multiplicity 1. This could also be done with cont_exp, but this
# procedure is faster.
cont_exp_m1:=proc(v,f,x,y,nnk)
local nk,a,result,machtx,r,rr;
	nk:=ceil(nnk*v[3])/v[3];
	machtx:=v[2];
	if nk<=machtx/v[3] then RETURN({v}) fi;
	result:=v[1]+a*x^v[2];
	# r:=g_expand(subs({x=x^v[3],y=result},f),v[4]);
	# r:=truncate(r/x^ldegree(r,x),v[3]*nk-machtx,x,v[4]);
	r:=truncate_subs(subs(x=x^v[3],f),x,y,result,v[3]*(nk+v[7]),v[4]);
	r:=expand(r/x^ldegree(r,x));
	while nk>machtx/v[3] do
		rr:=coeff(r,x,0);
		if degree(rr,a)<>1 then bug(`degree is not 1`) fi;
		rr:=g_normal(-coeff(rr,a,0)/coeff(rr,a,1));
		result:=g_expand(subs(a=rr+x*a,result),v[4]);
		machtx:=machtx+1;
		if nk>machtx/v[3] then
			r:=subs(a=rr+x*a,r);
			r:=g_expand(r,v[4],x);
			if ldegree(r,x)<>1 then bug(`degree is not 1`) fi;
			r:=truncate(r/x,v[3]*nk-machtx,x,v[4])
		fi
	od;
	{[subs(a=0,result),machtx,v[3],v[4],v[5],v[6],v[7]]}
end:

`help/text/puiseux`:=TEXT(
``,
`FUNCTION:  puiseux - determines the Puiseux expansions of an algebraic function`,
``,
`CALLING SEQUENCE: puiseux(a,x=v,n)`,
``,
`PARAMETERS:`,
`   a   - an algebraic function in RootOf form`,
`   x=v - gives the point around which the expansions are computed`,
`   n   - a real number`,
``,
`SYNOPSIS:`,
`- The number of Puiseux expansions is equal to the degree of the minimal`,
`  polynomial of a. The minimal polynomial is a polynomial over a field L(x).`,
`  The procedure puiseux determines the field L from a.`,
`  The groundfield L of the computation is taken the field Q extended with`, 
`  the algebraic numbers (which must be in RootOf form) that appear in the`, 
`  minimal polynomial of a. If more than 1 litteral appears in a, then these`, 
`  litterals are also elements of L, except the litteral on the left side of`, 
`  the second argument.`,
``,
`- If a number of Puiseux expansions are algebraically conjugated over L, then`, 
`  only one of these expansions is given.`,
``,
`- The Puiseux expansions are being computed modulo x^n, so if for instance`,
`  n=10, then the term x^(49/5) would be computed, but not the term x^10`,
`  The Puiseux expansions are being computed further than x^n if this would`,
`  be necessary to seperate all different branches. So if n=0, the computation`,
`  would stop precisely at the moment when all Puiseux expansions can be`,
`  distinguished from one another, that is when all computed expansions are`,
`  different.`,
``,
`EXAMPLES:`,
` alias(alpha=RootOf(x^3+7));`,
` f:=y^8+3*x*y^5+5*x^4+x^6*alpha;`,
` alias(beta=RootOf(f,y));`,
` puiseux(beta,x=0,0);`,
` puiseux(beta,x=0,9.5);`,
``,
`BUGS: The Maple "alias" function does not recognize an alias in terms of`,
`  another alias. Therefore, you must not use nested aliasses for algebraic`,
`  numbers.`,
``,
`EXAMPLE:`,
` alpha:=RootOf(x^3+7); # alias(alpha=RootOf(x^3+7)) would go wrong`,
` alias(beta=RootOf(x^4+alpha*x+7)); # Or beta:=RootOf(x^4+alpha*x+7)`,
` f:=y^5+alpha*x*y^4+beta*x^5*y^2+7*x^7;`,
` alias(delta=RootOf(f,y)); # Or delta:=RootOf(f,y)`,
` puiseux(delta,x=0,0);`,
``
):

###############################
# Integral basis computation  #
###############################

`help/text/integral_basis`:=TEXT(
``,
`FUNCTION:  integral_basis - the integral basis of an algebraic function field`,
``,
`CALLING SEQUENCE: integral_basis(a,x)`,
``,
`PARAMETERS:`,
`   a   - an algebraic function in one variable`,
`   x   - this variable (only needed if the number of litterals in a is greater`,
`         than 1)`,
``,
`SYNOPSIS:`,
`- This procedure computes an integral basis for an algebraic function that is`,
`  given as a RootOf. The method is based on Puiseux expansions. For more`,
`  information (a LaTeX file) send an e-mail to hoeij@sci.kun.nl`,
``,
`- If printlevel > 1, some information will be printed during the computation.`,
``,
`EXAMPLE:`,
` alias(alpha=RootOf(x^3+7));`,
` # Or: alpha:=RootOf(x^3+7); `,
` f:=y^8+3*x*y^5*alpha+5*x^4+x^6*alpha;`,
` alias(beta=RootOf(f,y));`,
` printlevel:=2;`,
` integral_basis(beta);`,
``,
`SEE ALSO:  maxorder, puiseux`
):

# There are 2 extra syntaxes now:
# integral_basis(f,x,y)
# and integral_basis(f,x,y,df) where df is a set of irreducible
# polynomials in x. In this case integral_basis will return a local
# integral basis for these factors in df.

# Comment: I think it would be faster first to compute all local (i.e. for
# only one irreducible factor of the discriminant) integral basis and then
# to mix these to one global integral basis.
# The code isn't particularly fast on very simple singularities (like a
# double point). Maybe it would be a good idea to have special code for just
# the case of a factor with multiplicity 2 or 3 of the discriminant.
integral_basis:=proc()
global trancendental_ext;
local alfa,f,x,y,extl,disc,df,basis,pl_transl,k,ext,nulp,f_translated,max_v7
 ,places,d,i,b,found_something,dummy,dummy2,q,power_of_k,equations
 ,values_basis_in_places,value_new_one,kk,max_power_k;
options remember;
	alfa:=eval(args[1]);
	if nargs=1 then x:=op(indets(alfa,`name`) minus {_Z}) else x:=args[2] fi;
	if not (type(alfa,algfun(rational)) and type(x,`name`)) then 
		ERROR(`Wrong arguments`)
	fi;
	if nargs>=3 then
		f:=args[1];
		x:=args[2];
		y:=args[3];
		alfa:=y;
		extl:=g_ext(f);
		f:=subs(g_conversion1,f)
	else
		extl:=g_ext(op(alfa));
		f:=subs(g_conversion1,op(alfa));
		for i from 1 to nops(extl) do
			f:=subs(extl[i]=b[i],f)
		od;
		f:=subs(_Z=y,f);
		for i from 1 to nops(extl) do
			f:=subs(b[i]=extl[i],f)
		od
	fi;
	f:=monic(f,y,extl,'q');
	trancendental_ext:=indets(f) minus
		{x,y,seq(`puiseux/rootof`.i,i=0..nops(g_conversion2))};
	q:=eval(q);
	 # f is the minimal polynomial of y over L(x), we made it monic
	if printlevel>1 then lprint(`Computing the discriminant ...`) fi;
if nargs<=3 then
	disc:=g_expand(discrim(f,y),extl);
	if printlevel>1 then
		lprint(`disc=`,disc);
		lprint();
		lprint(`Computing multiple factors of the discriminant  ...`)
	fi;
	df:=double_factors(disc,x,extl);
	if printlevel>1 then
		lprint(`Factors with multiplicity > 1 are:`,df);
		lprint()
	fi;
	 # df contains those factors k that appear more than once in the
	 # discriminant discrim(f,y)
else
	df:=args[4]
fi;
for k in df do
	power_of_k[k]:=1;
	# Later we will search for integrals of the shape (..)/k^power_of_k[k]
	nulp:=g_zero_of(k,x,'ext');
	ext:=eval(ext);
	ext:=[ext,op(extl)];
	 # the zero of k is now in nulp, if an algebraic extension is needed
	 # it is placed in ext
	f_translated:=g_expand(subs(x=x+nulp,f),ext);
	if printlevel > 1 then
		lprint(`Computing the Puiseux expansions for the factor`,k)
	fi;
	pl_transl:=[op(`puiseux/technical_answer`(f_translated,x,y,0,ext))];
	 # Now pl_transl contains the puiseux expansions, but
	 # translated x=x+nulp
	max_v7:=max(seq(dummy[7],dummy=pl_transl));
	 # this number is the maximum internal intersection multiplicity
	 # (however, the way we count it x^(1/2) and -(x^(1/2)) have
	 # intersection multiplicity 1/2), see also in the file puiseux the
	 # procedure cont_exp
	max_power_k[k]:=floor(max_v7);
	if printlevel >1 then
		lprint(`Maximum number of factors`,k,
		 `in the denominator is`,max_power_k[k]);
		 lprint()
	fi;
	pl_transl:=[seq(op(cont_exp_m1(dummy,f_translated,x,y,min(max_v7
	 -dummy[7]+dummy[2]/dummy[3],max_power_k[k]))),dummy=pl_transl)];
	   # Here I computed more terms of the puiseux expansions
	places[k]:=[seq([dummy[1],nulp,dummy[3],
	  # The following ext_needed(..) checks if we really need a certain
	  # algebraic extension. If we don't we get rid of it.
	  # In the previous versions of IntBasis we had the following line
	  # which was a bug.
	  # ext_needed(dummy[7],max_power_k[k],dummy[4],dummy[1],x)
	 ext_needed(dummy,max_v7,x)
	],dummy=pl_transl)];
	  # places[k][dummy] contains: (dummy=1..number of puiseux expansions)
	  # 1) The puiseux expansion (but translated x=x+nulp)
	  # 2) Contains the zero of k   (= this number nulp)
	  # 3) The least common multiple of the denominators of the powers
	  # of x, so x here must be interpreted as x^(1/d) if this lcm is d.
	  # If we take into account the translation it would be x^(1/d)-nulp.
	  # So in order to interpret a normal x we substitute
	  # subs(x=x^d+nulp,..)
	  # 4) The algebraic extension belonging to this puiseux expansion
	values_basis_in_places[k,1]:=[seq(1,dummy=pl_transl)]
	   # Because the first basis-element will be 1.
od;
i:=df;
for k in i do if max_power_k[k]=0 then
	df:=df minus {k}
fi od;
basis:=[1];
for d from 2 to degree(f,y) do
	if printlevel>1 then
		lprint(`Integral elements of degree less than`,d-1,`are:`);
		lprint(`the previous elements and`,basis[d-1]);
		lprint()
	fi;
	basis:=[op(basis),y*basis[d-1]];
	  # This basis is our first guess
	for k in df do
		  # Now we compute the values of this new basis-element in all
		  # places (that is, we substitute y=puiseux expansions)
		values_basis_in_places[k,d]:=[seq(truncate
		 (values_basis_in_places[k,d-1][dummy]*places[k][dummy][1]
		 ,places[k][dummy][3]*max_power_k[k],x
		 ,places[k][dummy][4]),dummy=1..nops(places[k]))]
	od;
	for k in df do
	  # using the values of this new basis-element in the places,
	  # we will see if we can find an integral, with a bigger denominator
	found_something:=true;
	while found_something and power_of_k[k]<=max_power_k[k] do
		for i from 1 to d-1 do b[i]:=evaln(b[i]) od;
		b[i]:=1;
		 # Now we compute the values of basis[1]*b[1]+..+basis[d]*b[d]
		 # in the places
		 # and we try to put an extra factor k in the denominator
		value_new_one:=[seq(truncate(sum('b[dummy2]
		 *values_basis_in_places[k,dummy2][dummy]',dummy2=1..d)
		 ,places[k][dummy][3]*power_of_k[k],x
		 ,places[k][dummy][4]),dummy=1..nops(places[k]))];
		 # All coefficients of powers of x less than power_of_k must
		 # be zero, in order for this new_one to be dividable by
		 # k^power_of_k[k]. So we find equations by taking the
		 # remainder
		equations:={seq(ext_to_coeffs(dummy,[places[k][1][2],extl]
		 ),dummy=value_new_one)};
		equations:={seq(coeffs(dummy,x),dummy=equations)};
		equations:=subs(g_conversion1,evala({g_solve(
		subs(g_conversion2,equations))}));
		 # Now we know what values b[1] .. b[d] must have
		if equations={} then
			found_something:=false
		else
			if degree_ext(places[k][1][2],extl)>1 then
				equations:=subs(places[k][1][2]=x,equations)
			fi;
			assign(op(equations));
			  # Now basis[1]*b[1]+..basis[d]*b[d] is dividable by k
			  # In the following for_do_od we compute the values of
			  # basis[1]*b[1]+..+basis[d]*b[d] in all places.
			for kk in df do
			 values_basis_in_places[kk,d]:=[seq(
			   truncate(sum('subs(x=x^places[kk][dummy
			   ][3]+places[kk][dummy][2],eval(b[dummy2]))*
			   values_basis_in_places[kk,dummy2][dummy]',dummy2=
			   1..d),places[kk][dummy][3]
			   *max_power_k[kk],x,places[kk][dummy][4])
			   ,dummy=1..nops(places[kk]))]
			od;
			 # Now we will put (basis[1]*b[1]+..+basis[d]*b[d])/k
			 # in the basis
			basis:=[seq(basis[dummy],dummy=1..d-1),normal(sum(
			'eval(b[dummy2])*basis[dummy2]',dummy2=1..d)/k)];
			  # Now we should divide values_basis_in_places[.. ,d]
			  # by k, but instead we multiply all other
			  # values_basis_in_places[.., <>d] by k, this will
			  # be done in the following nested for_do_od
			for i from 1 to d-1 do for kk in df do
			 values_basis_in_places[kk,i]:=[seq(
			   truncate(values_basis_in_places[kk,i][
			   dummy]*subs(x=x^places[kk][dummy][3]+places[kk][
			   dummy][2],k),places[kk][dummy][3]*max_power_k[kk]
			   ,x,places[kk][dummy][4]),dummy=1..nops(places[kk]))]
			od od;
			  # Now we increase power_of_k, so we will try to put
			  # more factors k in the denominator
			power_of_k[k]:=power_of_k[k]+1
		fi
	od od od;
	subs(g_conversion2,subs(y=alfa*q,basis))
end:

# This procedure finds out if we can miss a certain alg. extension.
ext_needed:=proc(v,max_v7,x)
local dummy;
RETURN(v[4]); # this procedure contained a bug, I fixed it this way because
# the time gain was minimal anyway.
end:

# This procedure gives those factors which appear more than once
double_factors:=proc(ff,x,ext)
local f,i,i2,result,v,vv,frem,j,ex;
	f:=numer(ff);ex:=ext;
	while nops(ex)>0 do
		f:=resultant(f,
		 subs(_Z=ex[1],subs(g_conversion1,op(subs(g_conversion2,ex[1]))))
		 ,ex[1]);
		ex:=[ex[2..nops(ex)]];
		if type(f,`^`) then f:=map(g_expand,f,ex)
		elif type(f,`*`) then
			f:=convert([seq(`if`(type(i,`^`),map(g_expand,i,ex),g_expand(i,ex)),i=[op(f)])],`*`)
		else
			f:=g_expand(f,ex)
		fi
	od;
	if type(f,`*`) then
		v:=[op(f)]
	else
		v:=[f]
	fi;
	vv:=1;
	for i in v do
		if type(i,`^`) then
			vv:=vv*op(i)[1]
		else
			vv:=vv*gcd(i,diff(i,x))
		fi
	od;
	f:=readlib(factors)(vv);
	v:={};
	for i in f[2] do if degree(i[1],x)>=1 then
		v:={op(v),i[1]}
	fi od;
	result:={};
	for i in v do
		if ext=[] then
			result:=result union {g_expand(i/lcoeff(i,x),ext)}
		else
			i2:=g_expand(i*i,ext);
			i2:=g_expand(i2*g_normal(1/lcoeff(i2,x)),ext);
			frem:=ff;
			while degree(frem,x)>=degree(i2,x) do
				frem:=g_expand(numer(frem-lcoeff(frem,x)*i2*x
				 ^(degree(frem,x)-degree(i2,x))),ext)
			od;
			if frem=0 then frem:=i2 fi;
			frem:=evala(Sqrfree(subs(g_conversion2,frem)))[2];
			f:=1;
			for j in frem do if j[2]>1 then
			    f:=f*evala(Factor(j[1],subs(g_conversion2,ext)))
			fi od;
			if type(f,`*`) then vv:={op(f)} else vv:={f} fi;
			for j in vv do if evala(Rem(subs(g_conversion2,i)
			 ,j,x))=0 then
				result:={op(result),j}
			fi od
		fi
	od;
	v:={};
	for i in result do
		if degree(i,x)>0 then v:=v union {i} fi
	od;
	subs(g_conversion1,v)
end:

##########################################
# genus and parametrizations computation #
##########################################

`help/text/genus`:=TEXT(
``,
`FUNCTION:  genus - determines genus of an algebraic curve`,
``,
`CALLING SEQUENCES:`,
`   genus(f,x,y)`,
`   genus(f,x,y,``parameter``,point)`,
`   genus(f,x,y,``parametrization``,t,point)`,
``,
`PARAMETERS:`,
`   f        - a polynomial in x and y describing the curve`,
`   x, y, t  - variables`,
`   point    - optional argument to prescribe the pole of the parameter`,
``,
`SYNOPSIS:`,
`- The genus of an irreducible algebraic curve is a nonnegative integer.`,
`  However, the procedure genus can give a negative outcome. This is only`,
`  possible if the given curve is reducible.`,
``,
`- An curve with genus 0 has a parameter. A parameter is a generator of the`, 
`  function field. It can be found using the property that a parameter has`, 
`  only one pole on the curve, with pole multiplicity 1. Given a point where`, 
`  the parameter must have a pole, the parameter is unique up to addition and`, 
`  multiplication by constants.`,
``,
`- When a regular point is known on the curve that does not use large algebraic`,
`  extensions on Q, it is best to specify this point. If the point is not`,
`  specified the algorithm will look for a regular point. If it needs algebraic`,
`  extensions to find this point, these extensions will also appear in the`,
`  result. A point can be given by a list [x,y,z] of the homogeneous`, 
`  coordinates of the point.`,
``,
`- genus(f,x,y,``parameter``) or genus(f,x,y,``parameter``,point) computes a`,
`  parameter. Note that f must be irreducible, otherwise the algorithm does`,
`  not return an answer. Default is no irreducibility checking, but you can`,
`  change this by giving an extra argument ``irr check```,
``,
`- genus(f,x,y,``parametrization``,t) or genus(f,x,y,``parametrization```,
`    ,t,point) gives a parametrization [x(t),y(t)] of the curve. This`,
`  means that x(t) and y(t) satisfy the equation given by f. It first`,
`  computes a parameter, and then expresses x and y as rational functions`,
`  of this parameter.`,
``,
`EXAMPLES:`,
`f:=761328152*x^6*z^4-5431439286*x^2*y^8+2494*x^2*z^8+228715574724*x^6*y^4+`,
` 9127158539954*x^10-15052058268*x^6*y^2*z^2+3212722859346*x^8*y^2-`,
` 134266087241*x^8*z^2-202172841*y^8*z^2-34263110700*x^4*y^6-6697080*y^6*z^4-`,
` 2042158*x^4*z^6-201803238*y^10+12024807786*x^4*y^4*z^2-128361096*x^4*y^2*z^4+`,
` 506101284*x^2*z^2*y^6+47970216*x^2*z^4*y^4+660492*x^2*z^6*y^2-`,
` z^10-474*z^8*y^2-84366*z^6*y^4;`,
`  # This f is a polynomial of degree 10 having a maximal number of cusps`,
`  # according to the Plucker formulas. It was found by Rob Koelman. It has`,
`  # 26 cusps and no other singularities, hence the genus is 10.`,
`f:=subs(x=1,f);`,
`genus(f,y,z);`,
``,
`f:=x^2+x*y+y^2;`,
`factor(f); # It is irreducible in Q[x,y]`,
`genus(f,x,y); # The result is impossible for an irreducible algebraic curve`,
`evala(AFactor(f)); # And we see that f is reducible`,
``,
`f:=(y+1)*(y-1)*(y+x)*(y-x)+(x-1)^2*(x+y+2)^2;`,
`f:=expand(f);`,
`genus(f,x,y,``parameter``);`,
`v:=genus(f,x,y,``parametrization``,t);`,
`evala(Normal(subs({x=v[1],y=v[2]},f))); # Checking the result`,
``,
`f:=-y-y^2-2*y^3-y^5-x*y-x*y^3-2*x^2-2*x^2*y^2+x^3+x^4*y;`,
`genus(f,x,y,``parametrization``,t,``irr check``);`,
``,
`f:=3610/907*x^2+14440/907*x+14440/907-4508/907*x^2*y-18032/907*y*x-16227/907*y+`,
`5406/907*x^2*y^2+16218/907*x*y^2+10812/907*y^2-2703/907*x^2*y^3-12626/907*x*y^3`,
`-13533/907*y^3+3610/907*x*y^4+3610/907*y^4+1805/907*y^5;`,
`genus(f,x,y,``parameter``,[-2,0,1]);`,
` # If we had not specified the point the result would be bigger`,
``,
`SEE ALSO: genus1`,
` For more examples of curves with genus 0 see the file IntBasis`
):

# f is a polynomial in x and y
# f must be irreducible over Qbar
# If f is reducible this procedure does not work properly.
genus:=proc()
	global  trancendental_ext;
	local f,x,y,d,t,p,v,result,i,point;
	options remember;
	if args[nargs]=`irr check` then
		f:=evala(AFactor(args[1]));
		if not type(f,`+`) then
			i:=2;
			if type(f,`*`) then
				i:=0;
				for t in [op(f)] do if indets(t) intersect
				 {args[2..3]} <> {} then i:=i+1
				fi od;
			fi;
			if i>1 then RETURN(args[1],`factored as`,f) fi;
		fi;
		RETURN(genus(args[1..nargs-1]))
	fi;
	if nargs<3 or (not type(args[2],name)) or (not type(args[3],name))
	or (not type(f,polynom(anything,[args[2..3]])))
	or (nargs>=4 and not member(args[4],{`parameter`,`parametrization`}))
	or (nargs>=4 and args[4]=`parametrization` and (nargs=4 or not
	 type(args[5],name)))
	or indets(args[1]) intersect {args[2..3]}={}
	then
		ERROR(`wrong number or type of arguments`)
	fi;
	f:=expand(args[1]);
	x:=args[2];
	y:=args[3];
	trancendental_ext:=indets(f) minus
		{x,y,seq(`puiseux/rootof`.i,i=0..nops(g_conversion2))};
	if nargs>3 and args[4]=`parametrization` then
		t:=args[5];
		p:=genus(args[1..3],`parameter`,args[6..nargs]);
		if p=`genus is not zero` then RETURN(p) fi;
		RETURN([compute_x(p,f,x,y,t),compute_x(p,f,y,x,t)])
	fi;
	d:=degree(f,{x,y});
	if coeff(f,y,d)=0 then
		if nargs=5 then
			v:=args[5];
			if nops(v)=3 then
				# v is a point
				if degree(f,x)=d then v:=[v[2],v[1],v[3]]
				else
					v:=[v[1]-v[2],v[2..3]]
				fi
			else
				# v[1] is a function with 1 pole in v[2],
				# where v[2]=`finite` or `infinity`
				if degree(f,x)=d then v:=subs({x=y,y=x},v)
				else
				 v:=[evala(Expand(subs(x=x+y,v[1]))),v[2]]
				fi
			fi
		else
			# No point, or function was specified yet.
			v:=NULL
		fi;
		if degree(f,x)=d then RETURN(subs({x=y,y=x},
		 genus(subs({x=y,y=x},f),args[2..min(4,nargs)],v)))
		else RETURN(subs(x=x-y,
		 genus(evala(Expand(subs(x=x+y,f))),args[2..min(4,nargs)],v)))
		fi
	fi;
	# Now coeff(f,y,d)<>0, so [0,1,0] is not a point on the curve anymore
	result:=(d-1)*(d-2)/2;
	v:=singularities(f,x,y);
	for i in v do
		result:=result-i[3]*degree_ext(i,f)
	od;
	if nargs=3 then RETURN(result) fi;
	if result<>0 then RETURN(`genus is not zero`) fi;
	if nargs>4 then
		point:=args[5]
	else
		point:=find_point(f,x,y)
	fi;
	if nops(point)=3 then
		# Now we compute a function that has a pole in point. If the
		# point is finite, this function has no other poles in finite
		# otherwise it has no other poles in infinity.
		# It may still have other poles, they will be removed by
		# the procedure function_with_one_pole
		if point[3]=0 then
			point:=[(homogeneous(evala(Quo(subs({t=0,x=1}
			 ,homogeneous(f,x,y,t,`polynom`)),y-point[2]/point[1],y))
			 ,y,t,x,`polynom`)),`infinity`];
		else
			point:=[evala(Quo(evala(Expand(subs(x=point[1]/point[3],f)))
			 ,y-point[2]/point[3],y))/(x-point[1]/point[3]),`finite`]
		fi
	fi;
	function_with_one_pole(f,x,y,point)
end:

# f is a genus 0 curve in x and y. p is the parameter, i.e. p generates
# the function field
# Output: a rational function in t that gives x as a function of p.
compute_x:=proc(p,f,x,y,t)
	local i,X,np,dp,R,ansatz,C,ar;
	if not member(y,indets(f)) then RETURN(g_solve(f)) fi;
	dp:=denom(p);
	np:=numer(p);
	X[1]:=-1;
	for i from 1 to 3 do
		R[i]:=0;
		while degree(R[i],t)<degree(f,y) do
			X[i]:=X[i]+1;
			ar:=evala(subs(x=X[i],[dp*t-np,f]));
			if member(t,indets(ar)) then
				R[i]:=evala(Resultant(op(ar),y))
			fi
		od;
		X[i+1]:=X[i]
	od;
	ansatz:=(-X[1]*C*R[2]+X[2]*R[1])/(R[1]-C*R[2]);
	ansatz:=g_normal(subs(g_solve({coeffs(collect(evala(Rem(collect(numer(
	 ansatz)-X[3]*denom(ansatz),t),R[3],t)),t),t)}),ansatz));
	if member(C,indets(ansatz)) then bug() fi;
	ansatz
end:
	
# The following procedure computes a parameter. A parameter is a function
# that has only 1 pole, with multiplicity 1. Then k(p)=k(x,y), p generates
# the function field. The pole will be in point. In fact, as the syntax
# is now, point is not a point, but a function with a pole in a certain point.
#
# point=[function,`finite`] or [function,`infinity`]. This function must
# have one pole in finite or infinity. This function is the start for our
# search for a parameter. We add an ansatz to this point. Since this function
# is already good in a part of the plane (either `finite` or `infinity`)
# we will change no poles in that part of the plane, and try to remove the
# poles in the other part.
#
# The point [0,1,0] must not be a point on the curve. The reason is that
# we have divided the projective plane in 2 affine parts, a line and
# a plane. Doing this, there remains 1 point, and therefore that point must
# not be on the curve.
#
# Note that the way we divide the plane in parts is not really relevant
# for the method. The fact that one part is only a line does not give much
# asymmetry in the algorithm. The important thing about these parts is that
# they are affine, and not projective. We need this to apply integral basis.
#
# We will compute an integral basis for the finite part of the plane, and a
# local integral basis for the line `infinity`. Using these integral basis
# we can find linear equations stating that an ansatz, or an ansatz + a
# given function is integral (i.e. no poles) in a part of the plane (the
# part `finite` and `infinity`.
function_with_one_pole:=proc(f,x,y,point)
	local d,j,z,f_infty,ib1,ib2,den1,den2,d1,d2,ansatz,a,f1,f2,equations,
	 ext,i,zero,point1,deg_inf;
	d:=degree(f,{x,y});
	f_infty:=subs(x=1,homogeneous(f,x,y,z,`polynom`));
	ib1:=integral_basis(f,x,y);
	ib2:=integral_basis(f_infty,z,y,{z});
	den1:=g_denom(ib1[d]);
	den2:=g_denom(ib2[d]);
	if point[2]=`infinity genus1` then
		# Syntax for genus1, the function point1 may not be changed
		# so we only determine the upper bound for the denominators now
		point1:=point[1];
		den1:=lcm(den1,denom(point1))
	elif point[2]=`infinity` then
		deg_inf:=degree(den1,x)-degree(point[1],{x,y})+1;
		point1:=x^deg_inf*point[1]/den1;
		if deg_inf<0 then den1:=evala(Expand(den1*x^(-deg_inf))) fi
	else
		point1:=point[1]
	fi;
	# den1 is an upperbound for denominators
	d1:=degree(den1,x);
	d2:=ldegree(den2,z);
	d2:=max(d2,degree(g_numer(point1),{x,y})-degree(g_denom(point1),{x,y}));
	# d2 is an upperbound for the denominators in infinity
	ansatz:=convert([seq(seq(a[i,j]*x^i*y^j,i=0..d1+d2-j)
	 ,j=0..d-1)],`+`)/den1;
	if point[2]=`finite` then
		f1:=ansatz;
		f2:=ansatz+point1
	else
		f1:=ansatz+point1;
		f2:=ansatz
	fi;
	# Now we solve the following set of linear equations:
	# {f1 has no poles in finite, and f2 has no poles in infinity}
	ext:=g_ext([f,f1,f2]);
	f1:=g_normal(f1);
	den1:=g_denom(f1);
	f1:=subs(g_conversion1,g_numer(f1));
	f2:=g_normal(subs(x=1,homogeneous(f2,x,y,z,`ratfunction`)));
	den2:=g_denom(f2);
	den2:=z^ldegree(den2,z);
	f2:=subs(g_conversion1,g_numer(f2));
	ib1:=subs(g_conversion1,[seq(g_expand(g_normal(i*den1),ext),i=ib1)]);
	ib2:=subs(g_conversion1,[seq(g_expand(g_normal(i*den2),ext),i=ib2)]);
	zero:=f1;
	for i from d-1 by -1 to 0 do while degree(coeff(zero,y,i),x)
	 >=degree(lcoeff(ib1[i+1],y),x) and coeff(zero,y,i)<>0 do
		zero:=g_expand(zero-ib1[i+1]*
		g_normal(lcoeff(coeff(zero,y,i),x)/lcoeff(lcoeff(ib1[i+1],y),x))
		*x^(degree(coeff(zero,y,i),x)-degree(lcoeff(ib1[i+1],y),x))
		,ext);
		if trancendental_ext<>{} then
			zero:=expand(collect(zero,[x,y],g_normal))
		fi
	od od;
	equations:={a[d1,0]=0,coeffs(zero,[x,y])};
	# These equations state that f1 is integral. We've added and 1 extra
	# linear condition, to make the solution unique. We still have to
	# add the equations for infinity:
	zero:=f2;
	for i from d-1 by -1 to 0 do while degree(coeff(zero,y,i),z)
	 >=degree(lcoeff(ib2[i+1],y),z) and coeff(zero,y,i)<>0 do
		zero:=g_expand(zero-ib2[i+1]*
		g_normal(lcoeff(coeff(zero,y,i),z)/lcoeff(lcoeff(ib2[i+1],y),z))
		*z^(degree(coeff(zero,y,i),z)-degree(lcoeff(ib2[i+1],y),z))
		,ext);
		if trancendental_ext<>{} then
			zero:=expand(collect(zero,[x,y],g_normal))
		fi
	od od;
	equations:=subs(g_conversion2,{op(equations),coeffs(zero,[y,z])});
	d:=g_normal(subs(g_solve(equations),ansatz+point1));
	if indets(d) minus indets([args]) <> {} then
		ERROR(`can not find parameter`,f=evala(AFactor(f)))
	fi;
	d
end:

`help/text/genus1`:=TEXT(
``,
`FUNCTION:  genus1 - compute an isomorphism for function fields with genus 1`,
``,
`CALLING SEQUENCE: genus1(f,x,y,x0,y0,opt)`,
``,
`PARAMETERS:`,
`   f          - a polynomial in x and y describing the curve of genus 1`,
`   x,y,x0,y0  - variables`,
`   opt        - (optional) a sequence of options`,
``,
`SYNOPSIS:`,
` - An algebraic function field Qbar(x)[y]/(f) where Qbar is the algebraic`,
`   closure of Q is isomorphic to the field Qbar(x0)[y0]/(f0) where f0 is of`,
`   the form y0^2 + square free polynomial in x0 of degree 3`,
`   if and only if the genus is equal to 1.`,
`   This procedure computes such an isomorphism.`,
``,
` - The output of genus is a list of 5 items containing:`,
`    The curve f0`,
`    The image of x0 under this isomorphism`,
`    The image of y0 under this isomorphism`,
`    The image of x under the inverse isomorphism`,
`    The image of y under the inverse isomorphism`,
``,
` - A regular point [x,y,z] on the curve can be specified as a 6'th argument.`,
`   Then x0 and y0 will have their poles in this point on the curve. If the`,
`   option ``no inverse`` is used then the inverse isomorphism is not`,
`   computed. The option ``Weierstrass`` results in a Weierstrass normal form`,
`   i.e. y0^2 - 4*x0^3 - a*x0 - b. If the option ``j invariant`` is specified`,
`   then only the j-invariant is given.`,
``,
`EXAMPLES:`,
`f:=x^4+y^4-2*x^3+x^2*y-2*y^3+x^2-x*y+y^2;`,
`v:=genus1(f,x,y,x0,y0);`,
`        3                      2      - 3 y + 2 x`,
`v := [x0  + 2/3 x0 + 1/108 + y0 , 1/3 -----------,`,
`                                           x`,
``,
`                 2    2    3      3        2`,
`      2 x y - 2 y  - x  + x  + 2 y  - 2 x y`,
`- 1/2 --------------------------------------,`,
`                     2`,
`                    x  (x - 1)`,
``,
`                                         2             3`,
`      - 36 y0 + 54 y0 x0 + 45 x0 - 108 x0  - 52 + 54 x0`,
`- 3/2 --------------------------------------------------,`,
`                               2         3        4`,
`            97 - 96 x0 + 216 x0  - 216 x0  + 81 x0`,
``,
`                                 2                  2         4               3`,
`    72 y0 - 216 y0 x0 + 162 y0 x0  - 246 x0 + 351 x0  + 162 x0  + 104 - 432 x0`,
`1/2 ---------------------------------------------------------------------------`,
`                                         2         3        4`,
`                      97 - 96 x0 + 216 x0  - 216 x0  + 81 x0`,
``,
`]`,
``,
`# check if the image of x and y still satisfy the relation f in`,
`# the field Qbar(x0)[y0]/(f0)`,
`subs({x=v[4],y=v[5]},f):`,
`evala(subs(y0=RootOf(v[1],y0),"));`,
`                                       0`,
``,
`# check if the image of x0 and y0 still satisfy the relation f0 in`,
`# the field Qbar(x)[y]/(f)`,
`subs({x0=v[2],y0=v[3]},v[1]):`,
`evala(subs(y=RootOf(f,y),"));`,
`                                       0`,
``,
`f:=y^5+4/3-23/3*y^2+11*y^3-17/3*y^4-16/3*x^2+16/3*x^3-4/3*x^4;`,
`genus1(f,x,y,x0,y0,``j invariant``);`,
`                                     1404928`,
`                                   - -------`,
`                                       171`,
``,
`f:=-207/82*x^4*y+y^3-64/41*x*y^3+621/82*x^3*y+135/82*x*y+4/41-521/82*x^2`,
`  *y-14/41*y-19/82*y^2-35/82*y^4-4/41*y^5-7/82*x-50/41*x^2+9/2*x^3-225/41`,
`  *x^4+90/41*x^5;`,
`genus1(f,x,y,x0,y0,``j invariant``);`,
`                             34406773155240934694401`,
`                           - -----------------------`,
`                              298818188938894417983`,
``
):

# Input: a polynomial f in x and y with genus 1, the variables x,y,x0 and y0
#  f must be irreducible over Qbar
# Optional 6'th argument: a regular point on the curve
# If the last argument is `no inverse` then the inverse of the isomorphism
#  is not computed (i.e. the image of x and y is not computed).
# If [0,1,0] is a point on the curve then:
#   - with the option `no inverse` the output is not necessarily
#     a polynomial in y.
#   - without the option `no inverse` it does become a polynomial, but this is
#     rather slow.
# Output: a list containing:
#  f0 such that Qbar(x)[y]/(f) is isomorphic with Qbar(x0)[y0]/(f0)
#  The image of x0 under this isomorphism
#  The image of y0
#  The image of x under the inverse isomorphism
#  The image of y
# f0 is of the form y0^2 + polynomial of degree 3 in x0^2
genus1:=proc(ff,x,y,x0,y0)
	global  trancendental_ext;
	local a,d,t,v,result,i,point,f,f0,x0v,y0v,ansatz,point2,point3;
	options remember;
	f:=expand(ff);
if not has(`no inverse`,{args}) then
	v:=genus1(args,`no inverse`);
	if v=`genus is not one` then RETURN(v) fi;
	f0:=v[1];
	if has(`j invariant`,{args}) then
		f0:=subs(y0=0,f0);
		f0:=f0/lcoeff(f0,x0);
		a:=coeff(f0,x0,1);
		d:=coeff(f0,x0,0);
		RETURN(evala(6912*a^3/(27*d^2+4*a^3)))
	fi;
	# Remove possible y's in the denominators (could be a costly operation)
	v:=subs(RootOf(f,y)=y,evala(subs(y=RootOf(f,y),[v[2..3]])));
	x0v:=v[1];
	y0v:=v[2];
	RETURN([f0,op(v),`genus1/compute_x`(x0v,y0v,x,y,f,f0,x0,y0),
	`genus1/compute_x`(x0v,y0v,y,x,f,f0,x0,y0)])
fi;
	# This set contains possible parameters appearing in f
	trancendental_ext:=indets(f) minus
		{x,y,seq(`puiseux/rootof`.i,i=0..nops(g_conversion2))};
	d:=degree(f,{x,y});
	if coeff(f,y,d)=0 then
		# [0,1,0] is a point on the curve, we want to avoid this
		if nargs>5 and type(args[6],list) then
			v:=args[6]; # user specified point on the curve
			v:=[v[1]-v[2],v[2..3]],args[7..nargs]
		else v:=args[6..nargs] # no point was specified
		fi;
		# apply transformation on f to remove the point [0,1,0]
		v:=genus1(evala(Expand(subs(x=x+y,f))),args[2..5],v);
		if nops(v)>1 then
			# Do the reverse transformation
			i:=subs(x=x-y,[v[2..3]]);
			if nops(v)>3 then
				v:=[v[1],op(i),v[4]+v[5],v[5]]
			else
				v:=[v[1],op(i)]
			fi
		fi;
		RETURN(v)
	fi;
	# Now coeff(f,y,d)<>0, so [0,1,0] is not a point on the curve
	result:=(d-1)*(d-2)/2;
	v:=singularities(f,x,y);
	# Compute the genus
	for i in v do
		result:=result-i[3]*degree_ext(i,f)
	od;
	if result<>1 then RETURN(`genus is not one`) fi;
	if nargs>5 and type(args[6],list) then
		point:=args[6]
	else
		point:=find_point(f,x,y) # Compute a regular point on the curve
	fi;
	# Now compute a function with a pole in that point:
	if point[3]=0 then # point in infinity
		point:=[(homogeneous(evala(Quo(subs({t=0,x=1}
		 ,homogeneous(f,x,y,t,`polynom`)),y-point[2]/point[1],y))
		 ,y,t,x,`polynom`))/x^(degree(f,y)-2),`infinity genus1`]
	else
		point:=[evala(Quo(evala(Expand(subs(x=point[1]/point[3],f)))
		 ,y-point[2]/point[3],y))/(x-point[1]/point[3]),`finite`]
	fi;
	# Now make it a double pole with an indeterminate as residue
	point2:=[evala(Normal(rem(expand(-point[1]^2+`genus/residue`
	 *point[1]),f,y))),point[2]];
	# Compute a function with this pole and no other poles:
	x0v:=function_with_one_pole(f,x,y,point2);
	point3:=[evala(Normal(rem(expand(-point[1]*x0v+`genus/residue`
	*point[1]),f,y))),point[2]];
	# Now compute a function with pole order 3:
	y0v:=function_with_one_pole(f,x,y,point3);
	# look for a relation f0 between x0v and y0v using an ansatz:
	i:=0;
	ansatz:=y0^2+x0^3+a[1]*y0+a[2]*x0*y0+a[3]+a[4]*x0+a[5]*x0^2;
	while has(ansatz,{a[1],a[2],a[3],a[4],a[5]}) do
		i:=i+1;
		# avoid certain values i for x
		while g_normal(subs(x=i,denom(x0v)*denom(y0v)))=0 do
			i:=i+1
		od;
		# find linear equations in a[1] .. a[5]
		ansatz:=subs(g_solve({coeffs(expand(rem(expand(
		 subs(x0=x0v,y0=y0v,x=i,ansatz)),expand(subs(x=i,f)),y)),y)})
		 ,ansatz)
	od; # repeat until all a[i] are determined (usually in the 1'st step)
	f0:=ansatz;
	# Now normalize f0 a little further:
	y0v:=g_normal(y0v+subs(x0=x0v,coeff(f0,y0,1))/2);
	# clear the coefficient of y0^1:
	f0:=collect(subs(y0=y0-coeff(f0,y0,1)/2,f0),
	 [x0,y0],distributed,g_normal);
	x0v:=g_normal(x0v+coeff(f0,x0,2)/3);
	# clear the coefficient of x0^2:
	f0:=collect(subs(x0=x0-coeff(f0,x0,2)/3,f0),
	 [x0,y0],distributed,g_normal);
	if has(`Weierstrass`,{args}) then # Weierstrass normal form
		x0v:=-x0v;
		y0v:=2*y0v;
		f0:=y0^2+4*subs(x0=-x0,y0=0,f0)
	fi;
	[f0,x0v,y0v]
end:

# Compute the characteristic polynomial of x over Qbar(x0v)
# f is a polynomial in x and y giving the algebraic relation between x and y
# x0v is an element of Qbar(x,y)
# We use the variable t for x0v.
`genus1/charpol`:=proc(x0v,x,y,f,t)
	local r,ansatz,a,n,j,k,i,aa;
	n:=degree(f,y);
#	Not necessarily	true: a[n,2]:=1;
	ansatz:=convert([seq(t^k*(a[k,2]*x^2+a[k,1]*x+a[k,0])
	 ,k=0..n)],`+`);
	aa:=indets(ansatz);
	j:=-1;
	i:=0;
while i<5 do
	j:=j+1;
	if evala(Gcd(op(subs(x=j,[denom(x0v),f]))))<>1 then next fi;
	r:=numer(t-subs(x=j,x0v));
	r:=collect(evala(Resultant(r,subs(x=j,f),y)),t,g_normal);
	r:=collect(r/lcoeff(r,t),t,g_normal);
	if degree(r,t)<>n then next fi;
	ansatz:=subs(solve({seq(subs(x=j,coeff(ansatz,t,k)-
	 coeff(r,t,k)*coeff(ansatz,t,n)),k=0..n-1)},aa),ansatz);
	i:=i+1
od;
	r:=expand(numer(g_normal(ansatz/lcoeff(ansatz,t))));
	collect(r/lcoeff(r,x),x,g_normal)
end:

# Express x as an expression in x0 and y0.
# f0 = algebraic relation between x0 and y0
# f = algebraic relation between x and y
# x0v = image of x0 in Qbar(x,y)
# y0v = image of x0 in Qbar(x,y)
`genus1/compute_x`:=proc(x0v,y0v,x,y,f,f0,x0,y0)
	local cp,result,v,i;
	cp:=`genus1/charpol`(x0v,x,y,f,x0);
	result:=-coeff(cp,x,1)/2;
	cp:=g_normal((x^2-subs(x=x+result,cp))/coeff(f0,y0,0));
if cp<>0 then
	# Now add y0*sqrt(cp) to result:
	v:=[seq(evala(Sqrfree(i))[2],i=[numer(cp),denom(cp)])];
	v:=g_normal(convert([seq(i[1]^(i[2]/2),i=v[1])],`*`)/
		    convert([seq(i[1]^(i[2]/2),i=v[2])],`*`));
	cp:=factors(numer(x^2+g_normal(cp/v^2))
	 ,indets([args],RootOf))[2];
	if has(cp,x0) then bug() fi;
	for i in cp do if has(i,x) and degree(i[1],x)=1 then
		v:=v*solve(i[1],x)*y0;
		result:=result+v,result-v;
		break
	fi od;
fi;
	i:=0;
	result:={result};
	while nops(result)>1 do for v in result do
		i:=i+1;
		# avoid certain values i for x
		while g_normal(subs(x=i,denom(x0v)*denom(y0v)))=0 do
			i:=i+1
		od;
		if g_normal(rem(expand(numer(subs(x0=x0v,y0=y0v,x=i,numer(v))
		 -i*subs(x0=x0v,y0=y0v,x=i,denom(v)))),subs(x=i,f),y))<>0 then
			result:=result minus {v};
			break
		fi
	od od;
	op(result)
end:

# Input: f is a polynomial in x and y
# Output is a list of the following lists:
# [ [x,y,z], multiplicity, contribution in the genus ]
singularities:=proc()
	local df,ext,extl,j, nulp,pl_transl,f_translated
	 ,f,x,y,z,n,k,ff,i,mult,contr,result;
	options remember;
	f:=args[1];
	x:=args[2];
	y:=args[3];
	if nargs=3 or (nargs=4 and args[4]=`points`) then
		ff:=[op(singularities(f,x,y,`finite`)),
		 op(singularities(f,x,y,`infinity`))];
		result:=NULL;
		for i in ff do if (nargs=3 and i[2]>1) or (nargs=4 and i[2]=1)
			then result:=result,i
		fi od;
		RETURN([result])
	fi;
n:=degree(f,{x,y});
if args[4]=`infinity` then
	k:=-1;
	while k=-1 or degree(subs(z=0,ff),x)<n do
		# Making sure that we may assume y=1
		k:=k+1;
		ff:=homogeneous(subs(y=y+k*x,f),x,y,z,`polynom`)
	od;
	RETURN([seq([[i[1][2],1+k*i[1][2],0],i[2],i[3]],
	 i=singularities(subs(y=1,ff),z,x,`around zero`,args[5..nargs]))])
fi;
	extl:=g_ext(f);
	f:=subs(g_conversion1,f);
	if args[4]=`around zero` then
		df:={x}
	else
		df:=double_factors(g_expand(discrim(f,y),extl),x,extl)
	fi;
	result:=NULL;
	for k in df do
		nulp:=g_zero_of(k,x,'ext'); # is the x value
		ext:=eval(ext);
		ext:=[ext,op(extl)];
		f_translated:=g_expand(subs(x=x+nulp,f),ext);
		pl_transl:=[op(`puiseux/technical_answer`(f_translated,x,y
		 ,0,ext))];
		ff:={seq(subs(x=0,i[1]),i=pl_transl)}; # the y values
		for i in ff do # i is y value
			mult:=0;
			contr:=0;
			for j in pl_transl do if subs(x=0,j[1])=i then
				mult:=mult+degree_ext(j[1],[i,op(ext)])*min(
				 max(ldegree(j[1]-subs(x=0,j[1]),x),1)/j[3],1);
				contr:=contr+degree_ext(j[1],[i,op(ext)])
				 *(j[7]-(j[3]-1)/j[3])/2
			fi od;
			if not (type(mult,integer) and type(contr,integer))
				then ERROR(`something is wrong`)
			fi;
			result:=result,[[nulp,i,1],mult,contr]
		od
	od;
	subs(g_conversion2,[result])
end:

# f is a polynomial or a rational function in x and y.
# Output: the corresponding homogeneous rational function or polynomial
homogeneous:=proc(ff,x,y,z,what)
local i,j,n,f;
	f:=expand(g_normal(ff));
	if what=`ratfunction` then
		i:=g_numer(f);
		j:=g_denom(f);
		n:=degree(i,[x,y,z])-degree(j,[x,y,z]);
		if n<0 then i:=i*z^(-n) elif n>0 then j:=j*z^n fi;
		RETURN(homogeneous(i,x,y,z,`polynom`)
		 /homogeneous(j,x,y,z,`polynom`))
	fi;
	n:=degree(f,{x,y});
	convert([seq(seq(coeff(coeff(f,x,i),y,j)*x^i*y^j*z^(n-i-j)
	,i=0..n-j),j=0..n)],`+`)
end:

# Searches a point in a small algebraic extension using the information
# gathered during the Puiseux series computation.
find_point:=proc(f,x,y)
local vp,i,m,j;
	vp:=singularities(f,x,y,`points`);
	m:=min(degree(f,{x,y}),seq(degree_ext(i,f),i=vp));
	if m=degree(f,{x,y}) and m>2 then
		m:=singularities(f,x,y);
		if nops(m)=1 and m[1][2]=degree(f,{x,y})-1 then
		# This is not a regular point, but it has mult n-1
		# and also a tangent line in the right direction
			m:=m[1][1];
			if m[3]=0 then RETURN([y-x*m[2]/m[1],`infinity`])
			else RETURN([(y-m[2])/(x-m[1]),`finite`])
			fi
		# a function with 1 pole is also allowed instead of a point
		fi;
		vp:={seq(i[1][1],i=m)};
		i:=0;
		while member(i,vp) do i:=i+1 od;
		vp:=readlib(factors)(subs(x=i,f));
		m:=min(seq(degree(j[1],y),j=vp[2]));
		for j in vp[2] do if degree(j[1],y)=m then
			RETURN([i,RootOf(j[1],y),1])
		fi od
	fi;
	for i in vp do if degree_ext(i,f)=m then RETURN(i[1]) fi od;
	bug()
end:

#########################
#   Singularity knots   #
#########################

# A large part of the code in this section written by Harm Derksen.

lprint(`For plotting type: ?plot_knot`);
macro(  element_of=`plot_knot/element_of`,
	equal_to_zero=`plot_knot/equal_to_zero`,
	all_substitutions=`plot_knot/all_substitutions`,
	curve=`plot_knot/curve`
):

# Test if subs(x^(1/d) = any d-th root of unity * x^(1/d),el) is
# an element of verz.
# x stands for x^(1/d) here
element_of:=proc(el,d,verz,x) local i,elz,j;
	# Remove the branches that don't go through 0:
	if coeff(el,x,0)<>0 then RETURN(true) fi;
	if d=1 then RETURN(false) fi;
	for i from 1 to d-1 do
		elz:=evala(evalf(subs(x=evalc(exp(i*2.0*Pi*I/d))*x,el)));
		for j in verz do if equal_to_zero(j[1]-elz,x) then
			RETURN(true)
		fi od
	od;
	false
end:

equal_to_zero:=proc(a,x) local som,i;
	som:=0.0;
	for i in {coeffs(a,x)} do som:=evalf(som+abs(i)^2) od;
	if som<0.000001 then RETURN(true) fi;
	false
end:

# All possible floating point substitutions for algebraic numbers.
# v: a list of RootOf's
all_substitutions:=proc(v) local i,a,b,result,j;
	if v=[] then RETURN({[]}) fi;
	a:=v[nops(v)];
	result:={seq([a=i],i=[fsolve(op(a),_Z,complex)])};
	b:=[v[1..nops(v)-1]];
	{seq(seq([op(j),op(i)],j=all_substitutions(subs(i,b))),i=result)}
end:

# f is a polynomial in x en y
plot_knot:=proc(f,x,y)
	local ext_f,pui,ext_pui,i,pui2,curven,true_or_not,t,opt,eps,cols;
	opt:={args[4..nargs]};
	if not member(numpoints,indets(opt)) then opt:=opt union {numpoints=150} fi;
	if not member(radius,indets(opt)) then opt:=opt union {radius=0.05} fi;
	if not member(tubepoints,indets(opt)) then opt:=opt union {tubepoints=5} fi;
	if not member(scaling,indets(opt)) then opt:=opt union {scaling=CONSTRAINED} fi;
	if not member(style,indets(opt)) then opt:=opt union {style=PATCH} fi;
	eps:=1;
	if member(epsilon,indets(opt)) then
		eps:=subs(opt,epsilon);
		opt:=opt minus {epsilon=eps}
	fi;
	cols:=[];
	if member(colours,indets(opt)) then
		cols:=subs(opt,colours);
		opt:=opt minus {colours=cols}
	fi;
	ext_f:=g_ext(f);
	pui:=`puiseux/technical_answer`(g_expand(subs(g_conversion1,f),ext_f)
	 ,x,y,0,ext_f);
	ext_f:=subs(g_conversion2,ext_f);
	ext_pui:=subs(g_conversion2,g_ext(pui));
	ext_pui:=[ext_pui[1..nops(ext_pui)-nops(ext_f)]];
	ext_f:=all_substitutions(ext_f);
	ext_f:=ext_f[1];
	ext_pui:=all_substitutions(subs(ext_f,ext_pui));
	pui:=subs(ext_f,subs(g_conversion2,pui));
	# Now replace the RootOf's by all corresponding floating point numbers
	pui2:={};
	for i in ext_pui do
		pui2:=pui2 union subs(i,pui)
	od;
	pui:=pui2;
	pui2:={};
	for i in pui do
		true_or_not:=element_of(i[1],i[3],pui2,x);
		if (not true_or_not) then pui2:=pui2 union {i} fi
	od;
	curven:={};
	for i in pui2 do curven:=curven union {curve(i[1],i[3],t,eps,x)} od;
	if printlevel>1 then print(`Number of branches:`,nops(curven)) fi;
	if cols=[] then
		plots[tubeplot](curven,t=0..2*Pi,op(opt))
	else
		curven:=[op(curven)];
		plots[display]([seq(plots[tubeplot](curven[i],t=0..2*Pi,op(opt),
		 color=cols[1+irem(i,nops(cols))]),i=1..nops(curven))])
	fi
end:

curve:=proc(functie_in_x,d,t,eps,x) local f,fre,fim,factor;
	f:=expand(subs(x=eps*(cos(t)+I*sin(t)),functie_in_x));
	fre:=coeff(f,I,0)/eps;fim:=coeff(f,I,1)/eps;
	factor:=1/(sqrt(1+fre^2+fim^2)-fre);
	[factor*cos(d*t),factor*sin(d*t),factor*fim];
end:

`help/text/plot_knot`:=TEXT(
``,
`FUNCTION:  plot_knot - make a tubeplot for a singularity knot`,
``,
`CALLING SEQUENCE: plot_knot(f,x,y,opt)`,
``,
`PARAMETERS:`,
`   f   - a polynomial in x and y with a singularity in x=0, y=0`,
`   x,y - variables`,
`   opt - (optional) a sequence of options`,
``,
`SYNOPSIS:`,
` - The polynomial f represents a curve in the complex plane C^2. By identifying`,
`   C^2 with R^4 the curve can be viewed as a 2 dimensional surface over the`,
`   real numbers. This procedure computes the intersection of this surface with`,
`   a sphere in R^4 with radius epsilon and center 0. The intersection consists`,
`   of a number of closed curves over the real numbers. After applying a`,
`   projection from the sphere (is 3 dimensional over R) to R^3 these curves can`,
`   be plotted by the tubeplot command in the plots package. Such a plot gives`,
`   information about the singularity of f in the point 0. See also:`,
`   E. Brieskorn, H. Kn"orrer: Ebene Algebraische Kurven, Birkh"auser 1981.`,
``,
` - If printlevel > 1 the number of curves will be printed to the screen.`,
``,
`OPTIONS:`,
`   epsilon:    the radius of the sphere. Default is set to 1. If there are`,
`               other singularities near 0 then a smaller number must be chosen.`,
`   colours:    Specifying a list of colours results in a plot where each`,
`               branch gets its own colour.`,
``,
`   Also the options for tubeplot can be used. The defaults for numpoints,`, 
`   radius, tubepoints, scaling and style have been set to the values 150,`, 
`   0.05, 5, CONSTRAINED and PATCH.`,
``,
`SEE ALSO: tubeplot`,
``,
`EXAMPLE:`,
` printlevel:=2;`,
` plot_knot(y^2-x^3,x,y);`,
` f:=(y^3-x^7)*(y^2-2*x^5);`,
` plot_knot(f,x,y,epsilon=0.8,radius=0.03,colours=[blue,red]);`,
` plot_knot(f+y^3,x,y,epsilon=0.8);`,
` plot_knot(f+y^3,y,x,epsilon=0.8); # Not exactly the same because the picture`,
`                                   # depends on the projection point`,
` f:=(y^3-x^7)*(y^2-2*x^5)*(y^2+2*x^5);`,
` plot_knot(f,y,x,epsilon=0.8,radius=0.03,colours=[blue,red,pink]);`,
` g:=(y^3-x^7)*(y^3-x^7+100*x^13)*(y^3-x^7-100*x^13);`,
`plot_knot(g,x,y,epsilon=0.8,radius=0.03,numpoints=250,colours=[blue,red,green]);`,
``
):

# save `IntBasis.m`;
# quit