/*
	extreme.c
	Written by Joseph O'Rourke
	orourke@cs.smith.edu
*/
#include	<stdio.h>
#include	<math.h>
#include <stdlib.h>
#define	ERROR	1
#define	X	0
#define	Y	1
typedef	enum { FALSE, TRUE }	bool;

#define DIM     2               /* Dimension of points */
typedef int     tPointi[DIM];   /* type integer point */
typedef double  tPointd[DIM];   /* type double point */
#define PMAX    1000            /* Max # of pts in polygon */
   
typedef tPointi tPolygoni[PMAX];/* type integer polygon */

void	SubVec( tPointi a, tPointi b, tPointi c );
int     Dot( tPointi a, tPointi b );
void	PrintPoly( int n, tPolygoni P, int labels[] );
void	PrintPoint( tPointi p );
int	Midway( int a, int b, int n );
int	Extreme( tPointi u, tPolygoni P, int n );
int	ReadPoly( tPolygoni P );

main()
{
	int	n, c;
	tPolygoni	P;
	tPointi	u;

	n = ReadPoly( P );
	u[0] = 0; u[1] = 1;
	c = Extreme( u, P, n );
	printf("Extreme returns %d\n", c);
	u[0] = 0; u[1] = -1;
	c = Extreme( u, P, n );
	printf("Extreme returns %d\n", c);
	u[0] = 1; u[1] = 0;
	c = Extreme( u, P, n );
	printf("Extreme returns %d\n", c);
	u[0] = -1; u[1] = 0;
	c = Extreme( u, P, n );
	printf("Extreme returns %d\n", c);
}

/*
	Returns index midway ccw from a to b, mod n.
*/
int	Midway( int a, int b, int n )
{
	if (a < b)
		return	( a + b ) / 2;
	else
		return	( ( a + b + n ) / 2 ) % n;
}

/*
	Returns index of a point extreme in direction u.
*/
int	Extreme( tPointi u, tPolygoni P, int n )
{
	int	a,a1, b;/* [a,b] includes extreme; a1=a+1. */
	int	c, c1;	/* index midway; c1 is c +- 1. */
	tPointi	A, C, B;/* edge vectors after a, after c, before c. */
	int	Adot, Cdot, Bdot;	/* dots with u. */
	int	y;	/* height difference: P[a] - P[c] in dir. u. */
	tPointi	ur;	/* u rotated cw by pi/2 */

	printf("\n==>Extreme: u = "); PrintPoint(u); putchar('\n');
	ur[0] = u[1];	ur[1] = -u[0];
	a = 0; b = 0;
	do {
		c = Midway( a, b, n );
		printf("Extreme: a,b,c=%d\t%d\t%d\n", a, b, c);

		/* Compute basic vectors and dots. */
		a1 = ( a + 1 ) % n;
		SubVec( P[a1], P[a], A );
		c1 = ( c + 1 ) % n;
		SubVec( P[c1], P[c], C );
		c1 = ( c + ( n-1 ) ) % n;
		SubVec( P[c], P[c1], B );
		printf("Extreme: A,C,B: ");
		PrintPoint( A ); putchar('\t');
		PrintPoint( C ); putchar('\t');
		PrintPoint( B ); putchar('\n');
		Adot = Dot(u,A);
		Cdot = Dot(u,C);
		Bdot = Dot(u,B);
		y = Dot(u,P[a]) - Dot(u,P[c]);
		printf("Extreme:\tAdot=%d\t", Adot);
		printf("Cdot=%d\t", Cdot);
		printf("Bdot=%d\t", Bdot);
		printf("y = %d\n", y);

		/* Termination conditions */
		/* If either A or C points left of u, then at extreme. */
		if ( (Adot == 0) && (Dot(ur,A) < 0) ) {
			printf("Extreme: A points left, so return a=%d\n",a);
			return	a;
		}
		if ( (Cdot == 0) && (Dot(ur,C) < 0) ) {
			printf("Extreme: C points left, so return c=%d\n",c);
			return	c;
		}
		/* From here on, can assume that zero dots put point on bot */

		if ( (Cdot < 0) && (Bdot > 0) ) {
			printf("Extreme: unique extreme c=%d\n", c);
			return	c;
		}
		if (a == c) {
			printf("Extreme: a=c=%d: choosing a or b=a+1\n", a);
			if (Adot > 0)
				return	b;
			else
				return	a;
		}

		/* Halving interval */
		if ( (Adot >= 0) && (Cdot <= 0) )
			b = c;		/* new: (a,c) */
		else if ( (Adot <= 0) && (Cdot >= 0) )
			a = c;		/* new: (c,b) */
		else if ( (Adot > 0) && (Cdot > 0) ) {
			if ( y > 0 )
				b = c;	/* new: (a,c) */
			else
				a = c;	/* new: (c,b) */
		}
		else if ( (Adot < 0) && (Cdot < 0) ) {
			if ( y < 0 )
				b = c;	/* new: (a,c) */
			else
				a = c;	/* new: (c,b) */
		}
		else {
			printf("Error in Extreme: shouldn't reach here\n");
			exit(1);
		}
	}
	while	( TRUE );
}
/*
        Puts a - b into c.
*/
void	SubVec( tPointi a, tPointi b, tPointi c )
{
        int     i;
        
        for( i = 0; i < DIM; i++ )
                c[i] = a[i] - b[i];
}

/*
        Returns the dot product of the two input vectors.
*/
int     Dot( tPointi a, tPointi b )
{
        int     i;
        int     sum;
	sum = 0;
/*
	printf("Dot: a, b =\n");
	PrintPoint(a); putchar('\n');
	PrintPoint(b); putchar('\n');
*/

        for( i = 0; i < DIM; i++ ) {
		/*printf("before: i=%d, sum=%d, a[i]=%d, b[i]=%d\n", 
			i, sum, a[i], b[i]);*/
                sum += a[i] * b[i];
		/*printf("after: i=%d, sum=%d, a[i]=%d, b[i]=%d\n", 
			i, sum, a[i], b[i]);*/
	}
	/*printf("Dot: returning %d\n", sum);*/
        
        return  sum;
}

/*
	Implements a = b, assignment of points/vectors.
	Assignment between arrays is not possible in C.
*/
void	PointAssign( tPointi a, tPointi b )
{
	int	i;

	for ( i = 0; i < DIM; i++ )
		a[i] = b[i];
}

void	PrintPoint( tPointi p )
{
	int	i;

	putchar('(');
	for ( i = 0; i < DIM; i++ ) {
		printf("%d", p[i]);
		if ( i != DIM-1 ) putchar(',');
	}
	putchar(')');
}
/*
	Reads in the coordinates of the vertices of a polygon from stdin,
	puts them into P, and returns n, the number of vertices.
	Formatting conventions: etc.
*/
int	ReadPoly( tPolygoni P )
{
	int	n = 0;

	printf("Polygon:\n");
	printf("  i   x   y\n");
	while ( (n < PMAX) && (scanf("%d %d",&P[n][0],&P[n][1]) != EOF) ) {
		printf("%3d%4d%4d\n", n, P[n][0], P[n][1]);
		++n;
	}
	if (n < PMAX)
		printf("n = %3d vertices read\n",n);
	else	printf("Error in read_poly:  too many points; max is %d\n", PMAX);
	putchar('\n');

	return	n;
}

void	PrintPoly( int n, tPolygoni P, int labels[] )
{
	int	i;

	printf("Polygon:\n");
	printf("  i   l   x   y\n");
	for( i = 0; i < n; i++ )
		printf("%3d%4d%4d%4d\n", i, labels[i], P[i][0], P[i][1]);
}