/* 
 * lp_base_case.c
 *
 * Copyright (c) 1990 Michael E. Hohmeyer,
 *       hohmeyer@miro.berkeley.edu,
 * Permission is granted to modify and re-distribute this code in any manner
 * as long as this notice is preserved.  All standard disclaimers apply.
 *
 */
#include "math.h"
#define not_zero(a) ((a) >= 2*EPS || (a) <= -2*EPS)
#define cross2(a,b) (a)[0]*(b)[1] - (a)[1]*(b)[0]
#define dot2(a,b) (a)[0]*(b)[0] + (a)[1]*(b)[1]
#include <math.h>
#include "lp.h"
#ifdef GCC
inline
#endif
unit2(a,b)
float a[], b[];
{
	float size;

	size = fsqrt(a[0]*a[0] + a[1]*a[1]);
	if(size < 2*EPS) return(1);
	b[0] = a[0]/size;
	b[1] = a[1]/size;
	return(0);
}
/*
 * return the minimum on the projective line
 *
 * halves --- half lines
 * m      --- terminal marker
 * n_vec  --- numerator funciton
 * d_vec  --- denominator function
 * opt    --- optimum 
 * next, prev  --- double linked list of indices
 */
lp_base_case(halves,m,n_vec,d_vec,opt, next, prev)
float halves[][2], n_vec[2], d_vec[2], opt[2];
int m;
int *next, *prev;
{
	float cw_vec[2], ccw_vec[2];
	float d_cw;
	int i, degen;
	int status;
	float ab;

/* find the feasible region of the line */
	status = wedge(halves,m,next,prev,cw_vec,ccw_vec,&degen);

	if(status==INFEASIBLE) return(status);
/* no non-trivial constraints one the plane: return the unconstrained
** optimum */
	if(status==UNBOUNDED) {
		return(lp_no_constraints(1,n_vec,d_vec,opt));
	}
        ab = ABS(cross2(n_vec,d_vec));
	if(ab < 2*EPS) {
		if(dot2(n_vec,n_vec) < 4*EPS*EPS || dot2(d_vec,d_vec) > 4*EPS*EPS) {
/* numerator is zero or numerator and denominator are linearly dependent */
			opt[0] = cw_vec[0];
			opt[1] = cw_vec[1];
			status = AMBIGUOUS;
		} else {
/* numerator is non-zero and denominator is zero 
** minimize linear functional on circle */
			if(!degen && cross2(cw_vec,n_vec) <= 0.0 &&
			cross2(n_vec,ccw_vec) <= 0.0 ) {
/* optimum is in interior of feasible region */
				opt[0] = -n_vec[0];
				opt[1] = -n_vec[1];
			} else if(dot2(n_vec,cw_vec) > dot2(n_vec,ccw_vec) ) {
/* optimum is at counter-clockwise boundary */
				opt[0] = ccw_vec[0];
				opt[1] = ccw_vec[1];
			} else {
/* optimum is at clockwise boundary */
				opt[0] = cw_vec[0];
				opt[1] = cw_vec[1];
			}
			status = MINIMUM;
		}
	} else {
/* niether numerator nor denominator is zero */
		lp_min_lin_rat(degen,cw_vec,ccw_vec,n_vec,d_vec,opt);
		status = MINIMUM;
	}
#ifdef CHECK
	for(i=0; i!=m; i=next[i]) {
		d_cw = dot2(opt,halves[i]);
		if(d_cw < -2*EPS) {
			printf("error at base level\n");
			exit(1);
		}
	}
#endif
	return(status);
}
lp_min_lin_rat(degen,cw_vec,ccw_vec,n_vec,d_vec,opt)
int degen;
float cw_vec[2], ccw_vec[2], n_vec[2], d_vec[2];
float opt[2];
{
	float d_cw, d_ccw, n_cw, n_ccw;

/* linear rational function case */
	d_cw = dot2(cw_vec,d_vec);
	d_ccw = dot2(ccw_vec,d_vec);
	n_cw = dot2(cw_vec,n_vec);
	n_ccw = dot2(ccw_vec,n_vec);
	if(degen) {
/* if degenerate simply compare values */
		if(n_cw/d_cw < n_ccw/d_ccw) {
			opt[0] = cw_vec[0];
			opt[1] = cw_vec[1];
		} else {
			opt[0] = ccw_vec[0];
			opt[1] = ccw_vec[1];
		}
/* check that the clock-wise and counter clockwise bounds are not near a poles */
	} else if(not_zero(d_cw) && not_zero(d_ccw)) {
/* the valid region does not contain a poles */
		if(d_cw*d_ccw > 0.0) {
/* find which end has the minimum value */
			if(n_cw/d_cw < n_ccw/d_ccw) {
				opt[0] = cw_vec[0];
				opt[1] = cw_vec[1];
			} else {
				opt[0] = ccw_vec[0];
				opt[1] = ccw_vec[1];
			}
		} else {
/* the valid region does contain a poles */
			if(d_cw > 0.0) {
				opt[0] = -d_vec[1];
				opt[1] = d_vec[0];
			} else {
				opt[0] = d_vec[1];
				opt[1] = -d_vec[0];
			}
		}
	} else if(not_zero(d_cw)) {
/* the counter clockwise bound is near a pole */
		if(n_ccw*d_cw > 0.0) {
/* counter clockwise bound is a positive pole */
			opt[0] = cw_vec[0];
			opt[1] = cw_vec[1];
		} else {
/* counter clockwise bound is a negative pole */
			opt[0] = ccw_vec[0];
			opt[1] = ccw_vec[1];
		}
	} else if(not_zero(d_ccw)) {
/* the clockwise bound is near a pole */
		if(n_cw*d_ccw > 2*EPS) {
/* clockwise bound is at a positive pole */
			opt[0] = ccw_vec[0];
			opt[1] = ccw_vec[1];
		} else {
/* clockwise bound is at a negative pole */
			 opt[0] = cw_vec[0];
			 opt[1] = cw_vec[1];
		} 
	} else {
/* both bounds are near poles */
		if(cross2(d_vec,n_vec) > 0.0) {
			opt[0] = cw_vec[0];
			opt[1] = cw_vec[1];
		} else {
			opt[0] = ccw_vec[0];
			opt[1] = ccw_vec[1];
		}
	}
}
wedge(halves,m,next,prev,cw_vec,ccw_vec,degen)
float halves[][2], cw_vec[2], ccw_vec[2];
int m, *next, *prev, *degen;
{
	int i;
	float d_cw, d_ccw;
	int offensive;

	*degen = 0;
	for(i=0;i!=m;i = next[i]) {
		if(!unit2(halves[i],ccw_vec)) {
/* clock-wise */
			cw_vec[0] = ccw_vec[1];
			cw_vec[1] = -ccw_vec[0];
/* counter-clockwise */
			ccw_vec[0] = -cw_vec[0];
			ccw_vec[1] = -cw_vec[1];
			break;
		}
	}
	if(i==m) return(UNBOUNDED);
	i = 0;
	while(i!=m) {
		offensive = 0;
		d_cw = dot2(cw_vec,halves[i]);
		d_ccw = dot2(ccw_vec,halves[i]);
		if(d_ccw >= 2*EPS) {
			if(d_cw <= -2*EPS) {
				cw_vec[0] = halves[i][1];
				cw_vec[1] = -halves[i][0];
				(void)unit2(cw_vec,cw_vec);
				offensive = 1;
			}
		} else if(d_cw >= 2*EPS) {
			if(d_ccw <= -2*EPS) {
				ccw_vec[0] = -halves[i][1];
				ccw_vec[1] = halves[i][0];
				(void)unit2(ccw_vec,ccw_vec);
				offensive = 1;
			}
		} else if(d_ccw <= -2*EPS && d_cw <= -2*EPS) {
			return(INFEASIBLE);
		} else if((d_cw <= -2*EPS) ||
			(d_ccw <= -2*EPS) ||
			(cross2(cw_vec,halves[i]) < 0.0)) {
/* degenerate */
			if(d_cw <= -2*EPS) {
				(void)unit2(ccw_vec,cw_vec);
			} else if(d_ccw <= -2*EPS) { 
				(void)unit2(cw_vec,ccw_vec);
			}
			*degen = 1;
			offensive = 1;
		}
/* place this offensive plane in second place */
		if(offensive) i = move_to_front(i,next,prev);
		i = next[i];
		if(*degen) break;
	}
	if(*degen) {
		while(i!=m) {
			d_cw = dot2(cw_vec,halves[i]);
			d_ccw = dot2(ccw_vec,halves[i]);
			if(d_cw < -2*EPS) {
				if(d_ccw < -2*EPS) {
					return(INFEASIBLE);
				} else {
					cw_vec[0] = ccw_vec[0];
					cw_vec[1] = ccw_vec[1];
				}
			} else if(d_ccw < -2*EPS) {
				ccw_vec[0] = cw_vec[0];
				ccw_vec[1] = cw_vec[1];
			}
			i = next[i];
		}
	}
	return(MINIMUM);
}