!<arch>
linprog.c       789421916   100   20    100644  7180      `
/* 
 * linprog.c
 *
 * Copyright (c) 1990 Michael E. Hohmeyer,
 *       hohmeyer@icemcfd.com
 * 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>
#include "lp.h"
#include "localmath.h"

int lp_no_constraints(int d,FLOAT n_vec[],FLOAT d_vec[],FLOAT opt[]);
int move_to_front(int i,int next[],int prev[]);

/* unitize a d+1 dimensional point */
lp_d_unit(int d, FLOAT a[], FLOAT b[]) {
	int i;
	FLOAT size;

	size = 0.0;
	for(i=0; i<=d; i++) 
		size += a[i]*a[i];
	if(size < (d+1)*EPS*EPS) return(1);
	size = 1.0/sqrt(size);
	for(i=0; i<=d; i++)
		b[i] = a[i]*size;
	return(0);
}
void vector_down(FLOAT elim_eqn[], int ivar, int idim,
		FLOAT old_vec[], FLOAT new_vec[]) {
	int i;
	FLOAT fac, ve, ee;
	ve = 0.0;
	ee = 0.0;
	for(i=0; i<=idim; i++) {
		ve += old_vec[i]*elim_eqn[i];
		ee += elim_eqn[i]*elim_eqn[i];
	}
	fac = ve/ee;
	for(i=0; i<ivar; i++) {
		new_vec[i] = old_vec[i] - elim_eqn[i]*fac;
	}
	for(i=ivar+1; i<=idim; i++) {
		new_vec[i-1] = old_vec[i] - elim_eqn[i]*fac;
	}
}
void plane_down(FLOAT elim_eqn[], int ivar, int idim,
	FLOAT old_plane[], FLOAT new_plane[]) {
	register FLOAT crit;
 	register int i;

	crit = old_plane[ivar]/elim_eqn[ivar];
	for(i=0; i<ivar; i++)  {
		new_plane[i] = old_plane[i] - elim_eqn[i]*crit;
	}
	for(i=ivar+1; i<=idim; i++)  {
		new_plane[i-1] = old_plane[i] - elim_eqn[i]*crit;
	}
}


#ifdef DOUBLE
dlinprog
#else
slinprog
#endif
(FLOAT halves[], /* halves  --- half spaces */
	int istart,     /* istart  --- should be zero
				 unless doing incremental algorithm */
	int m,  	/* m       --- terminal marker */
	FLOAT n_vec[], 	/* n_vec   --- numerator vector */
	FLOAT d_vec[], 	/* d_vec   --- denominator vector */
	int d, 		/* d       --- projective dimension */
	FLOAT opt[],	/* opt     --- optimum */
	FLOAT work[], 	/* work    --- work space (see below) */
	int next[], 	/* next    --- array of indices into halves */
	int prev[], 	/* prev    --- array of indices into halves */
	int max_size) 	/* max_size --- size of halves array */
/*
**
** half-spaces are in the form
** halves[i][0]*x[0] + halves[i][1]*x[1] + 
** ... + halves[i][d-1]*x[d-1] + halves[i][d]*x[d] >= 0
**
** coefficients should be normalized
** half-spaces should be in random order
** the order of the half spaces is 0, next[0] next[next[0]] ...
** and prev[next[i]] = i
**
** halves[max_size][d+1]
**
** the optimum has been computed for the half spaces
** 0 , next[0], next[next[0]] , ... , prev[istart]
** the next plane that needs to be tested is istart
**
** m is the index of the first plane that is NOT on the list
** i.e. m is the terminal marker for the linked list.
**
** the objective function is dot(x,nvec)/dot(x,dvec)
** if you want the program to solve standard d dimensional linear programming
** problems then n_vec = ( x0, x1, x2, ..., xd-1, 0)
** and           d_vec = (  0,  0,  0, ...,    0, 1)
** and halves[0] = (0, 0, ... , 1)
**
** work points to (max_size+3)*(d+2)*(d-1)/2 FLOAT space
*/
{
	int status;
	int i, j, k, imax;
	FLOAT *new_opt, *new_n_vec, *new_d_vec,  *new_halves, *new_work;
	FLOAT *plane_i;
	FLOAT val;

	if(d==1 && m!=0) {
		return(lp_base_case((FLOAT (*)[2])halves,m,n_vec,d_vec,opt,
			next,prev));
	} else {
		int d_vec_zero;
		val = 0.0;
		for(j=0; j<=d; j++) val += d_vec[j]*d_vec[j];
		d_vec_zero = (val < (d+1)*EPS*EPS);

/* find the unconstrained minimum */
		if(!istart) {
			status = lp_no_constraints(d,n_vec,d_vec,opt); 
		} else {
			status = MINIMUM;
		}
		if(m==0) return(status);
/* allocate memory for next level of recursion */
		new_opt = work;
		new_n_vec = new_opt + d;
		new_d_vec = new_n_vec + d;
		new_halves = new_d_vec + d;
		new_work = new_halves + max_size*d;
		for(i = istart; i!=m; i=next[i]) {
#ifdef CHECK
			if(i<0 || i>=max_size) {
				printf("index error\n");
				exit(1);
			}
#endif
/* if the optimum is not in half space i then project the problem
** onto that plane */
			plane_i = halves + i*(d+1);
/* determine if the optimum is on the correct side of plane_i */
			val = 0.0;
			for(j=0; j<=d; j++) val += opt[j]*plane_i[j];
			if(val<-(d+1)*EPS) {
/* find the largest of the coefficients to eliminate */
			    findimax(plane_i,d,&imax);
/* eliminate that variable */
			    if(i!=0) {
				FLOAT fac;
				fac = 1.0/plane_i[imax];
				for(j=0; j!=i; j=next[j]) {
					FLOAT *old_plane, *new_plane;
					int k;
					FLOAT crit;

					old_plane = halves + j*(d+1);
					new_plane = new_halves + j*d;
					crit = old_plane[imax]*fac;
					for(k=0; k<imax; k++)  {
						new_plane[k] = old_plane[k] - plane_i[k]*crit;
					}
					for(k=imax+1; k<=d; k++)  {
						new_plane[k-1] = old_plane[k] - plane_i[k]*crit;
					}
				}
			    }
/* project the objective function to lower dimension */
			    if(d_vec_zero) {
				vector_down(plane_i,imax,d,n_vec,new_n_vec);
				for(j=0; j<d; j++) new_d_vec[j] = 0.0;
			    } else {
			        plane_down(plane_i,imax,d,n_vec,new_n_vec);
			        plane_down(plane_i,imax,d,d_vec,new_d_vec);
			    }
/* solve sub problem */
			    status = linprog(new_halves,0,i,new_n_vec,
			    new_d_vec,d-1,new_opt,new_work,next,prev,max_size);
/* back substitution */
			    if(status!=INFEASIBLE) {
				    vector_up(plane_i,imax,d,new_opt,opt);
				{
/* in line code for unit */
				FLOAT size;
				size = 0.0;
				for(j=0; j<=d; j++) 
				    size += opt[j]*opt[j];
				size = 1.0/sqrt(size);
				for(j=0; j<=d; j++)
				    opt[j] *= size;
				}
			    } else {
				    return(status);
			    }
/* place this offensive plane in second place */
			    i = move_to_front(i,next,prev);
#ifdef CHECK
			    j=0;
			    while(1) {
/* check the validity of the result */
				val = 0.0;
				for(k=0; k<=d; k++) 
					val += opt[k]*halves[j*(d+1)+k];
				if(val <-(d+1)*EPS) {
				    printf("error\n");
				    exit(1);
				}
				if(j==i)break;
				j=next[j];
			    }
#endif
			} 
		}
 		return(status);
	}
}
/* returns the index of the plane that is in i's place */
move_to_front(int i,int next[],int prev[]) 
{
	int previ;
	if(i==0 || i == next[0]) return i;
	previ = prev[i];
/* remove i from it's current position */
	next[prev[i]] = next[i];
	prev[next[i]] = prev[i];
/* put i at the front */
	next[i] = next[0];
	prev[i] = 0;
	prev[next[i]] = i;
	next[0] = i;
	return(previ);
}
/* optimize the objective function when there are no contraints */
lp_no_constraints(int d,FLOAT n_vec[],FLOAT d_vec[],FLOAT opt[])
{
	int i;
	FLOAT n_dot_d, d_dot_d;

	n_dot_d = 0.0;
	d_dot_d = 0.0;
	for(i=0; i<=d; i++) {
		n_dot_d += n_vec[i]*d_vec[i];
		d_dot_d += d_vec[i]*d_vec[i];
	}
	if(d_dot_d < EPS*EPS) {
		d_dot_d = 1.0;
		n_dot_d = 0.0;
	}
	for(i=0; i<=d; i++) {
		opt[i] = -n_vec[i] + d_vec[i]*n_dot_d/d_dot_d;
	}
/* normalize the optimal point */
	if(lp_d_unit(d,opt,opt)) {
		opt[d] = 1.0;
		return(AMBIGUOUS);
	} else {
		return(MINIMUM);
	}
}
/* find the largest coefficient in a plane */
void findimax(FLOAT pln[],int idim,int *imax) {
	FLOAT rmax;
	int i;

	*imax = 0;
	rmax = ABS(pln[0]);
	for(i=1; i<=idim; i++) {
		FLOAT ab;
                ab = ABS(pln[i]);
		if(ab>rmax) {
			*imax = i;
			rmax = ab;
		}
	}
}

vector_up.c     789422187   100   20    100644  663       `
/* 
 * vector_up.c
 *
 * Copyright (c) 1990 Michael E. Hohmeyer,
 *       hohmeyer@icemcfd.com
 * 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 "tol.h"
#include "lp.h"


void vector_up(FLOAT equation[],int ivar,int idim,
		FLOAT low_vector[],FLOAT vector[])
{
	int i;

	vector[ivar] = 0.0;
	for(i=0; i<ivar; i++) {
		vector[i] = low_vector[i];
		vector[ivar] -= equation[i]*low_vector[i];
	}
	for(i=ivar+1; i<=idim; i++) {
		vector[i] = low_vector[i-1];
		vector[ivar] -= equation[i]*low_vector[i-1];
	}
	vector[ivar] /= equation[ivar];
}

lp_base_case.c  789422106   100   20    100644  6388      `
/* 
 * lp_base_case.c
 *
 * Copyright (c) 1990 Michael E. Hohmeyer,
 *       hohmeyer@icemcfd.com
 * 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>
#include "localmath.h"
#define not_zero(a) ((a) >= 2*EPS || (a) <= -2*EPS)

#include "lp.h"

int lp_no_constraints(int d,FLOAT n_vec[],FLOAT d_vec[],FLOAT opt[]);
/*
void lp_min_lin_rat(int degen, FLOAT cw_vec[2], 
		FLOAT ccw_vec[2], FLOAT n_vec[2], FLOAT d_vec[2], FLOAT opt[2]);
*/
int move_to_front(int i,int next[],int prev[]);


/*
int
unit2(FLOAT a[],FLOAT b[],FLOAT eps)
{
	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
 *
 */
lp_base_case(FLOAT halves[][2], 	/* halves --- half lines */
	int m, 				/* m      --- terminal marker */
	FLOAT n_vec[2],			/* n_vec  --- numerator funciton */
	FLOAT d_vec[2],			/* d_vec  --- denominator function */
	FLOAT opt[2],			/* opt    --- optimum  */
	int next[],
	int prev[])			/* next, prev  --- 
					double linked list of indices */
{
	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*EPS) {
		if(dot2(n_vec,n_vec) < 2*EPS*EPS ||
			 dot2(d_vec,d_vec) > 2*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(int degen,
		FLOAT cw_vec[2],
		FLOAT ccw_vec[2],
		FLOAT n_vec[2],
		FLOAT 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(FLOAT halves[][2],
	int m,
	int next[],
	int prev[],
	FLOAT cw_vec[],
	FLOAT ccw_vec[],
	int *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,EPS)) {
/* 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,EPS);
				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,EPS);
				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,EPS);
			} else if(d_ccw <= -2*EPS) { 
				(void)unit2(cw_vec,ccw_vec,EPS);
			}
			*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);
}
randperm.c      789422224   100   20    100644  449       `
/* 
 * randperm.c
 *
 * Copyright (c) 1990 Michael E. Hohmeyer,
 *       hohmeyer@icemcfd.com
 * Permission is granted to modify and re-distribute this code in any manner
 * as long as this notice is preserved.  All standard disclaimers apply.
 *       
 */
int rand();
randperm(int n,int perm[])
{
	int i, j, t;

	for(i=0; i<n; i++)
		perm[i] = i;
	for(i=0; i<n; i++) {
		j = rand()%(n-i)+i;
		t = perm[j];
		perm[j] = perm[i];
		perm[i] = t;
	}
}

randomize.c     664586038   100   20    100644  161       `
int rand();
void randomize(int n,int *perm)
{
	int i, j, t;

	for(i=0; i<n; i++) {
		j = rand()%(n-i)+i;
		t = perm[j];
		perm[j] = perm[i];
		perm[i] = t;
	}
}

unit2.c         719187711   100   20    100644  402       `
#include <math.h>

#ifdef DOUBLE
int
d_unit2(double a[],double b[],double eps)
{
	double size;
	size = sqrt(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);
}
#else
int
s_unit2(float a[],float b[],float eps)
{
	float size;
	size = sqrt(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);
}
#endif
README          781899599   100   20    100644  663       `
if you make "do_lp" and "randp" you can run the
following test.

randp 3 100 | do_lp

randp generates 100 random planes tangent to the hyperboloid
in 3 dimensions centered at the origin with z as an axis of
symmetry. Randp puts that on standard output. do_lp reads that and
feeds it to the routine linprog, which solves the problem.
If you look into do_lp.c you can observe how to call linprog.

Documentation is still rather sketchy so any questions can
be directed to hohmeyer@icemcfd.com

If you want to attempt to make this double precision, pay special
attention to the value EPS defined in tol.h which should be
the smallest number so that 1.0 + EPS  > 1.0

makefile        789428854   100   20    100644  1646      `
# DOUBLE = -DDOUBLE 
DOUBLE = 
CFLAGS = -c -O $(DOUBLE)
CFLAGS2 = -O $(DOUBLE)


#
# SGI
#
MACHINE = sgi
CC = cc
RANLIB = 

#
# HP
#
#MACHINE = hp
#CC = c89
#RANLIB = 

#
# SOL
#
#MACHINE = sol
#CC = cc
#RANLIB =

#
# SUN
#
#MACHINE = sun
#CC = acc
#RANLIB = ranlib Lib.sun.a

#
# IBM
#
#cc = c89 -c -O -I../ -D_XOPEN_SOURCE -DIBM
#MACHINE = ibm
#LIBS = -lgl -lm

# plane_down.c 
DCS =  linprog.c \
	vector_up.c \
	lp_base_case.c  \
	randperm.c \
	randomize.c \
	unit2.c


#plane_down.o 
DOS = linprog.o \
	vector_up.o \
	lp_base_case.o  \
	randperm.o \
	randomize.o \
	unit2.o

LIBA = Lib.$(MACHINE).a

$(LIBA):  $(DOS)
	rm -f $(LIBA); ar q $(LIBA) $(DOS); $(RANLIB)

llib-llinear.ln: $(DCS)
	lint -Clinear $(DCS)


lint: $(LIBA)
	lint -u -DCHECK $(DCS)

randperm.o: randperm.c
	$(CC) $(CFLAGS)  randperm.c 
lp_base_case.o: lp_base_case.c lp.h
	$(CC) $(CFLAGS)  lp_base_case.c 
linprog.o: linprog.c lp.h localmath.h
	$(CC) $(CFLAGS)  linprog.c 
#plane_down.o: plane_down.c
#	$(CC) $(CFLAGS)  plane_down.c
vector_up.o: vector_up.c
	$(CC) $(CFLAGS)  vector_up.c

randomize.o: randomize.c
	$(CC) $(CFLAGS)  randomize.c

do_lp: do_lp.c $(LIBA) lp.h unit2.o
	$(CC) $(CFLAGS2) -o do_lp do_lp.c unit2.o $(LIBA) -lm
lintdo_lp: llib-llinear.ln do_lp.c
	lint -u do_lp.c llib-llinear.ln
randp: randp.c
	$(CC) -g -o randp randp.c -lm
findeps: findeps.c
	$(CC) -o findeps findeps.c -lm

plane_gen: plane_gen.c
	$(CC) -g -o plane_gen plane_gen.c

clean:
	rm -f $(DOS) do_lp randp a.out

source.a: $(DCS) README makefile tol.h localmath.h do_lp.c lp.h randp.c
	rm -f source.a; \
	ar qv source.a $(DCS) README makefile tol.h localmath.h do_lp.c \
	lp.h randp.c
tol.h           695534610   100   20    100777  249       `
#ifndef TOL_H
#define TOL_H

#ifdef DOUBLE

/* this is the worst error we'll get doing a few additions
 * of numbers on the order of unity */
#define EPS		2.0e-16
#define FLOAT	double

#else

#define EPS 		1.0e-7
#define FLOAT	float

#endif

#endif

localmath.h     786262723   100   20    100644  4287      `
#ifndef LOCALMATH
#define LOCALMATH

#include "tol.h"

typedef FLOAT twoD[2];
#define ABS(a) ((a)>0.0 ? (a) : -(a))

/* line dot product line vs. point */
#define ldot(a,b) ((a)[0]*(b)[0] + (a)[1]*(b)[1] + (a)[2])

/*
 * line cross product a line from two points
 *  (az, ay, 1) x (bx, bx, 1) 
 */
#define lcross(a, b, c) {       (c)[0] = (a)[1] - (b)[1]; \
				(c)[1] = (b)[0] - (a)[0]; \
				(c)[2] = (a)[0]*(b)[1] - (a)[1]*(b)[0]; }

#define vsub2(a, b, c) 	{ (c)[0] = (a)[0] - (b)[0]; \
                        (c)[1] = (a)[1] - (b)[1]; }
#define vadd2(a, b, c)	{ (c)[0] = (a)[0] + (b)[0]; \
			(c)[1] = (a)[1] + (b)[1]; }
#define vcopy2(a, b) 	{ (b)[0] = (a)[0]; \
                        (b)[1] = (a)[1]; }

#define vmult2(s,a,b) 	{ (b)[0] = (s)*(a)[0]; \
			(b)[1] = (s)*(a)[1]; }
#define dot2(a,b) 	((a)[0]*(b)[0] + (a)[1]*(b)[1])
#define cross2(a,b) 	((a)[0]*(b)[1] - (a)[1]*(b)[0])
#define aprxlen2(a)	(ABS((a)[0]) + ABS((a)[1]))


#define midpoint(a, b, c)	{ (c)[0] = ((a)[0]  + (b)[0])*.5; \
				 (c)[1] = ((a)[1]  + (b)[1])*.5; \
				 (c)[2] = ((a)[2]  + (b)[2])*.5; }

#define dot(a,b) ((a)[0]*(b)[0] + (a)[1]*(b)[1] + (a)[2]*(b)[2])
#define min(a,b) ((a)<(b) ?(a):(b))
#define max(a,b) ((a)>(b) ?(a):(b))
#define square(a) ((a)*(a))
typedef FLOAT point[3];
typedef FLOAT vector[3];
typedef FLOAT triple[3];
typedef FLOAT plane[4];

#define vadd(a, b, c)	{ (c)[0] = (a)[0] + (b)[0]; \
			(c)[1] = (a)[1] + (b)[1]; \
			(c)[2] = (a)[2] + (b)[2]; }
#define vsub(a, b, c)	{ (c)[0] = (a)[0] - (b)[0]; \
			(c)[1] = (a)[1] - (b)[1]; \
			(c)[2] = (a)[2] - (b)[2]; }
#define vmult(s, a, b) 	{ (b)[0] = (s)*(a)[0]; \
			(b)[1] = (s)*(a)[1]; \
			(b)[2] =(s)*(a)[2]; }
#define size2(a) dot((a),(a))
#define vzero2(a)	{ (a)[0] = (a)[1] = 0.0; }


#define vlinc(s, a, t, b, c)  { (c)[0] = (s)*(a)[0] + (t)*(b)[0]; \
				(c)[1] = (s)*(a)[1] + (t)*(b)[1]; \
				(c)[2] = (s)*(a)[2] + (t)*(b)[2]; }

/* multiply and add */
#define vmadd(s, a, b, c)  { (c)[0] = (s)*(a)[0] + (b)[0]; \
				(c)[1] = (s)*(a)[1] + (b)[1]; \
				(c)[2] = (s)*(a)[2] + (b)[2]; }

#define vzero(a)	{ (a)[0] = (a)[1] = (a)[2] = 0.0; }

#define vcopy(a, b)	{ (b)[0] = (a)[0]; \
			  (b)[1] = (a)[1]; \
			  (b)[2] = (a)[2]; }
#define vneg(a,b) 	{ (b)[0] = -(a)[0]; \
			  (b)[1] = -(a)[1]; \
			  (b)[2] = -(a)[2]; }


#ifdef DOUBLE

#define unit2(a,b,eps)		d_unit2(a,b,eps)
#define triple_cross(a,b,c)	d_triple_cross(a,b,c)
#define cross(a,b,c)		d_cross(a,b,c)
#define unit(a,b)		d_unit(a,b)
#define solve3(a,x,b)		d_solve3(a,x,b)
#define gen_perp(a,b)		d_gen_perp(a,b)
#define distance(a,b)		d_distance(a,b)
#define vswap(a,b)		d_vswap(a,b)
#define Fsort(a,b)		d_sort(a,b)
#define Fmerge_sort(a,b,c)	d_merge_sort(a,b,c)

#else

#define unit2(a,b,eps)		s_unit2(a,b,eps)
#define triple_cross(a,b,c)	s_triple_cross(a,b,c)
#define cross(a,b,c)		s_cross(a,b,c)
#define unit(a,b)		s_unit(a,b)
#define solve3(a,x,b)		s_solve3(a,x,b)
#define gen_perp(a,b)		s_gen_perp(a,b)
#define distance(a,b)		s_distance(a,b)
#define vswap(a,b)		s_vswap(a,b)
#define Fsort(a,b)		s_sort(a,b)
#define Fmerge_sort(a,b,c)	s_merge_sort(a,b,c)

#endif


#ifdef __cplusplus
extern "C" {
#endif

void i_sort(int size, int *list);

#ifdef DOUBLE

int d_unit2(double *a, double *b, double eps);
double d_triple_cross(double *a, double *b, double *c);
void d_cross(double *a, double *b, double *c);
int d_unit(double *a, double *b); 
int d_solve3(double a[3][3], double *b, double *x);
void d_gen_perp(double *v, double *perp);
double d_distance(double *a, double *b);
void d_vswap(double *a, double *b);
void d_sort(int i, double *a);
void d_merge_sort(int size, float *list, float *temp);

#else 

int s_unit2(float *a, float *b, float eps);
float s_triple_cross(float *a, float *b, float *c);
void s_cross(float *a, float *b, float *c);
int s_unit(float *a, float *b); 
int s_solve3(float a[3][3], float *b, float *x);
void s_gen_perp(float *v, float *perp);
float s_distance(float *a, float *b);
void s_vswap(float *a, float *b);
void s_sort(int i, float *a);
void s_merge_sort(int size, float *list, float *temp);
int invert3(float from[3][3], float to[3][3]);

#endif


#ifdef __cplusplus
}
#endif

#define SWAP(type,a,b) {type temp = (a); (a) = (b); (b) = temp; }
#define ORDER(type, a, b)       if((a) > (b)) SWAP(type, a, b);

#endif

do_lp.c         694931840   100   20    100644  3957      `
#define REAL 0
#define PROJECTIVE 1

#include "lp.h"
char *malloc();
#include <stdio.h>


main(argc,argv)
int argc;
char *argv[];
{
	int *next, *prev;
	FLOAT *halves, *n_vec, *opt, *d_vec;
	int *perm, r;
	int m, d, i, j, k;
	int status, type;
	FLOAT *work, val;
	int repeat;

	if(argc>1) {
		repeat = atoi(argv[1]);
	} else {
		repeat = 1;
	}
	if(read_lp(&halves,&d,&m,&n_vec, &d_vec, &type)) {
		perm = (int*)malloc((unsigned)(m-1)*sizeof(int));
		next = (int*)malloc((unsigned)(m)*sizeof(int));
		prev = (int*)malloc((unsigned)(m)*sizeof(int));
		printf("dimension: %d, number of planes: %d, repeat: %d\n",
			d,m,repeat);
		opt = (FLOAT *)malloc((unsigned)(d+1)*sizeof(FLOAT));
		work = (FLOAT *)malloc((unsigned)(m+3)*(d+2)*(d-1)/2*sizeof(FLOAT));
		for(r=0; r<repeat; r++) {
/* randomize the input planes */
			randperm(m-1,perm);
/* previous to 0 should never be used */
			prev[0] = 1234567890;
/* link the zero position in at the beginning */
			next[0] = perm[0]+1;
			prev[perm[0]+1] = 0;
/* link the other planes */
			for(i=0; i<m-2; i++) {
				next[perm[i]+1] = perm[i+1]+1;
				prev[perm[i+1]+1] = perm[i]+1;
			}
/* flag the last plane */
			next[perm[m-2]+1] = m;
			status = linprog(halves,0,m,n_vec,d_vec,d,opt,work,
				next,prev,m);
		}
		switch(status) {
			case INFEASIBLE:
				(void)printf("no feasible solution\n");
				break;
			case MINIMUM:
				(void)printf("minimum attained at\n");
				break;
			case UNBOUNDED:
				(void)printf("region is unbounded: last vertex is\n");
				break;
			case AMBIGUOUS:
			(void)printf("region is bounded by plane orthogonal\n");
 			(void)printf("to minimization vector: one vertex is\n");
				break;
			default:
				(void)printf("unknown case returned from linprog\n");
		}
		if(status!=INFEASIBLE) {
			if(type == PROJECTIVE) {
				(void)printf("(");
				for(i=0; i<d; i++) {
					if(opt[d]==0.0) {
						(void)printf("%f",opt[i]);
					} else {
						(void)printf("%f",opt[i]/opt[d]);
					}
					if(i!=d-1)(void)printf(",");
				}
				(void)printf(")\n");
				if(opt[d]==0.0) (void)printf("at infinite\n");
				for(j=0; j!=m; j=next[j]) {
				    FLOAT val = 0.0;
			    for(k=0; k<=d; k++) val += opt[k]*halves[j*(d+1)+k];
                                    if(val <-d*EPS) {
                                        printf("error\n");
                                        exit(1);
                                    }
                                }
			} else {
				(void)printf("(");
				for(i=0; i<=d; i++) {
					(void)printf("%f",opt[i]);
					if(i!=d)(void)printf(",");
				}
				(void)printf(")\n");
			}
		}
		free((char *)halves);
		free((char *)next);
		free((char *)prev);
		free((char *)perm);
		free((char *)n_vec);
		free((char *)opt);
		free((char *)work);
	} else {
		(void)printf("parse error\n");
	}
}
#ifdef  DOUBLE
#define	CONVERSION	"%lf"
#else
#define CONVERSION	"%f"
#endif
read_lp(halves,d,m,n_vec,d_vec,itype)
FLOAT **halves, **n_vec, **d_vec;
int *d, *m, *itype;
{
	int i, j;
	char type[100];

	if(scanf("dimension: %d, number of planes: %d\n",d,m) !=2)return(0);
	if((i=scanf("%s",type)) !=1)return(0);
	*n_vec = (FLOAT *)malloc((unsigned)(*d+1)*sizeof(FLOAT));
	*d_vec = (FLOAT *)malloc((unsigned)(*d+1)*sizeof(FLOAT));
	*halves = (FLOAT *)malloc((unsigned)*m*(*d+1)*sizeof(FLOAT));
	for(i=0; i<=*d; i++) {
		(*d_vec)[i] = 0.0;
	}
	if(!strcmp(type,"projective")) {
		*itype = PROJECTIVE;
		for(i=0; i<*d; i++) {
			if(scanf(CONVERSION,(*n_vec)+i)!=1)  goto err;
		}
		(*d_vec)[*d] = 1.0;
		(*n_vec)[*d] = 0.0;
	} else if(!strcmp(type,"real")){
		*itype = REAL;
		for(i=0; i<=*d; i++) {
			if(scanf(CONVERSION,(*n_vec)+i)!=1)  goto err;
		}
	} else {
		goto err;
	}
	for(i=0; i<*m; i++) {
		for(j=0; j<=*d; j++) {
			if(scanf(CONVERSION,(*halves)+i*(*d+1)+j)!=1) {
				goto err;
			}
		}
		(void)lp_d_unit(*d,(*halves) + i*(*d+1),(*halves) + i*(*d+1));
	}
	return(1);
 err: 	free((char *)*d_vec);
	free((char *)*n_vec);
	free((char *)*halves);
	return(0);
}

lp.h            789422092   100   20    100644  1685      `
/* 
 * lp.h
 *
 * Copyright (c) 1990 Michael E. Hohmeyer, 
 *       hohmeyer@icemcfd.com
 * Permission is granted to modify and re-distribute this code in any manner
 * as long as this notice is preserved.  All standard disclaimers apply.
 *
 * version history
 *     1/22/1991 - original version
 *     1/5/1995  - fix bug in projection of degenerate objective
 *                 function to lower dimension
 */
/* status from lp_intersect or linprog */ 
#define INFEASIBLE 0
#define MINIMUM 1
#define UNBOUNDED 2
#define AMBIGUOUS 3

/* status from plane_down */
#define REDUNDANT 0
#define PROPER 1

#include "tol.h"

#ifdef DOUBLE
#define linprog(v,istart, n,num,den,dim,opt,work,next,prev,max_size)  \
dlinprog(v,istart, n,num,den,dim,opt,work,next,prev,max_size)
#else
#define linprog(v,istart, n,num,den,dim,opt,work,next,prev,max_size)  \
slinprog(v,istart, n,num,den,dim,opt,work,next,prev,max_size)
#endif

#ifdef __cplusplus
extern "C" {
#endif
    void randperm(int i, int *p);
    void randomize(int n, int *perm);
    int linprog(FLOAT *v, int istart,int n, FLOAT *num, FLOAT *den,
    int dim, FLOAT *opt, FLOAT *work, int *next, int *prev, int max_size);
    int lp_base_case(FLOAT halves[][2], int m, FLOAT n_vec[2],
        FLOAT d_vec[2], FLOAT opt[2], int *next, int *prev);
    int wedge(FLOAT halves[][2], int m, int *next, int *prev,
        FLOAT cw_vec[2], FLOAT ccw_vec[2], int *degen);
    void plane_down(FLOAT *elim_eqn, int ivar, int idim,
        FLOAT *old_plane, FLOAT *new_plane);
    void findimax(FLOAT *pl,int idim,int *imax);
    void vector_up(FLOAT *equation,int ivar,int idim,
        FLOAT *low_vector,FLOAT *vec);
#ifdef __cplusplus
}
#endif

randp.c         724809153   100   20    100644  662       `
#include <math.h>
#include <stdio.h>
#include "tol.h"

int rand();
double pow();

main(int argc,char *argv[])
{
	int i, j, d, m;
	FLOAT x, z, fac;

	if(argc!=3) {
		fprintf(stderr,"incorrect arguments\n");
		exit(1);
	}
	d = atoi(argv[1]);
	m = atoi(argv[2]);
	printf("dimension: %d, number of planes: %d\n",d,m+1);
	printf("projective\n");
	for(i=0; i<d-1; i++) printf("0 ");
	printf("1\n");
	for(i=0; i<d; i++) printf("0 ");
	printf("1\n");
	fac = 1.0/pow(2.0,14.0);
	for(i=0; i<m; i++) {
		z = 0.0;
		for(j=0; j<d-1; j++) {
			x = rand()*fac - 1.0;
			x = (fabs(x) < 1.0e-5) ? 1.0e-5 : x ;
			z += x*x;
			printf("%f ",-2.0*x);
		}
		printf("1 %f\n",z);
	}
}