# The 17 functions:
# Latest version: 1997-10-27
# (c) copyright 1997, 1998, 2001 Leroy J. Dickey
#

with(linalg) ; 


# procedure group 1: cross_point

cross_point2 	:=	proc(w,x,  y,z) ; 
	# find the coordinates of the point
	# where line WX meets line YZ. 
	crossprod ( w, x ) ; 
	crossprod ( y, z ) ; 
	crossprod ( " , "" ) ; 
	map ( normal , " ) ; 
	end:

cross_point 	:=	proc(w,x,  y,z) ; 
	cross_point2 (w,x, y,z);
	tidyS ( " ) ;
	end:


# procedure group 2: Pappus, gcd3

gcd3		:=	proc ( v ) ; 
	# find the gcd of the first three entries in a vector
	gcd ( v[1] , gcd ( v[2], v[3] ) ) ; 
	end:

Pappus		:=	proc ( a, b, c,  d, e, f ) ; 
# To find the Pappus line of the hexagon a e c d b f.
#		A	B	C
#		D	E	F
	    cross_point ( a,e,  b,d ) ; 
	    cross_point ( b,f,  c,e ) ; 
	    crossprod ( ", "" ) ; 
	    tidyN ( " ) ; 
	    tidyS ( " ) ; 
	end:


# procedure group 3: Steiner_point, point_on_line

Steiner_point	:=	proc ( a,b,c,  d,e,f ) ; 
	Pappus ( a, b, c, d, e, f ) ; 
	Pappus ( a, b, c, e, f, d ) ; 
	crossprod ( " , "" ) ; 
	map ( normal , " ) ; 
	tidyS ( " ) ; 
	end:

point_on_line 	:=	proc ( point, line ) ; 
	# test to see if the point is on the line.
	dotprod ( point , line ) ; 
	normal ( " ) ; 
	if ( " = 0 ) then
		RETURN ( `TRUE` ) ; 
	else
		RETURN ( `FALSE` ) ; 
	fi ; 
	end:
 

# procedure group 4: conic, point_on_conic, matform
 
conic 		:= 	proc ( P, Q, R, S, T ) 
	local a, b, c, f, g, h, f1, f2, f3, f4, f5, vv, ww ; 
	unassign ( 'a', 'b', 'c', 'f', 'g', 'h' ) ; 
  f1:=a*P[1]^2+b*P[2]^2+c*P[3]^2+2*f*P[2]*P[3]+2*g*P[3]*P[1]+2*h*P[1]*P[2] ; 
  f2:=a*Q[1]^2+b*Q[2]^2+c*Q[3]^2+2*f*Q[2]*Q[3]+2*g*Q[3]*Q[1]+2*h*Q[1]*Q[2] ; 
  f3:=a*R[1]^2+b*R[2]^2+c*R[3]^2+2*f*R[2]*R[3]+2*g*R[3]*R[1]+2*h*R[1]*R[2] ; 
  f4:=a*S[1]^2+b*S[2]^2+c*S[3]^2+2*f*S[2]*S[3]+2*g*S[3]*S[1]+2*h*S[1]*S[2] ; 
  f5:=a*T[1]^2+b*T[2]^2+c*T[3]^2+2*f*T[2]*T[3]+2*g*T[3]*T[1]+2*h*T[1]*T[2] ; 
	solve ( { f1, f2, f3, f4, f5 }, { a, b, c, f, g, h } ) ; 
	assign ( " ) ;  
	vv := [ a , b, c, f, g, h ] ;  
	map ( denom , " ) ; 
	lcm ( "[1] , "[2], "[3], "[4], "[5], " [6] ) ; 
	scalarmul ( vv , " ) ; 
	ww := map ( normal , " ) ; 
	gcd3 ( [ ww[1], ww[2], ww[3] ] ) ; 
	gcd3 ( [ ww[4], ww[5], ww[6] ] ) ; 
	gcd ( "", " ) ; 
	scalarmul ( ww , 1/ " ) ; 
	map ( normal , " ) ; 
	matform ( " ) ;
	end:

point_on_conic 	:=	proc ( point, matrix ) 
	local zz ; 
	# test is the point is on the conic
	dotprod ( point , multiply ( matrix, point ) ) ; 
	zz := simplify ( " ) ; 
	if ( " = 0 ) then
		RETURN ( `TRUE` ) ; 
	else
		RETURN ( zz ) ; 
	fi ; 
	end:

matform		:=	proc ( v ) ; 
	# given vector [ a, b, c, f, g, h ] ,
	# produce matrix [ a, h, g ] , [ h, b, f ] , [g , f, c ]
 	matrix ( 3, 3, [v[1], v[6], v[5],
		   	v[6], v[2], v[4],
		   	v[5], v[4], v[3] ] ) ; 
	end:


# procedure group 5: Steiner_conic, new_conic. 

Steiner_conic 	:=	proc ( A, B, C, D, E ) ; 
	# the "even" Steiner conic 
	#	A  	B	C
	#	    D       E    
	conic ( D, E, 	cross_point ( A, D,  C, E ) , 
			cross_point ( B, D,  A, E ) ,
			cross_point ( C, D,  B, E ) ) ; 
	RETURN ( " ) ; 
	end:

new_conic	:= 	proc ( A, B, C,   D, E ) 
	local p1, p2, p3, p4, p5, p6 ; 
	p1 := cross_point ( A, D,  B, E ) ; 
	p2 := cross_point ( B, E,  C, D ) ; 
	p3 := cross_point ( C, D,  A, E ) ; 
	p4 := cross_point ( A, E,  B, D ) ; 
	p5 := cross_point ( B, D,  C, E ) ; 
	# p6 := cross_point ( C, E,  A, D ) ; 
	# Take first 5
	conic ( p1, p2, p3, p4, p5 ) ; 
	RETURN ( " ) ; 
	end:


# procedure group 6: tidy, line_meet_conic

tidyD		:=	proc (v) ; 
	# clean up the denominator
	map ( denom, v ) ; 
	lcm ( "[1], "[2], "[3] ) ; 
	scalarmul ( v, " ) ; 
	end:

tidyN		:=	proc (v) ; 
	# clean up the numerator 
	gcd3 ( v ) ; 
	scalarmul ( v, 1/" ) ; 
	map ( normal, " ) ; 
	end:

tidyS		:=	proc (v) ; 
	# clean up the sign of the first non-zero coordinate
	scalarmul ( v, sign ( v [3] ) ) ; 
	scalarmul ( ", sign ( " [2] ) ) ; 
	scalarmul ( ", sign ( " [1] ) ) ; 
	end:

tidy		:=	proc (v) ; 
	tidyD (v) ; 
	tidyN (") ; 
	tidyS (") ; 
	end:

line_meet_conic	:= 	proc ( p,q, M ) 
	# test:  find the points where the line meets the conic
	local r, r1, r2, u, u1, u2, fun, s ; 
	
	dotprod ( p, multiply ( M, p ) ) ; 
	if  ( " = 0 ) then
		dotprod ( q, multiply ( M, q ) ) ; 
		if ( " = 0 ) then
			# two solutions are P and Q 
			matrix ( 2, 3 , [ tidyD(p), tidyD(q) ] ) ; 
			RETURN ( " ) ; 
		else
			# use the form  P + x Q 
	  	        r:=[ p[1]+u*q[1],p[2]+u*q[2],p[3]+u*q[3] ] ; 
			fun := dotprod ( r, multiply ( M, r ) ) ; 
			s := [ solve ( fun ,  u ) ] ; 
			u1 := s [1] ; 
			u2 := s [2] ; 
	  	        r1:=[p[1]+u1*q[1],p[2]+u1*q[2],p[3]+u1*q[3]] ; 
	  	        r2:=[p[1]+u2*q[1],p[2]+u2*q[2],p[3]+u2*q[3]] ; 
			r1 := tidyD ( r1 ) ; 
			r2 := tidyD ( r2 ) ; 
			matrix ( 2, 3 , [ r1, r2 ] ) ; 
			RETURN ( " ) ; 
		fi ; 
	else
		# use the form  Q+xP 
		r:=[u*p[1]+q[1], u*p[2]+q[2], u*p[3]+q[3] ] ; 
		fun := dotprod ( r, multiply ( M, r ) ) ; 
		s := [ solve ( fun ,  u ) ] ; 
		u1 := s [1] ; 
		u2 := s [2] ; 
		r1:=[u1*p[1]+q[1], u1*p[2]+q[2], u1*p[3]+q[3] ] ; 
		r1 := tidyD ( r1 ) ; 
		r2:=[u2*p[1]+q[1], u2*p[2]+q[2], u2*p[3]+q[3] ] ; 
		r2 := tidyD ( r2 ) ; 
		matrix ( 2, 3 , [ r1, r2 ] ) ; 
		RETURN ( " ) ; 
	fi ; 
	end:

row		:=	proc ( r , M ) ; 
	# give a row of a matrix having three columns.
	[ M[r,1], M[r,2], M[r,3] ] ; 
	end: