int *bestfreq;
#include <math.h>
#include <stdio.h>

#define N 200

struct point {
	double x;
	double y;
	double phi;
	char *name;
};
struct point testpts[20];
struct point *xarray;
struct point zero = {
	0., 0., 0., "zero"};

struct way {
	double value;
	double t,u,v,w;
	char *name;
	struct point *from, *to;
};

#define REVERSE	1
#define REFLECT 2

struct way D();
int debug;
char *calloc();
char *strdup();
char *getenv();
double mod2pi();


main(ac, av)
char **av;
{
	struct way **dist;
	struct point *p;
	struct point *q;
	int i, j, k, n = N;

	tabdub();
	if(ac == 3) {
		srand(atoi(av[2]));
	}
	if(ac >= 2) {
		n = atoi(av[1]);
	}
	for(i=0; i<20; i++) {
		char buf[100];

		randpoint( &testpts[i]);
		sprintf(buf, "r%d", i+1);
		testpts[i].name = strdup(buf);
	}
	dochecks(20, &testpts[0]);
	if(av[1][0] == '-')  {
		wx();
		return(0);
	}

	p = xarray = (struct point *)calloc((unsigned)n, sizeof(*p));

	for(i=0; i<n; p++,i++) {
		char buf[100];

		randpoint( p);
		sprintf(buf, "p%d", i+1);
		p->name = strdup(buf);
	}


	dist = (struct way **)calloc((unsigned)n, sizeof (struct way *));
	for(i=0; i<n; i++)
		dist[i] = (struct way *)calloc((unsigned)n, sizeof(struct way));
	for(i=0; i<n; i++)
		for(j=0; j<n; j++) {
			p = &xarray[i];
			q = &xarray[j];
			if(i != j)dfunc(p, q, &dist[i][j]);
		}
#define trip(i,j) &dist[i][j]
#define d(i,j) dist[i][j].value
	for(i=0; i<n; i++) {
		for(j=i+1; j<n; j++) {
			for(k=j+1; k<n; k++) {
				if(d(i,j) + d(j,k) +.000001 < d(i,k) ) tritest(trip(i,j), trip(j,k), trip(i,k));
				if(d(j,k) + d(k,i) +.000001 < d(j,i) ) tritest(trip(j,k), trip(k,i), trip(j,i));
				if(d(k,i) + d(i,j) +.000001 < d(k,j) ) tritest(trip(k,i), trip(i,j), trip(k,j));
			}
		}
	}
	report();
	return 0;
}
#undef d
#undef trip

tritest(wab, wbc, wac)
struct way *wab, *wbc, *wac;
{
	struct point *a, *b, *c;

	if(wac->value <= wab->value+wbc->value) return;

	a = wab->from;
	b = wab->to;
	c = wac->to;

	if(a != wac->from) {
		printf("ABSURD 1\n"); 
		exit(1);
	}
	if(b != wbc->from) {
		printf("ABSURD 2\n"); 
		exit(1);
	}
	if(c != wbc->to) {
		printf("ABSURD 3\n"); 
		exit(1);
	}

	printf("d(%s,%s)=%g > d(%s,%s)+d(%s,%s) = %g + %g = %g\n",
	    a->name, c->name, wac->value,
	    a->name, b->name, b->name, c->name,
	    wab->value, wbc->value, wab->value + wbc->value);

	pshow(a);
	pshow(b);
	pshow(c);

	wshow(wab);
	wshow(wbc);
	wshow(wac);

	if(getenv("DEBUG")) {
		printf("%s->%s\n", c->name, b->name); 
		dfunc(c, b, (struct way *)0);
		printf("%s->%s\n", a->name, b->name); 
		dfunc(a, b, (struct way *)0);
		printf("%s->%s\n", a->name, c->name); 
		dfunc(a, c, (struct way *)0);
	}
	printf("\n\n");
	fflush(stdout);
}
randpoint(p)
struct point *p;
{
	gauss(&p->x, &p->y);
	if(frand()<.5) {
		p->x *= .1;
		p->y *= .1;
	}
	p->phi = mod2pi(PI2*frand());
}
pshow(p)
struct point *p;
{
	printf("%s:	[%10.6g %10.6g %10.6g]\n", p->name, p->x, p->y, p->phi);
}
wshow(p)
struct way *p;
{
	printf("%10.6g %s: t=%g u=%g v=%g w=%g",
	    p->value, p->name, p->t, p->u, p->v, p->w);
	printf("\t%s->%s\n",
	    p->from->name, p->to->name);
}
gauss(xp, yp)
double *xp, *yp;
{
	double v1, v2, s;

	do {
		v1 = 2*frand()-1;
		v2 = 2*frand()-1;
		s = v1*v1 + v2*v2;
	} while (s == 0 || s >= 1);

	s = sqrt(-2*log(s)/s);

	*xp = v1*s;
	*yp = v2*s;
}

struct way  D(p0, p1)
struct point *p0, *p1;
{
	double r, theta;
	double x, y, phi;
	double xr, yr;
	double xv[4], yv[4], phiv[4];
	struct way *wbest, val[4];
	int *fbest, i;

	rec2polar(p1->x-p0->x, p1->y-p0->y, &r, &theta);
	x = r*cos(theta-p0->phi);
	y = r*sin(theta-p0->phi);
	phi = p1->phi-p0->phi;

	theta += PI;
	/* rec2polar(p0->x-p1->x, p0->y-p1->y, &r, &theta); */
	xr = r*cos(theta-p1->phi);
	yr = r*sin(theta-p1->phi);

	xv[0] = xv[REFLECT] = x;
	xv[REVERSE] = xv[REVERSE|REFLECT] = xr;
	yv[0] = y;
	yv[REFLECT] = -y;
	yv[REVERSE] = yr;
	yv[REVERSE|REFLECT] = -yr;
	phiv[0] = phi;
	phiv[REFLECT] = -phi;
	phiv[REVERSE] = -phi;
	phiv[REVERSE|REFLECT] = phi;

	wbest = 0;
	for(i=0; i<4; i++) {
		val[i].from = p0;
		val[i].to = p1;
		if(L(xv[i], yv[i], phiv[i], &val[i], i) == 0) continue;
		if(wbest == 0 || val[i].value < wbest->value) {
			fbest = bestfreq;
			wbest = &val[i];
		}
	}
	if(wbest == 0) {
		printf("XXXXXX\n");
		exit(1);
	}
	bestfreq = fbest;
	i = wbest - val;

	if(i & REVERSE) {
		double t;

		t = wbest->t;
		wbest->t = wbest->w;
		wbest->w = t;

		t = wbest->u;
		wbest->u = wbest->v;
		wbest->v = t;

		wbest->t = - wbest->t;
		wbest->u = - wbest->u;
		wbest->v = - wbest->v;
		wbest->w = - wbest->w;
	}

	checkit(wbest, 0);

	return *wbest;
}
checkit(wp, trips)
struct way *wp;
{
	double norm();
	struct point *p0, *p1;
	struct point q, q2;

	p0 = wp->from;
	p1 = wp->to;
	q.name = "drive temp";
	q2.name = "drive temp 2";


	drive(wp->name[0], 	p0, 	&q, 	wp->t, 	trips);
	drive(wp->name[1],       &q, &q2, wp->u, trips);
	drive(wp->name[2],      &q2,  &q, wp->v, trips);
	drive(wp->name[3],       &q, &q2, wp->w, trips);

	if(trips == 1)return;
	if(norm(&q2, p1)  > .001) {
		printf("NONO %s\n", wp->name);
		pshow(p0);
		pshow(p1);
		checkit(wp, 1);
		fflush(stdout);
	}
}
double norm (x,y)
struct point *x, *y;
{
	double s, d;

	s = 0;
	d = x->x - y->x ; 
	s += d*d;
	d = x->y - y->y ; 
	s += d*d;
	d = mod2pi(x->phi - y->phi) ; 
	s += d*d;

	return sqrt(s);
}
drive(c, xold, xnew, t, flag)
struct point *xold, *xnew;
double t;
{
	switch(c){
	case 's':
		xnew->x  = xold->x + t*cos(xold->phi);
		xnew->y  = xold->y + t*sin(xold->phi);
		xnew->phi = xold->phi;
		break;
	case 'r':
		xnew->x  = xold->x - sin(xold->phi - t) + sin(xold->phi);
		xnew->y  = xold->y + cos(xold->phi - t) - cos(xold->phi);
		xnew->phi = mod2pi(xold->phi - t);
		break;
	case 'l':
		xnew->x  = xold->x + sin(xold->phi + t) - sin(xold->phi);
		xnew->y  = xold->y - cos(xold->phi + t) + cos(xold->phi);
		xnew->phi = mod2pi(xold->phi + t);
		break;
	default:
		*xnew = *xold;
		return;
		break;
	}
	if(flag) {
		printf("%c(%g) sends ", c, t);
		pshow(xold);
		printf("to ");
		pshow(xnew);
	}
}
dfunc(p1, p2, q)
struct point *p1, *p2;
struct way *q;
{
	if(q == 0) {
		debug = 1;
		(void)D(p1, p2);
		debug = 0;
	} else {
		*q = D(p1,p2);
		(*bestfreq)++;
	}
}


#define DECLARE(func) func(x,y,phi,wp)double x,y,phi;struct way *wp;
#define RET(_t,_u,_v,_w) {wp->t=_t;\
			wp->u=_u;\
			wp->v=_v;\
			wp->w=_w;\
			wp->value=fabs(_t)+fabs(_u)+fabs(_v)+fabs(_w);\
			return 1;}
int lrsla(), lrslb();
int lrsl3();
int lrsl4();
int lrsl5();
int lrsr1();
int lrsr2();
int lrsr4();
int lrsr5();
int lrs();
int lrlr();
int lrl(), lsr1(), lsr2(), lsl();

struct tab {
	char *s;
	int (*f)();
	char *ss[4];
	int freq;
} tab[] = {
	{"lrs ", lrs, 0},
	{"lrsla", lrsla, 0},
	{"lrslb", lrslb, 0},
	{"lrsl3", lrsl3, 0},
	{"lrsl4", lrsl4, 0},
	{"lrsl5", lrsl5, 0},
	{"lrsr1", lrsr1, 0},
	{"lrsr2", lrsr2, 0},
	{"lrsr4", lrsr4, 0},
	{"lrsr5", lrsr5, 0},
	{"lrlr", lrlr, 0},
	{"lrl ", lrl, 0},
	{"lsr +", lsr1, 0},
	{"lsr -", lsr2, 0},
	{"lsl ", lsl, 0},
	{(char *)0, 0}
};

tabdub() {
	int i;
	char name[50], *strsup();
	struct tab *p;

	for(p=tab; p->s; p++) {
		for(i=0; i<4; i++) {
			sprintf(name, "%s(%d.%d)", p->s, p-tab, i);
			p->ss[i] = strsup(name, i);
		}
	}
}
char *strsup(s, code)
char *s;
{
	char *t;

	t = strdup(s);
	if(REFLECT&code)
		for(s=t; s<t+4; s++) {
			switch(*s) {
			case 'l': 
				*s = 'r';	
				break;
			case 'r': 
				*s = 'l';	
				break;
			}
		}
	if(REVERSE&code) {
		char c;

		c = t[0]; 
		t[0] =  t[3]; 
		t[3] = c;
		c = t[1]; 
		t[1] =  t[2]; 
		t[2] = c;
	}
	return t;
}
dochecks(n, pt)
struct point *pt;
{
	struct way w;
	struct tab *p;
	int i;

	printf("first, some simple checks\n");

	for(i=0; i<n; pt++,i++) {
		w.from = &zero;
		w.to = pt;

		for(p=tab; p->s; p++) {
			w.name = p->ss[0];
			if((*p->f)(pt->x, pt->y, pt->phi, &w))
				checkit(&w, 0);
		}
	}
	printf("checks done\n");
	fflush(stdout);
}
L(x, y, phi, wp, code)
double x, y, phi;
struct way *wp;
{
	struct way xxx;
	double f, d;
	struct tab *p, *pbest;


	d = x*x + y*y + 1000;
	pbest = 0;

	xxx = *wp;
	for(p=tab; p->s; p++) {
		if((*p->f)(x,y,phi, &xxx)) {
			f = xxx.value;

			if(debug) wshow(&xxx);

			if(pbest == 0 || f<d) {
				pbest = p;
				d = f;
				*wp = xxx;
			}
		}
	}
	if(pbest == 0) return 0;
	wp->name = pbest->ss[code];
	bestfreq = &pbest->freq;
	return 1;
}
report() {
	struct tab *p;

	for(p=tab; p->s; p++)
		if(p->freq)
			printf("%s used %d times\n", p->s, p->freq);
}
double _trash;
#define TRASH &_trash

double mod2pi(t)
double t;
{
	while (t < -PI) t += PI2;
	while (t >= PI) t -= PI2;
	return t;
}
rec2polar(x, y, pr, ptheta)
double x, y, *pr, *ptheta;
{
	*pr = hypot(x,y);
	*ptheta = mod2pi(atan2(y,x));
}
DECLARE(lsl) {
	double t,u;
	double xi, eta;

	xi = x - sin(phi);
	eta = y-1+cos(phi);
	rec2polar(xi,eta, &u, &t);
	RET(mod2pi(t), u, mod2pi(phi-t), 0.);
}
DECLARE(lsr1) {
	double t,u;
	double r, theta;
	double xi, eta;

	xi = x + sin(phi);
	eta = y-1-cos(phi);
	r = xi*xi + eta*eta;
	if(r < 4) return 0;
	u = sqrt(r-4);
	rec2polar(u, 2., TRASH, &theta);
	rec2polar(xi,eta, TRASH, &t);
	t += theta;
	RET(mod2pi(t), u, mod2pi(t-phi), 0.);
}
DECLARE(lsr2) {
	double t,u;
	double r, theta;
	double xi, eta;

	xi = x + sin(phi);
	eta = y-1-cos(phi);
	r = xi*xi + eta*eta;
	if(r < 4) return 0;
	u = -sqrt(r-4);
	rec2polar(u, 2., TRASH, &theta);
	rec2polar(xi,eta, TRASH, &t);
	t += theta;
	RET(mod2pi(t), u, mod2pi(t-phi), 0.);
}
DECLARE(lrs) {
	double t1, u1, v1, L1;
	double t2, u2, v2, L2;
	double zz;

	zz = (x*sin(phi) - (y-1)*cos(phi) + 1)/2;
	if(fabs(zz) >= 1) return 0;

	u1 = acos(zz);
	u2 = -u1;

	t1 = mod2pi(u1+phi);
	t2 = mod2pi(u2+phi);

	v1 = x*cos(phi)+(y-1)*sin(phi)-2*sin(u1);
	v2 = x*cos(phi)+(y-1)*sin(phi)-2*sin(u2);

	L1 = fabs(t1)+fabs(u1)+fabs(v1);
	L2 = fabs(t2)+fabs(u2)+fabs(v2);

	if(L1<L2) {
		RET(t1,u1,v1,0.);
	} else {
		RET(t2,u2,v2,0.);
	}
}
DECLARE(lrl)
{
	double xx, yy;
	double R, theta, u1, u2, as;
	double v1, v2, w1, w2;

	rec2polar(x-sin(phi), y-1+cos(phi), &R, &theta);

	if(R > 4.) return 0;

	as = asin(R/4);
	u1 = mod2pi(as + theta);
	u2 = mod2pi(PI-as+theta);
	if(u1>u2){
		double z;
		z = u1;
		u1 = u2;
		u2 = z;
	}
	xx = 2-sin(u1)*(x-sin(phi))+(y-1+cos(phi))*cos(u1);
	yy =   cos(u1)*(x-sin(phi))+(y-1+cos(phi))*sin(u1);
	rec2polar(xx, yy, TRASH, &v1);
	w1 = mod2pi(phi+v1-u1);

	xx = 2-sin(u2)*(x-sin(phi))+(y-1+cos(phi))*cos(u2);
	yy =   cos(u2)*(x-sin(phi))+(y-1+cos(phi))*sin(u2);
	rec2polar(xx, yy, TRASH, &v2);
	w2 = mod2pi(phi+v2-u2);

	if(fabs(u1)+fabs(v1)+fabs(w1) < fabs(u2)+fabs(v2)+fabs(w2)) {
		RET(u1, v1, w1, 0.);
	} else {
		RET(u2, v2, w2, 0.);
	}
}



DECLARE(lrlr)
{
	double u1, u2;
	double v1, v2;
	double w1, w2;
	double R, theta;
	double L1, L2;

	rec2polar(x+sin(phi), y-1-cos(phi), &R, &theta);
	if(fabs((20-R*R)/16) > 1.) return 0;

	v1 = acos((20-R*R)/16);
	v2 = -v1;

	rec2polar(
	    -4*sin(theta)+2*sin(theta+v1), 
	    4*cos(theta)-2*cos(theta+v1), 
	    TRASH, &u1);
	w1 = mod2pi(u1 - phi);

	rec2polar(
	    -4*sin(theta)+2*sin(theta+v2), 
	    4*cos(theta)-2*cos(theta+v2), 
	    TRASH, &u2);
	w2 = mod2pi(u2 - phi);

	L1 = fabs(u1)+2*fabs(v1)+fabs(w1);
	L2 = fabs(u2)+2*fabs(v2)+fabs(w2);

	if(L1 < L2) {
		RET(u1, v1, v1, w1);
	} else {
		RET(u2, v2, v2, w2);
	}
}
/*
 * -+++
 */
DECLARE(lrsla) {
	double t,u,v;
	double RR, xi, eta;

	xi = x - sin(phi);
	eta = y-1+cos(phi);
	RR = xi * xi + eta * eta;

	if(RR>4)
		v = 2 - sqrt(RR-4);
	else
		return 0;

	u = - PI/2;

	rec2polar((v-2)*eta - 2*xi, -xi*(v-2) -2*eta, TRASH, &t);

	RET(mod2pi(t),mod2pi(u),v,mod2pi(phi+u-t));
}
/*
 * +++-
 */
DECLARE(lrslb) {
	double t, u, v;
	double A, B, RR, xi, eta;


	xi = x - sin(phi);
	eta = y-1+cos(phi);
	RR = xi * xi + eta * eta;
	if(RR<16) return 0;
	v = sqrt(RR-16);
	u = 2*atan(2/v);

	A = v * sin(u) + 2 * (1-cos(u) );
	B = v * cos(u) + 2 * sin(u);

	rec2polar( (xi * B - eta * A),  (xi * A + eta * B), TRASH, &t);

	RET(mod2pi(t), mod2pi(u),v, mod2pi(phi+u-t));
}

/*
 * ++-+
 */
DECLARE(lrsl3) {
	double q, R, xi, eta, A, B;

	xi = x - sin(phi);
	eta = y-1+cos(phi);
	R = xi * xi + eta * eta;
	R = sqrt(R);
	q = (-R*R + R * sqrt(R*R+64) )/32;
	if(q<=0)return 0;
	q = sqrt(q);
	wp->u = 2*asin(q);
	wp->v = -4 * tan(wp->u/2);
	A = wp->v * sin(wp->u) + 2 * (1-cos(wp->u) );
	B = wp->v * cos(wp->u) + 2 * sin(wp->u);
	if((xi*B-eta*A) == 0) return 0;

	rec2polar(xi*B - eta*A, xi*A+eta*B, TRASH, &wp->t);
	wp->w = mod2pi(phi + wp->u - wp->t);
	wp->value = fabs(wp->t)+fabs(wp->u)+fabs(wp->v)+fabs(wp->w);
	return 1;
}
/*
 * ++--
 */
DECLARE(lrsl4) {
	double t,u,v;
	double eta, xi, A, B;

	xi = x - sin(phi);
	eta = y-1+cos(phi);

	v = -sqrt(xi*xi + eta*eta);
	u = 2*atan(-v/2);

	A = v * sin(u) + 2 * (1-cos(u) );
	B = v * cos(u) + 2 * sin(u);

	rec2polar(xi*B - eta*A, xi*A+eta*B, TRASH, &t);

	RET(t,u,v,mod2pi(phi+u-t));
}
/*
 * +-++
 */
DECLARE(lrsl5) {
	double u1, t1, v1, w1, L1;
	double u2, t2, v2, w2, L2;
	double eta, xi, R, A, B;

	xi = x - sin(phi);
	eta = y-1+cos(phi);
	R = sqrt(xi*xi+eta*eta);

	u1 = (2+R)/2;
	if(fabs(u1) >= 1) {
		u1 = -acos(1/u1);
		v1 = -2*tan(u1);
		A = v1 * sin(u1) + 2 * (1-cos(u1) );
		B = v1 * cos(u1) + 2 * sin(u1);
		rec2polar(xi*B - eta*A, xi*A+eta*B, TRASH, &t1);
		w1 = mod2pi(phi+u1-t1);
		L1 = fabs(t1)+fabs(u1)+fabs(v1)+fabs(w1);
	}
	else L1 = -1;
	u2 = (2-R)/2;
	if(fabs(u2) >= 1) {
		u2 = -acos(1/u2);
		v2 = -2*tan(u2);
		A = v2 * sin(u2) + 2 * (1-cos(u2) );
		B = v2 * cos(u2) + 2 * sin(u2);
		rec2polar(xi*B - eta*A, xi*A+eta*B, TRASH, &t2);
		w2 = mod2pi(phi+u2-t2);
		L2 = fabs(t2)+fabs(u2)+fabs(v2)+fabs(w2);
	}
	else L2 = -1;
	if(L1 < 0 && L2 < 0) return 0;
	if(L1<0) RET(t2,u2,v2,w2);
	if(L2<0) RET(t1,u1,v1,w1);
	if(L1<L2) { 
		RET(t1,u1,v1,w1);
	}
	else { 
		RET(t2,u2,v2,w2); 
	}
}

/*
 * 38. lrsr +---
 */
DECLARE(lrsr1) {
	double xi, eta, R;
	double t,u,v,w;

	xi = x+sin(phi);
	eta = y-1-cos(phi);
	R = sqrt(xi*xi + eta*eta);

	u = -PI/2;
	v = 2 - R;
	if(v == 2) return 0;
	rec2polar( eta/(v-2),  -xi/(v-2), TRASH, &t);
	w = mod2pi(t-u-phi);
	RET(t,u,v,w);
}
/*
 * 36. lrsr +-+-
 */
DECLARE(lrsr2) {
	double xi, eta, R;
	double t,u,v,w;

	xi = x+sin(phi);
	eta = y-1-cos(phi);
	R = sqrt(xi*xi + eta*eta);

	u = -PI/2;
	v = 2 + R;
	if(v == 2) return 0;
	rec2polar( eta/(v-2),  -xi/(v-2), TRASH, &t);
	w = mod2pi(t-u-phi);
	RET(t,u,v,w);
}
/*
 * 32. lrsr +++- (minus root)
 */
DECLARE(lrsr4) {
	double xi, eta, R;
	double t,u,v,w;
	double A, B;

	xi = x+sin(phi);
	eta = y-1-cos(phi);
	R = sqrt(xi*xi + eta*eta);

	u = 2*acos(sqrt( (R*R+16 - R*sqrt(R*R+32) )/32) );
	v = -2*tan (u/2);
	A = v*sin(u)+2;
	B = v*cos(u);
	rec2polar(xi*B-eta*A, xi*A+eta*B, TRASH, &t);

	w = mod2pi(t-u-phi);
	RET(t,u,v,w);
}
/*
 * 33. lrsr +-++ (cubic case)
 */
DECLARE(lrsr5) {
	double sq(), cb(), cbrt();
	double y1, z, p, q, p3, p32, q2, qq, alph, beta;
	double xi, eta, R;
	double t,u,v,w;
	double A, B;

	xi = x+sin(phi);
	eta = y-1-cos(phi);
	R = sqrt(xi*xi + eta*eta);

	p = -16/3. - R*R;
	q = 304./27 - 4*R*R/3.;
	p3 = cb(p/3);
	p32 = sqrt(fabs(p3) );
	q2 = sq(q/2);
	qq = p3+q2;

	if(p == 0) {
		y1 = cbrt(q);
		z = y1 + 4./3;
		if(z < 2) return 0;
	}
	else if(qq < 0) {
		alph = acos(-q/(2*p32) );
		y1 = 2 * sqrt(-p/3) * cos(alph/3);
 if(fabs(cb(y1)+p*y1+q) > .000001) printf("%f %f %f %f %f,TROUBLE0\n",x,y,phi,y1,qq);
/* if(2*cbrt(-p/3)*cos( (PI+alph) /3) >= 2./3) printf("%f %f %f %f %f,TROUBLE1\n",x,y,phi,y1,qq);
 if(2*cbrt(-p/3)*cos( (-PI+alph) /3) >= 2./3) printf("%f %f %f %f %f,TROUBLE2\n",x,y,phi,y1,qq);*/
		z = y1 + 4./3;
		if(z < 2) return 0;
	}
	if(qq >= 0 && p > 0) {
		beta = atan(2*p32/q);
		alph = atan( cbrt( tan(beta/2) ) );
		y1 = -2 * sqrt( p/3) / tan(2*alph);
 if(fabs(cb(y1)+p*y1+q) > .000001) printf("%f %f %f %f %f,TROUBLE3\n",x,y,phi,y1,qq);
		z = y1 + 4./3;
		if(z < 2) return 0;
	}
	if(qq >= 0 && p < 0) {
		beta = asin(2*p32/q);
		alph = atan( cbrt( tan(beta/2) ) );
		y1 = -2 * sqrt(-p/3) / sin(2*alph);
 if(fabs(cb(y1)+p*y1+q) > .000001) printf("%f %f %f %f %f,TROUBLE4\n",x,y,phi,y1,qq);
		z = y1 + 4./3;
		if(z < 2) return 0;
	}
	v = sqrt(z*z-4);
	u = -atan(v/2);
	A = v*sin(u)+2;
	B = v*cos(u);
	rec2polar(xi*B-eta*A, xi*A+eta*B, TRASH, &t);

	w = mod2pi(t-u-phi);
	RET(t,u,v,w);
}

wx() {
	double fin(), func();
	struct point A, C, B;
	struct way AB, AC, BC;
	double dbest;
	int n = 0;

	A.x = fin("A.x");
	A.y = fin("A.y");
	A.phi = fin("A.phi");

	C.x = fin("C.x");
	C.y = fin("C.y");
	C.phi = fin("C.phi");

	A.name = "A";
	B.name = "B";
	C.name = "C";

	dfunc(&A, &C, &AC);

	pshow(&A);
	pshow(&C);
	wshow(&AC);

	dbest = AC.value;

	for(n=0;;n++) {
		if(n%1000 == 0) {
			printf("n=%d\n", n);
			fflush(stdout);
		}
		randpoint(&B);

		dfunc(&A, &B, &AB);
		dfunc(&B, &C, &BC);
		if(dbest > AB.value + BC.value) {
			printf("n=%d\n", n);
			dbest = AB.value + BC.value;
			tritest(&AB, &BC, &AC);
		}
	}
}
double fin(s)
char *s;
{
	double x;
	do 
		printf("%s: ", s); 
	while(scanf("%lf", &x) != 1);
	return x;
}
double sq(x)
double x;
{
	return x*x;
}
double cb(x)
double x;
{
	return x*x*x;
}
double cbrt(x)
double x;
{
	if(x >  .000001) return  exp( log( x)/3.);
	if(x < -.000001) return -exp( log(-x)/3.);
	return 0;
}