#define mark_edge_done(v,e) ((e) ->type|= (((v) ==get_vertex(e,0) ) ?0x200:0x400) ) 
#define edge_marked_done(v,e) ((e) ->type&(((v) ==get_vertex(e,0) ) ?0x200:0x400) )  \
 \

#define slower_and_finer 1 \

#define min(x,y) (((x) > (y) ) ?y:x) 
#define max(x,y) (((x) > (y) ) ?x:y) 
#define in_MorseExtract 1 \

#define is_xdeadend(s) ((s) ->type&0x1000) 
#define make_xdeadend(s) ((s) ->type|= 0x1000)  \

#define olist_next_free(l) ((l) ->free>=0) ?(l) ->free:list_count((l) ->entries) 
#define set_value(t,f) ((t) ->type|= 0x80,(t) ->h= f) 
#define smallest_vertex_in_edge(e) (vertex_compare(get_vertex(e,0) ,get_vertex(e,1) ) <0? \
get_vertex(e,0) :get_vertex(e,1) ) 
#define largest_vertex_in_edge(e) (vertex_compare(get_vertex(e,0) ,get_vertex(e,1) ) <0? \
get_vertex(e,1) :get_vertex(e,0) ) 
#define unmake_critical(s) (s) ->type&= 0xff87
#define make_critical(s) ((s) ->type&= 0xff87,(s) ->type|= 8,plist_push(crit[dimension(s) ],s) ) 
#define is_persistent(t) ((t) ->type&0x20) 
#define make_persistent(t) ((t) ->type|= 0x20) 
#define unmake_persistent(t) ((t) ->type&= 0xffdf) 
/*1:*/


/*18:*/


#include <stdio.h> 
#include "MorseExtract.h"

/*:18*/


/*19:*/


void abort_message(char*s);



/*:19*/


/*2:*/


list*crit[4];
long num_simp[4];
vertex*vertexlist;
edge*elist;
triangle*flist;
tetrahedron*tlist;


/*:2*/


/*84:*/


void list_initialize(list**l,unsigned int sz)
{
*l= (list*)malloc(sizeof(list));
if((*l)==NULL)abort_message("Out of memory");
(*l)->size= sz;
(*l)->padded_size= (sz+sizeof(long)-1)/sizeof(long);
(*l)->length= 0;
(*l)->body= NULL;
(*l)->body_length= 0;
}


/*:84*//*85:*/



void list_abandon(list**l)
{
if((*l)->body_length> 0)free((*l)->body);
free(*l);
*l= NULL;
}



/*:85*//*86:*/



void*list_push(list*l,void*q)
{
long*p;
if(l->body_length<=l->length)
{
l->body_length= l->length+1000;
l->body= (long*)realloc(l->body,l->body_length*l->padded_size*sizeof(long));
if(l->body==NULL)abort_message("out of memory");
}

p= list_entry(l,l->length);

memcpy(p,q,l->size);
l->length++;
return p;
}

void plist_push(list*l,void*q)
{
void*p= q;
list_push(l,&p);
}


/*:86*//*87:*/



void list_pop(list*l,void*q)
{
long*p;

if(l->length==0)abort_message("pop from empty list");

l->length--;
p= list_entry(l,l->length);
memcpy(q,p,l->size);
}

void*plist_pop(list*l)
{
void*p;

list_pop(l,&p);
return p;
}


/*:87*//*89:*/


void*list_read(list*l)
{
void*p,*pp;

if(l->read_index>=l->length)
{

l->length= l->read_index= l->read_delete_index;
return(NULL);
}

pp= list_entry(l,l->read_delete_index);


if(l->read_index!=l->read_delete_index)
{
p= list_entry(l,l->read_index);
memcpy(pp,p,l->size);
}

l->read_delete_index++;
l->read_index++;
return pp;
}

void*plist_read(list*l)
{
void**r;
r= list_read(l);
if(r!=NULL)return*r;
else return NULL;
}

/*:89*//*90:*/


void*list_read_insert(list*l,void*p)
{
if(l->read_index==l->read_delete_index)
list_push(l,p);
else
{
memcpy(list_entry(l,l->read_delete_index),p,l->size);
l->read_delete_index++;
}
}


/*:90*//*93:*/


void olist_initialize(olist**l,int sz)
{
*l= (olist*)malloc(sizeof(olist));
if((*l)==NULL)abort_message("Out of memory");
list_initialize(&((*l)->keys),sizeof(struct olist_key));
list_initialize(&((*l)->entries),sz);
(*l)->top= -1;
(*l)->free= -1;
(*l)->sz= sz;
}

/*:93*//*94:*/


void olist_abandon(olist**l)
{
list_abandon(&((*l)->keys));
list_abandon(&((*l)->entries));
free(*l);
*l= NULL;
}

/*:94*//*95:*/


void olist_clear(olist*l)
{
list_clear(l->keys);
list_clear(l->entries);
l->top= -1;
l->free= -1;
}


/*:95*//*96:*/


long olist_min(olist*l,void*p)
{
struct olist_key*q,*r;
long min_index;


min_index= l->top;
if(min_index<0)abort_message("min of empty list");
q= (struct olist_key*)list_entry(l->keys,min_index);
if(q->low>=0)
{
while(q->low>=0)
{
r= q;
min_index= q->low;
q= (struct olist_key*)list_entry(l->keys,min_index);
}
if(q->eq>=0)
{
r= q;
min_index= q->eq;
q= (struct olist_key*)list_entry(l->keys,min_index);
r->eq= q->eq;
}
else r->low= q->high;
}
else if(q->eq>=0)
{
r= q;
min_index= q->eq;
q= (struct olist_key*)list_entry(l->keys,min_index);
r->eq= q->eq;
}
else l->top= q->high;

memcpy(p,list_entry(l->entries,min_index),l->sz);

q->low= l->free;
l->free= min_index;

return q->k;
}


/*:96*//*97:*/



long olist_add(olist*l,long m,void*p)
{
struct olist_key*q,*r;
long b,*last;


b= olist_next_free(l);
/*99:*/


{
struct olist_key rr;

rr.high= rr.low= rr.eq= -1;
rr.k= m;

if(l->free>=0)
{
memcpy(list_entry(l->entries,b),p,l->sz);
r= (struct olist_key*)list_entry(l->keys,b);
l->free= r->low;
memcpy(r,&rr,sizeof(struct olist_key));
}
else
{
list_push(l->entries,p);
r= (struct olist_key*)list_push(l->keys,&rr);
}
}

/*:99*/



last= &(l->top);

while((*last)>=0)
{
q= (struct olist_key*)list_entry(l->keys,*last);
if(m<q->k)last= &(q->low);
else if(m> q->k)last= &(q->high);
else
{
r->eq= *last;
r->low= q->low;
r->high= q->high;
q->low= -1;
q->high= -1;
break;
}
}
*last= b;
return b;
}

/*:97*//*98:*/



long olist_find_add(olist*l,long m,void*p,int*flag)
{
struct olist_key*q,*r;
long b,*last;

last= &(l->top);

while((*last)>=0)
{
q= (struct olist_key*)list_entry(l->keys,*last);
if(m<q->k)last= &(q->low);
else if(m> q->k)last= &(q->high);
else
{
if(flag!=NULL)*flag= 1;
return*last;
}
}
if(p!=NULL)
{
b= *last= olist_next_free(l);
/*99:*/


{
struct olist_key rr;

rr.high= rr.low= rr.eq= -1;
rr.k= m;

if(l->free>=0)
{
memcpy(list_entry(l->entries,b),p,l->sz);
r= (struct olist_key*)list_entry(l->keys,b);
l->free= r->low;
memcpy(r,&rr,sizeof(struct olist_key));
}
else
{
list_push(l->entries,p);
r= (struct olist_key*)list_push(l->keys,&rr);
}
}

/*:99*/


}
else b= -1;

if(flag!=NULL)*flag= 0;
return b;
}


/*:98*/


/*102:*/



void get_triangle_vertices(triangle*t,vertex*vlist[3])

{
vertex*v;
edge*e;

e= get_edge(t,0);

get_edge_vertices(e,vlist);

e= get_edge(t,1);

v= get_vertex(e,1);
if(v!=vlist[0]&&v!=vlist[1])vlist[2]= v;
else vlist[2]= get_vertex(e,0);
}

void get_tetrahedron_vertices(tetrahedron*t,vertex*vlist[4])

{
vertex*vlistp[3];
int i;


get_triangle_vertices(get_face(t,0),vlist);
get_triangle_vertices(get_face(t,1),vlistp);

for(i= 0;i<3;i++)
{
if(vlistp[i]!=vlist[0]&&vlistp[i]!=vlist[1]&&vlistp[i]!=vlist[2])
{
vlist[3]= vlistp[i];
break;
}
}
}

/*:102*//*103:*/


int edge_orient(triangle*t,int i)
{
vertex*vl[2],*vlp[2];

if(i==0)return 1;
get_edge_vertices(get_edge(t,0),vl);
get_edge_vertices(get_edge(t,i),vlp);
if(vlp[0]==vl[1]||vlp[1]==vl[0])return 1;
return-1;
}


/*:103*//*104:*/


long vertex_compare(vertex*v,vertex*w)
{
if(value(v)<value(w))return-1;
else if(value(v)> value(w))return 1;
else return(long)v-(long)w;
}

/*:104*//*105:*/



int smallest_vertex(vertex*v[2],int n)
{
int i,j;
for(i= 0,j= 1;j<n;j++)
{
if(vertex_compare(v[j],v[i])<0)i= j;
}
return i;
}

int largest_vertex(vertex*v[2],int n)
{
int i,j;
for(i= 0,j= 1;j<n;j++)
{
if(vertex_compare(v[j],v[i])> 0)i= j;
}
return i;
}

/*:105*//*106:*/



int is_in_lower_Star(tetrahedron*t,vertex*v)
{
vertex*vlist[4];

get_tetrahedron_vertices(t,vlist);

return v==vlist[largest_vertex(vlist,4)];
}




/*:106*//*108:*/


int triangle_in_tetrahedron(triangle*s,tetrahedron*t)
{
int i;

for(i= 0;i<4;i++)
{
if(s==get_face(t,i))return i;
}
return-1;
}


/*:108*//*109:*/


edge*other_edge(vertex*v,edge*f,triangle*t)

{
int i;
edge*e;


for(i= 0;i<3;i++)
{
e= get_edge(t,i);
if(e==f)continue;
if(get_vertex(e,0)==v||get_vertex(e,1)==v)return e;
}

return NULL;
}

triangle*
other_face(edge*e,triangle*f,tetrahedron*t)

{
int i;
triangle*s;

for(i= 0;i<4;i++)
{
s= get_face(t,i);
if(f!=s&&edge_in_triangle(e,s)>=0)return s;
}
return NULL;
}



/*:109*//*110:*/



void clean_crit(void)
{
tetrahedron*t;
triangle*f;
edge*e;
vertex*v;
list*l;

l= crit[0];
list_read_init(l);
while((v= (vertex*)plist_read(l))!=NULL)
{
if(!is_critical(v))list_read_delete(l);
}
l= crit[1];
list_read_init(l);
while((e= (edge*)plist_read(l))!=NULL)
{
if(!is_critical(e))list_read_delete(l);
}
l= crit[2];
list_read_init(l);
while((f= (triangle*)plist_read(l))!=NULL)
{
if(!is_critical(f))list_read_delete(l);
}
l= crit[3];
list_read_init(l);
while((t= (tetrahedron*)plist_read(l))!=NULL)
{
if(!is_critical(t))list_read_delete(l);
}


}


/*:110*//*112:*/


void pair01(vertex*s0,edge*s1)
{
int n;

s0->type&= 0xFF87;
s1->type&= 0xFF87;

s0->links[0]= s1;
n= vertex_in_edge(s0,s1);
s1->type|= (16+(n<<5));
}

void pair12(edge*s0,triangle*s1)
{
int n;

s0->type&= 0xFF87;
s1->type&= 0xFF87;

s0->links[2]= s1;
n= edge_in_triangle(s0,s1);
s1->type|= (16+(n<<5));
}

void pair23(triangle*s0,tetrahedron*s1)
{
int n;

s0->type&= 0xFF87;
s1->type&= 0xFF87;

if(s0->links[4]==s1)s0->type|= 0x20;
n= triangle_in_tetrahedron(s0,s1);
s1->type|= (16+(n<<5));
}

/*:112*//*113:*/


int max_value(void*s,int d)
{
tetrahedron*te;
triangle*t;
edge*e;
int m0,m1;

if(s==NULL)return-0x8000;

switch(d)
{
case 3:
te= (tetrahedron*)s;
m0= value(get_face(te,0));
m1= value(get_face(te,1));
break;
case 2:
t= (triangle*)s;
m0= value(get_edge(t,0));
m1= value(get_edge(t,1));
break;
case 1:
e= (edge*)s;
m0= value(get_vertex(e,0));
m1= value(get_vertex(e,1));
break;
default:
abort_message("max value of simplex of unknown dimension");
break;
}
return max(m0,m1);
}

/*:113*//*114:*/


int min_value(void*s,int d)
{
tetrahedron*te;
triangle*t;
edge*e;
int m0,m1;

if(s==NULL)return 0x7fff;
switch(d)
{
case 3:
te= (tetrahedron*)s;
m0= min_value(get_face(te,0),2);
m1= min_value(get_face(te,1),2);
break;
case 2:
t= (triangle*)s;
m0= min_value(get_edge(t,0),1);
m1= min_value(get_edge(t,1),1);
break;
case 1:
e= (edge*)s;
m0= value(get_vertex(e,0));
m1= value(get_vertex(e,1));
break;
default:
abort_message("min value of simplex of unknown dimension");
break;
}
return min(m0,m1);
}


/*:114*/


/*115:*/


void abort_message(char*s)
{
printf("Error: %s",s);
exit(1);
}





/*:115*/


/*117:*/


vertex*new_vertex(long id,short val,short flag)
{
vertex*v;

v= id2vertex(id);

v->type= (flag)?0x84:0x80;
v->links[0]= NULL;
v->h= val;

return v;
}

/*:117*//*118:*/


edge*new_edge(long id,long vids[2],short flag)
{
edge*e;
short flagp;
vertex*v0,*v1;

e= id2edge(id);
v0= id2vertex(vids[0]);
v1= id2vertex(vids[1]);
e->links[0]= v0;
e->links[1]= v1;
e->links[2]= NULL;


e->h= max(v0->h,v1->h);

flagp= (flag&&is_in_K(v0)&&is_in_K(v1));


e->type= (flagp)?5:1;

if(flag> 0&&flagp==0)
abort_message("vertices of edge in K must be in K");

if(v0->links[0]==NULL)v0->links[0]= e;
if(v1->links[0]==NULL)v1->links[0]= e;

return e;
}

/*:118*//*119:*/


triangle*new_triangle(long id,long eids[3],short flag)
{
triangle*t;
vertex*vl[3],*v[2];
edge*e0,*e1,*e2;
short flagp;

e0= id2edge(eids[0]);
e1= id2edge(eids[1]);
e2= id2edge(eids[2]);

get_edge_vertices(e0,vl);
get_edge_vertices(e1,v);
if(vertex_in_edge(v[0],e0)<0)vl[2]= v[0];
else if(vertex_in_edge(v[1],e0)<0)vl[2]= v[1];
else abort_message("new triangle has more than three vertices");

get_edge_vertices(e2,v);
if(v[0]!=vl[0]&&v[0]!=vl[1]&&v[0]!=vl[2])
abort_message("new triangle has more than three vertices-1");
if(v[1]!=vl[0]&&v[1]!=vl[1]&&v[1]!=vl[2])
abort_message("new triangle has more than three vertices-2");


t= id2triangle(id);
t->links[0]= e0;
t->links[1]= e1;
t->links[2]= e2;

t->h= max(e0->h,e1->h);

flagp= (flag&&is_in_K(e0)&&is_in_K(e1)&&is_in_K(e2));

t->type= (flagp)?6:2;
t->links[3]= NULL;
t->links[4]= NULL;

if(e0->links[2]==NULL)e0->links[2]= t;
if(e1->links[2]==NULL)e1->links[2]= t;
if(e2->links[2]==NULL)e2->links[2]= t;

if(flag> 0&&flagp==0)
abort_message("edges of triangle in K must be in K");
return t;
}

/*:119*//*120:*/


tetrahedron*new_tetrahedron(long id,long fids[4],short flag)
{
tetrahedron*t;
vertex*vl[4],*vlp[3];
int i;
short flagp;
triangle*f0,*f1,*f2,*f3;

f0= id2triangle(fids[0]);
f1= id2triangle(fids[1]);
f2= id2triangle(fids[2]);
f3= id2triangle(fids[3]);

get_triangle_vertices(f0,vl);
get_triangle_vertices(f1,vlp);
for(i= 0;i<3;i++)
{
if(vlp[i]!=vl[0]&&vlp[i]!=vl[1]&&vlp[i]!=vl[2])
{
vl[3]= vlp[i];
break;
}
}
/*121:*/


for(i= 0;i<3;i++)
{
if(vlp[i]!=vl[0]&&vlp[i]!=vl[1]&&vlp[i]!=vl[2]&&vlp[i]!=vl[3])
abort_message("new tetrahedron has more than four vertices");
}






/*:121*/


get_triangle_vertices(f2,vlp);
/*121:*/


for(i= 0;i<3;i++)
{
if(vlp[i]!=vl[0]&&vlp[i]!=vl[1]&&vlp[i]!=vl[2]&&vlp[i]!=vl[3])
abort_message("new tetrahedron has more than four vertices");
}






/*:121*/


get_triangle_vertices(f3,vlp);
/*121:*/


for(i= 0;i<3;i++)
{
if(vlp[i]!=vl[0]&&vlp[i]!=vl[1]&&vlp[i]!=vl[2]&&vlp[i]!=vl[3])
abort_message("new tetrahedron has more than four vertices");
}






/*:121*/




t= id2tetrahedron(id);


t->links[0]= f0;
t->links[1]= f1;
t->links[2]= f2;
t->links[3]= f3;

t->h= max(f0->h,f1->h);

flagp= (flag&&is_in_K(f0)&&is_in_K(f1)&&is_in_K(f2)&&is_in_K(f3));

t->type= (flagp)?7:3;

if(f0->links[3]==NULL)f0->links[3]= t;
else if(f0->links[4]==NULL)f0->links[4]= t;
else abort_message("a triangle can be in only two tetrahedra");
if(f1->links[3]==NULL)f1->links[3]= t;
else if(f1->links[4]==NULL)f1->links[4]= t;
else abort_message("a triangle can be in only two tetrahedra-1");
if(f2->links[3]==NULL)f2->links[3]= t;
else if(f2->links[4]==NULL)f2->links[4]= t;
else abort_message("a triangle can be in only two tetrahedra-2");
if(f3->links[3]==NULL)f3->links[3]= t;
else if(f3->links[4]==NULL)f3->links[4]= t;
else abort_message("a triangle can be in only two tetrahedra-3");

if(flag> 0&&flagp==0)
abort_message("faces of tetrahedron in K must be in K");
return t;
}

/*:120*//*123:*/


void
read_in_complex(FILE*df)
{
long i,n[4];
int res;

/*124:*/


{
int val;
int vertex_in_K,c;

res= fscanf(df,"%ld\n",n);
if(res!=1)abort_message("bad number of vertices");


num_simp[0]= n[0];
vertexlist= (vertex*)malloc(n[0]*sizeof(vertex));
if(vertexlist==NULL)abort_message("Out of memory");

for(i= 0;i<n[0];i++)
{
c= getc(df);
if(c=='x')vertex_in_K= 0;
else
{
vertex_in_K= 1;
ungetc(c,df);
}
res= fscanf(df,"%ld\n",&val);
if(res!=1)abort_message("bad vertex");
new_vertex(i,val,vertex_in_K);
}
}

/*:124*/


/*125:*/


{
long ids[2];
int edge_in_K,c;

res= fscanf(df,"%ld\n",n+1);
if(res!=1)abort_message("bad number of edges");
num_simp[1]= n[1];

elist= (edge*)malloc(n[1]*sizeof(edge));
if(elist==NULL)abort_message("Out of memory");

for(i= 0;i<n[1];i++)
{
c= getc(df);
if(c=='x')edge_in_K= 0;
else
{
edge_in_K= -1;
ungetc(c,df);
}
res= fscanf(df,"%ld %ld\n",ids,ids+1);
if(res!=2)abort_message("bad edge");
new_edge(i,ids,edge_in_K);
}
}

/*:125*/


/*126:*/


{
long ids[3];
int face_in_K,c;

res= fscanf(df,"%ld\n",n+2);
if(res!=1)abort_message("bad number of triangles");
num_simp[2]= n[2];

flist= (triangle*)malloc(n[2]*sizeof(triangle));
if(flist==NULL)abort_message("Out of memory");

for(i= 0;i<n[2];i++)
{
c= getc(df);
if(c=='x')face_in_K= 0;
else
{
face_in_K= -1;
ungetc(c,df);
}
res= fscanf(df,"%ld %ld %ld\n",ids,ids+1,ids+2);
if(res!=3)abort_message("bad triangle");
new_triangle(i,ids,face_in_K);
}
}

/*:126*/


/*127:*/


{
long ids[4];
int tetrahedron_in_K,c;

res= fscanf(df,"%ld\n",n+3);
if(res!=1)abort_message("bad number of tetrahedra");
num_simp[3]= n[3];

tlist= (tetrahedron*)malloc(n[3]*sizeof(tetrahedron));
if(tlist==NULL)abort_message("Out of memory");

for(i= 0;i<n[3];i++)
{
c= getc(df);
if(c=='x')tetrahedron_in_K= 0;
else
{
tetrahedron_in_K= -1;
ungetc(c,df);
}
res= fscanf(df,"%ld %ld %ld %ld\n",ids,ids+1,ids+2,ids+3);
if(res!=4)abort_message("bad tetrahedron");
new_tetrahedron(i,ids,tetrahedron_in_K);
}
}

/*:127*/


}

/*:123*/


/*51:*/


vertex*FindGrad01(vertex*u,int m)
{
vertex*v;
edge*e;

v= u;
while(!is_critical(v)&&value(v)> m)
{
e= r01(v);
v= other_vertex_in_edge(v,e);
}
if(value(v)<=m)return NULL;
return v;
}

/*:51*//*52:*/


tetrahedron*FindGrad23(tetrahedron*tau,int m)

{
tetrahedron*t;
triangle*s;

if(tau==NULL)return NULL;
t= tau;
while(is_in_K(t)&&value(t)<m)
{
if(is_critical(t))return t;
s= r32(t);
if(is_deadend(s))
{
make_deadend(t);
return NULL;
}
t= other_coface(s,t);
if(is_deadend(t))
{
make_deadend(s);
return NULL;
}
}
if(!is_in_K(t)&&t!=tau)make_deadend(s);
return NULL;
}

/*:52*//*53:*/


edge*FindGradPaths12(triangle*sigma,int p,list*grad_path,int flags)
{
static olist*graph;
static list*to_do;
static list*crits;
static int first= 0;
struct edge_graph
{
long up;
long count;
}r,*q;
edge*e;
triangle*t;
long m;

/*54:*/


if(first==0)
{
olist_initialize(&graph,sizeof(struct edge_graph));
list_initialize(&crits,sizeof(long));
list_initialize(&to_do,sizeof(long));
first= 1;
}
else
{
olist_clear(graph);
list_clear(crits);
list_clear(to_do);
}


/*:54*/



m= -1;
e= NULL;
t= sigma;
/*55:*/


{
int i;
long n;
int flag;
edge*ep;
int livecount;

livecount= 0;

for(i= 0;i<3;i++)
{
ep= get_edge(t,i);
if(e!=ep&&((is_paired_up(ep)&&!is_deadend(ep))||is_critical(ep)))
{
livecount++;

if(value(ep)<=value(sigma)-p)continue;

r.up= m;
r.count= 0;

n= olist_find_add(graph,(long)ep,&r,&flag);
q= (struct edge_graph*)olist_entry(graph,n);
q->count++;

if(!flag)
{
if(is_critical(ep))list_push(crits,&n);
else list_push(to_do,&n);
}
}
}
if(livecount==0&&e!=NULL)make_deadend(e);
}

/*:55*/



while(!list_is_empty(to_do))
{
list_pop(to_do,&m);
e= (edge*)olist_get_key(graph,m);
t= r12(e);
/*55:*/


{
int i;
long n;
int flag;
edge*ep;
int livecount;

livecount= 0;

for(i= 0;i<3;i++)
{
ep= get_edge(t,i);
if(e!=ep&&((is_paired_up(ep)&&!is_deadend(ep))||is_critical(ep)))
{
livecount++;

if(value(ep)<=value(sigma)-p)continue;

r.up= m;
r.count= 0;

n= olist_find_add(graph,(long)ep,&r,&flag);
q= (struct edge_graph*)olist_entry(graph,n);
q->count++;

if(!flag)
{
if(is_critical(ep))list_push(crits,&n);
else list_push(to_do,&n);
}
}
}
if(livecount==0&&e!=NULL)make_deadend(e);
}

/*:55*/


}

if(flags&1)/*82:*/


{
edge*ee;

e= NULL;
while(!list_is_empty(crits))
{
list_pop(crits,&m);
ee= (edge*)olist_get_key(graph,m);
if(is_persistent(ee))
{
e= ee;
break;
}
}
}

/*:82*/


else/*56:*/


{
int best_value,val,bestm;

bestm= -1;

while(!list_is_empty(crits))
{
list_pop(crits,&m);
q= (struct edge_graph*)olist_entry(graph,m);
val= value((edge*)olist_get_key(graph,m));
if(q->count==1&&(bestm<0||val> best_value))
/*57:*/


{

while(q->up>=0)
{
q= (struct edge_graph*)olist_entry(graph,q->up);
if(q->count> 1)break;
}
if(q->count==1)
{
best_value= val;
bestm= m;
}
}

/*:57*/


}
if(bestm>=0)e= (edge*)olist_get_key(graph,bestm);
else e= NULL;
if(e!=NULL&&grad_path!=NULL)/*58:*/


{
list_clear(grad_path);

m= bestm;

while(m>=0)
{
plist_push(grad_path,(edge*)olist_get_key(graph,m));
m= ((struct edge_graph*)olist_entry(graph,m))->up;
}
}


/*:58*/



}

/*:56*/



return e;
}

/*:53*//*59:*/


int
splitrejoin(tetrahedron*t)
{
int i,j;
edge*e,*ep;
triangle*s;

for(i= 0;i<4;i++)
{
s= get_face(t,i);
if(is_paired_down(s))
{
e= r21(s);
for(j= 0;j<3;j++)
{
ep= get_edge(s,j);
if(ep!=e
&&(!is_paired_up(ep)||triangle_in_tetrahedron(r12(ep),t)<0))break;
}
if(j==3)return i;
}
}
return-1;
}

/*:59*//*60:*/


triangle*unsplitrejoin(tetrahedron*t,int n)
{
triangle*s,*sp;
edge*e;

s= get_face(t,n);
e= r21(s);
sp= other_face(e,s,t);

if(!is_in_K(t))return NULL;
if(is_critical(sp))return NULL;
if(is_paired_down(sp))
abort_message("tetrahedron surrounded by 2->1 triangles");
pair23(s,t);
pair12(e,sp);
return sp;
}

/*:60*//*62:*/



long
find_all_grad12_paths(triangle*sigma,olist*edges,int options)
{
int kk= -1;
edge*ep;
long todo= -2;
long this;
long crit= -1;
struct grad12_struct*p;
triangle*f;

olist_clear(edges);

do
{
this= todo;
if(this>=0)
{
p= (struct grad12_struct*)olist_entry(edges,this);
ep= id2edge(olist_get_key(edges,this));
kk= p->flags&3;
if(kk)
{
todo= p->links[0];
p->links[0]= -1;
}
else
{
todo= p->links[1];
p->links[1]= -1;
}
f= r12(ep);
}
else f= sigma;
/*63:*/


{
int k;
edge*e;
int livecount= 0;

for(k= 0;k<3;k++)
{
if(k==kk)continue;
e= get_edge(f,k);
/*64:*/


{
struct grad12_struct r,*q;
triangle*t;
int flag;
int j;
long id;

if(is_paired_up(e)&&(!is_deadend(e)||(options&3)==0))
{
livecount++;
id= olist_find_add(edges,edge_id(e),&r,&flag);
q= (struct grad12_struct*)olist_entry(edges,id);
if(flag)
{
p->links[k]= q->links[q->flags&3];
q->links[q->flags&3]= this;
}
else
{
q->links[1]= q->links[2]= -1;
t= r12(e);
q->flags= edge_in_triangle(e,t);
q->links[q->flags]= this;
if(q->flags)q->links[0]= todo;
else q->links[1]= todo;
todo= id;
}
}
else if(is_critical(e))
{
livecount++;
id= olist_find_add(edges,edge_id(e),&r,&flag);
q= (struct grad12_struct*)olist_entry(edges,id);
if(flag)
{
p->links[k]= q->links[0];
q->links[0]= this;
}
else
{
q->links[0]= this;
q->links[1]= crit;
q->flags= 3;
crit= id;
}
}
}

/*:64*/


}
if(livecount==0&&(options&3)&&this>=0)
{
make_deadend(f);
make_deadend(ep);
}
}


/*:63*/


}while(todo>=0);

if((options&3)==2)
/*65:*/


{
long j;
struct grad12_struct*p;
edge*e;
list*todo_list;

list_initialize(&todo_list,sizeof(long));

todo= crit;

while(todo>=0)
{
/*66:*/


{
long next,id;
struct grad12_struct*q;
int i;

list_clear(todo_list);
list_push(todo_list,&todo);

while(!list_is_empty(todo_list))
{
list_pop(todo_list,&next);
p= (struct grad12_struct*)olist_entry(edges,next);

if(p->flags&4)continue;
p->flags|= 4;

i= p->flags&3;
if(i==3)id= p->links[0];
else id= p->links[i];

e= id2edge(olist_get_key(edges,next));

while(id>=0)
{
q= (struct grad12_struct*)olist_entry(edges,id);

ep= id2edge(olist_get_key(edges,id));

if((q->flags&4)==0)
{
list_push(todo_list,&id);
}
i= edge_in_triangle(e,r12(ep));
if(i<0)abort_message("error");
id= q->links[i];
}
}
}


/*:66*/


p= (struct grad12_struct*)olist_entry(edges,todo);
todo= p->links[1];
}

list_abandon(&todo_list);

for(j= 0;j<olist_count(edges);j++)
{
p= (struct grad12_struct*)olist_entry(edges,j);
if(p->flags&4)continue;
e= id2edge(olist_get_key(edges,j));
make_deadend(e);
make_deadend(r12(e));
}
}

/*:65*/



if(options&4)
/*67:*/


{
list*todo_list;
long id= -1,*idp;
struct grad12_struct*q;
int error_check;

list_initialize(&todo_list,sizeof(long));
list_push(todo_list,&id);

while(!list_is_empty(todo_list))
{
list_read_init(todo_list);
error_check= 1;
while((idp= (long*)list_read(todo_list))!=NULL)
{
if(*idp<0)
{
list_read_delete(todo_list);
kk= -1;
f= sigma;
}
else
{
p= (struct grad12_struct*)olist_entry(edges,*idp);
if(p->flags&24)continue;

/*68:*/


{
edge*e;
int i;

ep= id2edge(olist_get_key(edges,*idp));
kk= p->flags&3;
if(kk==3)id= p->links[0];
else id= p->links[kk];

while(id>=0)
{
q= (struct grad12_struct*)olist_entry(edges,id);
if((q->flags&8)==0)break;
e= id2edge(olist_get_key(edges,id));
i= edge_in_triangle(ep,r12(e));
id= q->links[i];
}

if(id<0)
/*71:*/


{
edge*e;
int i;
long ct= 0;
long orient;
triangle*t;

p->flags&= 0xffffffef;
p->flags|= 8;
list_read_delete(todo_list);

kk= p->flags&3;
if(kk==3)id= p->links[0];
else id= p->links[kk];

while(id>=0)
{
q= (struct grad12_struct*)olist_entry(edges,id);
e= id2edge(olist_get_key(edges,id));
t= r12(e);
i= edge_in_triangle(ep,t);
id= q->links[i];
ct+= edge_orient(t,i)*q->count;
}
if(id==-2)
{
ct+= edge_orient(sigma,edge_in_triangle(ep,sigma));
}
if(is_critical(ep))p->count= ct;
else p->count= -ct*edge_orient(r12(ep),p->flags&3);
}

/*:71*/


else p->flags|= 16;

if(id>=0||is_critical(ep))f= NULL;
else f= r12(ep);
}

/*:68*/


}

if(f!=NULL)
{
error_check= 0;
/*69:*/


{
int k;
edge*e;

for(k= 0;k<3;k++)
{
if(k==kk)continue;
e= get_edge(f,k);
/*70:*/


{
int flag;
long qid;

if(is_critical(e)||(is_paired_up(e)&&(!is_deadend(e)||(options&3)==0)))
{
qid= olist_find_add(edges,edge_id(e),NULL,&flag);
if(!flag)abort_message("count error 2");
q= (struct grad12_struct*)olist_entry(edges,qid);
if((q->flags&32)==0)
{
list_read_insert(todo_list,&qid);
q->flags|= 32;
}
else q->flags&= 0xffffffef;
}
}

/*:70*/


}
}

/*:69*/


}
}
if(error_check&&!list_is_empty(todo_list))abort_message("count error");
}

list_abandon(&todo_list);
}

/*:67*/



return crit;
}

/*:62*//*72:*/


long
find_all_backward_grad12_paths(edge*tau,olist*triangles,list*crits,int options)
{
int kk= -1;
triangle*f;
edge*ep;
long todo= -1;
long this;
struct grad12_struct*p;
long start= -1;
list*todo_list;

olist_clear(triangles);
if(crits!=NULL)list_clear(crits);
if((options&3)==2)list_initialize(&todo_list,sizeof(long));

do
{
this= todo;
if(this>=0)
{
p= (struct grad12_struct*)olist_entry(triangles,this);
f= id2triangle(olist_get_key(triangles,this));
ep= r21(f);
kk= p->flags&3;
todo= p->links[kk];
p->links[kk]= -1;
}
else ep= tau;
/*73:*/


{
tetrahedron*t;
int link_count= 0;

f= ep->links[2];
t= coface(f,0);
do
{
if(is_in_K(f)&&(is_critical(f)||
((!is_xdeadend(f)||(options&3)==0)&&
is_paired_down(f)&&ep!=r21(f))))
{
link_count++;
/*74:*/


{
struct grad12_struct r,*q;
long id;
int flag;
int k;

id= olist_find_add(triangles,triangle_id(f),&r,&flag);
q= (struct grad12_struct*)olist_entry(triangles,id);
k= edge_in_triangle(ep,f);
if(this>=0)
p= (struct grad12_struct*)olist_entry(triangles,this);
if(flag)
{
if(q->links[k]>=0)abort_message("backgrad error 1");
if(this<0)abort_message("backgrad error 2");
q->links[k]= p->links[kk];
p->links[kk]= id;
}
else
{
q->links[0]= q->links[1]= q->links[2]= -2;
if(this>=0)
{
q->links[k]= p->links[kk];
p->links[kk]= id;
}
else
{
q->links[k]= start;
start= id;
}
if(is_critical(f))
{
q->flags= 3;
if(crits!=NULL)list_push(crits,&id);
if((options&3)==2)list_push(todo_list,&id);
}
else
{
q->flags= edge_in_triangle(r21(f),f);
q->links[q->flags]= todo;
todo= id;
}
}
}

/*:74*/


}
t= other_coface(f,t);
f= other_face(ep,f,t);
}while(f!=ep->links[2]);
if(link_count==0&&(options&3)&&this>=0)
{
make_xdeadend(ep);
make_xdeadend(r12(ep));
}
}

/*:73*/


}while(todo>=0);
if((options&3)==2)
/*75:*/


{
long id;
int k;
edge*e;
struct grad12_struct r,*q;
triangle*f;
int flag;

while(!list_is_empty(todo_list))
{
list_pop(todo_list,&id);
p= (struct grad12_struct*)olist_entry(triangles,id);
if(p->flags&4)continue;
p->flags|= 4;
kk= p->flags&3;
f= id2triangle(olist_get_key(triangles,id));
for(k= 0;k<3;k++)
{
if(k==kk||p->links[k]<-1)continue;
e= get_edge(f,k);
if(is_paired_up(e)&&!is_deadend(e))
{
if(is_xdeadend(e))abort_message("back grad error 3");
id= olist_find_add(triangles,triangle_id(r12(e)),&r,&flag);
q= (struct grad12_struct*)olist_entry(triangles,id);
if(flag==0)
{
q->flags= edge_in_triangle(e,r12(e));
printf("back grad puzzle\n");
}
if((q->flags&4)==0)list_push(todo_list,&id);
}
}
}
list_abandon(&todo_list);

for(id= 0;id<olist_count(triangles);id++)
{
p= (struct grad12_struct*)olist_entry(triangles,id);
if(p->flags&4)continue;
f= id2triangle(olist_get_key(triangles,id));
make_deadend(f);
make_deadend(r21(f));
}
}

/*:75*/


return start;
}

/*:72*/


/*47:*/


void LocalCancel12(vertex*v,edge*sigma,edge*sigmap,triangle*tau)
{
edge*u;
triangle*w,*wp;

u= sigmap;
w= tau;
unmake_critical(tau);
while(!is_critical(u))
{
wp= r12(u);
pair12(u,w);
w= wp;
u= other_edge(v,u,w);
}
unmake_critical(sigma);
pair12(sigma,w);
}

/*:47*//*48:*/


void Cancel23(triangle*sigma,int n,tetrahedron*tau)
{
triangle*s,*sp;
tetrahedron*t;

s= sigma;
t= coface(s,n);
unmake_critical(sigma);

while(!is_critical(t))
{
sp= r32(t);
pair23(s,t);
s= sp;
t= other_coface(s,t);
}
unmake_critical(tau);
pair23(s,tau);
}

/*:48*//*49:*/




void Cancel01(vertex*v,int n,edge*kappa)
{
vertex*u;
edge*e,*ep;

u= get_vertex(kappa,n);

e= kappa;
unmake_critical(kappa);

while(!is_critical(u))
{
ep= r01(u);
pair01(u,e);
e= ep;
u= other_vertex_in_edge(u,e);
}
unmake_critical(v);
pair01(v,e);
}


/*:49*//*50:*/


void Cancel12(edge*sigma,triangle*kappa,list*grad_path,list*changed)
{
triangle*t,*nextt;
edge*e;

unmake_critical(kappa);
unmake_critical(sigma);

t= kappa;

do
{
e= plist_pop(grad_path);
if(e==sigma)break;
nextt= r12(e);
pair12(e,t);
plist_push(changed,t);
t= nextt;
}while(1);
pair12(e,t);
plist_push(changed,t);
}




/*:50*/


/*21:*/


void
LocalCancel(vertex*v,list*lcrit2)
{
triangle*t,*s;
tetrahedron*te[2];
edge*e[2],*ep[2];
int j,i,in;


while(!list_is_empty(lcrit2))
{
t= plist_pop(lcrit2);
/*22:*/



for(i= 0,j= 0;j<3;j++)
{
if(vertex_in_edge(v,get_edge(t,j))>=0)
{
e[i++]= get_edge(t,j);
}
}

/*:22*/

;
/*23:*/



for(i= 0;i<2;i++)
{
ep[i]= e[i];
while(!is_critical(ep[i]))
{
if(is_paired_down(ep[i]))
{
ep[i]= NULL;
break;
}
s= r12(ep[i]);
ep[i]= other_edge(v,ep[i],s);
}
}

/*:23*/

;

if(ep[0]!=ep[1])
{
/*25:*/


if(ep[0]==NULL)LocalCancel12(v,ep[1],e[1],t);
else if(ep[1]==NULL)LocalCancel12(v,ep[0],e[0],t);
else if(min_value(ep[0],1)> min_value(ep[1],1))
LocalCancel12(v,ep[0],e[0],t);
else LocalCancel12(v,ep[1],e[1],t);

/*:25*/

;
}
else
{

for(i= 0;i<2;i++)
{
/*24:*/


te[i]= coface(t,i);
s= t;
while(is_in_K(te[i]))
{
in= is_in_lower_Star(te[i],v);
if(!in)break;
if(is_critical(te[i]))break;
s= r32(te[i]);
te[i]= other_coface(s,te[i]);
}
if(!is_in_K(te[i])||!in)
te[i]= NULL;

/*:24*/

;
}
if(te[0]!=te[1])
{
/*26:*/


if(te[0]==NULL)Cancel23(t,1,te[1]);
else if(te[1]==NULL)Cancel23(t,0,te[0]);
else if(min_value(te[1],3)> min_value(te[0],3))
Cancel23(t,0,te[0]);
else Cancel23(t,1,te[1]);







/*:26*/

;
}
}
}
}


/*:21*//*27:*/



void ExtractCancel1(int p)
{
olist*c0;
olist*c1;
struct ccrit1
{
edge*s;
long c[2];
}s;
int i;
vertex*v[2];
int thisk;

olist_initialize(&c0,sizeof(long));
olist_initialize(&c1,sizeof(struct ccrit1));

/*28:*/


{
edge*e;

list_read_init(crit[1]);
while((e= (edge*)plist_read(crit[1]))!=NULL)
{
if(!is_critical(e))list_read_delete(crit[1]);
else
{

v[0]= FindGrad01(get_vertex(e,0),value(e)-p);
v[1]= FindGrad01(get_vertex(e,1),value(e)-p);
if(v[0]!=v[1])
{
/*29:*/


{
long link;

link= -1;
s.s= e;
s.c[0]= (v[0]==NULL)?-1:olist_find_add(c0,(long)v[0],&link,NULL);
s.c[1]= (v[1]==NULL)?-1:olist_find_add(c0,(long)v[1],&link,NULL);
i= /*34:*/


(v[0]==NULL)?1:
(v[1]==NULL)?0:
(value(v[0])> value(v[1]))?0:
(value(v[0])<value(v[1]))?1:
is_critical(v[0])?1:0




/*:34*/

;
olist_add(c1,value(e)-value(v[i]),&s);
}

/*:29*/

;
}
}
}
}

/*:28*/

;

while(!olist_is_empty(c1))
{
thisk= olist_min(c1,&s);
/*30:*/


v[0]= (s.c[0]<0)?NULL:(vertex*)olist_get_key(c0,s.c[0]);
v[1]= (s.c[1]<0)?NULL:(vertex*)olist_get_key(c0,s.c[1]);


/*:30*/


i= /*34:*/


(v[0]==NULL)?1:
(v[1]==NULL)?0:
(value(v[0])> value(v[1]))?0:
(value(v[0])<value(v[1]))?1:
is_critical(v[0])?1:0




/*:34*/

;
if(is_critical(v[i])&&thisk==value(s.s)-value(v[i]))/*31:*/



{
Cancel01(v[i],i,s.s);
*((long*)olist_entry(c0,s.c[i]))= s.c[1-i];

}

/*:31*/


else/*32:*/


{
int m;

for(i= 0;i<2;i++)
{
if(s.c[i]<0)continue;
/*33:*/


{
long a;

for(a= s.c[i];a>=0;)
{
v[i]= (vertex*)olist_get_key(c0,a);
if(is_critical(v[i]))break;
a= *((long*)olist_entry(c0,a));
}
*((long*)olist_entry(c0,s.c[i]))= a;
if(value(v[i])<=value(s.s)-p)
{
v[i]= NULL;
s.c[i]= -1;
}
else s.c[i]= a;
}

/*:33*/


}


if(s.c[1]!=s.c[0])
{
i= /*34:*/


(v[0]==NULL)?1:
(v[1]==NULL)?0:
(value(v[0])> value(v[1]))?0:
(value(v[0])<value(v[1]))?1:
is_critical(v[0])?1:0




/*:34*/

;
m= value(s.s)-value(v[i]);

if(m<p)olist_add(c1,m,&s);
}
}

/*:32*/


}
olist_abandon(&c0);
olist_abandon(&c1);
}

/*:27*//*35:*/



void ExtractCancel3(int p)
{
olist*c3;
olist*c2;
struct ccrit2
{
triangle*s;
long c[2];
}r;
int i;
tetrahedron*t[2];
int thisk;

olist_initialize(&c2,sizeof(struct ccrit2));
olist_initialize(&c3,sizeof(long));

/*36:*/


{
triangle*s;

list_read_init(crit[2]);
while((s= (triangle*)plist_read(crit[2]))!=NULL)
{
if(!is_critical(s))list_read_delete(crit[2]);
else
{

t[0]= FindGrad23(coface(s,0),value(s)+p);
t[1]= FindGrad23(coface(s,1),value(s)+p);
if(t[0]!=t[1])
/*37:*/


{
long link;

link= -1;
r.s= s;
r.c[0]= (t[0]==NULL)?-1:olist_find_add(c3,(long)t[0],&link,NULL);
r.c[1]= (t[1]==NULL)?-1:olist_find_add(c3,(long)t[1],&link,NULL);

i= /*42:*/


(t[0]==NULL)?1:
(t[1]==NULL)?0:
(value(t[0])> value(t[1]))?1:
(value(t[0])<value(t[1]))?0:
is_critical(t[0])?1:0





/*:42*/

;
olist_add(c2,value(t[i])-value(s),&r);
}

/*:37*/


}
}
}

/*:36*/

;
while(!olist_is_empty(c2))
{
thisk= olist_min(c2,&r);
/*38:*/


t[0]= (r.c[0]<0)?NULL:(tetrahedron*)olist_get_key(c3,r.c[0]);
t[1]= (r.c[1]<0)?NULL:(tetrahedron*)olist_get_key(c3,r.c[1]);

/*:38*/


i= /*42:*/


(t[0]==NULL)?1:
(t[1]==NULL)?0:
(value(t[0])> value(t[1]))?1:
(value(t[0])<value(t[1]))?0:
is_critical(t[0])?1:0





/*:42*/

;
if(is_critical(t[i])&&thisk==value(t[i])-value(r.s))/*39:*/


{
Cancel23(r.s,i,t[i]);
*((long*)olist_entry(c3,r.c[i]))= r.c[1-i];
}

/*:39*/


else/*40:*/


{
int m;

for(i= 0;i<2;i++)
{
if(r.c[i]<0)continue;
/*41:*/


{
long a;

for(a= r.c[i];a>=0;)
{
t[i]= (tetrahedron*)olist_get_key(c3,a);
if(is_critical(t[i]))break;
a= *((long*)olist_entry(c3,a));
}
*((long*)olist_entry(c3,r.c[i]))= a;
if(value(t[i])>=value(r.s)+p)
{
t[i]= NULL;
r.c[i]= -1;
}
else r.c[i]= a;
}

/*:41*/


}


if(r.c[1]!=r.c[0])
{
i= /*42:*/


(t[0]==NULL)?1:
(t[1]==NULL)?0:
(value(t[0])> value(t[1]))?1:
(value(t[0])<value(t[1]))?0:
is_critical(t[0])?1:0





/*:42*/

;

m= value(t[i])-value(r.s);

if(m<p)olist_add(c2,m,&r);
}
}

/*:40*/


}
olist_abandon(&c2);
olist_abandon(&c3);
}

/*:35*//*43:*/


void ExtractCancel2(int p)
{
list*changed;
list*grad_path;
olist*goodpairs;
int lastp;

list_initialize(&changed,sizeof(triangle*));
list_initialize(&grad_path,sizeof(edge*));
olist_initialize(&goodpairs,sizeof(triangle*));
lastp= 0;
do
{
/*44:*/


{
triangle*t;
edge*e;

list_read_init(crit[2]);
while((t= (triangle*)plist_read(crit[2]))!=NULL)
{
if(!is_critical(t))list_read_delete(crit[2]);
else
{
e= FindGradPaths12(t,p,NULL,0);
if(e!=NULL)
olist_add(goodpairs,value(t)-value(e),&t);
}
}
}




/*:44*/


if(olist_is_empty(goodpairs))break;
/*45:*/


{
triangle*t;
edge*e;
int bp,thisp;

while(!olist_is_empty(goodpairs))
{
bp= olist_min(goodpairs,&t);
e= FindGradPaths12(t,p,grad_path,0);
if(e==NULL)continue;
thisp= value(t)-value(e);
if(thisp<=bp)
{
list_clear(changed);
Cancel12(e,t,grad_path,changed);
/*46:*/


{
int i,j,k;
tetrahedron*te;

k= 0;
while(!list_is_empty(changed)&&k<10000)
{
t= plist_pop(changed);
if(t==NULL)continue;
for(i= 0;i<2;i++)
{
te= coface(t,i);
if(is_in_K(te)&&!is_critical(te))
{
j= splitrejoin(te);
if(j>=0)
{
plist_push(changed,unsplitrejoin(te,j));
k++;
}
}
}
}
if(!list_is_empty(changed))printf("May be infinite loop in split-rejoin");
}




/*:46*/

;
#ifdef verbose
if(thisp<lastp)printf("\n persistence out of order: %d > %d\n",lastp,thisp);
#endif
lastp= thisp;
}
else olist_add(goodpairs,thisp,&t);
}
}






/*:45*/


}while(1);
list_abandon(&changed);
list_abandon(&grad_path);
olist_abandon(&goodpairs);
}

/*:43*/


/*3:*/


void Extract(int p)
{
long i;
long t1,t2,t3,t4,t5,t6;

t1= clock();
for(i= 0;i<4;i++)list_initialize(crit+i,sizeof(void*));

/*4:*/


{
edge*beste;
vertex*v;
int is_lower_Star_empty;
long vid;
list*lcrit2;
list*edges_todo;

list_initialize(&lcrit2,sizeof(triangle*));
list_initialize(&edges_todo,sizeof(edge*));

for(vid= 0;vid<num_simp[0];vid++)
{
v= id2vertex(vid);
if(!is_in_K(v))continue;
beste= NULL;

/*5:*/


{
edge*e;
vertex*w;
int val,bestval,flag;

val= value(v);

list_clear(edges_todo);

plist_push(edges_todo,v->links[0]);

is_lower_Star_empty= 1;
while(!list_is_empty(edges_todo))
{
/*6:*/


e= plist_pop(edges_todo);
if(edge_marked_done(v,e))continue;
mark_edge_done(v,e);


/*:6*/


/*7:*/


flag= 0;
if(is_in_K(e))
{
w= smallest_vertex_in_edge(e);

if(w!=v)
{
set_value(e,val);
is_lower_Star_empty= 0;
flag= 1;
}
}


/*:7*/



/*8:*/


{
triangle*s,*bests,*sp;
tetrahedron*t;
edge*ep;
int state= 0;
int bestsval;
int first_in_lower_Star;
vertex*min_vertex;

#ifdef slower_and_finer 
triangle*bests_in_run,*first_bests_in_run;
tetrahedron*bestt_in_run,*first_bestt_in_run;
int bestsval_in_run,first_bestsval;
#else
vertex*first_min_vertex;
#endif 

s= sp= e->links[2];
t= coface(s,0);

do
{
ep= other_edge(v,e,s);
if(!edge_marked_done(v,ep))plist_push(edges_todo,ep);
#ifdef slower_and_finer 
if(flag)/*14:*/


{
int in_lower_Star_of_e;

/*10:*/


in_lower_Star_of_e= 0;
if(is_in_K(s))
{
vertex*vlist[3];

get_triangle_vertices(s,vlist);
min_vertex= vlist[smallest_vertex(vlist,3)];
in_lower_Star_of_e= (vertex_in_edge(min_vertex,e)<0);
if(in_lower_Star_of_e)set_value(s,val);
}

/*:10*/

;

switch(state)
{
case 0:
if(in_lower_Star_of_e)
{
first_in_lower_Star= (s==sp&&is_in_K(t));

if(first_in_lower_Star)state= 1;
else state= 3;

bests= NULL;
bests_in_run= s;
bestt_in_run= t;
bestsval_in_run= value(min_vertex);
}
break;
case 1:
case 3:
if(!in_lower_Star_of_e||!is_in_K(t))
{
/*15:*/


if(bests==NULL)
{
if(state==1)
{
first_bests_in_run= bests_in_run;
first_bestt_in_run= bestt_in_run;
first_bestsval= bestsval_in_run;
}
else
{
bests= bests_in_run;
bestsval= bestsval_in_run;
}
}
else if(bestsval_in_run<bestval)
{
make_critical(bests);
plist_push(lcrit2,bests);
bests= bests_in_run;
bestsval= bestsval_in_run;
}
else
{
make_critical(bests_in_run);
plist_push(lcrit2,bests_in_run);
}

/*:15*/


}
if(!in_lower_Star_of_e)
{
state= 2;
break;
}
else if(is_in_K(t))
{
set_value(t,val);

if(value(min_vertex)<bestsval_in_run)
/*16:*/


{
triangle*ss;
tetrahedron*tt;

do
{
tt= other_coface(bests_in_run,bestt_in_run);
ss= r32(tt);
pair23(bests_in_run,tt);
bests_in_run= ss;
bestt_in_run= tt;
}while(bestt_in_run!=t);
bests_in_run= s;
bestsval_in_run= value(min_vertex);
}

/*:16*/


else pair23(s,t);
break;
}

case 2:
if(in_lower_Star_of_e)
{
bests_in_run= s;
bestt_in_run= t;
bestsval_in_run= value(min_vertex);
state= 3;
}
break;
}
}

/*:14*/


#else
if(flag)/*9:*/


{
int in_lower_Star_of_e;

/*10:*/


in_lower_Star_of_e= 0;
if(is_in_K(s))
{
vertex*vlist[3];

get_triangle_vertices(s,vlist);
min_vertex= vlist[smallest_vertex(vlist,3)];
in_lower_Star_of_e= (vertex_in_edge(min_vertex,e)<0);
if(in_lower_Star_of_e)set_value(s,val);
}

/*:10*/

;

switch(state)
{
case 0:
if(in_lower_Star_of_e)
{
first_in_lower_Star= (s==sp&&is_in_K(t));

if(first_in_lower_Star)
{
state= 1;
bests= NULL;
first_min_vertex= min_vertex;
}
else
{
state= 3;
bests= s;
bestsval= value(min_vertex);
}
}
break;
case 1:
case 3:
if(!in_lower_Star_of_e)
{
state= 2;
break;
}
else if(is_in_K(t))
{
set_value(t,val);
pair23(s,t);
break;
}

case 2:
if(in_lower_Star_of_e)
{
/*11:*/


if(bests==NULL)
{
bests= s;
bestsval= value(min_vertex);
}
else if(value(min_vertex)<bestval)
{
make_critical(bests);
plist_push(lcrit2,bests);
bests= s;
bestsval= value(min_vertex);
}
else
{
make_critical(s);
plist_push(lcrit2,s);
}


/*:11*/


state= 3;
}
break;
}
}

/*:9*/


#endif 
t= other_coface(s,t);
s= other_face(e,s,t);
}while(s!=sp);
#ifdef slower_and_finer 
if(flag)/*17:*/


{
triangle*sss;

switch(state)
{
case 0:
/*13:*/


if(beste==NULL)
{
beste= e;
bestval= value(w);
}
else if(value(w)<bestval)
{
make_critical(beste);
beste= e;
bestval= value(w);
}
else make_critical(e);


/*:13*/


break;
case 1:
make_critical(t);
set_value(t,val);
pair12(e,bests_in_run);
break;
case 3:
if(first_in_lower_Star)
{
set_value(t,val);
if(first_bestsval> bestsval_in_run)
{
sss= bests_in_run;
bestt_in_run= other_coface(first_bests_in_run,first_bestt_in_run);
bests_in_run= first_bests_in_run;
}
else sss= first_bests_in_run;
/*16:*/


{
triangle*ss;
tetrahedron*tt;

do
{
tt= other_coface(bests_in_run,bestt_in_run);
ss= r32(tt);
pair23(bests_in_run,tt);
bests_in_run= ss;
bestt_in_run= tt;
}while(bestt_in_run!=t);
bests_in_run= s;
bestsval_in_run= value(min_vertex);
}

/*:16*/


bests_in_run= sss;
}
/*15:*/


if(bests==NULL)
{
if(state==1)
{
first_bests_in_run= bests_in_run;
first_bestt_in_run= bestt_in_run;
first_bestsval= bestsval_in_run;
}
else
{
bests= bests_in_run;
bestsval= bestsval_in_run;
}
}
else if(bestsval_in_run<bestval)
{
make_critical(bests);
plist_push(lcrit2,bests);
bests= bests_in_run;
bestsval= bestsval_in_run;
}
else
{
make_critical(bests_in_run);
plist_push(lcrit2,bests_in_run);
}

/*:15*/


pair12(e,bests);
break;
case 2:
if(first_in_lower_Star)
{
bests_in_run= first_bests_in_run;
bestsval_in_run= first_bestsval;
/*15:*/


if(bests==NULL)
{
if(state==1)
{
first_bests_in_run= bests_in_run;
first_bestt_in_run= bestt_in_run;
first_bestsval= bestsval_in_run;
}
else
{
bests= bests_in_run;
bestsval= bestsval_in_run;
}
}
else if(bestsval_in_run<bestval)
{
make_critical(bests);
plist_push(lcrit2,bests);
bests= bests_in_run;
bestsval= bestsval_in_run;
}
else
{
make_critical(bests_in_run);
plist_push(lcrit2,bests_in_run);
}

/*:15*/


}
pair12(e,bests);
break;
}
}





/*:17*/


#else
if(flag)/*12:*/


switch(state)
{
case 0:
/*13:*/


if(beste==NULL)
{
beste= e;
bestval= value(w);
}
else if(value(w)<bestval)
{
make_critical(beste);
beste= e;
bestval= value(w);
}
else make_critical(e);


/*:13*/


break;
case 1:
make_critical(t);
set_value(t,val);
pair12(e,s);
break;
case 3:
if(first_in_lower_Star)
{
pair23(s,t);
set_value(t,val);
}
pair12(e,bests);
break;
case 2:
if(first_in_lower_Star)
{
min_vertex= first_min_vertex;
/*11:*/


if(bests==NULL)
{
bests= s;
bestsval= value(min_vertex);
}
else if(value(min_vertex)<bestval)
{
make_critical(bests);
plist_push(lcrit2,bests);
bests= s;
bestsval= value(min_vertex);
}
else
{
make_critical(s);
plist_push(lcrit2,s);
}


/*:11*/


}
pair12(e,bests);
break;
}





/*:12*/


#endif 
}

/*:8*/


}
}


/*:5*/

;
if(is_lower_Star_empty)make_critical(v);
else
{
pair01(v,beste);
LocalCancel(v,lcrit2);
}
}
list_abandon(&lcrit2);
list_abandon(&edges_todo);
}


/*:4*/

;
t2= clock();


ExtractCancel1(p);
t3= clock();


ExtractCancel3(p);
t4= clock();


ExtractCancel2(p);
t5= clock();

#ifdef verbose

printf(
"\nExtract Raw = %d, Cancel1 = %d, Cancel2 = %d, Cancel3 = %d\n",
t2-t1,t3-t2,t5-t4,t4-t3);
#endif
}



/*:3*/



/*:1*/