1/* Copyright (c) 2000, 2010 Oracle and/or its affiliates. All rights reserved.
2 Copyright (C) 2011 Monty Program Ab.
3
4 This program is free software; you can redistribute it and/or modify
5 it under the terms of the GNU General Public License as published by
6 the Free Software Foundation; version 2 of the License.
7
8 This program is distributed in the hope that it will be useful,
9 but WITHOUT ANY WARRANTY; without even the implied warranty of
10 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 GNU General Public License for more details.
12
13 You should have received a copy of the GNU General Public License
14 along with this program; if not, write to the Free Software
15 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02111-1301 USA */
16
17
18#include "mariadb.h"
19#include <my_sys.h>
20#include <m_string.h>
21
22#ifdef HAVE_SPATIAL
23
24#include "gcalc_slicescan.h"
25
26
27#define PH_DATA_OFFSET 8
28#define coord_to_float(d) ((double) d)
29#define coord_eq(a, b) (a == b)
30
31typedef int (*sc_compare_func)(const void*, const void*);
32
33#define LS_LIST_ITEM Gcalc_dyn_list::Item
34#define LS_COMPARE_FUNC_DECL sc_compare_func compare,
35#define LS_COMPARE_FUNC_CALL(list_el1, list_el2) (*compare)(list_el1, list_el2)
36#define LS_NEXT(A) (A)->next
37#define LS_SET_NEXT(A,val) (A)->next= val
38#define LS_P_NEXT(A) &(A)->next
39#define LS_NAME sort_list
40#define LS_SCOPE static
41#define LS_STRUCT_NAME sort_list_stack_struct
42#include "plistsort.c"
43
44
45#define GCALC_COORD_MINUS 0x80000000
46#define FIRST_DIGIT(d) ((d) & 0x7FFFFFFF)
47#define GCALC_SIGN(d) ((d) & 0x80000000)
48
49static Gcalc_scan_iterator::point *eq_sp(const Gcalc_heap::Info *pi)
50{
51 GCALC_DBUG_ASSERT(pi->type == Gcalc_heap::nt_eq_node);
52 return (Gcalc_scan_iterator::point *) pi->node.eq.data;
53}
54
55
56static Gcalc_scan_iterator::intersection_info *i_data(const Gcalc_heap::Info *pi)
57{
58 GCALC_DBUG_ASSERT(pi->type == Gcalc_heap::nt_intersection);
59 return (Gcalc_scan_iterator::intersection_info *) pi->node.intersection.data;
60}
61
62
63#ifndef GCALC_DBUG_OFF
64
65int gcalc_step_counter= 0;
66
67void GCALC_DBUG_CHECK_COUNTER()
68{
69 if (++gcalc_step_counter == 0)
70 GCALC_DBUG_PRINT(("step_counter_0"));
71 else
72 GCALC_DBUG_PRINT(("%d step_counter", gcalc_step_counter));
73}
74
75
76const char *gcalc_ev_name(int ev)
77{
78 switch (ev)
79 {
80 case scev_none:
81 return "n";
82 case scev_thread:
83 return "t";
84 case scev_two_threads:
85 return "tt";
86 case scev_end:
87 return "e";
88 case scev_two_ends:
89 return "ee";
90 case scev_intersection:
91 return "i";
92 case scev_point:
93 return "p";
94 case scev_single_point:
95 return "sp";
96 default:;
97 };
98 GCALC_DBUG_ASSERT(0);
99 return "unk";
100}
101
102
103static int gcalc_pi_str(char *str, const Gcalc_heap::Info *pi, const char *postfix)
104{
105 return sprintf(str, "%s %d %d | %s %d %d%s",
106 GCALC_SIGN(pi->node.shape.ix[0]) ? "-":"", FIRST_DIGIT(pi->node.shape.ix[0]),pi->node.shape.ix[1],
107 GCALC_SIGN(pi->node.shape.iy[0]) ? "-":"", FIRST_DIGIT(pi->node.shape.iy[0]),pi->node.shape.iy[1],
108 postfix);
109
110}
111
112
113static void GCALC_DBUG_PRINT_PI(const Gcalc_heap::Info *pi)
114{
115 char buf[128];
116 int n_buf;
117 if (pi->type == Gcalc_heap::nt_intersection)
118 {
119 const Gcalc_scan_iterator::intersection_info *ic= i_data(pi);
120
121 GCALC_DBUG_PRINT(("intersection point %d %d",
122 ic->edge_a->thread, ic->edge_b->thread));
123 return;
124 }
125 if (pi->type == Gcalc_heap::nt_eq_node)
126 {
127 const Gcalc_scan_iterator::point *e= eq_sp(pi);
128 GCALC_DBUG_PRINT(("eq point %d", e->thread));
129 return;
130 }
131 n_buf= gcalc_pi_str(buf, pi, "");
132 buf[n_buf]= 0;
133 GCALC_DBUG_PRINT(("%s", buf));
134}
135
136
137static void GCALC_DBUG_PRINT_SLICE(const char *header,
138 const Gcalc_scan_iterator::point *slice)
139{
140 size_t nbuf;
141 char buf[1024];
142 nbuf= strlen(header);
143 strcpy(buf, header);
144 for (; slice; slice= slice->get_next())
145 {
146 size_t lnbuf= nbuf;
147 lnbuf+= sprintf(buf + lnbuf, "%d\t", slice->thread);
148 lnbuf+= sprintf(buf + lnbuf, "%s\t", gcalc_ev_name(slice->event));
149
150 lnbuf+= gcalc_pi_str(buf + lnbuf, slice->pi, "\t");
151 if (slice->is_bottom())
152 lnbuf+= sprintf(buf+lnbuf, "bt\t");
153 else
154 lnbuf+= gcalc_pi_str(buf+lnbuf, slice->next_pi, "\t");
155 buf[lnbuf]= 0;
156 GCALC_DBUG_PRINT(("%s", buf));
157 }
158}
159
160
161#else
162#define GCALC_DBUG_CHECK_COUNTER() do { } while(0)
163#define GCALC_DBUG_PRINT_PI(pi) do { } while(0)
164#define GCALC_DBUG_PRINT_SLICE(a, b) do { } while(0)
165#define GCALC_DBUG_PRINT_INTERSECTIONS(a) do { } while(0)
166#define GCALC_DBUG_PRINT_STATE(a) do { } while(0)
167#endif /*GCALC_DBUG_OFF*/
168
169
170Gcalc_dyn_list::Gcalc_dyn_list(size_t blk_size, size_t sizeof_item):
171 m_blk_size(blk_size - ALLOC_ROOT_MIN_BLOCK_SIZE),
172 m_sizeof_item(ALIGN_SIZE(sizeof_item)),
173 m_points_per_blk((uint)((m_blk_size - PH_DATA_OFFSET) / m_sizeof_item)),
174 m_blk_hook(&m_first_blk),
175 m_free(NULL),
176 m_keep(NULL)
177{}
178
179
180void Gcalc_dyn_list::format_blk(void* block)
181{
182 Item *pi_end, *cur_pi, *first_pi;
183 GCALC_DBUG_ASSERT(m_free == NULL);
184 first_pi= cur_pi= (Item *)(((char *)block) + PH_DATA_OFFSET);
185 pi_end= ptr_add(first_pi, m_points_per_blk - 1);
186 do {
187 cur_pi= cur_pi->next= ptr_add(cur_pi, 1);
188 } while (cur_pi<pi_end);
189 cur_pi->next= m_free;
190 m_free= first_pi;
191}
192
193
194Gcalc_dyn_list::Item *Gcalc_dyn_list::alloc_new_blk()
195{
196 void *new_block= my_malloc(m_blk_size, MYF(MY_WME));
197 if (!new_block)
198 return NULL;
199 *m_blk_hook= new_block;
200 m_blk_hook= (void**)new_block;
201 format_blk(new_block);
202 return new_item();
203}
204
205
206static void free_blk_list(void *list)
207{
208 void *next_blk;
209 while (list)
210 {
211 next_blk= *((void **)list);
212 my_free(list);
213 list= next_blk;
214 }
215}
216
217
218void Gcalc_dyn_list::cleanup()
219{
220 *m_blk_hook= NULL;
221 free_blk_list(m_first_blk);
222 m_first_blk= NULL;
223 m_blk_hook= &m_first_blk;
224 m_free= NULL;
225}
226
227
228Gcalc_dyn_list::~Gcalc_dyn_list()
229{
230 cleanup();
231}
232
233
234void Gcalc_dyn_list::reset()
235{
236 *m_blk_hook= NULL;
237 if (m_first_blk)
238 {
239 free_blk_list(*((void **)m_first_blk));
240 m_blk_hook= (void**)m_first_blk;
241 m_free= NULL;
242 format_blk(m_first_blk);
243 }
244}
245
246
247/* Internal coordinate operations implementations */
248
249void gcalc_set_zero(Gcalc_internal_coord *d, int d_len)
250{
251 do
252 {
253 d[--d_len]= 0;
254 } while (d_len);
255}
256
257
258int gcalc_is_zero(const Gcalc_internal_coord *d, int d_len)
259{
260 do
261 {
262 if (d[--d_len] != 0)
263 return 0;
264 } while (d_len);
265 return 1;
266}
267
268
269#ifdef GCALC_CHECK_WITH_FLOAT
270static double *gcalc_coord_extent= NULL;
271
272long double gcalc_get_double(const Gcalc_internal_coord *d, int d_len)
273{
274 int n= 1;
275 long double res= (long double) FIRST_DIGIT(d[0]);
276 do
277 {
278 res*= (long double) GCALC_DIG_BASE;
279 res+= (long double) d[n];
280 } while(++n < d_len);
281
282 n= 0;
283 do
284 {
285 if ((n & 1) && gcalc_coord_extent)
286 res/= *gcalc_coord_extent;
287 } while(++n < d_len);
288
289 if (GCALC_SIGN(d[0]))
290 res*= -1.0;
291 return res;
292}
293#endif /*GCALC_CHECK_WITH_FLOAT*/
294
295
296static void do_add(Gcalc_internal_coord *result, int result_len,
297 const Gcalc_internal_coord *a,
298 const Gcalc_internal_coord *b)
299{
300 int n_digit= result_len-1;
301 gcalc_digit_t carry= 0;
302
303 do
304 {
305 if ((result[n_digit]=
306 a[n_digit] + b[n_digit] + carry) >= GCALC_DIG_BASE)
307 {
308 carry= 1;
309 result[n_digit]-= GCALC_DIG_BASE;
310 }
311 else
312 carry= 0;
313 } while (--n_digit);
314
315 result[0]= (a[0] + FIRST_DIGIT(b[0]) + carry);
316
317 GCALC_DBUG_ASSERT(FIRST_DIGIT(result[0]) < GCALC_DIG_BASE);
318}
319
320
321static void do_sub(Gcalc_internal_coord *result, int result_len,
322 const Gcalc_internal_coord *a,
323 const Gcalc_internal_coord *b)
324{
325 int n_digit= result_len-1;
326 gcalc_digit_t carry= 0;
327 gcalc_digit_t cur_b, cur_a;
328
329 do
330 {
331 cur_b= b[n_digit] + carry;
332 cur_a= a[n_digit];
333 if (cur_a < cur_b)
334 {
335 carry= 1;
336 result[n_digit]= (GCALC_DIG_BASE - cur_b) + cur_a;
337 }
338 else
339 {
340 carry= 0;
341 result[n_digit]= cur_a - cur_b;
342 }
343 } while (--n_digit);
344
345
346 result[0]= a[0] - FIRST_DIGIT(b[0]) - carry;
347
348 GCALC_DBUG_ASSERT(FIRST_DIGIT(a[0]) >= FIRST_DIGIT(b[0]) + carry);
349 GCALC_DBUG_ASSERT(!gcalc_is_zero(result, result_len));
350}
351/*
352static void do_sub(Gcalc_internal_coord *result, int result_len,
353 const Gcalc_internal_coord *a,
354 const Gcalc_internal_coord *b)
355{
356 int n_digit= result_len-1;
357 gcalc_digit_t carry= 0;
358
359 do
360 {
361 if ((result[n_digit]= a[n_digit] - b[n_digit] - carry) < 0)
362 {
363 carry= 1;
364 result[n_digit]+= GCALC_DIG_BASE;
365 }
366 else
367 carry= 0;
368 } while (--n_digit);
369
370
371 result[0]= a[0] - FIRST_DIGIT(b[0]) - carry;
372
373 GCALC_DBUG_ASSERT(FIRST_DIGIT(a[0]) - FIRST_DIGIT(b[0]) - carry >= 0);
374 GCALC_DBUG_ASSERT(!gcalc_is_zero(result, result_len));
375}
376*/
377
378static int do_cmp(const Gcalc_internal_coord *a,
379 const Gcalc_internal_coord *b, int len)
380{
381 int n_digit= 1;
382
383 if ((FIRST_DIGIT(a[0]) != FIRST_DIGIT(b[0])))
384 return FIRST_DIGIT(a[0]) > FIRST_DIGIT(b[0]) ? 1 : -1;
385
386 do
387 {
388 if ((a[n_digit] != b[n_digit]))
389 return a[n_digit] > b[n_digit] ? 1 : -1;
390 } while (++n_digit < len);
391
392 return 0;
393}
394
395
396#ifdef GCALC_CHECK_WITH_FLOAT
397static int de_weak_check(long double a, long double b, long double ex)
398{
399 long double d= a - b;
400 if (d < ex && d > -ex)
401 return 1;
402
403 d/= fabsl(a) + fabsl(b);
404 if (d < ex && d > -ex)
405 return 1;
406 return 0;
407}
408
409static int de_check(long double a, long double b)
410{
411 return de_weak_check(a, b, (long double) 1e-9);
412}
413#endif /*GCALC_CHECK_WITH_FLOAT*/
414
415
416void gcalc_mul_coord(Gcalc_internal_coord *result, int result_len,
417 const Gcalc_internal_coord *a, int a_len,
418 const Gcalc_internal_coord *b, int b_len)
419{
420 GCALC_DBUG_ASSERT(result_len == a_len + b_len);
421 GCALC_DBUG_ASSERT(a_len >= b_len);
422 int n_a, n_b, n_res;
423 gcalc_digit_t carry= 0;
424
425 gcalc_set_zero(result, result_len);
426
427 n_a= a_len - 1;
428 do
429 {
430 gcalc_coord2 cur_a= n_a ? a[n_a] : FIRST_DIGIT(a[0]);
431 n_b= b_len - 1;
432 do
433 {
434 gcalc_coord2 cur_b= n_b ? b[n_b] : FIRST_DIGIT(b[0]);
435 gcalc_coord2 mul= cur_a * cur_b + carry + result[n_a + n_b + 1];
436 result[n_a + n_b + 1]= mul % GCALC_DIG_BASE;
437 carry= (gcalc_digit_t) (mul / (gcalc_coord2) GCALC_DIG_BASE);
438 } while (n_b--);
439 if (carry)
440 {
441 for (n_res= n_a; (result[n_res]+= carry) >= GCALC_DIG_BASE;
442 n_res--)
443 {
444 result[n_res]-= GCALC_DIG_BASE;
445 carry= 1;
446 }
447 carry= 0;
448 }
449 } while (n_a--);
450 if (!gcalc_is_zero(result, result_len))
451 result[0]|= GCALC_SIGN(a[0] ^ b[0]);
452#ifdef GCALC_CHECK_WITH_FLOAT
453 GCALC_DBUG_ASSERT(de_check(gcalc_get_double(a, a_len) *
454 gcalc_get_double(b, b_len),
455 gcalc_get_double(result, result_len)));
456#endif /*GCALC_CHECK_WITH_FLOAT*/
457}
458
459
460inline void gcalc_mul_coord1(Gcalc_coord1 result,
461 const Gcalc_coord1 a, const Gcalc_coord1 b)
462{
463 return gcalc_mul_coord(result, GCALC_COORD_BASE2,
464 a, GCALC_COORD_BASE, b, GCALC_COORD_BASE);
465}
466
467
468void gcalc_add_coord(Gcalc_internal_coord *result, int result_len,
469 const Gcalc_internal_coord *a,
470 const Gcalc_internal_coord *b)
471{
472 if (GCALC_SIGN(a[0]) == GCALC_SIGN(b[0]))
473 do_add(result, result_len, a, b);
474 else
475 {
476 int cmp_res= do_cmp(a, b, result_len);
477 if (cmp_res == 0)
478 gcalc_set_zero(result, result_len);
479 else if (cmp_res > 0)
480 do_sub(result, result_len, a, b);
481 else
482 do_sub(result, result_len, b, a);
483 }
484#ifdef GCALC_CHECK_WITH_FLOAT
485 GCALC_DBUG_ASSERT(de_check(gcalc_get_double(a, result_len) +
486 gcalc_get_double(b, result_len),
487 gcalc_get_double(result, result_len)));
488#endif /*GCALC_CHECK_WITH_FLOAT*/
489}
490
491
492void gcalc_sub_coord(Gcalc_internal_coord *result, int result_len,
493 const Gcalc_internal_coord *a,
494 const Gcalc_internal_coord *b)
495{
496 if (GCALC_SIGN(a[0] ^ b[0]))
497 do_add(result, result_len, a, b);
498 else
499 {
500 int cmp_res= do_cmp(a, b, result_len);
501 if (cmp_res == 0)
502 gcalc_set_zero(result, result_len);
503 else if (cmp_res > 0)
504 do_sub(result, result_len, a, b);
505 else
506 {
507 do_sub(result, result_len, b, a);
508 result[0]^= GCALC_COORD_MINUS;
509 }
510 }
511#ifdef GCALC_CHECK_WITH_FLOAT
512 GCALC_DBUG_ASSERT(de_check(gcalc_get_double(a, result_len) -
513 gcalc_get_double(b, result_len),
514 gcalc_get_double(result, result_len)));
515#endif /*GCALC_CHECK_WITH_FLOAT*/
516}
517
518
519inline void gcalc_sub_coord1(Gcalc_coord1 result,
520 const Gcalc_coord1 a, const Gcalc_coord1 b)
521{
522 return gcalc_sub_coord(result, GCALC_COORD_BASE, a, b);
523}
524
525
526int gcalc_cmp_coord(const Gcalc_internal_coord *a,
527 const Gcalc_internal_coord *b, int len)
528{
529 int n_digit= 0;
530 int result= 0;
531
532 do
533 {
534 if (a[n_digit] == b[n_digit])
535 {
536 n_digit++;
537 continue;
538 }
539 if (a[n_digit] > b[n_digit])
540 result= GCALC_SIGN(a[0]) ? -1 : 1;
541 else
542 result= GCALC_SIGN(b[0]) ? 1 : -1;
543 break;
544
545 } while (n_digit < len);
546
547#ifdef GCALC_CHECK_WITH_FLOAT
548 if (result == 0)
549 GCALC_DBUG_ASSERT(de_check(gcalc_get_double(a, len),
550 gcalc_get_double(b, len)));
551 else if (result == 1)
552 GCALC_DBUG_ASSERT(de_check(gcalc_get_double(a, len),
553 gcalc_get_double(b, len)) ||
554 gcalc_get_double(a, len) > gcalc_get_double(b, len));
555 else
556 GCALC_DBUG_ASSERT(de_check(gcalc_get_double(a, len),
557 gcalc_get_double(b, len)) ||
558 gcalc_get_double(a, len) < gcalc_get_double(b, len));
559#endif /*GCALC_CHECK_WITH_FLOAT*/
560 return result;
561}
562
563
564#define gcalc_cmp_coord1(a, b) gcalc_cmp_coord(a, b, GCALC_COORD_BASE)
565
566int gcalc_set_double(Gcalc_internal_coord *c, double d, double ext)
567{
568 int sign;
569 double ds= d * ext;
570 if ((sign= ds < 0))
571 ds= -ds;
572 c[0]= (gcalc_digit_t) (ds / (double) GCALC_DIG_BASE);
573 c[1]= (gcalc_digit_t) (ds - ((double) c[0]) * (double) GCALC_DIG_BASE);
574 if (c[1] >= GCALC_DIG_BASE)
575 {
576 c[1]= 0;
577 c[0]++;
578 }
579 if (sign && (c[0] | c[1]))
580 c[0]|= GCALC_COORD_MINUS;
581#ifdef GCALC_CHECK_WITH_FLOAT
582 GCALC_DBUG_ASSERT(de_check(d, gcalc_get_double(c, 2)));
583#endif /*GCALC_CHECK_WITH_FLOAT*/
584 return 0;
585}
586
587
588typedef gcalc_digit_t Gcalc_coord4[GCALC_COORD_BASE*4];
589typedef gcalc_digit_t Gcalc_coord5[GCALC_COORD_BASE*5];
590
591
592void Gcalc_scan_iterator::intersection_info::do_calc_t()
593{
594 Gcalc_coord1 a2_a1x, a2_a1y;
595 Gcalc_coord2 x1y2, x2y1;
596
597 gcalc_sub_coord1(a2_a1x, edge_b->pi->node.shape.ix, edge_a->pi->node.shape.ix);
598 gcalc_sub_coord1(a2_a1y, edge_b->pi->node.shape.iy, edge_a->pi->node.shape.iy);
599
600 GCALC_DBUG_ASSERT(!gcalc_is_zero(edge_a->dy, GCALC_COORD_BASE) ||
601 !gcalc_is_zero(edge_b->dy, GCALC_COORD_BASE));
602
603 gcalc_mul_coord1(x1y2, edge_a->dx, edge_b->dy);
604 gcalc_mul_coord1(x2y1, edge_a->dy, edge_b->dx);
605 gcalc_sub_coord(t_b, GCALC_COORD_BASE2, x1y2, x2y1);
606
607
608 gcalc_mul_coord1(x1y2, a2_a1x, edge_b->dy);
609 gcalc_mul_coord1(x2y1, a2_a1y, edge_b->dx);
610 gcalc_sub_coord(t_a, GCALC_COORD_BASE2, x1y2, x2y1);
611 t_calculated= 1;
612}
613
614
615void Gcalc_scan_iterator::intersection_info::do_calc_y()
616{
617 GCALC_DBUG_ASSERT(t_calculated);
618
619 Gcalc_coord3 a_tb, b_ta;
620
621 gcalc_mul_coord(a_tb, GCALC_COORD_BASE3,
622 t_b, GCALC_COORD_BASE2, edge_a->pi->node.shape.iy, GCALC_COORD_BASE);
623 gcalc_mul_coord(b_ta, GCALC_COORD_BASE3,
624 t_a, GCALC_COORD_BASE2, edge_a->dy, GCALC_COORD_BASE);
625
626 gcalc_add_coord(y_exp, GCALC_COORD_BASE3, a_tb, b_ta);
627 y_calculated= 1;
628}
629
630
631void Gcalc_scan_iterator::intersection_info::do_calc_x()
632{
633 GCALC_DBUG_ASSERT(t_calculated);
634
635 Gcalc_coord3 a_tb, b_ta;
636
637 gcalc_mul_coord(a_tb, GCALC_COORD_BASE3,
638 t_b, GCALC_COORD_BASE2, edge_a->pi->node.shape.ix, GCALC_COORD_BASE);
639 gcalc_mul_coord(b_ta, GCALC_COORD_BASE3,
640 t_a, GCALC_COORD_BASE2, edge_a->dx, GCALC_COORD_BASE);
641
642 gcalc_add_coord(x_exp, GCALC_COORD_BASE3, a_tb, b_ta);
643 x_calculated= 1;
644}
645
646
647static int cmp_node_isc(const Gcalc_heap::Info *node,
648 const Gcalc_heap::Info *isc)
649{
650 GCALC_DBUG_ASSERT(node->type == Gcalc_heap::nt_shape_node);
651 Gcalc_scan_iterator::intersection_info *inf= i_data(isc);
652 Gcalc_coord3 exp;
653 int result;
654
655 inf->calc_t();
656 inf->calc_y_exp();
657
658 gcalc_mul_coord(exp, GCALC_COORD_BASE3,
659 inf->t_b, GCALC_COORD_BASE2, node->node.shape.iy, GCALC_COORD_BASE);
660
661 result= gcalc_cmp_coord(exp, inf->y_exp, GCALC_COORD_BASE3);
662#ifdef GCALC_CHECK_WITH_FLOAT
663 long double int_x, int_y;
664 isc->calc_xy_ld(&int_x, &int_y);
665 if (result < 0)
666 {
667 if (!de_check(int_y, node->node.shape.y) && node->node.shape.y > int_y)
668 GCALC_DBUG_PRINT(("floatcheck cmp_nod_iscy %g < %LG", node->node.shape.y, int_y));
669 }
670 else if (result > 0)
671 {
672 if (!de_check(int_y, node->node.shape.y) && node->node.shape.y < int_y)
673 GCALC_DBUG_PRINT(("floatcheck cmp_nod_iscy %g > %LG", node->node.shape.y, int_y));
674 }
675 else
676 {
677 if (!de_check(int_y, node->node.shape.y))
678 GCALC_DBUG_PRINT(("floatcheck cmp_nod_iscy %g == %LG", node->node.shape.y, int_y));
679 }
680#endif /*GCALC_CHECK_WITH_FLOAT*/
681 if (result)
682 goto exit;
683
684
685 inf->calc_x_exp();
686 gcalc_mul_coord(exp, GCALC_COORD_BASE3,
687 inf->t_b, GCALC_COORD_BASE2, node->node.shape.ix, GCALC_COORD_BASE);
688
689 result= gcalc_cmp_coord(exp, inf->x_exp, GCALC_COORD_BASE3);
690#ifdef GCALC_CHECK_WITH_FLOAT
691 if (result < 0)
692 {
693 if (!de_check(int_x, node->node.shape.x) && node->node.shape.x > int_x)
694 GCALC_DBUG_PRINT(("floatcheck cmp_nod_iscx failed %g < %LG",
695 node->node.shape.x, int_x));
696 }
697 else if (result > 0)
698 {
699 if (!de_check(int_x, node->node.shape.x) && node->node.shape.x < int_x)
700 GCALC_DBUG_PRINT(("floatcheck cmp_nod_iscx failed %g > %LG",
701 node->node.shape.x, int_x));
702 }
703 else
704 {
705 if (!de_check(int_x, node->node.shape.x))
706 GCALC_DBUG_PRINT(("floatcheck cmp_nod_iscx failed %g == %LG",
707 node->node.shape.x, int_x));
708 }
709#endif /*GCALC_CHECK_WITH_FLOAT*/
710exit:
711 return result;
712}
713
714
715static int cmp_intersections(const Gcalc_heap::Info *i1,
716 const Gcalc_heap::Info *i2)
717{
718 Gcalc_scan_iterator::intersection_info *inf1= i_data(i1);
719 Gcalc_scan_iterator::intersection_info *inf2= i_data(i2);
720 Gcalc_coord5 exp_a, exp_b;
721 int result;
722
723 inf1->calc_t();
724 inf2->calc_t();
725
726 inf1->calc_y_exp();
727 inf2->calc_y_exp();
728
729 gcalc_mul_coord(exp_a, GCALC_COORD_BASE5,
730 inf1->y_exp, GCALC_COORD_BASE3, inf2->t_b, GCALC_COORD_BASE2);
731 gcalc_mul_coord(exp_b, GCALC_COORD_BASE5,
732 inf2->y_exp, GCALC_COORD_BASE3, inf1->t_b, GCALC_COORD_BASE2);
733
734 result= gcalc_cmp_coord(exp_a, exp_b, GCALC_COORD_BASE5);
735#ifdef GCALC_CHECK_WITH_FLOAT
736 long double x1, y1, x2, y2;
737 i1->calc_xy_ld(&x1, &y1);
738 i2->calc_xy_ld(&x2, &y2);
739
740 if (result < 0)
741 {
742 if (!de_check(y1, y2) && y2 > y1)
743 GCALC_DBUG_PRINT(("floatcheck cmp_intersections_y failed %LG < %LG",
744 y2, y1));
745 }
746 else if (result > 0)
747 {
748 if (!de_check(y1, y2) && y2 < y1)
749 GCALC_DBUG_PRINT(("floatcheck cmp_intersections_y failed %LG > %LG",
750 y2, y1));
751 }
752 else
753 {
754 if (!de_check(y1, y2))
755 GCALC_DBUG_PRINT(("floatcheck cmp_intersections_y failed %LG == %LG",
756 y2, y1));
757 }
758#endif /*GCALC_CHECK_WITH_FLOAT*/
759
760 if (result != 0)
761 return result;
762
763
764 inf1->calc_x_exp();
765 inf2->calc_x_exp();
766 gcalc_mul_coord(exp_a, GCALC_COORD_BASE5,
767 inf1->x_exp, GCALC_COORD_BASE3, inf2->t_b, GCALC_COORD_BASE2);
768 gcalc_mul_coord(exp_b, GCALC_COORD_BASE5,
769 inf2->x_exp, GCALC_COORD_BASE3, inf1->t_b, GCALC_COORD_BASE2);
770
771 result= gcalc_cmp_coord(exp_a, exp_b, GCALC_COORD_BASE5);
772#ifdef GCALC_CHECK_WITH_FLOAT
773 if (result < 0)
774 {
775 if (!de_check(x1, x2) && x2 > x1)
776 GCALC_DBUG_PRINT(("floatcheck cmp_intersectionsx failed %LG < %LG",
777 x2, x1));
778 }
779 else if (result > 0)
780 {
781 if (!de_check(x1, x2) && x2 < x1)
782 GCALC_DBUG_PRINT(("floatcheck cmp_intersectionsx failed %LG > %LG",
783 x2, x1));
784 }
785 else
786 {
787 if (!de_check(x1, x2))
788 GCALC_DBUG_PRINT(("floatcheck cmp_intersectionsx failed %LG == %LG",
789 x2, x1));
790 }
791#endif /*GCALC_CHECK_WITH_FLOAT*/
792 return result;
793}
794/* Internal coordinates implementation end */
795
796
797#define GCALC_SCALE_1 1e18
798
799static double find_scale(double extent)
800{
801 double scale= 1e-2;
802 while (scale < extent)
803 scale*= (double ) 10;
804 return GCALC_SCALE_1 / scale / 10;
805}
806
807
808void Gcalc_heap::set_extent(double xmin, double xmax, double ymin, double ymax)
809{
810 xmin= fabs(xmin);
811 xmax= fabs(xmax);
812 ymin= fabs(ymin);
813 ymax= fabs(ymax);
814
815 if (xmax < xmin)
816 xmax= xmin;
817 if (ymax < ymin)
818 ymax= ymin;
819
820 coord_extent= xmax > ymax ? xmax : ymax;
821 coord_extent= find_scale(coord_extent);
822#ifdef GCALC_CHECK_WITH_FLOAT
823 gcalc_coord_extent= &coord_extent;
824#endif /*GCALC_CHECK_WITH_FLOAT*/
825}
826
827
828void Gcalc_heap::free_point_info(Gcalc_heap::Info *i,
829 Gcalc_dyn_list::Item **i_hook)
830{
831 if (m_hook == &i->next)
832 m_hook= i_hook;
833 *i_hook= i->next;
834 free_item(i);
835 m_n_points--;
836}
837
838
839Gcalc_heap::Info *Gcalc_heap::new_point_info(double x, double y,
840 gcalc_shape_info shape)
841{
842 Info *result= (Info *)new_item();
843 if (!result)
844 return NULL;
845 *m_hook= result;
846 m_hook= &result->next;
847 result->node.shape.x= x;
848 result->node.shape.y= y;
849 result->node.shape.shape= shape;
850 result->node.shape.top_node= 1;
851 result->type= nt_shape_node;
852 gcalc_set_double(result->node.shape.ix, x, coord_extent);
853 gcalc_set_double(result->node.shape.iy, y, coord_extent);
854
855 m_n_points++;
856 return result;
857}
858
859
860static Gcalc_heap::Info *new_intersection(
861 Gcalc_heap *heap, Gcalc_scan_iterator::intersection_info *ii)
862{
863 Gcalc_heap::Info *isc= (Gcalc_heap::Info *)heap->new_item();
864 if (!isc)
865 return 0;
866 isc->type= Gcalc_heap::nt_intersection;
867 isc->node.intersection.p1= ii->edge_a->pi;
868 isc->node.intersection.p2= ii->edge_a->next_pi;
869 isc->node.intersection.p3= ii->edge_b->pi;
870 isc->node.intersection.p4= ii->edge_b->next_pi;
871 isc->node.intersection.data= ii;
872 return isc;
873}
874
875
876static Gcalc_heap::Info *new_eq_point(
877 Gcalc_heap *heap, const Gcalc_heap::Info *p,
878 Gcalc_scan_iterator::point *edge)
879{
880 Gcalc_heap::Info *eqp= (Gcalc_heap::Info *)heap->new_item();
881 if (!eqp)
882 return 0;
883 eqp->type= Gcalc_heap::nt_eq_node;
884 eqp->node.eq.node= p;
885 eqp->node.eq.data= edge;
886 return eqp;
887}
888
889
890void Gcalc_heap::Info::calc_xy(double *x, double *y) const
891{
892 double b0_x= node.intersection.p2->node.shape.x - node.intersection.p1->node.shape.x;
893 double b0_y= node.intersection.p2->node.shape.y - node.intersection.p1->node.shape.y;
894 double b1_x= node.intersection.p4->node.shape.x - node.intersection.p3->node.shape.x;
895 double b1_y= node.intersection.p4->node.shape.y - node.intersection.p3->node.shape.y;
896 double b0xb1= b0_x * b1_y - b0_y * b1_x;
897 double t= (node.intersection.p3->node.shape.x - node.intersection.p1->node.shape.x) * b1_y - (node.intersection.p3->node.shape.y - node.intersection.p1->node.shape.y) * b1_x;
898
899 t/= b0xb1;
900
901 *x= node.intersection.p1->node.shape.x + b0_x * t;
902 *y= node.intersection.p1->node.shape.y + b0_y * t;
903}
904
905
906#ifdef GCALC_CHECK_WITH_FLOAT
907void Gcalc_heap::Info::calc_xy_ld(long double *x, long double *y) const
908{
909 long double b0_x= ((long double) p2->node.shape.x) - p1->node.shape.x;
910 long double b0_y= ((long double) p2->node.shape.y) - p1->node.shape.y;
911 long double b1_x= ((long double) p4->node.shape.x) - p3->node.shape.x;
912 long double b1_y= ((long double) p4->node.shape.y) - p3->node.shape.y;
913 long double b0xb1= b0_x * b1_y - b0_y * b1_x;
914 long double ax= ((long double) p3->node.shape.x) - p1->node.shape.x;
915 long double ay= ((long double) p3->node.shape.y) - p1->node.shape.y;
916 long double t_a= ax * b1_y - ay * b1_x;
917 long double hx= (b0xb1 * (long double) p1->node.shape.x + b0_x * t_a);
918 long double hy= (b0xb1 * (long double) p1->node.shape.y + b0_y * t_a);
919
920 if (fabs(b0xb1) < 1e-15)
921 {
922 *x= p1->node.shape.x;
923 *y= p1->node.shape.y;
924 return;
925 }
926
927 *x= hx/b0xb1;
928 *y= hy/b0xb1;
929}
930#endif /*GCALC_CHECK_WITH_FLOAT*/
931
932
933static int cmp_point_info(const Gcalc_heap::Info *i0,
934 const Gcalc_heap::Info *i1)
935{
936 int cmp_y= gcalc_cmp_coord1(i0->node.shape.iy, i1->node.shape.iy);
937 if (cmp_y)
938 return cmp_y;
939 return gcalc_cmp_coord1(i0->node.shape.ix, i1->node.shape.ix);
940}
941
942
943static inline void trim_node(Gcalc_heap::Info *node, Gcalc_heap::Info *prev_node)
944{
945 if (!node)
946 return;
947 node->node.shape.top_node= 0;
948 GCALC_DBUG_ASSERT((node->node.shape.left == prev_node) || (node->node.shape.right == prev_node));
949 if (node->node.shape.left == prev_node)
950 node->node.shape.left= node->node.shape.right;
951 node->node.shape.right= NULL;
952 GCALC_DBUG_ASSERT(cmp_point_info(node, prev_node));
953}
954
955
956static int compare_point_info(const void *e0, const void *e1)
957{
958 const Gcalc_heap::Info *i0= (const Gcalc_heap::Info *)e0;
959 const Gcalc_heap::Info *i1= (const Gcalc_heap::Info *)e1;
960 return cmp_point_info(i0, i1) > 0;
961}
962
963
964void Gcalc_heap::prepare_operation()
965{
966 Info *cur;
967 GCALC_DBUG_ASSERT(m_hook);
968 *m_hook= NULL;
969 m_hook= NULL; /* just to check it's not called twice */
970 m_first= sort_list(compare_point_info, m_first, m_n_points);
971
972 /* TODO - move this to the 'normal_scan' loop */
973 for (cur= get_first(); cur; cur= cur->get_next())
974 {
975 trim_node(cur->node.shape.left, cur);
976 trim_node(cur->node.shape.right, cur);
977 }
978}
979
980
981void Gcalc_heap::reset()
982{
983 if (m_n_points)
984 {
985 free_list(m_first);
986 m_n_points= 0;
987 }
988 m_hook= &m_first;
989}
990
991
992int Gcalc_shape_transporter::int_single_point(gcalc_shape_info Info,
993 double x, double y)
994{
995 Gcalc_heap::Info *point= m_heap->new_point_info(x, y, Info);
996 if (!point)
997 return 1;
998 point->node.shape.left= point->node.shape.right= 0;
999 return 0;
1000}
1001
1002
1003int Gcalc_shape_transporter::int_add_point(gcalc_shape_info Info,
1004 double x, double y)
1005{
1006 Gcalc_heap::Info *point;
1007 Gcalc_dyn_list::Item **hook;
1008
1009 hook= m_heap->get_cur_hook();
1010
1011 if (!(point= m_heap->new_point_info(x, y, Info)))
1012 return 1;
1013 if (m_first)
1014 {
1015 if (cmp_point_info(m_prev, point) == 0)
1016 {
1017 /* Coinciding points, do nothing */
1018 m_heap->free_point_info(point, hook);
1019 return 0;
1020 }
1021 GCALC_DBUG_ASSERT(!m_prev || m_prev->node.shape.x != x || m_prev->node.shape.y != y);
1022 m_prev->node.shape.left= point;
1023 point->node.shape.right= m_prev;
1024 }
1025 else
1026 m_first= point;
1027 m_prev= point;
1028 m_prev_hook= hook;
1029 return 0;
1030}
1031
1032
1033void Gcalc_shape_transporter::int_complete()
1034{
1035 GCALC_DBUG_ASSERT(m_shape_started == 1 || m_shape_started == 3);
1036
1037 if (!m_first)
1038 return;
1039
1040 /* simple point */
1041 if (m_first == m_prev)
1042 {
1043 m_first->node.shape.right= m_first->node.shape.left= NULL;
1044 return;
1045 }
1046
1047 /* line */
1048 if (m_shape_started == 1)
1049 {
1050 m_first->node.shape.right= NULL;
1051 m_prev->node.shape.left= m_prev->node.shape.right;
1052 m_prev->node.shape.right= NULL;
1053 return;
1054 }
1055
1056 /* polygon */
1057 if (cmp_point_info(m_first, m_prev) == 0)
1058 {
1059 /* Coinciding points, remove the last one from the list */
1060 m_prev->node.shape.right->node.shape.left= m_first;
1061 m_first->node.shape.right= m_prev->node.shape.right;
1062 m_heap->free_point_info(m_prev, m_prev_hook);
1063 }
1064 else
1065 {
1066 GCALC_DBUG_ASSERT(m_prev->node.shape.x != m_first->node.shape.x || m_prev->node.shape.y != m_first->node.shape.y);
1067 m_first->node.shape.right= m_prev;
1068 m_prev->node.shape.left= m_first;
1069 }
1070}
1071
1072
1073inline void calc_dx_dy(Gcalc_scan_iterator::point *p)
1074{
1075 gcalc_sub_coord1(p->dx, p->next_pi->node.shape.ix, p->pi->node.shape.ix);
1076 gcalc_sub_coord1(p->dy, p->next_pi->node.shape.iy, p->pi->node.shape.iy);
1077 if (GCALC_SIGN(p->dx[0]))
1078 {
1079 p->l_border= &p->next_pi->node.shape.ix;
1080 p->r_border= &p->pi->node.shape.ix;
1081 }
1082 else
1083 {
1084 p->r_border= &p->next_pi->node.shape.ix;
1085 p->l_border= &p->pi->node.shape.ix;
1086 }
1087}
1088
1089
1090Gcalc_scan_iterator::Gcalc_scan_iterator(size_t blk_size) :
1091 Gcalc_dyn_list(blk_size, sizeof(point) > sizeof(intersection_info) ?
1092 sizeof(point) :
1093 sizeof(intersection_info))
1094{
1095 state.slice= NULL;
1096 m_bottom_points= NULL;
1097 m_bottom_hook= &m_bottom_points;
1098}
1099
1100
1101void Gcalc_scan_iterator::init(Gcalc_heap *points)
1102{
1103 GCALC_DBUG_ASSERT(points->ready());
1104 GCALC_DBUG_ASSERT(!state.slice);
1105
1106 if (!(m_cur_pi= points->get_first()))
1107 return;
1108 m_heap= points;
1109 state.event_position_hook= &state.slice;
1110 state.event_end= NULL;
1111#ifndef GCALC_DBUG_OFF
1112 m_cur_thread= 0;
1113#endif /*GCALC_DBUG_OFF*/
1114 GCALC_SET_TERMINATED(killed, 0);
1115}
1116
1117void Gcalc_scan_iterator::reset()
1118{
1119 state.slice= NULL;
1120 m_bottom_points= NULL;
1121 m_bottom_hook= &m_bottom_points;
1122 Gcalc_dyn_list::reset();
1123}
1124
1125
1126int Gcalc_scan_iterator::point::cmp_dx_dy(const Gcalc_coord1 dx_a,
1127 const Gcalc_coord1 dy_a,
1128 const Gcalc_coord1 dx_b,
1129 const Gcalc_coord1 dy_b)
1130{
1131 Gcalc_coord2 dx_a_dy_b;
1132 Gcalc_coord2 dy_a_dx_b;
1133 gcalc_mul_coord1(dx_a_dy_b, dx_a, dy_b);
1134 gcalc_mul_coord1(dy_a_dx_b, dy_a, dx_b);
1135
1136 return gcalc_cmp_coord(dx_a_dy_b, dy_a_dx_b, GCALC_COORD_BASE2);
1137}
1138
1139
1140int Gcalc_scan_iterator::point::cmp_dx_dy(const Gcalc_heap::Info *p1,
1141 const Gcalc_heap::Info *p2,
1142 const Gcalc_heap::Info *p3,
1143 const Gcalc_heap::Info *p4)
1144{
1145 Gcalc_coord1 dx_a, dy_a, dx_b, dy_b;
1146 gcalc_sub_coord1(dx_a, p2->node.shape.ix, p1->node.shape.ix);
1147 gcalc_sub_coord1(dy_a, p2->node.shape.iy, p1->node.shape.iy);
1148 gcalc_sub_coord1(dx_b, p4->node.shape.ix, p3->node.shape.ix);
1149 gcalc_sub_coord1(dy_b, p4->node.shape.iy, p3->node.shape.iy);
1150 return cmp_dx_dy(dx_a, dy_a, dx_b, dy_b);
1151}
1152
1153
1154int Gcalc_scan_iterator::point::cmp_dx_dy(const point *p) const
1155{
1156 GCALC_DBUG_ASSERT(!is_bottom());
1157 return cmp_dx_dy(dx, dy, p->dx, p->dy);
1158}
1159
1160
1161#ifdef GCALC_CHECK_WITH_FLOAT
1162void Gcalc_scan_iterator::point::calc_x(long double *x, long double y,
1163 long double ix) const
1164{
1165 long double ddy= gcalc_get_double(dy, GCALC_COORD_BASE);
1166 if (fabsl(ddy) < (long double) 1e-20)
1167 {
1168 *x= ix;
1169 }
1170 else
1171 *x= (ddy * (long double) pi->node.shape.x + gcalc_get_double(dx, GCALC_COORD_BASE) *
1172 (y - pi->node.shape.y)) / ddy;
1173}
1174#endif /*GCALC_CHECK_WITH_FLOAT*/
1175
1176
1177static int compare_events(const void *e0, const void *e1)
1178{
1179 const Gcalc_scan_iterator::point *p0= (const Gcalc_scan_iterator::point *)e0;
1180 const Gcalc_scan_iterator::point *p1= (const Gcalc_scan_iterator::point *)e1;
1181 return p0->cmp_dx_dy(p1) > 0;
1182}
1183
1184
1185int Gcalc_scan_iterator::arrange_event(int do_sorting, int n_intersections)
1186{
1187 int ev_counter;
1188 point *sp;
1189 point **sp_hook;
1190
1191 ev_counter= 0;
1192
1193 *m_bottom_hook= NULL;
1194 for (sp= m_bottom_points; sp; sp= sp->get_next())
1195 sp->ev_next= sp->get_next();
1196
1197 for (sp= state.slice, sp_hook= &state.slice;
1198 sp; sp_hook= sp->next_ptr(), sp= sp->get_next())
1199 {
1200 if (sp->event)
1201 {
1202 state.event_position_hook= sp_hook;
1203 break;
1204 }
1205 }
1206
1207 for (sp= *(sp_hook= state.event_position_hook);
1208 sp && sp->event; sp_hook= sp->next_ptr(), sp= sp->get_next())
1209 {
1210 ev_counter++;
1211 if (sp->get_next() && sp->get_next()->event)
1212 sp->ev_next= sp->get_next();
1213 else
1214 sp->ev_next= m_bottom_points;
1215 }
1216
1217#ifndef GCALC_DBUG_OFF
1218 {
1219 point *cur_p= sp;
1220 for (; cur_p; cur_p= cur_p->get_next())
1221 GCALC_DBUG_ASSERT(!cur_p->event);
1222 }
1223#endif /*GCALC_DBUG_OFF*/
1224
1225 state.event_end= sp;
1226
1227 if (ev_counter == 2 && n_intersections == 1)
1228 {
1229 /* If we had only intersection, just swap the two points. */
1230 sp= *state.event_position_hook;
1231 *state.event_position_hook= sp->get_next();
1232 sp->next= (*state.event_position_hook)->next;
1233 (*state.event_position_hook)->next= sp;
1234
1235 /* The list of the events should be restored. */
1236 (*state.event_position_hook)->ev_next= sp;
1237 sp->ev_next= m_bottom_points;
1238 }
1239 else if (ev_counter == 2 && get_events()->event == scev_two_threads)
1240 {
1241 /* Do nothing. */
1242 }
1243 else if (ev_counter > 1 && do_sorting)
1244 {
1245 point *cur_p;
1246 *sp_hook= NULL;
1247 sp= (point *) sort_list(compare_events, *state.event_position_hook,
1248 ev_counter);
1249 /* Find last item in the list, it's changed after the sorting. */
1250 for (cur_p= sp->get_next(); cur_p->get_next();
1251 cur_p= cur_p->get_next())
1252 {}
1253 cur_p->next= state.event_end;
1254 *state.event_position_hook= sp;
1255 /* The list of the events should be restored. */
1256 for (; sp && sp->event; sp= sp->get_next())
1257 {
1258 if (sp->get_next() && sp->get_next()->event)
1259 sp->ev_next= sp->get_next();
1260 else
1261 sp->ev_next= m_bottom_points;
1262 }
1263 }
1264
1265#ifndef GCALC_DBUG_OFF
1266 {
1267 const event_point *ev= get_events();
1268 for (; ev && ev->get_next(); ev= ev->get_next())
1269 {
1270 if (ev->is_bottom() || ev->get_next()->is_bottom())
1271 break;
1272 GCALC_DBUG_ASSERT(ev->cmp_dx_dy(ev->get_next()) <= 0);
1273 }
1274 }
1275#endif /*GCALC_DBUG_OFF*/
1276 return 0;
1277}
1278
1279
1280int Gcalc_heap::Info::equal_pi(const Info *pi) const
1281{
1282 if (type == nt_intersection)
1283 return node.intersection.equal;
1284 if (pi->type == nt_eq_node)
1285 return 1;
1286 if (type == nt_eq_node || pi->type == nt_intersection)
1287 return 0;
1288 return cmp_point_info(this, pi) == 0;
1289}
1290
1291int Gcalc_scan_iterator::step()
1292{
1293 int result= 0;
1294 int do_sorting= 0;
1295 int n_intersections= 0;
1296 point *sp;
1297 GCALC_DBUG_ENTER("Gcalc_scan_iterator::step");
1298 GCALC_DBUG_ASSERT(more_points());
1299
1300 if (GCALC_TERMINATED(killed))
1301 GCALC_DBUG_RETURN(0xFFFF);
1302
1303 /* Clear the old event marks. */
1304 if (m_bottom_points)
1305 {
1306 free_list((Gcalc_dyn_list::Item **) &m_bottom_points,
1307 (Gcalc_dyn_list::Item **) m_bottom_hook);
1308 m_bottom_points= NULL;
1309 m_bottom_hook= &m_bottom_points;
1310 }
1311 for (sp= *state.event_position_hook;
1312 sp != state.event_end; sp= sp->get_next())
1313 sp->event= scev_none;
1314
1315//#ifndef GCALC_DBUG_OFF
1316 state.event_position_hook= NULL;
1317 state.pi= NULL;
1318//#endif /*GCALC_DBUG_OFF*/
1319
1320 do
1321 {
1322#ifndef GCALC_DBUG_OFF
1323 if (m_cur_pi->type == Gcalc_heap::nt_intersection &&
1324 m_cur_pi->get_next()->type == Gcalc_heap::nt_intersection &&
1325 m_cur_pi->node.intersection.equal)
1326 GCALC_DBUG_ASSERT(cmp_intersections(m_cur_pi, m_cur_pi->get_next()) == 0);
1327#endif /*GCALC_DBUG_OFF*/
1328 GCALC_DBUG_CHECK_COUNTER();
1329 GCALC_DBUG_PRINT_SLICE("step:", state.slice);
1330 GCALC_DBUG_PRINT_PI(m_cur_pi);
1331 if (m_cur_pi->type == Gcalc_heap::nt_shape_node)
1332 {
1333 if (m_cur_pi->is_top())
1334 {
1335 result= insert_top_node();
1336 if (!m_cur_pi->is_bottom())
1337 do_sorting++;
1338 }
1339 else if (m_cur_pi->is_bottom())
1340 remove_bottom_node();
1341 else
1342 {
1343 do_sorting++;
1344 result= node_scan();
1345 }
1346 if (result)
1347 GCALC_DBUG_RETURN(result);
1348 state.pi= m_cur_pi;
1349 }
1350 else if (m_cur_pi->type == Gcalc_heap::nt_eq_node)
1351 {
1352 do_sorting++;
1353 eq_scan();
1354 }
1355 else
1356 {
1357 /* nt_intersection */
1358 do_sorting++;
1359 n_intersections++;
1360 intersection_scan();
1361 if (!state.pi || state.pi->type == Gcalc_heap::nt_intersection)
1362 state.pi= m_cur_pi;
1363 }
1364
1365 m_cur_pi= m_cur_pi->get_next();
1366 } while (m_cur_pi && state.pi->equal_pi(m_cur_pi));
1367
1368 GCALC_DBUG_RETURN(arrange_event(do_sorting, n_intersections));
1369}
1370
1371
1372static int node_on_right(const Gcalc_heap::Info *node,
1373 const Gcalc_heap::Info *edge_a, const Gcalc_heap::Info *edge_b)
1374{
1375 Gcalc_coord1 a_x, a_y;
1376 Gcalc_coord1 b_x, b_y;
1377 Gcalc_coord2 ax_by, ay_bx;
1378 int result;
1379
1380 gcalc_sub_coord1(a_x, node->node.shape.ix, edge_a->node.shape.ix);
1381 gcalc_sub_coord1(a_y, node->node.shape.iy, edge_a->node.shape.iy);
1382 gcalc_sub_coord1(b_x, edge_b->node.shape.ix, edge_a->node.shape.ix);
1383 gcalc_sub_coord1(b_y, edge_b->node.shape.iy, edge_a->node.shape.iy);
1384 gcalc_mul_coord1(ax_by, a_x, b_y);
1385 gcalc_mul_coord1(ay_bx, a_y, b_x);
1386 result= gcalc_cmp_coord(ax_by, ay_bx, GCALC_COORD_BASE2);
1387#ifdef GCALC_CHECK_WITH_FLOAT
1388 {
1389 long double dx= gcalc_get_double(edge_b->node.shape.ix, GCALC_COORD_BASE) -
1390 gcalc_get_double(edge_a->node.shape.ix, GCALC_COORD_BASE);
1391 long double dy= gcalc_get_double(edge_b->node.shape.iy, GCALC_COORD_BASE) -
1392 gcalc_get_double(edge_a->node.shape.iy, GCALC_COORD_BASE);
1393 long double ax= gcalc_get_double(node->node.shape.ix, GCALC_COORD_BASE) -
1394 gcalc_get_double(edge_a->node.shape.ix, GCALC_COORD_BASE);
1395 long double ay= gcalc_get_double(node->node.shape.iy, GCALC_COORD_BASE) -
1396 gcalc_get_double(edge_a->node.shape.iy, GCALC_COORD_BASE);
1397 long double d= ax * dy - ay * dx;
1398 if (result == 0)
1399 GCALC_DBUG_ASSERT(de_check(d, 0.0));
1400 else if (result < 0)
1401 GCALC_DBUG_ASSERT(de_check(d, 0.0) || d < 0);
1402 else
1403 GCALC_DBUG_ASSERT(de_check(d, 0.0) || d > 0);
1404 }
1405#endif /*GCALC_CHECK_WITH_FLOAT*/
1406 return result;
1407}
1408
1409
1410static int cmp_tops(const Gcalc_heap::Info *top_node,
1411 const Gcalc_heap::Info *edge_a, const Gcalc_heap::Info *edge_b)
1412{
1413 int cmp_res_a, cmp_res_b;
1414
1415 cmp_res_a= gcalc_cmp_coord1(edge_a->node.shape.ix, top_node->node.shape.ix);
1416 cmp_res_b= gcalc_cmp_coord1(edge_b->node.shape.ix, top_node->node.shape.ix);
1417
1418 if (cmp_res_a <= 0 && cmp_res_b > 0)
1419 return -1;
1420 if (cmp_res_b <= 0 && cmp_res_a > 0)
1421 return 1;
1422 if (cmp_res_a == 0 && cmp_res_b == 0)
1423 return 0;
1424
1425 return node_on_right(edge_a, top_node, edge_b);
1426}
1427
1428
1429int Gcalc_scan_iterator::insert_top_node()
1430{
1431 point *sp= state.slice;
1432 point **prev_hook= &state.slice;
1433 point *sp1= NULL;
1434 point *sp0= new_slice_point();
1435 int cmp_res;
1436
1437 GCALC_DBUG_ENTER("Gcalc_scan_iterator::insert_top_node");
1438 if (!sp0)
1439 GCALC_DBUG_RETURN(1);
1440 sp0->pi= m_cur_pi;
1441 sp0->next_pi= m_cur_pi->node.shape.left;
1442#ifndef GCALC_DBUG_OFF
1443 sp0->thread= m_cur_thread++;
1444#endif /*GCALC_DBUG_OFF*/
1445 if (m_cur_pi->node.shape.left)
1446 {
1447 calc_dx_dy(sp0);
1448 if (m_cur_pi->node.shape.right)
1449 {
1450 if (!(sp1= new_slice_point()))
1451 GCALC_DBUG_RETURN(1);
1452 sp1->event= sp0->event= scev_two_threads;
1453 sp1->pi= m_cur_pi;
1454 sp1->next_pi= m_cur_pi->node.shape.right;
1455#ifndef GCALC_DBUG_OFF
1456 sp1->thread= m_cur_thread++;
1457#endif /*GCALC_DBUG_OFF*/
1458 calc_dx_dy(sp1);
1459 /* We have two threads so should decide which one will be first */
1460 cmp_res= cmp_tops(m_cur_pi, m_cur_pi->node.shape.left, m_cur_pi->node.shape.right);
1461 if (cmp_res > 0)
1462 {
1463 point *tmp= sp0;
1464 sp0= sp1;
1465 sp1= tmp;
1466 }
1467 else if (cmp_res == 0)
1468 {
1469 /* Exactly same direction of the edges. */
1470 cmp_res= gcalc_cmp_coord1(m_cur_pi->node.shape.left->node.shape.iy, m_cur_pi->node.shape.right->node.shape.iy);
1471 if (cmp_res != 0)
1472 {
1473 if (cmp_res < 0)
1474 {
1475 if (add_eq_node(sp0->next_pi, sp1))
1476 GCALC_DBUG_RETURN(1);
1477 }
1478 else
1479 {
1480 if (add_eq_node(sp1->next_pi, sp0))
1481 GCALC_DBUG_RETURN(1);
1482 }
1483 }
1484 else
1485 {
1486 cmp_res= gcalc_cmp_coord1(m_cur_pi->node.shape.left->node.shape.ix, m_cur_pi->node.shape.right->node.shape.ix);
1487 if (cmp_res != 0)
1488 {
1489 if (cmp_res < 0)
1490 {
1491 if (add_eq_node(sp0->next_pi, sp1))
1492 GCALC_DBUG_RETURN(1);
1493 }
1494 else
1495 {
1496 if (add_eq_node(sp1->next_pi, sp0))
1497 GCALC_DBUG_RETURN(1);
1498 }
1499 }
1500 }
1501 }
1502 }
1503 else
1504 sp0->event= scev_thread;
1505 }
1506 else
1507 sp0->event= scev_single_point;
1508
1509
1510 /* Check if we already have an event - then we'll place the node there */
1511 for (; sp && !sp->event; prev_hook= sp->next_ptr(), sp=sp->get_next())
1512 {}
1513 if (!sp)
1514 {
1515 sp= state.slice;
1516 prev_hook= &state.slice;
1517 /* We need to find the place to insert. */
1518 for (; sp; prev_hook= sp->next_ptr(), sp=sp->get_next())
1519 {
1520 if (sp->event || gcalc_cmp_coord1(*sp->r_border, m_cur_pi->node.shape.ix) < 0)
1521 continue;
1522 cmp_res= node_on_right(m_cur_pi, sp->pi, sp->next_pi);
1523 if (cmp_res == 0)
1524 {
1525 /* The top node lies on the edge. */
1526 /* Nodes of that edge will be handled in other places. */
1527 sp->event= scev_intersection;
1528 }
1529 else if (cmp_res < 0)
1530 break;
1531 }
1532 }
1533
1534 if (sp0->event == scev_single_point)
1535 {
1536 /* Add single point to the bottom list. */
1537 *m_bottom_hook= sp0;
1538 m_bottom_hook= sp0->next_ptr();
1539 state.event_position_hook= prev_hook;
1540 }
1541 else
1542 {
1543 *prev_hook= sp0;
1544 sp0->next= sp;
1545 if (add_events_for_node(sp0))
1546 GCALC_DBUG_RETURN(1);
1547
1548 if (sp0->event == scev_two_threads)
1549 {
1550 *prev_hook= sp1;
1551 sp1->next= sp;
1552 if (add_events_for_node(sp1))
1553 GCALC_DBUG_RETURN(1);
1554
1555 sp0->next= sp1;
1556 *prev_hook= sp0;
1557 }
1558 }
1559
1560 GCALC_DBUG_RETURN(0);
1561}
1562
1563
1564void Gcalc_scan_iterator::remove_bottom_node()
1565{
1566 point *sp= state.slice;
1567 point **sp_hook= &state.slice;
1568 point *first_bottom_point= NULL;
1569
1570 GCALC_DBUG_ENTER("Gcalc_scan_iterator::remove_bottom_node");
1571 for (; sp; sp= sp->get_next())
1572 {
1573 if (sp->next_pi == m_cur_pi)
1574 {
1575 *sp_hook= sp->get_next();
1576 sp->pi= m_cur_pi;
1577 sp->next_pi= NULL;
1578 if (first_bottom_point)
1579 {
1580 first_bottom_point->event= sp->event= scev_two_ends;
1581 break;
1582 }
1583 first_bottom_point= sp;
1584 sp->event= scev_end;
1585 state.event_position_hook= sp_hook;
1586 }
1587 else
1588 sp_hook= sp->next_ptr();
1589 }
1590 GCALC_DBUG_ASSERT(first_bottom_point);
1591 *m_bottom_hook= first_bottom_point;
1592 m_bottom_hook= first_bottom_point->next_ptr();
1593 if (sp)
1594 {
1595 *m_bottom_hook= sp;
1596 m_bottom_hook= sp->next_ptr();
1597 }
1598
1599 GCALC_DBUG_VOID_RETURN;
1600}
1601
1602
1603int Gcalc_scan_iterator::add_events_for_node(point *sp_node)
1604{
1605 point *sp= state.slice;
1606 int cur_pi_r, sp_pi_r;
1607
1608 GCALC_DBUG_ENTER("Gcalc_scan_iterator::add_events_for_node");
1609
1610 /* Scan to the event point. */
1611 for (; sp != sp_node; sp= sp->get_next())
1612 {
1613 GCALC_DBUG_ASSERT(!sp->is_bottom());
1614 GCALC_DBUG_PRINT(("left cut_edge %d", sp->thread));
1615 if (sp->next_pi == sp_node->next_pi ||
1616 gcalc_cmp_coord1(*sp->r_border, *sp_node->l_border) < 0)
1617 continue;
1618 sp_pi_r= node_on_right(sp->next_pi, sp_node->pi, sp_node->next_pi);
1619 if (sp_pi_r < 0)
1620 continue;
1621 cur_pi_r= node_on_right(sp_node->next_pi, sp->pi, sp->next_pi);
1622 if (cur_pi_r > 0)
1623 continue;
1624 if (cur_pi_r == 0 && sp_pi_r == 0)
1625 {
1626 int cmp_res= cmp_point_info(sp->next_pi, sp_node->next_pi);
1627 if (cmp_res > 0)
1628 {
1629 if (add_eq_node(sp_node->next_pi, sp))
1630 GCALC_DBUG_RETURN(1);
1631 }
1632 else if (cmp_res < 0)
1633 {
1634 if (add_eq_node(sp->next_pi, sp_node))
1635 GCALC_DBUG_RETURN(1);
1636 }
1637 continue;
1638 }
1639
1640 if (cur_pi_r == 0)
1641 {
1642 if (add_eq_node(sp_node->next_pi, sp))
1643 GCALC_DBUG_RETURN(1);
1644 continue;
1645 }
1646 else if (sp_pi_r == 0)
1647 {
1648 if (add_eq_node(sp->next_pi, sp_node))
1649 GCALC_DBUG_RETURN(1);
1650 continue;
1651 }
1652
1653 if (sp->event)
1654 {
1655#ifndef GCALC_DBUG_OFF
1656 cur_pi_r= node_on_right(sp_node->pi, sp->pi, sp->next_pi);
1657 GCALC_DBUG_ASSERT(cur_pi_r == 0);
1658#endif /*GCALC_DBUG_OFF*/
1659 continue;
1660 }
1661 cur_pi_r= node_on_right(sp_node->pi, sp->pi, sp->next_pi);
1662 GCALC_DBUG_ASSERT(cur_pi_r >= 0);
1663 //GCALC_DBUG_ASSERT(cur_pi_r > 0); /* Is it ever violated? */
1664 if (cur_pi_r > 0 && add_intersection(sp, sp_node, m_cur_pi))
1665 GCALC_DBUG_RETURN(1);
1666 }
1667
1668 /* Scan to the end of the slice */
1669 sp= sp->get_next();
1670
1671 for (; sp; sp= sp->get_next())
1672 {
1673 GCALC_DBUG_ASSERT(!sp->is_bottom());
1674 GCALC_DBUG_PRINT(("right cut_edge %d", sp->thread));
1675 if (sp->next_pi == sp_node->next_pi ||
1676 gcalc_cmp_coord1(*sp_node->r_border, *sp->l_border) < 0)
1677 continue;
1678 sp_pi_r= node_on_right(sp->next_pi, sp_node->pi, sp_node->next_pi);
1679 if (sp_pi_r > 0)
1680 continue;
1681 cur_pi_r= node_on_right(sp_node->next_pi, sp->pi, sp->next_pi);
1682 if (cur_pi_r < 0)
1683 continue;
1684 if (cur_pi_r == 0 && sp_pi_r == 0)
1685 {
1686 int cmp_res= cmp_point_info(sp->next_pi, sp_node->next_pi);
1687 if (cmp_res > 0)
1688 {
1689 if (add_eq_node(sp_node->next_pi, sp))
1690 GCALC_DBUG_RETURN(1);
1691 }
1692 else if (cmp_res < 0)
1693 {
1694 if (add_eq_node(sp->next_pi, sp_node))
1695 GCALC_DBUG_RETURN(1);
1696 }
1697 continue;
1698 }
1699 if (cur_pi_r == 0)
1700 {
1701 if (add_eq_node(sp_node->next_pi, sp))
1702 GCALC_DBUG_RETURN(1);
1703 continue;
1704 }
1705 else if (sp_pi_r == 0)
1706 {
1707 if (add_eq_node(sp->next_pi, sp_node))
1708 GCALC_DBUG_RETURN(1);
1709 continue;
1710 }
1711
1712 if (sp->event)
1713 {
1714#ifndef GCALC_DBUG_OFF
1715 cur_pi_r= node_on_right(sp_node->pi, sp->pi, sp->next_pi);
1716 GCALC_DBUG_ASSERT(cur_pi_r == 0);
1717#endif /*GCALC_DBUG_OFF*/
1718 continue;
1719 }
1720 cur_pi_r= node_on_right(sp_node->pi, sp->pi, sp->next_pi);
1721 GCALC_DBUG_ASSERT(cur_pi_r <= 0);
1722 //GCALC_DBUG_ASSERT(cur_pi_r < 0); /* Is it ever violated? */
1723 if (cur_pi_r < 0 && add_intersection(sp_node, sp, m_cur_pi))
1724 GCALC_DBUG_RETURN(1);
1725 }
1726
1727 GCALC_DBUG_RETURN(0);
1728}
1729
1730
1731int Gcalc_scan_iterator::node_scan()
1732{
1733 point *sp= state.slice;
1734 Gcalc_heap::Info *cur_pi= m_cur_pi;
1735
1736 GCALC_DBUG_ENTER("Gcalc_scan_iterator::node_scan");
1737
1738 /* Scan to the event point. */
1739 /* Can be avoided if we add link to the sp to the Info. */
1740 for (; sp->next_pi != cur_pi; sp= sp->get_next())
1741 {}
1742
1743 GCALC_DBUG_PRINT(("node for %d", sp->thread));
1744 /* Handle the point itself. */
1745 sp->pi= cur_pi;
1746 sp->next_pi= cur_pi->node.shape.left;
1747 sp->event= scev_point;
1748 calc_dx_dy(sp);
1749
1750 GCALC_DBUG_RETURN(add_events_for_node(sp));
1751}
1752
1753
1754void Gcalc_scan_iterator::eq_scan()
1755{
1756 point *sp= eq_sp(m_cur_pi);
1757 GCALC_DBUG_ENTER("Gcalc_scan_iterator::eq_scan");
1758
1759#ifndef GCALC_DBUG_OFF
1760 {
1761 point *cur_p= state.slice;
1762 for (; cur_p && cur_p != sp; cur_p= cur_p->get_next())
1763 {}
1764 GCALC_DBUG_ASSERT(cur_p);
1765 }
1766#endif /*GCALC_DBUG_OFF*/
1767 if (!sp->event)
1768 {
1769 sp->event= scev_intersection;
1770 sp->ev_pi= m_cur_pi;
1771 }
1772
1773 GCALC_DBUG_VOID_RETURN;
1774}
1775
1776
1777void Gcalc_scan_iterator::intersection_scan()
1778{
1779 intersection_info *ii= i_data(m_cur_pi);
1780 GCALC_DBUG_ENTER("Gcalc_scan_iterator::intersection_scan");
1781
1782#ifndef GCALC_DBUG_OFF
1783 {
1784 point *sp= state.slice;
1785 for (; sp && sp != ii->edge_a; sp= sp->get_next())
1786 {}
1787 GCALC_DBUG_ASSERT(sp);
1788 for (; sp && sp != ii->edge_b; sp= sp->get_next())
1789 {}
1790 GCALC_DBUG_ASSERT(sp);
1791 }
1792#endif /*GCALC_DBUG_OFF*/
1793
1794 ii->edge_a->event= ii->edge_b->event= scev_intersection;
1795 ii->edge_a->ev_pi= ii->edge_b->ev_pi= m_cur_pi;
1796 free_item(ii);
1797 m_cur_pi->node.intersection.data= NULL;
1798
1799 GCALC_DBUG_VOID_RETURN;
1800}
1801
1802
1803int Gcalc_scan_iterator::add_intersection(point *sp_a, point *sp_b,
1804 Gcalc_heap::Info *pi_from)
1805{
1806 Gcalc_heap::Info *ii;
1807 intersection_info *i_calc;
1808 int cmp_res;
1809 int skip_next= 0;
1810
1811 GCALC_DBUG_ENTER("Gcalc_scan_iterator::add_intersection");
1812 if (!(i_calc= new_intersection_info(sp_a, sp_b)) ||
1813 !(ii= new_intersection(m_heap, i_calc)))
1814 GCALC_DBUG_RETURN(1);
1815
1816 ii->node.intersection.equal= 0;
1817
1818 for (;
1819 pi_from->get_next() != sp_a->next_pi &&
1820 pi_from->get_next() != sp_b->next_pi;
1821 pi_from= pi_from->get_next())
1822 {
1823 Gcalc_heap::Info *cur= pi_from->get_next();
1824 if (skip_next)
1825 {
1826 if (cur->type == Gcalc_heap::nt_intersection)
1827 skip_next= cur->node.intersection.equal;
1828 else
1829 skip_next= 0;
1830 continue;
1831 }
1832 if (cur->type == Gcalc_heap::nt_intersection)
1833 {
1834 cmp_res= cmp_intersections(cur, ii);
1835 skip_next= cur->node.intersection.equal;
1836 }
1837 else if (cur->type == Gcalc_heap::nt_eq_node)
1838 continue;
1839 else
1840 cmp_res= cmp_node_isc(cur, ii);
1841 if (cmp_res == 0)
1842 {
1843 ii->node.intersection.equal= 1;
1844 break;
1845 }
1846 else if (cmp_res > 0)
1847 break;
1848 }
1849
1850 /* Intersection inserted before the equal point. */
1851 ii->next= pi_from->get_next();
1852 pi_from->next= ii;
1853
1854 GCALC_DBUG_RETURN(0);
1855}
1856
1857
1858int Gcalc_scan_iterator::add_eq_node(Gcalc_heap::Info *node, point *sp)
1859{
1860 Gcalc_heap::Info *en;
1861
1862 GCALC_DBUG_ENTER("Gcalc_scan_iterator::add_intersection");
1863 en= new_eq_point(m_heap, node, sp);
1864 if (!en)
1865 GCALC_DBUG_RETURN(1);
1866
1867 /* eq_node iserted after teh equal point. */
1868 en->next= node->get_next();
1869 node->next= en;
1870
1871 GCALC_DBUG_RETURN(0);
1872}
1873
1874
1875void calc_t(Gcalc_coord2 t_a, Gcalc_coord2 t_b,
1876 Gcalc_coord1 dxa, Gcalc_coord1 dxb,
1877 const Gcalc_heap::Info *p1, const Gcalc_heap::Info *p2,
1878 const Gcalc_heap::Info *p3, const Gcalc_heap::Info *p4)
1879{
1880 Gcalc_coord1 a2_a1x, a2_a1y;
1881 Gcalc_coord2 x1y2, x2y1;
1882 Gcalc_coord1 dya, dyb;
1883
1884 gcalc_sub_coord1(a2_a1x, p3->node.shape.ix, p1->node.shape.ix);
1885 gcalc_sub_coord1(a2_a1y, p3->node.shape.iy, p1->node.shape.iy);
1886
1887 gcalc_sub_coord1(dxa, p2->node.shape.ix, p1->node.shape.ix);
1888 gcalc_sub_coord1(dya, p2->node.shape.iy, p1->node.shape.iy);
1889 gcalc_sub_coord1(dxb, p4->node.shape.ix, p3->node.shape.ix);
1890 gcalc_sub_coord1(dyb, p4->node.shape.iy, p3->node.shape.iy);
1891
1892 gcalc_mul_coord1(x1y2, dxa, dyb);
1893 gcalc_mul_coord1(x2y1, dya, dxb);
1894 gcalc_sub_coord(t_b, GCALC_COORD_BASE2, x1y2, x2y1);
1895
1896
1897 gcalc_mul_coord1(x1y2, a2_a1x, dyb);
1898 gcalc_mul_coord1(x2y1, a2_a1y, dxb);
1899 gcalc_sub_coord(t_a, GCALC_COORD_BASE2, x1y2, x2y1);
1900}
1901
1902
1903double Gcalc_scan_iterator::get_y() const
1904{
1905 if (state.pi->type == Gcalc_heap::nt_intersection)
1906 {
1907 Gcalc_coord1 dxa, dya;
1908 Gcalc_coord2 t_a, t_b;
1909 Gcalc_coord3 a_tb, b_ta, y_exp;
1910 calc_t(t_a, t_b, dxa, dya,
1911 state.pi->node.intersection.p1, state.pi->node.intersection.p2, state.pi->node.intersection.p3, state.pi->node.intersection.p4);
1912
1913
1914 gcalc_mul_coord(a_tb, GCALC_COORD_BASE3,
1915 t_b, GCALC_COORD_BASE2, state.pi->node.intersection.p1->node.shape.iy, GCALC_COORD_BASE);
1916 gcalc_mul_coord(b_ta, GCALC_COORD_BASE3,
1917 t_a, GCALC_COORD_BASE2, dya, GCALC_COORD_BASE);
1918
1919 gcalc_add_coord(y_exp, GCALC_COORD_BASE3, a_tb, b_ta);
1920
1921 return (get_pure_double(y_exp, GCALC_COORD_BASE3) /
1922 get_pure_double(t_b, GCALC_COORD_BASE2)) / m_heap->coord_extent;
1923 }
1924 else
1925 return state.pi->node.shape.y;
1926}
1927
1928
1929double Gcalc_scan_iterator::get_event_x() const
1930{
1931 if (state.pi->type == Gcalc_heap::nt_intersection)
1932 {
1933 Gcalc_coord1 dxa, dya;
1934 Gcalc_coord2 t_a, t_b;
1935 Gcalc_coord3 a_tb, b_ta, x_exp;
1936 calc_t(t_a, t_b, dxa, dya,
1937 state.pi->node.intersection.p1, state.pi->node.intersection.p2, state.pi->node.intersection.p3, state.pi->node.intersection.p4);
1938
1939
1940 gcalc_mul_coord(a_tb, GCALC_COORD_BASE3,
1941 t_b, GCALC_COORD_BASE2, state.pi->node.intersection.p1->node.shape.ix, GCALC_COORD_BASE);
1942 gcalc_mul_coord(b_ta, GCALC_COORD_BASE3,
1943 t_a, GCALC_COORD_BASE2, dxa, GCALC_COORD_BASE);
1944
1945 gcalc_add_coord(x_exp, GCALC_COORD_BASE3, a_tb, b_ta);
1946
1947 return (get_pure_double(x_exp, GCALC_COORD_BASE3) /
1948 get_pure_double(t_b, GCALC_COORD_BASE2)) / m_heap->coord_extent;
1949 }
1950 else
1951 return state.pi->node.shape.x;
1952}
1953
1954double Gcalc_scan_iterator::get_h() const
1955{
1956 double cur_y= get_y();
1957 double next_y;
1958 if (state.pi->type == Gcalc_heap::nt_intersection)
1959 {
1960 double x;
1961 state.pi->calc_xy(&x, &next_y);
1962 }
1963 else
1964 next_y= state.pi->next ? state.pi->get_next()->node.shape.y : 0.0;
1965 return next_y - cur_y;
1966}
1967
1968
1969double Gcalc_scan_iterator::get_sp_x(const point *sp) const
1970{
1971 double dy;
1972 if (sp->event & (scev_end | scev_two_ends | scev_point))
1973 return sp->pi->node.shape.x;
1974 dy= sp->next_pi->node.shape.y - sp->pi->node.shape.y;
1975 if (fabs(dy) < 1e-12)
1976 return sp->pi->node.shape.x;
1977 return sp->pi->node.shape.x + (sp->next_pi->node.shape.x - sp->pi->node.shape.x) * dy;
1978}
1979
1980
1981double Gcalc_scan_iterator::get_pure_double(const Gcalc_internal_coord *d,
1982 int d_len)
1983{
1984 int n= 1;
1985 long double res= (long double) FIRST_DIGIT(d[0]);
1986 do
1987 {
1988 res*= (long double) GCALC_DIG_BASE;
1989 res+= (long double) d[n];
1990 } while(++n < d_len);
1991
1992 if (GCALC_SIGN(d[0]))
1993 res*= -1.0;
1994 return res;
1995}
1996
1997
1998#endif /* HAVE_SPATIAL */
1999