[shkola] Re: Malko poqsneniq

> Za tezi nenadareni s prorocheski sposobnosti (i ne chetqshti edin
> drug maillist) malko razasneniq za posledniq post na Ljudo. Dolu e
> citiran mail na Ivo.
Az sledq wxprosniqt drug mailing list, no neshto ne shwanah za kakxw 
tochno posleden post na Ludo stawa wxpros.
Anyway, prashtam wi source-a. Mnogo se sxmnqwam che na nqkoi shte mu 
trqbwa wxobshte da go otwarq ili che sled kato go otwori shte razbere 
kakwo sxm mislil w 4 chasa prez noshtta, no nishto.

Otkirh edin interesen test, koito cql den se opitwam da si go obqsnq.
Programata mi wadi weren radius, no razlichni koordinati na centxra. Ne 
moga da si go obqsnq. Ako namiram prawilniqt radius i s krxga s centxr, 
kojto izwezhdam kato reshenie, moga da opisha tochkite. Znachi towa 
sxshto e reshenie! Ot druga strana ne moga da si obqsnq kakxw shte e 
tozi mnogoxgxlnik, kojto da mozhesh da go opishesh po dwa nachina s 
okrxzhnost s minimalen radius! Qsno e che ako ima nqkakwa greshka tq 
trqbwa da e w men.

Ako nqkoi izmisli neshto, da se obazhda!
-- 
Ivaylo Riskov <ivaylo_riskov@xxxxxxx>

"How much wood would a woodchuck chuck, 
if a woodchuck would chuck wood?"

#include <stdio.h>
#include <math.h>

typedef struct {
        double x;
        double y;
        double dist;
        double angle;
} SPoint;

SPoint pt[100024];
int n;
int hull[100024];
int ptr;
int x, y;
double best;
double ox, oy;

int cmp(double x, double y) {
        if (x - y > 1e-8)
                return 1;
        if (y - x > 1e-8)
                return -1;
        return 0;
}

int fcmp(const void *e1, const void *e2) {
        register SPoint *s1, *s2;
        int t;
        
        s1 = (SPoint *)e1;
        s2 = (SPoint *)e2;
        t = cmp(s1->angle, s2->angle); 
        if (t)
                return t;
        return cmp(s2->dist, s1->dist);
}

double get_dist(double dx, double dy) {
        return sqrt(dx*dx + dy*dy);
}

int is_right_oriented(SPoint *a, SPoint *b, SPoint *c) {
        double x;
        
        x = a->x*(b->y - c->y) + a->y*(c->x - b->x) + b->x*c->y - b->y*c->x;
        return cmp(x, 0.0) < 0;
}

void convex_hull() {
        int first, i;
        double angle;
        
        first = 0;
        for (i = 1; i < n; i++)
                if (cmp(pt[first].y, pt[i].y) > 0 ||
                   (cmp(pt[first].y, pt[i].y) == 0 && 
                        cmp(pt[first].x, pt[i].x) > 0))
                   first = i;

        pt[first].angle = -1.0;
        for (i = 0; i < n; i++)
                if (i != first) {
                        angle = atan2(pt[i].y - pt[first].y, pt[i].x - 
pt[first].x);
                        if (angle < 0)
                                angle += 4*acos(0);
                        pt[i].angle = angle;
                        pt[i].dist = get_dist(pt[i].y - pt[first].y, pt[i].x - 
pt[first].x);
                }
        
        qsort(pt, n, sizeof(SPoint), fcmp);

        hull[0] = 0;
        hull[1] = 1;
//      hull[2] = 2;
        ptr = 2;
        for (i = 2; i < n; i++) 
                if (cmp(pt[i].angle, pt[i - 1].angle)) {
                        while (is_right_oriented(pt + hull[ptr - 2], 
                                                                         pt + 
hull[ptr - 1], pt + i))
                                ptr--;
                        hull[ptr++] = i;
                }
        hull[ptr] = hull[0];
}

void get_farthest() {
        double d, t;
        int i, m;

        x = 0; y = 1;
        for (i = 2; i < ptr; i++)
                if (cmp(pt[ hull[i] ].dist, pt[ hull[y] ].dist) > 0)
                        y = i;
        best = pt[ hull[y] ].dist;

        m = y;
        for (i = 1; i < ptr; i++) {
                d = get_dist(pt[hull[i]].x - pt[hull[m]].x, 
                                         pt[hull[i]].y - pt[hull[m]].y);
                while (1) {
                        t = get_dist(pt[hull[i]].x - pt[hull[m + 1]].x, 
                                                 pt[hull[i]].y - pt[hull[m + 
1]].y);
                        if (cmp(d, t) > 0)
                                break;
                        d = t;
                        m++;
                        if (m == ptr)
                                m = 0;
                }
                if (cmp(d, best) > 0)
                        best = d, x = i, y = m;
        }
}

double calc_angle(SPoint *p) {
        double a, b;

        a = get_dist(pt[hull[x]].x - p->x, pt[hull[x]].y - p->y);
        b = get_dist(pt[hull[y]].x - p->x, pt[hull[y]].y - p->y);
        return acos((a*a + b*b - best*best)/(2.0*a*b));
}

void get_best() {
        int i;
        double angle, t;
        double r;
                
        angle = acos(0);
        for (i = 0; i < ptr; i++) 
                if (i != x && i != y) {
                        t = calc_angle(pt + hull[i]);
                        if (cmp(t, angle) < 0)
                                angle = t;
                }
        r = best/(2.0*sin(angle));
        angle = acos(0) - angle;
        angle += atan2(pt[hull[x]].y - pt[hull[y]].y, 
                                   pt[hull[x]].x - pt[hull[y]].x);
        ox = r*cos(angle) + pt[hull[y]].x;
        oy = r*sin(angle) + pt[hull[y]].y;
        printf("%.10lf\n", r);
        printf("%.10lf %.10lf\n", ox, oy);
}

int main() {
        int i;
        
        scanf("%i", &n);
        for (i = 0; i < n; i++)
                scanf("%lf %lf", &pt[i].x, &pt[i].y);

        if (n == 2) {
                printf("%.10lf\n", get_dist(pt[0].x - pt[1].x, pt[0].y - 
pt[1].y));
                printf("%.10lf %.10lf\n", (pt[0].x + pt[1].x)/2.0, 
                           (pt[0].y + pt[1].y)/2.0);
                return 0;
        }

        convex_hull();
        
        get_farthest();

        get_best();
        return 0;
}

Other related posts: