Knots := module()
	export KnotDiagram, AlexanderPolynomial, DowkerNotation;
	uses GraphTheory, LinearAlgebra, ArrayTools;
	local KnotCoordinates, orient, inter;
	option package;
 
	orient := (p1, p2, p3) -> signum(Determinant(Matrix(3, 3, [p1, p2, p3], 'fill' = 1)));
 
  inter := proc(p, q)
		local u := p[2] - p[1], v := q[2] - q[1];
		(p[1] - q[1]).(v.v*u - u.v*v)/((u.v)^2 - (u.u)*(v.v))*u + p[1]
	end proc;

	KnotCoordinates := proc(dow)
		local succ, ndp := nops(dow), edges, emb, v1, v2;
		succ := i -> `if`(i = 2*ndp, 1, i + 1);
		
		edges := `union`(seq(
			proc()
				v1, v2 := 2*i - 1, abs(dow[i]);
				{{v1, -succ(v1)}, {v2, -succ(v2)},
				 {-v1, -v2}, {-v1, v2}, {v1, -v2}, {v1, v2}}
			end proc(),
			i = 1..ndp));
		emb := DrawPlanar(Graph([seq(op([-i, i]), i = 1..2*ndp)], edges), 'internal');
		Transpose(Reshape(Array(map(Vector, emb)), 2, 2*ndp))
	end proc;

	KnotDiagram := proc(knot)
		local succ, dow, ndp, vcoord, edges, dpcoord, p1, p2, q1, q2, v1, v2;
		
		if op(0, knot) = 'TorusKnot' or op(0, knot) = 'PretzelKnot' then
			return `if`([op(knot)]::list(integer),
				KnotDiagram(DowkerNotation(knot)),
				'KnotDiagram(knot)')
		end if;

		if not knot::list then
			return 'KnotDiagram(knot)'
		end if;
		
		dow := knot; ndp := nops(dow); vcoord := KnotCoordinates(dow);
		succ := i -> `if`(i = 2*ndp, 1, i + 1);
		
		edges := `union`(seq(
			proc()
				v1, v2 := 2*i - 1, abs(dow[i]);
				p1, q1, p2, q2 := op(convert(
					vcoord[`if`(dow[i] > 0, [v1, v2], [v2, v1])], list));
				dpcoord := inter([p1, p2], [q1, q2]);
				{[vcoord[v1, 2], vcoord[succ(v1), 1]],
 				 [vcoord[v2, 2], vcoord[succ(v2), 1]],
				 [p1, p2],
				 [q1, max(1 - .025/Norm(dpcoord - q1), 0)*(dpcoord - q1) + q1],
				 [min(.025/Norm(q2 - dpcoord), 1)*(q2 - dpcoord) + dpcoord, q2]}
			end proc(),
			i = 1..ndp));
			
		PLOT(op(map(vv -> POLYGONS(map(v -> convert(v, list), vv),
						STYLE(LINE), COLOR(RGB, 0, 0, 1)),
					edges)),
			SCALING(CONSTRAINED), AXESSTYLE(NONE))
	end proc;

	AlexanderPolynomial := proc(knot, t)
		local dow, ndp, vcoord, poly, mat, arcs, dps, s, cur, v1, v2, p, q, r, i;

		if op(0, knot) = 'TorusKnot' then
			p, q := op(knot);
			return `if`([p, q]::list(integer),
				quo((t^(p*q) - 1)*(t - 1), (t^p - 1)*(t^q - 1), t),
				(t^(p*q) - 1)*(t - 1)/(t^p - 1)/(t^q - 1))
		end if;
		
		if op(0, knot) = 'PretzelKnot' then
			p, q, r := op(knot);
			poly := signum(p*q + p*r + q*r + 1)*
				((p*q + p*r + q*r)*(t^2 - 2*t + 1) + (t^2 + 2*t + 1))/4;
			return `if`([p, q, r]::list(integer),
				`if`([p, q, r]::list(odd),
					`if`(p*q + p*r + q*r <> -1, poly, 1),
					AlexanderPolynomial(DowkerNotation(knot), t)),
				piecewise(
					`and`(p::odd, q::odd, r::odd), piecewise(p*q + p*r + q*r <> -1, poly, 1),
					'AlexanderPolynomial(knot, t)'))			
		end if;
		
		if not knot::list then
			return 'AlexanderPolynomial(knot, t)'
		end if;
		
		dow := knot; ndp := nops(dow); vcoord := KnotCoordinates(dow);
		dps := Array(1..2*ndp);
		for i from 1 to ndp do
			dps[[2*i - 1, abs(dow[i])]] := Array(`if`(dow[i] > 0, [1, -1], [-1, 1]))
		end do;
		
		cur := 1;
		arcs := [seq(
			proc()
				if dps[i] > 0 then [cur, cur]
					elif cur = ndp
						then cur := 1; [ndp, cur]
					else cur := cur + 1; [cur - 1, cur]
				end if
			end proc(),
			i = 1..2*ndp)];

		mat := Matrix(ndp, ndp);
		for i from 1 to ndp do
			v1, v2 := 2*i - 1, abs(dow[i]);
			if dow[i] < 0 then v1, v2 := v2, v1 end if;
			s := orient(vcoord[v1, 1], vcoord[v1, 2], vcoord[v2, 2]);
			zip(proc(j, e) mat[i, j] := mat[i, j] + e end proc,
				[arcs[v1, 1], arcs[v2, 1], arcs[v2, 2]],
				`if`(s > 0, [-t + 1, -1, t], [-t + 1, t, -1]))
		end do;

		poly := Determinant(mat[1..-2, 1..-2]);
		expand(signum(lcoeff(poly, t))*poly/t^ldegree(poly, t))
	end proc;
	
	DowkerNotation := proc(knot)
		local dow, dps, sp, sq, sr, p, q, r, s;
		
		if op(0, knot) = 'TorusKnot' then
			p, q := op(knot);
			if not [p, q]::list(integer) then
				return 'DowkerNotation(knot)'
			end if;
			
			dow := [seq(seq([i*q - j*p, i*q + j*p], j = 1 .. q - 1), i = 0 .. p - 1)];
			dow := dow mod 2*p*q;
			dow := subs(zip(`=`, sort(map(op, dow)), [seq(1..2*p*(q - 1))]), dow);
			dps := [seq(seq(
				(2*(dow[i*(q - 1) + j, 1] mod 2) - 1)*(2*(i mod 2) - 1)*signum(j*p mod 2*q - q),
				j = 1 .. q - 1), i = 0 .. p - 1)];
			dow := zip((dp, s) -> [1, s]*~`if`(dp[1]::odd, dp, dp[[2, 1]]), dow, dps);
			return sort(dow)[1..-1, 2]
		end if;
		
		if op(0, knot) = 'PretzelKnot' then
			sp, sq, sr := op(knot);
			if not [sp, sq, sr]::list(integer) then
				return 'DowkerNotation(knot)'
			end if;
			
			p, q, r := op(abs~([sp, sq, sr]));
			s := p + q + r;			
			if [p, q, r]::list(odd) then
				return [signum(sp)*seq(s + p .. s + 1, -2),
								signum(sq)*seq(2*s - r - 1 .. s + p + 2, -2),
								signum(sr)*seq(2*s .. 2*s - r + 1, -2),
								signum(sp)*seq(p - 1 .. 2, -2),
								signum(sq)*seq(s - r .. p + 1, -2),
								signum(sr)*seq(s - 1 .. s - r + 2, -2)]
			end if;
			
			if q::even then
				return DowkerNotation(PretzelKnot(sq, sr, sp))
			elif r::even then
				return DowkerNotation(PretzelKnot(sr, sp, sq))
			end if;
			return [signum(sp)*seq(s + p .. s + 2, -2),
							signum(sq)*seq(2*s - q + 1 .. 2*s, 2),
							signum(sr)*seq(s + p + 2 .. 2*s - q - 1, 2),
							signum(sp)*seq(p .. 2, -2),
							signum(sr)*seq(s - r + 1 .. s, 2),
							signum(sq)*seq(p + 2 .. s - r - 1, 2)]
		end if;
		
		`if`(knot::list, knot, 'DowkerNotation(knot)')
	end proc
end module;