#include "basic.h"
#include <math.h>
#include <string.h>	/* for memset() */

int inserted(int segnum, int whichpt);
int math_logstar_n(int n);
int math_N(int n, int h);

static int q_idx;
static int tr_idx;


/* Return a new node to be added into the query tree */
static int newnode()
{
    if (q_idx < QSIZE)
        return q_idx++;
    else
    {
        fprintf(stderr, "newnode: Query-table overflow\n");
        return -1;
    }
}

/* Return a free trapezoid */
static int newtrap()
{
    if (tr_idx < TRSIZE)
    {
        tr[tr_idx].lseg = -1;
        tr[tr_idx].rseg = -1;
        tr[tr_idx].state = ST_VALID;
        return tr_idx++;
    }
    else
    {
        fprintf(stderr, "newtrap: Trapezoid-table overflow\n");
        return -1;
    }
}


/* Return the maximum of the two points into the yval structure */
static int _max(yval, v0, v1)
point_t *yval;
point_t *v0;
point_t *v1;
{
    if (v0->y > v1->y + C_EPS)
        *yval = *v0;
        else if (FP_EQUAL(v0->y, v1->y))
        {
            if (v0->x > v1->x + C_EPS)
                *yval = *v0;
            else
                *yval = *v1;
        }
        else
            *yval = *v1;
            
            return 0;
}


/* Return the minimum of the two points into the yval structure */
static int _min(yval, v0, v1)
point_t *yval;
point_t *v0;
point_t *v1;
{
    if (v0->y < v1->y - C_EPS)
        *yval = *v0;
        else if (FP_EQUAL(v0->y, v1->y))
        {
            if (v0->x < v1->x)
                *yval = *v0;
            else
                *yval = *v1;
        }
        else
            *yval = *v1;
            
            return 0;
}


int _greater_than(point_t* v0, point_t* v1)
{
    if (v0->y > v1->y + C_EPS)
        return TRUE;
    else if (v0->y < v1->y - C_EPS)
        return FALSE;
    else
        return (v0->x > v1->x);
}


int _equal_to(point_t* v0, point_t* v1)
{
    return (FP_EQUAL(v0->y, v1->y) && FP_EQUAL(v0->x, v1->x));
}

int _greater_than_equal_to(point_t* v0, point_t* v1)
{
    if (v0->y > v1->y + C_EPS)
        return TRUE;
    else if (v0->y < v1->y - C_EPS)
        return FALSE;
    else
        return (v0->x >= v1->x);
}

int _less_than(point_t* v0, point_t* v1)
{
    if (v0->y < v1->y - C_EPS)
        return TRUE;
    else if (v0->y > v1->y + C_EPS)
        return FALSE;
    else
        return (v0->x < v1->x);
}


/* Initilialise the query structure (Q) and the trapezoid table (T)
 * when the first segment is added to start the trapezoidation
 */
int init_query_structure(segnum)
int segnum;
{
    int i1, i2, i3, i4, i5, i6, i7, root;
    int t1, t2, t3, t4;
    segment_t *s = &seg[segnum];
    
    memset((void *)tr, 0, sizeof(tr));
    memset((void *)qs, 0, sizeof(qs));
    
    i1 = newnode();
    qs[i1].nodetype = T_Y;
    _max(&qs[i1].yval, &s->v0, &s->v1); /* root */
    root = i1;
    
    qs[i1].right = i2 = newnode();
    qs[i2].nodetype = T_SINK;
    qs[i2].parent = i1;
    
    qs[i1].left = i3 = newnode();
    qs[i3].nodetype = T_Y;
    _min(&qs[i3].yval, &s->v0, &s->v1); /* root */
    qs[i3].parent = i1;
    
    qs[i3].left = i4 = newnode();
    qs[i4].nodetype = T_SINK;
    qs[i4].parent = i3;
    
    qs[i3].right = i5 = newnode();
    qs[i5].nodetype = T_X;
    qs[i5].segnum = segnum;
    qs[i5].parent = i3;
    
    qs[i5].left = i6 = newnode();
    qs[i6].nodetype = T_SINK;
    qs[i6].parent = i5;
    
    qs[i5].right = i7 = newnode();
    qs[i7].nodetype = T_SINK;
    qs[i7].parent = i5;
    
    t1 = newtrap();		/* middle left */
    t2 = newtrap();		/* middle right */
    t3 = newtrap();		/* bottom-most */
    t4 = newtrap();		/* topmost */
    
    tr[t1].hi = tr[t2].hi = tr[t4].lo = qs[i1].yval;
    tr[t1].lo = tr[t2].lo = tr[t3].hi = qs[i3].yval;
    tr[t4].hi.y = (double) (INFINITY);
    tr[t4].hi.x = (double) (INFINITY);
    tr[t3].lo.y = (double) -1* (INFINITY);
    tr[t3].lo.x = (double) -1* (INFINITY);
    tr[t1].rseg = tr[t2].lseg = segnum;
    tr[t1].u0 = tr[t2].u0 = t4;
    tr[t1].d0 = tr[t2].d0 = t3;
    tr[t4].d0 = tr[t3].u0 = t1;
    tr[t4].d1 = tr[t3].u1 = t2;
    
    tr[t1].sink = i6;
    tr[t2].sink = i7;
    tr[t3].sink = i4;
    tr[t4].sink = i2;
    
    tr[t1].state = tr[t2].state = ST_VALID;
    tr[t3].state = tr[t4].state = ST_VALID;
    
    qs[i2].trnum = t4;
    qs[i4].trnum = t3;
    qs[i6].trnum = t1;
    qs[i7].trnum = t2;
    
    s->is_inserted = TRUE;
    return root;
}


/* Retun TRUE if the vertex v is to the left of line segment no.
 * segnum
 */

static int is_left_of(segnum, v)
int segnum;
point_t *v;
{
    segment_t *s = &seg[segnum];
    double area;
    
    if (_greater_than(&s->v1, &s->v0)) /* seg. going upwards */
    {
        if (FP_EQUAL(s->v1.y, v->y))
        {
            if (v->x < s->v1.x)
                area = 1.0;
            else
                area = -1.0;
        }
        else if (FP_EQUAL(s->v0.y, v->y))
        {
            if (v->x < s->v0.x)
                area = 1.0;
            else
                area = -1.0;
        }
        else
            area = CROSS(s->v0, s->v1, (*v));
    }
    else				/* v0 > v1 */
    {
        if (FP_EQUAL(s->v1.y, v->y))
        {
            if (v->x < s->v1.x)
                area = 1.0;
            else
                area = -1.0;
        }
        else if (FP_EQUAL(s->v0.y, v->y))
        {
            if (v->x < s->v0.x)
                area = 1.0;
            else
                area = -1.0;
        }
        else
            area = CROSS(s->v1, s->v0, (*v));
    }
    
    if (area > 0.0)
        return TRUE;
    else
        return FALSE;
}


int is_collinear(segnum, v, is_swapped)
int segnum;
point_t *v;
int is_swapped;
{
    int n;
    
    /* First check if the endpoint is already inserted */
    if (!is_swapped)
        n = MODULO_NEXT(segnum + 1, global.nseg);
        else
            if ((n = segnum - 1) == 0)
                n = 1;
                
                return seg[n].is_inserted;
}


/* This is query routine which determines which trapezoid does the
 * point v lie in. The return value is the trapezoid number
 */

int locate_endpoint(point_t* v, point_t* vo, int r)
{
    node_t *rptr = &qs[r];
    
    switch (rptr->nodetype)
    {
        case T_SINK:
            return rptr->trnum;
            
        case T_Y:
            if (_greater_than(v, &rptr->yval)) /* above */
                return locate_endpoint(v, vo, rptr->right);
            else if (_equal_to(v, &rptr->yval)) /* the point is already */
            {			          /* inserted. */
                if (_greater_than(vo, &rptr->yval)) /* above */
                    return locate_endpoint(v, vo, rptr->right);
                else
                    return locate_endpoint(v, vo, rptr->left); /* below */
            }
            else
                return locate_endpoint(v, vo, rptr->left); /* below */
            
        case T_X:
            if (_equal_to(v, &seg[rptr->segnum].v0) ||
                _equal_to(v, &seg[rptr->segnum].v1))
            {
                if (FP_EQUAL(v->y, vo->y)) /* horizontal segment */
                {
                    if (vo->x < v->x)
                        return locate_endpoint(v, vo, rptr->left); /* left */
                    else
                        return locate_endpoint(v, vo, rptr->right); /* right */
                }
                
                else if (is_left_of(rptr->segnum, vo))
                    return locate_endpoint(v, vo, rptr->left); /* left */
                else
                    return locate_endpoint(v, vo, rptr->right); /* right */
            }
            else if (is_left_of(rptr->segnum, v))
                return locate_endpoint(v, vo, rptr->left); /* left */
            else
                return locate_endpoint(v, vo, rptr->right); /* right */
            
        default:
            fprintf(stderr, "Haggu !!!!!\n");
            return 0;
    }
}


/* Thread in the segment into the existing trapezoidation. The
 * limiting trapezoids are given by tfirst and tlast (which are the
 * trapezoids containing the two endpoints of the segment
 */

int merge_trapezoids(segnum, tfirst, tlast, side)
int segnum;
int tfirst;
int tlast;
int side;
{
    int t, tnext, cond;
    int ptnext;
    
    /* First merge polys on the LHS */
    t = tfirst;
    while ((t > 0) && _greater_than_equal_to(&tr[t].lo, &tr[tlast].lo))
    {
        if (side == S_LEFT)
            cond = ((((tnext = tr[t].d0) > 0) && (tr[tnext].rseg == segnum)) ||
                    (((tnext = tr[t].d1) > 0) && (tr[tnext].rseg == segnum)));
        else
            cond = ((((tnext = tr[t].d0) > 0) && (tr[tnext].lseg == segnum)) ||
                    (((tnext = tr[t].d1) > 0) && (tr[tnext].lseg == segnum)));
        
        if (cond)
        {
            if ((tr[t].lseg == tr[tnext].lseg) &&
                (tr[t].rseg == tr[tnext].rseg)) /* good neighbours */
            {			              /* merge them */
                /* Use the upper node as the new node i.e. t */
                
                ptnext = qs[tr[tnext].sink].parent;
                
                if (qs[ptnext].left == tr[tnext].sink)
                    qs[ptnext].left = tr[t].sink;
                else
                    qs[ptnext].right = tr[t].sink;	/* redirect parent */
                
                
                /* Change the upper neighbours of the lower trapezoids */
                
                if ((tr[t].d0 = tr[tnext].d0) > 0) {
                    if (tr[tr[t].d0].u0 == tnext)
                        tr[tr[t].d0].u0 = t;
                    else if (tr[tr[t].d0].u1 == tnext)
                        tr[tr[t].d0].u1 = t;
                }
                
                if ((tr[t].d1 = tr[tnext].d1) > 0) {
                    if (tr[tr[t].d1].u0 == tnext)
                        tr[tr[t].d1].u0 = t;
                    else if (tr[tr[t].d1].u1 == tnext)
                        tr[tr[t].d1].u1 = t;
                }
                
                tr[t].lo = tr[tnext].lo;
                tr[tnext].state = ST_INVALID; /* invalidate the lower */
				            /* trapezium */
            }
            else		    /* not good neighbours */
                t = tnext;
        }
        else		    /* do not satisfy the outer if */
            t = tnext;
        
    } /* end-while */
    
    return 0;
}


/* Add in the new segment into the trapezoidation and update Q and T
 * structures
 */
int add_segment(int segnum)
{
    segment_t s;
    int tu, tl, sk, tfirst, tlast, tnext;
    int tfirstr, tlastr, tfirstl, tlastl;
    int i1, i2, t, tn;
    point_t vper, tpt;
    int tritop = 0, tribot = 0, is_swapped = 0;
    int tmptriseg;
    
    s = seg[segnum];
    if (_greater_than(&s.v1, &s.v0)) /* Get higher vertex in v0 */
    {
        int tmp;
        tpt = s.v0;
        s.v0 = s.v1;
        s.v1 = tpt;
        tmp = s.root0;
        s.root0 = s.root1;
        s.root1 = tmp;
        is_swapped = TRUE;
    }
    
    if ((is_swapped) ? !inserted(segnum, LASTPT) :
        !inserted(segnum, FIRSTPT))     /* insert v0 in the tree */
    {
        int tmp_d;
        
        tu = locate_endpoint(&s.v0, &s.v1, s.root0);
        tl = newtrap();		/* tl is the new lower trapezoid */
        tr[tl].state = ST_VALID;
        tr[tl] = tr[tu];
        tr[tu].lo.y = tr[tl].hi.y = s.v0.y;
        tr[tu].lo.x = tr[tl].hi.x = s.v0.x;
        tr[tu].d0 = tl;
        tr[tu].d1 = 0;
        tr[tl].u0 = tu;
        tr[tl].u1 = 0;
        
        if (((tmp_d = tr[tl].d0) > 0) && (tr[tmp_d].u0 == tu))
            tr[tmp_d].u0 = tl;
        if (((tmp_d = tr[tl].d0) > 0) && (tr[tmp_d].u1 == tu))
            tr[tmp_d].u1 = tl;
        
        if (((tmp_d = tr[tl].d1) > 0) && (tr[tmp_d].u0 == tu))
            tr[tmp_d].u0 = tl;
        if (((tmp_d = tr[tl].d1) > 0) && (tr[tmp_d].u1 == tu))
            tr[tmp_d].u1 = tl;
        
        /* Now update the query structure and obtain the sinks for the */
        /* two trapezoids */
        
        i1 = newnode();		/* Upper trapezoid sink */
        i2 = newnode();		/* Lower trapezoid sink */
        sk = tr[tu].sink;
        
        qs[sk].nodetype = T_Y;
        qs[sk].yval = s.v0;
        qs[sk].segnum = segnum;	/* not really reqd ... maybe later */
        qs[sk].left = i2;
        qs[sk].right = i1;
        
        qs[i1].nodetype = T_SINK;
        qs[i1].trnum = tu;
        qs[i1].parent = sk;
        
        qs[i2].nodetype = T_SINK;
        qs[i2].trnum = tl;
        qs[i2].parent = sk;
        
        tr[tu].sink = i1;
        tr[tl].sink = i2;
        tfirst = tl;
    }
    else				/* v0 already present */
    {       /* Get the topmost intersecting trapezoid */
        vper.x = s.v0.x + EPS * (s.v1.x - s.v0.x);
        vper.y = s.v0.y + EPS * (s.v1.y - s.v0.y);
        tfirst = locate_endpoint(&s.v0, &s.v1, s.root0);
        tritop = 1;
    }
    
    
    if ((is_swapped) ? !inserted(segnum, FIRSTPT) :
        !inserted(segnum, LASTPT))     /* insert v1 in the tree */
    {
        int tmp_d;
        
        tu = locate_endpoint(&s.v1, &s.v0, s.root1);
        
        tl = newtrap();		/* tl is the new lower trapezoid */
        tr[tl].state = ST_VALID;
        tr[tl] = tr[tu];
        tr[tu].lo.y = tr[tl].hi.y = s.v1.y;
        tr[tu].lo.x = tr[tl].hi.x = s.v1.x;
        tr[tu].d0 = tl;
        tr[tu].d1 = 0;
        tr[tl].u0 = tu;
        tr[tl].u1 = 0;
        
        if (((tmp_d = tr[tl].d0) > 0) && (tr[tmp_d].u0 == tu))
            tr[tmp_d].u0 = tl;
        if (((tmp_d = tr[tl].d0) > 0) && (tr[tmp_d].u1 == tu))
            tr[tmp_d].u1 = tl;
        
        if (((tmp_d = tr[tl].d1) > 0) && (tr[tmp_d].u0 == tu))
            tr[tmp_d].u0 = tl;
        if (((tmp_d = tr[tl].d1) > 0) && (tr[tmp_d].u1 == tu))
            tr[tmp_d].u1 = tl;
        
        /* Now update the query structure and obtain the sinks for the */
        /* two trapezoids */
        
        i1 = newnode();		/* Upper trapezoid sink */
        i2 = newnode();		/* Lower trapezoid sink */
        sk = tr[tu].sink;
        
        qs[sk].nodetype = T_Y;
        qs[sk].yval = s.v1;
        qs[sk].segnum = segnum;	/* not really reqd ... maybe later */
        qs[sk].left = i2;
        qs[sk].right = i1;
        
        qs[i1].nodetype = T_SINK;
        qs[i1].trnum = tu;
        qs[i1].parent = sk;
        
        qs[i2].nodetype = T_SINK;
        qs[i2].trnum = tl;
        qs[i2].parent = sk;
        
        tr[tu].sink = i1;
        tr[tl].sink = i2;
        tlast = tu;
    }
    else				/* v1 already present */
    {       /* Get the lowermost intersecting trapezoid */
        vper.x = s.v1.x + EPS * (s.v0.x - s.v1.x);
        vper.y = s.v1.y + EPS * (s.v0.y - s.v1.y);
        tlast = locate_endpoint(&s.v1, &s.v0, s.root1);
        tribot = 1;
    }
    
    /* Thread the segment into the query tree creating a new X-node */
    /* First, split all the trapezoids which are intersected by s into */
    /* two */
    
    t = tfirst;			/* topmost trapezoid */
    
    while ((t > 0) &&
           _greater_than_equal_to(&tr[t].lo, &tr[tlast].lo))
				/* traverse from top to bot */
    {
        int t_sav, tn_sav;
        sk = tr[t].sink;
        i1 = newnode();		/* left trapezoid sink */
        i2 = newnode();		/* right trapezoid sink */
        
        qs[sk].nodetype = T_X;
        qs[sk].segnum = segnum;
        qs[sk].left = i1;
        qs[sk].right = i2;
        
        qs[i1].nodetype = T_SINK;	/* left trapezoid (use existing one) */
        qs[i1].trnum = t;
        qs[i1].parent = sk;
        
        qs[i2].nodetype = T_SINK;	/* right trapezoid (allocate new) */
        qs[i2].trnum = tn = newtrap();
        tr[tn].state = ST_VALID;
        qs[i2].parent = sk;
        
        if (t == tfirst)
            tfirstr = tn;
        if (_equal_to(&tr[t].lo, &tr[tlast].lo))
            tlastr = tn;
        
        tr[tn] = tr[t];
        tr[t].sink = i1;
        tr[tn].sink = i2;
        t_sav = t;
        tn_sav = tn;
        
        /* error */
        
        if ((tr[t].d0 <= 0) && (tr[t].d1 <= 0)) /* case cannot arise */
        {
            fprintf(stderr, "add_segment: error\n");
            break;
        }
        
        /* only one trapezoid below. partition t into two and make the */
        /* two resulting trapezoids t and tn as the upper neighbours of */
        /* the sole lower trapezoid */
        
        else if ((tr[t].d0 > 0) && (tr[t].d1 <= 0))
        {			/* Only one trapezoid below */
            if ((tr[t].u0 > 0) && (tr[t].u1 > 0))
            {			/* continuation of a chain from abv. */
                if (tr[t].usave > 0) /* three upper neighbours */
                {
                    if (tr[t].uside == S_LEFT)
                    {
                        tr[tn].u0 = tr[t].u1;
                        tr[t].u1 = -1;
                        tr[tn].u1 = tr[t].usave;
                        
                        tr[tr[t].u0].d0 = t;
                        tr[tr[tn].u0].d0 = tn;
                        tr[tr[tn].u1].d0 = tn;
                    }
                    else		/* intersects in the right */
                    {
                        tr[tn].u1 = -1;
                        tr[tn].u0 = tr[t].u1;
                        tr[t].u1 = tr[t].u0;
                        tr[t].u0 = tr[t].usave;
                        
                        tr[tr[t].u0].d0 = t;
                        tr[tr[t].u1].d0 = t;
                        tr[tr[tn].u0].d0 = tn;
                    }
                    
                    tr[t].usave = tr[tn].usave = 0;
                }
                else		/* No usave.... simple case */
                {
                    tr[tn].u0 = tr[t].u1;
                    tr[t].u1 = tr[tn].u1 = -1;
                    tr[tr[tn].u0].d0 = tn;
                }
            }
            else
            {			/* fresh seg. or upward cusp */
                int tmp_u = tr[t].u0;
                int td0, td1;
                if (((td0 = tr[tmp_u].d0) > 0) &&
                    ((td1 = tr[tmp_u].d1) > 0))
                {		/* upward cusp */
                    if ((tr[td0].rseg > 0) &&
                        !is_left_of(tr[td0].rseg, &s.v1))
                    {
                        tr[t].u0 = tr[t].u1 = tr[tn].u1 = -1;
                        tr[tr[tn].u0].d1 = tn;
                    }
                    else		/* cusp going leftwards */
                    {
                        tr[tn].u0 = tr[tn].u1 = tr[t].u1 = -1;
                        tr[tr[t].u0].d0 = t;
                    }
                }
                else		/* fresh segment */
                {
                    tr[tr[t].u0].d0 = t;
                    tr[tr[t].u0].d1 = tn;
                }
            }
            
            if (FP_EQUAL(tr[t].lo.y, tr[tlast].lo.y) &&
                FP_EQUAL(tr[t].lo.x, tr[tlast].lo.x) && tribot)
            {		/* bottom forms a triangle */
                
                if (is_swapped)
                {
                    tmptriseg = segnum - 1;
                    if (tmptriseg == 0)
                        tmptriseg = global.nseg;
                }
                else
                    tmptriseg = MODULO_NEXT(segnum + 1, global.nseg);
                
                if ((tmptriseg > 0) && is_left_of(tmptriseg, &s.v0))
                {
                    /* L-R downward cusp */
                    tr[tr[t].d0].u0 = t;
                    tr[tn].d0 = tr[tn].d1 = -1;
                }
                else
                {
                    /* R-L downward cusp */
                    tr[tr[tn].d0].u1 = tn;
                    tr[t].d0 = tr[t].d1 = -1;
                }
            }
            else
            {
                if ((tr[tr[t].d0].u0 > 0) && (tr[tr[t].d0].u1 > 0))
                {
                    if (tr[tr[t].d0].u0 == t) /* passes thru LHS */
                    {
                        tr[tr[t].d0].usave = tr[tr[t].d0].u1;
                        tr[tr[t].d0].uside = S_LEFT;
                    }
                    else
                    {
                        tr[tr[t].d0].usave = tr[tr[t].d0].u0;
                        tr[tr[t].d0].uside = S_RIGHT;
                    }
                }
                tr[tr[t].d0].u0 = t;
                tr[tr[t].d0].u1 = tn;
            }
            
            t = tr[t].d0;
        }
        
        
        else if ((tr[t].d0 <= 0) && (tr[t].d1 > 0))
        {			/* Only one trapezoid below */
            if ((tr[t].u0 > 0) && (tr[t].u1 > 0))
            {			/* continuation of a chain from abv. */
                if (tr[t].usave > 0) /* three upper neighbours */
                {
                    if (tr[t].uside == S_LEFT)
                    {
                        tr[tn].u0 = tr[t].u1;
                        tr[t].u1 = -1;
                        tr[tn].u1 = tr[t].usave;
                        
                        tr[tr[t].u0].d0 = t;
                        tr[tr[tn].u0].d0 = tn;
                        tr[tr[tn].u1].d0 = tn;
                    }
                    else		/* intersects in the right */
                    {
                        tr[tn].u1 = -1;
                        tr[tn].u0 = tr[t].u1;
                        tr[t].u1 = tr[t].u0;
                        tr[t].u0 = tr[t].usave;
                        
                        tr[tr[t].u0].d0 = t;
                        tr[tr[t].u1].d0 = t;
                        tr[tr[tn].u0].d0 = tn;
                    }
                    
                    tr[t].usave = tr[tn].usave = 0;
                }
                else		/* No usave.... simple case */
                {
                    tr[tn].u0 = tr[t].u1;
                    tr[t].u1 = tr[tn].u1 = -1;
                    tr[tr[tn].u0].d0 = tn;
                }
            }
            else
            {			/* fresh seg. or upward cusp */
                int tmp_u = tr[t].u0;
                int td0, td1;
                if (((td0 = tr[tmp_u].d0) > 0) &&
                    ((td1 = tr[tmp_u].d1) > 0))
                {		/* upward cusp */
                    if ((tr[td0].rseg > 0) &&
                        !is_left_of(tr[td0].rseg, &s.v1))
                    {
                        tr[t].u0 = tr[t].u1 = tr[tn].u1 = -1;
                        tr[tr[tn].u0].d1 = tn;
                    }
                    else
                    {
                        tr[tn].u0 = tr[tn].u1 = tr[t].u1 = -1;
                        tr[tr[t].u0].d0 = t;
                    }
                }
                else		/* fresh segment */
                {
                    tr[tr[t].u0].d0 = t;
                    tr[tr[t].u0].d1 = tn;
                }
            }
            
            if (FP_EQUAL(tr[t].lo.y, tr[tlast].lo.y) &&
                FP_EQUAL(tr[t].lo.x, tr[tlast].lo.x) && tribot)
            {		/* bottom forms a triangle */
                int tmpseg;
                if (is_swapped)
                {
                    tmpseg = segnum - 1;
                    if (tmpseg == 0)
                        tmpseg = global.nseg;
                }
                else
                    tmpseg = MODULO_NEXT(segnum + 1, global.nseg);
                
                if ((tmpseg > 0) && is_left_of(tmpseg, &s.v0))
                {
                    /* L-R downward cusp */
                    tr[tr[t].d1].u0 = t;
                    tr[tn].d0 = tr[tn].d1 = -1;
                }
                else
                {
                    /* R-L downward cusp */
                    tr[tr[tn].d1].u1 = tn;
                    tr[t].d0 = tr[t].d1 = -1;
                }
            }
            else
            {
                if ((tr[tr[t].d1].u0 > 0) && (tr[tr[t].d1].u1 > 0))
                {
                    if (tr[tr[t].d1].u0 == t) /* passes thru LHS */
                    {
                        tr[tr[t].d1].usave = tr[tr[t].d1].u1;
                        tr[tr[t].d1].uside = S_LEFT;
                    }
                    else
                    {
                        tr[tr[t].d1].usave = tr[tr[t].d1].u0;
                        tr[tr[t].d1].uside = S_RIGHT;
                    }
                }
                tr[tr[t].d1].u0 = t;
                tr[tr[t].d1].u1 = tn;
            }
            
            t = tr[t].d1;
        }
        
        /* two trapezoids below. Find out which one is intersected by */
        /* this segment and proceed down that one */
        
        else
        {
            double y0, yt;
            point_t tmppt;
            int i_d0;
            
            i_d0 = FALSE;
            if (FP_EQUAL(tr[t].lo.y, s.v0.y))
            {
                if (tr[t].lo.x > s.v0.x)
                    i_d0 = TRUE;
            }
            else
            {
                tmppt.y = y0 = tr[t].lo.y;
                yt = (y0 - s.v0.y)/(s.v1.y - s.v0.y);
                tmppt.x = s.v0.x + yt * (s.v1.x - s.v0.x);
                
                if (_less_than(&tmppt, &tr[t].lo))
                    i_d0 = TRUE;
            }
            
            /* check continuity from the top so that the lower-neighbour */
            /* values are properly filled for the upper trapezoid */
            
            if ((tr[t].u0 > 0) && (tr[t].u1 > 0))
            {			/* continuation of a chain from abv. */
                if (tr[t].usave > 0) /* three upper neighbours */
                {
                    if (tr[t].uside == S_LEFT)
                    {
                        tr[tn].u0 = tr[t].u1;
                        tr[t].u1 = -1;
                        tr[tn].u1 = tr[t].usave;
                        
                        tr[tr[t].u0].d0 = t;
                        tr[tr[tn].u0].d0 = tn;
                        tr[tr[tn].u1].d0 = tn;
                    }
                    else		/* intersects in the right */
                    {
                        tr[tn].u1 = -1;
                        tr[tn].u0 = tr[t].u1;
                        tr[t].u1 = tr[t].u0;
                        tr[t].u0 = tr[t].usave;
                        
                        tr[tr[t].u0].d0 = t;
                        tr[tr[t].u1].d0 = t;
                        tr[tr[tn].u0].d0 = tn;		      
                    }
                    
                    tr[t].usave = tr[tn].usave = 0;
                }
                else		/* No usave.... simple case */
                {
                    tr[tn].u0 = tr[t].u1;
                    tr[tn].u1 = -1;
                    tr[t].u1 = -1;
                    tr[tr[tn].u0].d0 = tn;
                }
            }
            else 
            {			/* fresh seg. or upward cusp */
                int tmp_u = tr[t].u0;
                int td0, td1;
                if (((td0 = tr[tmp_u].d0) > 0) && 
                    ((td1 = tr[tmp_u].d1) > 0))
                {		/* upward cusp */
                    if ((tr[td0].rseg > 0) &&
                        !is_left_of(tr[td0].rseg, &s.v1))
                    {
                        tr[t].u0 = tr[t].u1 = tr[tn].u1 = -1;
                        tr[tr[tn].u0].d1 = tn;
                    }
                    else 
                    {
                        tr[tn].u0 = tr[tn].u1 = tr[t].u1 = -1;
                        tr[tr[t].u0].d0 = t;
                    }
                }
                else		/* fresh segment */
                {
                    tr[tr[t].u0].d0 = t;
                    tr[tr[t].u0].d1 = tn;
                }
            }
            
            if (FP_EQUAL(tr[t].lo.y, tr[tlast].lo.y) && 
                FP_EQUAL(tr[t].lo.x, tr[tlast].lo.x) && tribot)
            {
                /* this case arises only at the lowest trapezoid.. i.e.
                 tlast, if the lower endpoint of the segment is
                 already inserted in the structure */
                
                tr[tr[t].d0].u0 = t;
                tr[tr[t].d0].u1 = -1;
                tr[tr[t].d1].u0 = tn;
                tr[tr[t].d1].u1 = -1;
                
                tr[tn].d0 = tr[t].d1;
                tr[t].d1 = tr[tn].d1 = -1;
                
                tnext = tr[t].d1;	      
            }
            else if (i_d0)
            /* intersecting d0 */
            {
                tr[tr[t].d0].u0 = t;
                tr[tr[t].d0].u1 = tn;
                tr[tr[t].d1].u0 = tn;
                tr[tr[t].d1].u1 = -1;
                
                /* new code to determine the bottom neighbours of the */
                /* newly partitioned trapezoid */
                
                tr[t].d1 = -1;
                
                tnext = tr[t].d0;
            }
            else			/* intersecting d1 */
            {
                tr[tr[t].d0].u0 = t;
                tr[tr[t].d0].u1 = -1;
                tr[tr[t].d1].u0 = t;
                tr[tr[t].d1].u1 = tn;
                
                /* new code to determine the bottom neighbours of the */
                /* newly partitioned trapezoid */
                
                tr[tn].d0 = tr[t].d1;
                tr[tn].d1 = -1;
                
                tnext = tr[t].d1;
            }	    
            
            t = tnext;
        }
        
        tr[t_sav].rseg = tr[tn_sav].lseg  = segnum;
    } /* end-while */
    
    /* Now combine those trapezoids which share common segments. We can */
    /* use the pointers to the parent to connect these together. This */
    /* works only because all these new trapezoids have been formed */
    /* due to splitting by the segment, and hence have only one parent */
    
    tfirstl = tfirst; 
    tlastl = tlast;
    merge_trapezoids(segnum, tfirstl, tlastl, S_LEFT);
    merge_trapezoids(segnum, tfirstr, tlastr, S_RIGHT);
    
    seg[segnum].is_inserted = TRUE;
    return 0;
}


/* Update the roots stored for each of the endpoints of the segment.
 * This is done to speed up the location-query for the endpoint when
 * the segment is inserted into the trapezoidation subsequently
 */
static int find_new_roots(int segnum)
{
    segment_t *s = &seg[segnum];
    point_t vper;
    
    if (s->is_inserted)
        return 0;
    
    vper.x = s->v0.x + EPS * (s->v1.x - s->v0.x);
    vper.y = s->v0.y + EPS * (s->v1.y - s->v0.y);
    s->root0 = locate_endpoint(&s->v0, &s->v1, s->root0);
    s->root0 = tr[s->root0].sink;
    
    vper.x = s->v1.x + EPS * (s->v0.x - s->v1.x);
    vper.y = s->v1.y + EPS * (s->v0.y - s->v1.y);
    s->root1 = locate_endpoint(&s->v1, &s->v0, s->root1);
    s->root1 = tr[s->root1].sink;  
    return 0;
}

int choose_segment();

/* Main routine to perform trapezoidation */
void construct_trapezoids(int nseg, segment_t seg[])
{
    register int i;
    int root, h;
    
    q_idx = tr_idx = 1;
    /* Add the first segment and get the query structure and trapezoid */
    /* list initialised */
    root = init_query_structure(choose_segment());
    
#ifdef SIMPLE			/* no randomization */
    
    for (i = 1; i <= nseg; i++)
        seg[i].root0 = seg[i].root1 = root;
    
    for (i = 2; i <= nseg; i++)
        add_segment(choose_segment());
    
#else
    
    for (i = 1; i <= nseg; i++)
        seg[i].root0 = seg[i].root1 = root;
    
    for (h = 1; h <= math_logstar_n(nseg); h++)
    {
        for (i = math_N(nseg, h -1) + 1; i <= math_N(nseg, h); i++)
            add_segment(choose_segment());
        
        /* Find a new root for each of the segment endpoints */
        for (i = 1; i <= nseg; i++)
            find_new_roots(i);
    }
    
    for (i = math_N(nseg, math_logstar_n(nseg)) + 1; i <= nseg; i++)
        add_segment(choose_segment());
    
#endif
    
}