Назад к лучшим решениям Status: AC, problem ZELC, contest ZEL07. By werewolf (Werewolf), 2007-03-10 03:57:25.
  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <string.h>
  4. #include <math.h>
  5.  
  6. #define PREC 0.001
  7. #define MAX_NODES 3004
  8. #define MAX_TRIANGLES (2*MAX_NODES+1)
  9. #define MAX_EDGES (2*MAX_TRIANGLES +1)
  10.  
  11. #define CACHE_INCREASE_COEF 12
  12. #define MAX_CACHE_SIZE 64
  13.  
  14. typedef struct _COORD
  15. {
  16. 	float x;
  17. 	float y;
  18. } COORD, *PCOORD;
  19.  
  20. typedef struct _NODE
  21. {
  22. 	COORD crd;
  23. 	int id;
  24. } NODE, *PNODE;
  25.  
  26.  
  27. typedef struct _TRIANGLE
  28. {
  29. 	NODE * n1, * n2, * n3;
  30. 	struct _TRIANGLE * t1, *t2, *t3;
  31. 	float d1, d2, d3;
  32. 	char usedEdges;
  33. } TRIANGLE, *PTRIANGLE;
  34.  
  35. typedef struct _REBUILD_ITEM
  36. {
  37. 	PTRIANGLE toProcess;
  38. 	PTRIANGLE toExclude;
  39. }REBUILD_ITEM, *PREBUILD_ITEM;
  40.  
  41. typedef struct _EDGE
  42. {
  43. 	float d;
  44. 	int n1;
  45. 	int n2;	
  46. 	PTRIANGLE trg1, trg2;
  47. }EDGE, *PEDGE;
  48.  
  49. typedef struct _TRG_SQUARE
  50. {
  51. 	float sqr;
  52. 	PTRIANGLE trg;	
  53. }TRG_SQUARE, *PTRG_SQUARE;
  54.  
  55. struct NORMALIZE
  56. {
  57. 	float minx, maxx, miny, maxy;
  58. 	float divx;
  59. 	float divy;
  60. } norm;
  61.  
  62. PTRIANGLE cache[MAX_CACHE_SIZE][MAX_CACHE_SIZE];
  63. int cache_size = 2;
  64.  
  65. REBUILD_ITEM rebuild_queue[MAX_TRIANGLES];
  66. int start_queue = 0, end_queue = 0;
  67.  
  68. #define vect_abs_2(a) ((a).x*(a).x + (a).y*(a).y)
  69. #define segm_abs_2(a,b) (((a).x-(b).x)*((a).x-(b).x) + ((a).y-(b).y)*((a).y-(b).y))
  70. #define dots_equal(a, b) ((a).x==(b).x && (a).y==(b).y)
  71.  
  72. TRIANGLE triangles[MAX_TRIANGLES];
  73. PTRIANGLE pnexttrg = triangles;
  74. int nTriangles = 0;
  75.  
  76. NODE nodes[MAX_NODES+MAX_NODES-2];
  77. int nNodes = 0;
  78. int nSteinerNodes = 0;
  79.  
  80. EDGE edges[MAX_EDGES];
  81. int nEdges = 0;
  82.  
  83.  
  84. typedef void (*MSTCALLBACK)(PEDGE pArc);
  85.  
  86. __inline void addEdge(int node1, int node2, float d, PTRIANGLE trg1, PTRIANGLE trg2)
  87. {
  88. 	edges[nEdges].trg1 = trg1;
  89. 	edges[nEdges].trg2 = trg2;
  90. 	edges[nEdges].d = d;
  91. 	edges[nEdges].n1 = node1;
  92. 	edges[nEdges++].n2 = node2;
  93. }
  94. __inline PTRIANGLE addTriangle(PNODE a, PNODE b, PNODE c)
  95. {
  96. 	pnexttrg->n1 = a;
  97. 	pnexttrg->n2 = b;
  98. 	pnexttrg->n3 = c;
  99. 	pnexttrg->d1 = pnexttrg->d2 = pnexttrg->d3 = 0;
  100. 	pnexttrg->usedEdges = 0;
  101. 	nTriangles++;
  102. 	return pnexttrg++;
  103. }
  104. __inline void pushToRebuild(PTRIANGLE trg, PTRIANGLE excl)
  105. {
  106. 	int i;
  107. 	for(i = start_queue; i < end_queue; i++)
  108. 		if(trg == rebuild_queue[i].toProcess) return;
  109.  
  110. 	rebuild_queue[end_queue].toProcess = trg;
  111. 	rebuild_queue[end_queue++].toExclude = excl;
  112. }
  113. void increaseCache()
  114. {
  115. 	int i, j;
  116. 	for(i = cache_size-1; i >= 0; i--)
  117. 	{
  118. 		for(j = cache_size-1; j >= 0; j--)
  119. 		{
  120. 			cache[2*i][2*j] = cache[2*i][2*j+1] =
  121. 			cache[2*i+1][2*j] = cache[2*i+1][2*j+1] = cache[i][j];
  122. 		}
  123. 	}
  124. 	cache_size*=2;
  125. }
  126. __inline PTRIANGLE popToRebuild(PTRIANGLE * pexcl)
  127. {
  128. 	PTRIANGLE ret;
  129. 	if(start_queue == end_queue) return 0;
  130. 	ret = rebuild_queue[start_queue].toProcess;
  131. 	*pexcl = rebuild_queue[start_queue++].toExclude;
  132.  
  133. 	if(start_queue == end_queue) start_queue = end_queue = 0;
  134. 	return ret;
  135. }
  136. enum {LEFT,  RIGHT,  BEYOND,  BEHIND, BETWEEN, START, END};
  137. int dotRelation(PCOORD dot, PCOORD start, PCOORD end)
  138. {
  139. 	COORD a = {end->x - start->x, end->y - start->y};
  140. 	COORD b = {dot->x - start->x, dot->y - start->y};
  141. 	
  142. 	float sa = a.x * b.y - b.x * a.y;
  143. 	
  144. 	if(sa > PREC)
  145. 	  return LEFT;
  146. 	if(sa < -PREC)
  147. 	  return RIGHT;
  148. 	if((a.x * b.x < 0.0) || (a.y * b.y < 0.0))
  149. 	  return BEHIND;
  150. 	if(vect_abs_2(a) < vect_abs_2(b))
  151. 	  return BEYOND;
  152. 	if(dots_equal(*start, *dot))
  153. 	  return START;
  154. 	if(dots_equal(*end, *dot))
  155. 	  return END;
  156. 	return BETWEEN;
  157. }
  158.  
  159. __inline int whereIsDot(PCOORD dot, PTRIANGLE trg)
  160. {
  161. 	if(dotRelation (dot, &trg->n1->crd, &trg->n2->crd) == LEFT)
  162. 		return 1;
  163. 	if(dotRelation (dot, &trg->n2->crd, &trg->n3->crd) == LEFT)
  164. 		return 2;
  165. 	if(dotRelation (dot, &trg->n3->crd, &trg->n1->crd) == LEFT)
  166. 		return 3;
  167. 	return 0;
  168. }
  169. PTRIANGLE findTriangle(PCOORD dot)
  170. {
  171. 	int a, b;
  172. 	PTRIANGLE tmp;
  173. 	a = (int)floor(((dot->x - norm.minx)*cache_size)/norm.divx);
  174. 	b = (int)floor(((dot->y - norm.miny)*cache_size)/norm.divy);
  175. 	tmp = cache[a][b];
  176.  
  177. 	for(;;)
  178. 	{
  179. 		switch(whereIsDot(dot, tmp))
  180. 		{
  181. 		case 0:
  182. 			break;
  183. 		case 1:
  184. 			tmp = tmp->t1;
  185. 			continue;
  186. 		case 2:
  187. 			tmp = tmp->t2;
  188. 			continue;
  189. 		case 3:
  190. 			tmp = tmp->t3;
  191. 			continue;
  192. 		}
  193. 		break;
  194. 	}
  195. 	cache[a][b] = tmp;
  196. 	return tmp;
  197. }
  198. __inline void rotateNodesFor(PTRIANGLE trg)
  199. {
  200. 	PNODE tmpn;
  201. 	PTRIANGLE tmpt;
  202.  
  203. 	tmpn = trg->n3;
  204. 	trg->n3 = trg->n2;
  205. 	trg->n2 = trg->n1;
  206. 	trg->n1 = tmpn;
  207.  
  208. 	tmpt = trg->t3;
  209. 	trg->t3 = trg->t2;
  210. 	trg->t2 = trg->t1;
  211. 	trg->t1 = tmpt;
  212. }
  213. __inline void rotateNodesRev(PTRIANGLE trg)
  214. {
  215. 	PNODE tmpn;
  216. 	PTRIANGLE tmpt;
  217.  
  218. 	tmpn = trg->n1;
  219. 	trg->n1 = trg->n2;
  220. 	trg->n2 = trg->n3;
  221. 	trg->n3 = tmpn;
  222.  
  223. 	tmpt = trg->t1;
  224. 	trg->t1 = trg->t2;
  225. 	trg->t2 = trg->t3;
  226. 	trg->t3 = tmpt;
  227. }
  228. __inline void updateTriangleRef(PTRIANGLE trg, PTRIANGLE ref, PTRIANGLE upd)
  229. {
  230. 	if( trg->t1 == ref )
  231. 		trg->t1 = upd;
  232. 	else if( trg->t2 == ref )
  233. 		trg->t2 = upd;
  234. 	else if( trg->t3 == ref )
  235. 		trg->t3 = upd;
  236. }
  237. int checkDelaunay(PCOORD dot, PTRIANGLE trg)
  238. {
  239. 	// the 1-2 edge must be between the dot and node3
  240. 	float sa, sb;
  241. 	PCOORD c0 = dot, c1 = &trg->n2->crd, c2 = &trg->n3->crd, c3 = &trg->n1->crd;
  242.  
  243. 	sa = (c1->x - c0->x)*(c3->x - c0->x) + (c1->y - c0->y)*(c3->y - c0->y);
  244. 	sb = (c1->x - c2->x)*(c3->x - c2->x) + (c1->y - c2->y)*(c3->y - c2->y);
  245. 	
  246. 	if(sa < 0 && sb < 0) return 0;
  247. 	if(sa >= 0 && sb >= 0) return 1;
  248.  
  249. 	return ((((c3->x - c0->x)*(c1->y - c0->y) - (c1->x - c0->x)*(c3->y - c0->y)) * sb +
  250. 		   ((c1->x - c2->x)*(c3->y - c2->y) - (c3->x - c2->x)*(c1->y - c2->y)) * sa) >= 0);
  251. }
  252. void rebuildTriangles()
  253. {
  254. 	int nCheck;
  255. 	PTRIANGLE trg1, trg2, excl;
  256. 	while(trg1 = popToRebuild(&excl))
  257. 	{
  258. 		nCheck = 3;
  259. 		do 
  260. 		{
  261. 			trg2 = trg1->t1;
  262. 			if(!trg2 || trg2 == excl) continue;
  263. 			while(trg2->t1 != trg1) rotateNodesFor(trg2);
  264. 			if(checkDelaunay(&trg2->n3->crd, trg1)) continue;
  265.  
  266. 			// rebuild triangles
  267. 			trg1->n2 = trg2->n3;
  268. 			trg2->n2 = trg1->n3;
  269.  
  270. 			if(trg1->t2) updateTriangleRef(trg1->t2, trg1, trg2);
  271. 			if(trg2->t2) updateTriangleRef(trg2->t2, trg2, trg1);
  272.  
  273. 			trg1->t1 = trg2->t2;
  274. 			trg2->t1 = trg1->t2;
  275.  
  276. 			trg1->t2 = trg2;
  277. 			trg2->t2 = trg1;
  278.  
  279. 			pushToRebuild(trg2, trg1);
  280. 			nCheck = 3;
  281. 			rotateNodesRev(trg1);
  282. 		} while(rotateNodesRev(trg1), --nCheck);
  283. 	}
  284. }
  285. void splitTriangle(PNODE node, PTRIANGLE trg)
  286. {
  287. 	int r1, r2, r3;
  288. 	PTRIANGLE trg1 = 0, trg2 = 0;
  289. 	r1 = dotRelation (&node->crd, &trg->n1->crd, &trg->n2->crd);
  290. 	r2 = dotRelation (&node->crd, &trg->n2->crd, &trg->n3->crd);
  291. 	r3 = dotRelation (&node->crd, &trg->n3->crd, &trg->n1->crd);
  292. 	if(r1 == RIGHT && r2 == RIGHT && r3 == RIGHT)
  293. 	{
  294. 		// the dot in the triangle
  295. 		trg1 = addTriangle(node, trg->n1, trg->n2);
  296. 		trg2 = addTriangle(node, trg->n3, trg->n1);
  297. 		trg->n1 = node;
  298.  
  299. 		if(trg->t1) updateTriangleRef(trg->t1, trg, trg1);
  300. 		if(trg->t3) updateTriangleRef(trg->t3, trg, trg2);
  301.  
  302. 		trg1->t1 = trg2;
  303. 		trg1->t2 = trg->t1;
  304. 		trg1->t3 = trg;
  305.  
  306. 		trg2->t1 = trg;
  307. 		trg2->t2 = trg->t3;
  308. 		trg2->t3 = trg1;
  309.  
  310. 		trg->t1 = trg1;
  311. 		trg->t3 = trg2;
  312. 		pushToRebuild(trg, trg2);
  313. 		pushToRebuild(trg1, trg);
  314. 		pushToRebuild(trg2, trg1);
  315. 	}else
  316. 	{
  317. 		if(r1 == BETWEEN) {}
  318. 		else if(r2 == BETWEEN)
  319. 			rotateNodesRev(trg);
  320. 		else if(r3 == BETWEEN)
  321. 			rotateNodesFor(trg);
  322. 		else
  323. 			return;
  324.  
  325. 		if(trg->t1)
  326. 			while(trg->t1->t1 != trg) rotateNodesFor(trg->t1);
  327.  
  328. 		// the dot on the edge 1-2 now
  329. 		trg1 = addTriangle(node, trg->n3, trg->n1);
  330. 		trg->n1 = node;
  331. 		if(trg->t1)
  332. 		{
  333. 			trg2 = addTriangle(node, trg->t1->n3, trg->t1->n1);
  334. 			trg->t1->n1 = node;
  335. 		}
  336.  
  337. 		trg1->t1 = trg;
  338. 		trg1->t2 = trg->t3;
  339. 		trg1->t3 = trg->t1;
  340.  
  341. 		if(trg->t1)
  342. 		{
  343. 			trg2->t1 = trg->t1;
  344. 			trg2->t2 = trg->t1->t3;
  345. 			trg2->t3 = trg;
  346. 			
  347. 			if(trg->t1->t3) updateTriangleRef(trg->t1->t3, trg->t1, trg2);
  348. 			trg->t1->t1 = trg1;
  349. 			trg->t1->t3 = trg2;
  350.  
  351. 			pushToRebuild(trg->t1, trg2);
  352. 			pushToRebuild(trg2, trg);
  353. 		}
  354.  
  355. 		pushToRebuild(trg, trg1);
  356. 		pushToRebuild(trg1, trg->t1);
  357.  
  358. 		if(trg->t3) updateTriangleRef(trg->t3, trg, trg1);
  359. 		trg->t1 = trg2;
  360. 		trg->t3 = trg1;
  361.  
  362. 	}
  363. 	rebuildTriangles();
  364. }
  365. __inline int checkEqualNodes(PNODE node, PNODE node1)
  366. {
  367. 	if(dots_equal(node->crd, node1->crd))
  368. 	{
  369. 		addEdge(node->id, node1->id, segm_abs_2(node->crd, node1->crd), 0, 0);
  370. 		return 1;
  371. 	}
  372. 	return 0;	
  373. }
  374. void delaunayTriangulation()
  375. {
  376. 	int i;
  377. 	PTRIANGLE trg1, trg2;
  378.  
  379. 	cache_size = 2;
  380. 	start_queue = 0;
  381. 	end_queue = 0;
  382. 	nEdges = 0;
  383. 	pnexttrg = triangles;
  384. 	nTriangles = 0;
  385. 	
  386. 	nodes[nNodes].crd.x = norm.minx;
  387. 	nodes[nNodes].crd.y = norm.miny;
  388. 	nodes[nNodes].id = nNodes;
  389.  
  390. 	nodes[nNodes+1].crd.x = norm.maxx;
  391. 	nodes[nNodes+1].crd.y = norm.miny;
  392. 	nodes[nNodes+1].id = nNodes+1;
  393.  
  394. 	nodes[nNodes+2].crd.x = norm.maxx;
  395. 	nodes[nNodes+2].crd.y = norm.maxy;
  396. 	nodes[nNodes+2].id = nNodes+2;
  397.  
  398. 	nodes[nNodes+3].crd.x = norm.minx;
  399. 	nodes[nNodes+3].crd.y = norm.maxy;
  400. 	nodes[nNodes+3].id = nNodes+3;
  401. 	trg1 = addTriangle(nodes+nNodes, nodes+nNodes+3, nodes+nNodes+1);
  402. 	trg2 = addTriangle(nodes+nNodes+1, nodes+nNodes+3, nodes+nNodes+2);
  403.  
  404. 	cache[0][0] = cache[0][1]  = trg1;
  405. 	cache[1][0] = cache[1][1]  = trg2;
  406.  
  407. 	trg1->t1 = trg1->t3 = 0;
  408. 	trg1->t2 = trg2;
  409. 	trg2->t2 = trg2->t3 = 0;
  410. 	trg2->t1 = trg1;
  411. 	
  412. 	for(i = 0; i< nNodes; i++)
  413. 	{
  414. 		trg1 = findTriangle(&nodes[i].crd);
  415. 		if(checkEqualNodes(trg1->n1, nodes+i) || 
  416. 			checkEqualNodes(trg1->n2, nodes+i) || 
  417. 			checkEqualNodes(trg1->n3, nodes+i))
  418. 			continue;
  419.  
  420. 		splitTriangle(nodes+i, trg1);
  421.  
  422. 		if(i >= CACHE_INCREASE_COEF*cache_size*cache_size)
  423. 			increaseCache();
  424. 			
  425. 	}
  426. }
  427. __inline int getEdgeMask(PTRIANGLE trg, int startNodeId)
  428. {
  429. 	if(trg->n1->id == startNodeId) return 1;
  430. 	if(trg->n2->id == startNodeId) return 2;
  431. 	if(trg->n3->id == startNodeId) return 4;
  432. 	return 0;
  433. }
  434.  
  435. void extractEdges()
  436. {
  437. 	char edgesmap[MAX_TRIANGLES] = {0};
  438. 	int i = 0;
  439. 	PTRIANGLE trg;
  440.  
  441. 	for(i = 0; i < nTriangles; i++)
  442. 	{
  443. 		trg = triangles + i;
  444. 		
  445. 		trg->d1 = segm_abs_2(trg->n1->crd, trg->n2->crd); 
  446. 		trg->d2 = segm_abs_2(trg->n2->crd, trg->n3->crd); 
  447. 		trg->d3 = segm_abs_2(trg->n3->crd, trg->n1->crd); 
  448.  
  449. 		if(trg->n1->id < nNodes && trg->n2->id < nNodes && trg->n3->id < nNodes)
  450. 		{
  451. 			if(trg->d1 > trg->d2)
  452. 			{
  453. 				if(trg->d1 > trg->d3)
  454. 				{
  455. 					if(trg->t1) edgesmap[trg->t1-triangles] |= getEdgeMask(trg->t1, trg->n2->id);
  456. 					edgesmap[i] |= 1;
  457. 				}
  458. 				else
  459. 				{
  460. 					if(trg->t3) edgesmap[trg->t3-triangles] |= getEdgeMask(trg->t3, trg->n1->id);
  461. 					edgesmap[i] |= 4;
  462. 				}
  463.  
  464.  
  465. 			}else
  466. 			{
  467. 				if(trg->d2 > trg->d3)
  468. 				{
  469. 					if(trg->t2) edgesmap[trg->t2-triangles] |= getEdgeMask(trg->t2, trg->n3->id);
  470. 					edgesmap[i] |= 2;
  471. 				}
  472. 				else
  473. 				{
  474. 					if(trg->t3) edgesmap[trg->t3-triangles] |= getEdgeMask(trg->t3, trg->n1->id);
  475. 					edgesmap[i] |= 4;
  476. 				}
  477. 			}
  478. 		}else
  479. 		{
  480. 			if(trg->n1->id >= nNodes)
  481. 				edgesmap[i] |= 5;
  482. 			if(trg->n2->id >= nNodes)
  483. 				edgesmap[i] |= 3;
  484. 			if(trg->n3->id >= nNodes)
  485. 				edgesmap[i] |= 6;
  486. 		}
  487.  
  488. 		if((edgesmap[i] & 1) == 0)
  489. 		{
  490. 			addEdge(trg->n1->id, trg->n2->id, trg->d1, trg, trg->t1);
  491. 			if(trg->t1) edgesmap[trg->t1-triangles] |= getEdgeMask(trg->t1, trg->n2->id);
  492. 		}
  493. 		
  494. 		if((edgesmap[i] & 2)==0)
  495. 		{
  496. 			addEdge(trg->n2->id, trg->n3->id, trg->d2, trg, trg->t2);
  497. 			if(trg->t2) edgesmap[trg->t2-triangles] |= getEdgeMask(trg->t2, trg->n3->id);
  498.  
  499. 		}
  500. 		
  501. 		if((edgesmap[i] & 4)==0)
  502. 		{
  503. 			addEdge(trg->n3->id, trg->n1->id, trg->d3, trg, trg->t3);
  504. 			if(trg->t3) edgesmap[trg->t3-triangles] |= getEdgeMask(trg->t3, trg->n1->id);
  505. 		}
  506. 		
  507. 	}
  508.  
  509. }
  510. __inline int floatcmp( const void *arg1, const void *arg2 )
  511. {
  512. 	return (*(float*)arg1)<(*(float*)arg2)?-1:(*(float*)arg1)==(*(float*)arg2)?0:1;
  513. }
  514. void findMST(MSTCALLBACK callback)
  515. {
  516. 	int i, j, q = 0;
  517. 	short selected[MAX_NODES];
  518.  
  519. 	for(i = 0; i < nNodes; i++)
  520. 		selected[i] = i;
  521.  
  522. 	qsort(edges, nEdges, sizeof(EDGE), floatcmp);
  523.  
  524. 	for(i = 0; i < nEdges; i++)
  525. 	{
  526. 		if(selected[edges[i].n1] != selected[edges[i].n2])
  527. 		{
  528. 			int s1 = selected[edges[i].n1],
  529. 				s2 = selected[edges[i].n2];
  530.  
  531. 			for(j = 0; j < nNodes; j++)
  532. 				if(selected[j] == s2) selected[j] = s1;
  533.  
  534. 			callback(edges + i);
  535.  
  536. 			q++;
  537. 			if(q == nNodes - 1)
  538. 				break;
  539. 		}
  540. 	}
  541. }
  542. #define two_edges(a) ((a)>4||(a)==3)
  543. #define triangle_square(a,b,c) ((a).x*((b).y-(c).y)+(b).x*((c).y-(a).y)+(c).x*((a).y-(b).y))/2
  544. #define sqrt32 (float)0.86602540378
  545. __inline int floatrcmp( const void *arg1, const void *arg2 )
  546. {
  547. 	return (*(float*)arg1)>(*(float*)arg2)?-1:(*(float*)arg1)==(*(float*)arg2)?0:1;
  548. }
  549. int getSteinerPoint(PCOORD crd1, PCOORD crd2, PCOORD crd3, float d1, float d2, float d3, float sqr, PCOORD spoint)
  550. {
  551. 	COORD P, O, M;
  552. 	double R2, a1, a2, b1, b2, a, b, c, t;
  553. 	double e1, e2, e3;
  554. 	if(d3-d1-d2 > 0.4*sqrt(d1*d2))
  555. 		return 0;
  556.  
  557. 	P.y = crd1->x-crd3->x+crd3->y;
  558. 	P.x = crd3->y-crd1->y+crd3->x;
  559.  
  560. 	if(dotRelation(&P, crd3, crd1) != LEFT)
  561. 	{
  562. 		P.y=-(P.y-crd3->y);
  563. 		P.x=-(P.x-crd3->x);
  564. 		
  565. 	}else
  566. 	{
  567. 		P.y-=crd3->y;
  568. 		P.x-=crd3->x;
  569. 	}
  570. 	P.x*=sqrt32;
  571. 	P.y*=sqrt32;
  572. 	M.x = (crd1->x+crd3->x)/2;
  573. 	M.y = (crd1->y+crd3->y)/2;
  574. 	O.x = P.x/3+M.x;
  575. 	O.y = P.y/3+M.y;
  576. 	P.x += M.x;
  577. 	P.y += M.y;
  578. 	R2 = d3/3;
  579. 	
  580. 	a1 = crd2->x - P.x;
  581. 	a2 = crd2->y - P.y;
  582. 	b1 = P.x - O.x;
  583. 	b2 = P.y - O.y;
  584. 	a = sqrt(a1*a1+a2*a2);
  585. 	a1 /= a;
  586. 	a2 /= a;
  587. 	b = a1*b1+a2*b2;
  588. 	c = b1*b1+b2*b2-R2;
  589. 	
  590. 	if(b>0)
  591. 		t = (-b-sqrt(b*b-c));
  592. 	else
  593. 		t = (-b+sqrt(b*b-c));
  594.  
  595. 	spoint->x = P.x + a1*t;
  596. 	spoint->y = P.y + a2*t;
  597.  
  598. 	return 1;
  599. }
  600. void findSMT()
  601. {
  602. 	TRG_SQUARE squares[MAX_TRIANGLES];
  603. 	PTRIANGLE trg;
  604. 	int i, nSquares = 0, bExists;
  605.  
  606. 	for(i =0; i < nEdges; i++)
  607. 		if(edges[i].d != 0) break;
  608.  
  609. 	nEdges = i;
  610. 	nSteinerNodes = nNodes;
  611.  
  612. 	for(i = 0; i < nTriangles; i++ )
  613. 	{
  614. 		trg = triangles+i;
  615. 		if( two_edges(trg->usedEdges) )
  616. 		{
  617. 			squares[nSquares].trg = trg;
  618. 			squares[nSquares++].sqr = triangle_square(trg->n3->crd, trg->n2->crd, trg->n1->crd);
  619. 		}
  620. 	}
  621. 	qsort(squares, nSquares, sizeof(TRG_SQUARE), floatrcmp);
  622.  
  623. 	for(i = 0; i < nSquares; i++)
  624. 	{
  625. 		trg = squares[i].trg;
  626. 		if( !two_edges(trg->usedEdges) )
  627. 			continue;
  628.  
  629. 		if(trg->usedEdges == 3)
  630. 			bExists = getSteinerPoint(&trg->n1->crd, &trg->n2->crd, &trg->n3->crd, trg->d1, trg->d2, trg->d3, squares[i].sqr, &nodes[nSteinerNodes].crd);
  631. 		else if(trg->usedEdges == 5)
  632. 			bExists = getSteinerPoint(&trg->n3->crd, &trg->n1->crd, &trg->n2->crd, trg->d3, trg->d1, trg->d2, squares[i].sqr, &nodes[nSteinerNodes].crd);
  633. 		else if(trg->usedEdges == 6)
  634. 			bExists = getSteinerPoint(&trg->n2->crd, &trg->n3->crd, &trg->n1->crd, trg->d2, trg->d3, trg->d1, squares[i].sqr, &nodes[nSteinerNodes].crd);
  635.  
  636. 		if(bExists)
  637. 		{
  638. 			if(trg->t1) trg->t1->usedEdges &= ~getEdgeMask(trg->t1, trg->n2->id);
  639. 			if(trg->t2) trg->t2->usedEdges &= ~getEdgeMask(trg->t2, trg->n3->id);
  640. 			if(trg->t3) trg->t3->usedEdges &= ~getEdgeMask(trg->t3, trg->n1->id);
  641. 			nodes[nSteinerNodes].id = nSteinerNodes;
  642.  
  643. 			addEdge(nSteinerNodes, trg->n1->id, 0, 0, 0);
  644. 			addEdge(nSteinerNodes, trg->n2->id, 0, 0, 0);
  645. 			addEdge(nSteinerNodes, trg->n3->id, 0, 0, 0);
  646. 			trg->usedEdges = 0;
  647. 			nSteinerNodes++;
  648. 		}
  649. 	}
  650.  
  651. 	for(i = 0; i < nTriangles; i++)
  652. 	{
  653. 		trg = triangles + i;
  654. 		
  655. 		if(trg->usedEdges & 1)
  656. 		{
  657. 			addEdge(trg->n1->id, trg->n2->id, 0, 0, 0);
  658. 			if(trg->t1) trg->t1->usedEdges &= ~getEdgeMask(trg->t1, trg->n2->id);
  659. 		}
  660. 		
  661. 		if(trg->usedEdges & 2)
  662. 		{
  663. 			addEdge(trg->n2->id, trg->n3->id, 0, 0, 0);
  664. 			if(trg->t2) trg->t2->usedEdges &= ~getEdgeMask(trg->t2, trg->n3->id);
  665. 		}
  666. 		
  667. 		if(trg->usedEdges & 4)
  668. 		{
  669. 			addEdge(trg->n3->id, trg->n1->id, 0, 0, 0);
  670. 			if(trg->t3) trg->t3->usedEdges &= ~getEdgeMask(trg->t3, trg->n1->id);
  671. 		}
  672.  
  673. 	}
  674.  
  675. }
  676.  
  677. void readData()
  678. {
  679. 	int i;
  680. 	scanf("%d", &nNodes);
  681. 	norm.minx = norm.miny = 100000.0;
  682. 	norm.maxx = norm.maxy = 0.0;
  683.  
  684. 	for(i = 0; i < nNodes; i++)
  685. 	{
  686. 		scanf("%f %f", &nodes[i].crd.x, &nodes[i].crd.y);
  687. 		if(nodes[i].crd.x < norm.minx)
  688. 			norm.minx = nodes[i].crd.x;
  689. 		if(nodes[i].crd.x > norm.maxx)
  690. 			norm.maxx = nodes[i].crd.x;
  691. 		if(nodes[i].crd.y < norm.miny)
  692. 			norm.miny = nodes[i].crd.y;
  693. 		if(nodes[i].crd.y > norm.maxy)
  694. 			norm.maxy = nodes[i].crd.y;
  695. 		nodes[i].id = i;
  696. 	}
  697. 	norm.maxx += 1;
  698. 	norm.maxy += 1;
  699. 	norm.minx -= 1;
  700. 	norm.miny -= 1;
  701. 	
  702. 	norm.divx = (norm.maxx - norm.minx);
  703. 	norm.divy = (norm.maxy - norm.miny);
  704. }
  705.  
  706.  
  707. void mstCallback(PEDGE pEdge)
  708. {
  709. 	if(pEdge->trg1) pEdge->trg1->usedEdges |= getEdgeMask(pEdge->trg1, pEdge->n1);
  710. 	if(pEdge->trg2) pEdge->trg2->usedEdges |= getEdgeMask(pEdge->trg2, pEdge->n2);
  711. }
  712. void writeData()
  713. {
  714. 	int i;
  715. 	printf("%d\n", nSteinerNodes - nNodes);
  716. 	for(i = nNodes; i < nSteinerNodes; i++)
  717. 	{
  718. 		printf("%f %f\n", nodes[i].crd.x, nodes[i].crd.y);
  719. 	}
  720. 	printf("%d\n", nEdges);
  721. 	for(i = 0 ; i < nEdges; i++)
  722. 	{
  723. 		printf("%d %d\n", edges[i].n1, edges[i].n2);
  724. 	}
  725. }
  726.  
  727.  
  728. int main()
  729. {
  730. 	int t;
  731.  
  732. 	scanf("%d", &t);
  733. 	while(t--)
  734. 	{
  735. 		readData();
  736. 		delaunayTriangulation();
  737. 		extractEdges();
  738. 		findMST(mstCallback);
  739. 		findSMT();
  740. 		writeData();
  741. 	}
  742. 	return 0;
  743. }