diff options
Diffstat (limited to '')
| -rw-r--r-- | sql/gcalc_slicescan.cc | 2011 |
1 files changed, 2011 insertions, 0 deletions
diff --git a/sql/gcalc_slicescan.cc b/sql/gcalc_slicescan.cc new file mode 100644 index 00000000..93561f5f --- /dev/null +++ b/sql/gcalc_slicescan.cc @@ -0,0 +1,2011 @@ +/* Copyright (c) 2000, 2010 Oracle and/or its affiliates. All rights reserved. + Copyright (C) 2011 Monty Program Ab. + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; version 2 of the License. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program; if not, write to the Free Software + Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1335 USA */ + + +#include "mariadb.h" +#include <my_sys.h> +#include <m_string.h> + +#ifdef HAVE_SPATIAL + +#include "gcalc_slicescan.h" + + +#define PH_DATA_OFFSET 8 +#define coord_to_float(d) ((double) d) +#define coord_eq(a, b) (a == b) + +typedef int (*sc_compare_func)(const void*, const void*); + +#define LS_LIST_ITEM Gcalc_dyn_list::Item +#define LS_COMPARE_FUNC_DECL sc_compare_func compare, +#define LS_COMPARE_FUNC_CALL(list_el1, list_el2) (*compare)(list_el1, list_el2) +#define LS_NEXT(A) (A)->next +#define LS_SET_NEXT(A,val) (A)->next= val +#define LS_P_NEXT(A) &(A)->next +#define LS_NAME sort_list +#define LS_SCOPE static +#define LS_STRUCT_NAME sort_list_stack_struct +#include "plistsort.c" + + +#define GCALC_COORD_MINUS 0x80000000 +#define FIRST_DIGIT(d) ((d) & 0x7FFFFFFF) +#define GCALC_SIGN(d) ((d) & 0x80000000) + +static Gcalc_scan_iterator::point *eq_sp(const Gcalc_heap::Info *pi) +{ + GCALC_DBUG_ASSERT(pi->type == Gcalc_heap::nt_eq_node); + return (Gcalc_scan_iterator::point *) pi->node.eq.data; +} + + +static Gcalc_scan_iterator::intersection_info *i_data(const Gcalc_heap::Info *pi) +{ + GCALC_DBUG_ASSERT(pi->type == Gcalc_heap::nt_intersection); + return (Gcalc_scan_iterator::intersection_info *) pi->node.intersection.data; +} + + +#ifndef GCALC_DBUG_OFF + +int gcalc_step_counter= 0; + +void GCALC_DBUG_CHECK_COUNTER() +{ + if (++gcalc_step_counter == 0) + GCALC_DBUG_PRINT(("step_counter_0")); + else + GCALC_DBUG_PRINT(("%d step_counter", gcalc_step_counter)); +} + + +const char *gcalc_ev_name(int ev) +{ + switch (ev) + { + case scev_none: + return "n"; + case scev_thread: + return "t"; + case scev_two_threads: + return "tt"; + case scev_end: + return "e"; + case scev_two_ends: + return "ee"; + case scev_intersection: + return "i"; + case scev_point: + return "p"; + case scev_single_point: + return "sp"; + default:; + }; + GCALC_DBUG_ASSERT(0); + return "unk"; +} + + +static int gcalc_pi_str(char *str, const Gcalc_heap::Info *pi, const char *postfix) +{ + return sprintf(str, "%s %d %d | %s %d %d%s", + GCALC_SIGN(pi->node.shape.ix[0]) ? "-":"", FIRST_DIGIT(pi->node.shape.ix[0]),pi->node.shape.ix[1], + GCALC_SIGN(pi->node.shape.iy[0]) ? "-":"", FIRST_DIGIT(pi->node.shape.iy[0]),pi->node.shape.iy[1], + postfix); + +} + + +static void GCALC_DBUG_PRINT_PI(const Gcalc_heap::Info *pi) +{ + char buf[128]; + int n_buf; + if (pi->type == Gcalc_heap::nt_intersection) + { + const Gcalc_scan_iterator::intersection_info *ic= i_data(pi); + + GCALC_DBUG_PRINT(("intersection point %d %d", + ic->edge_a->thread, ic->edge_b->thread)); + return; + } + if (pi->type == Gcalc_heap::nt_eq_node) + { + const Gcalc_scan_iterator::point *e= eq_sp(pi); + GCALC_DBUG_PRINT(("eq point %d", e->thread)); + return; + } + n_buf= gcalc_pi_str(buf, pi, ""); + buf[n_buf]= 0; + GCALC_DBUG_PRINT(("%s", buf)); +} + + +static void GCALC_DBUG_PRINT_SLICE(const char *header, + const Gcalc_scan_iterator::point *slice) +{ + size_t nbuf; + char buf[1024]; + nbuf= strlen(header); + strcpy(buf, header); + for (; slice; slice= slice->get_next()) + { + size_t lnbuf= nbuf; + lnbuf+= sprintf(buf + lnbuf, "%d\t", slice->thread); + lnbuf+= sprintf(buf + lnbuf, "%s\t", gcalc_ev_name(slice->event)); + + lnbuf+= gcalc_pi_str(buf + lnbuf, slice->pi, "\t"); + if (slice->is_bottom()) + lnbuf+= sprintf(buf+lnbuf, "bt\t"); + else + lnbuf+= gcalc_pi_str(buf+lnbuf, slice->next_pi, "\t"); + buf[lnbuf]= 0; + GCALC_DBUG_PRINT(("%s", buf)); + } +} + + +#else +#define GCALC_DBUG_CHECK_COUNTER() do { } while(0) +#define GCALC_DBUG_PRINT_PI(pi) do { } while(0) +#define GCALC_DBUG_PRINT_SLICE(a, b) do { } while(0) +#define GCALC_DBUG_PRINT_INTERSECTIONS(a) do { } while(0) +#define GCALC_DBUG_PRINT_STATE(a) do { } while(0) +#endif /*GCALC_DBUG_OFF*/ + + +Gcalc_dyn_list::Gcalc_dyn_list(size_t blk_size, size_t sizeof_item): + m_blk_size(blk_size - ALLOC_ROOT_MIN_BLOCK_SIZE), + m_sizeof_item(ALIGN_SIZE(sizeof_item)), + m_points_per_blk((uint)((m_blk_size - PH_DATA_OFFSET) / m_sizeof_item)), + m_blk_hook(&m_first_blk), + m_free(NULL), + m_keep(NULL) +{} + + +Gcalc_dyn_list::Gcalc_dyn_list(const Gcalc_dyn_list &dl) +{ + m_blk_size= dl.m_blk_size; + m_sizeof_item= dl.m_sizeof_item; + m_points_per_blk= dl.m_points_per_blk; + m_blk_hook= &m_first_blk; + m_free= NULL; + m_keep= NULL; +} + + +void Gcalc_dyn_list::format_blk(void* block) +{ + Item *pi_end, *cur_pi, *first_pi; + GCALC_DBUG_ASSERT(m_free == NULL); + first_pi= cur_pi= (Item *)(((char *)block) + PH_DATA_OFFSET); + pi_end= ptr_add(first_pi, m_points_per_blk - 1); + do { + cur_pi= cur_pi->next= ptr_add(cur_pi, 1); + } while (cur_pi<pi_end); + cur_pi->next= m_free; + m_free= first_pi; +} + + +Gcalc_dyn_list::Item *Gcalc_dyn_list::alloc_new_blk() +{ + void *new_block= my_malloc(PSI_INSTRUMENT_ME, m_blk_size, MYF(MY_WME)); + if (!new_block) + return NULL; + *m_blk_hook= new_block; + m_blk_hook= (void**)new_block; + format_blk(new_block); + return new_item(); +} + + +static void free_blk_list(void *list) +{ + void *next_blk; + while (list) + { + next_blk= *((void **)list); + my_free(list); + list= next_blk; + } +} + + +void Gcalc_dyn_list::cleanup() +{ + *m_blk_hook= NULL; + free_blk_list(m_first_blk); + m_first_blk= NULL; + m_blk_hook= &m_first_blk; + m_free= NULL; +} + + +Gcalc_dyn_list::~Gcalc_dyn_list() +{ + cleanup(); +} + + +void Gcalc_dyn_list::reset() +{ + *m_blk_hook= NULL; + if (m_first_blk) + { + free_blk_list(*((void **)m_first_blk)); + m_blk_hook= (void**)m_first_blk; + m_free= NULL; + format_blk(m_first_blk); + } +} + + +/* Internal coordinate operations implementations */ + +void gcalc_set_zero(Gcalc_internal_coord *d, int d_len) +{ + do + { + d[--d_len]= 0; + } while (d_len); +} + + +int gcalc_is_zero(const Gcalc_internal_coord *d, int d_len) +{ + do + { + if (d[--d_len] != 0) + return 0; + } while (d_len); + return 1; +} + + +#ifdef GCALC_CHECK_WITH_FLOAT +static double *gcalc_coord_extent= NULL; + +long double gcalc_get_double(const Gcalc_internal_coord *d, int d_len) +{ + int n= 1; + long double res= (long double) FIRST_DIGIT(d[0]); + do + { + res*= (long double) GCALC_DIG_BASE; + res+= (long double) d[n]; + } while(++n < d_len); + + n= 0; + do + { + if ((n & 1) && gcalc_coord_extent) + res/= *gcalc_coord_extent; + } while(++n < d_len); + + if (GCALC_SIGN(d[0])) + res*= -1.0; + return res; +} +#endif /*GCALC_CHECK_WITH_FLOAT*/ + + +static void do_add(Gcalc_internal_coord *result, int result_len, + const Gcalc_internal_coord *a, + const Gcalc_internal_coord *b) +{ + int n_digit= result_len-1; + gcalc_digit_t carry= 0; + + do + { + if ((result[n_digit]= + a[n_digit] + b[n_digit] + carry) >= GCALC_DIG_BASE) + { + carry= 1; + result[n_digit]-= GCALC_DIG_BASE; + } + else + carry= 0; + } while (--n_digit); + + result[0]= (a[0] + FIRST_DIGIT(b[0]) + carry); + + GCALC_DBUG_ASSERT(FIRST_DIGIT(result[0]) < GCALC_DIG_BASE); +} + + +static void do_sub(Gcalc_internal_coord *result, int result_len, + const Gcalc_internal_coord *a, + const Gcalc_internal_coord *b) +{ + int n_digit= result_len-1; + gcalc_digit_t carry= 0; + gcalc_digit_t cur_b, cur_a; + + do + { + cur_b= b[n_digit] + carry; + cur_a= a[n_digit]; + if (cur_a < cur_b) + { + carry= 1; + result[n_digit]= (GCALC_DIG_BASE - cur_b) + cur_a; + } + else + { + carry= 0; + result[n_digit]= cur_a - cur_b; + } + } while (--n_digit); + + + result[0]= a[0] - FIRST_DIGIT(b[0]) - carry; + + GCALC_DBUG_ASSERT(FIRST_DIGIT(a[0]) >= FIRST_DIGIT(b[0]) + carry); + GCALC_DBUG_ASSERT(!gcalc_is_zero(result, result_len)); +} +/* +static void do_sub(Gcalc_internal_coord *result, int result_len, + const Gcalc_internal_coord *a, + const Gcalc_internal_coord *b) +{ + int n_digit= result_len-1; + gcalc_digit_t carry= 0; + + do + { + if ((result[n_digit]= a[n_digit] - b[n_digit] - carry) < 0) + { + carry= 1; + result[n_digit]+= GCALC_DIG_BASE; + } + else + carry= 0; + } while (--n_digit); + + + result[0]= a[0] - FIRST_DIGIT(b[0]) - carry; + + GCALC_DBUG_ASSERT(FIRST_DIGIT(a[0]) - FIRST_DIGIT(b[0]) - carry >= 0); + GCALC_DBUG_ASSERT(!gcalc_is_zero(result, result_len)); +} +*/ + +static int do_cmp(const Gcalc_internal_coord *a, + const Gcalc_internal_coord *b, int len) +{ + int n_digit= 1; + + if ((FIRST_DIGIT(a[0]) != FIRST_DIGIT(b[0]))) + return FIRST_DIGIT(a[0]) > FIRST_DIGIT(b[0]) ? 1 : -1; + + do + { + if ((a[n_digit] != b[n_digit])) + return a[n_digit] > b[n_digit] ? 1 : -1; + } while (++n_digit < len); + + return 0; +} + + +#ifdef GCALC_CHECK_WITH_FLOAT +static int de_weak_check(long double a, long double b, long double ex) +{ + long double d= a - b; + if (d < ex && d > -ex) + return 1; + + d/= fabsl(a) + fabsl(b); + if (d < ex && d > -ex) + return 1; + return 0; +} + +static int de_check(long double a, long double b) +{ + return de_weak_check(a, b, (long double) 1e-9); +} +#endif /*GCALC_CHECK_WITH_FLOAT*/ + + +void gcalc_mul_coord(Gcalc_internal_coord *result, int result_len, + const Gcalc_internal_coord *a, int a_len, + const Gcalc_internal_coord *b, int b_len) +{ + GCALC_DBUG_ASSERT(result_len == a_len + b_len); + GCALC_DBUG_ASSERT(a_len >= b_len); + int n_a, n_b, n_res; + gcalc_digit_t carry= 0; + + gcalc_set_zero(result, result_len); + + n_a= a_len - 1; + do + { + gcalc_coord2 cur_a= n_a ? a[n_a] : FIRST_DIGIT(a[0]); + n_b= b_len - 1; + do + { + gcalc_coord2 cur_b= n_b ? b[n_b] : FIRST_DIGIT(b[0]); + gcalc_coord2 mul= cur_a * cur_b + carry + result[n_a + n_b + 1]; + result[n_a + n_b + 1]= mul % GCALC_DIG_BASE; + carry= (gcalc_digit_t) (mul / (gcalc_coord2) GCALC_DIG_BASE); + } while (n_b--); + if (carry) + { + for (n_res= n_a; (result[n_res]+= carry) >= GCALC_DIG_BASE; + n_res--) + { + result[n_res]-= GCALC_DIG_BASE; + carry= 1; + } + carry= 0; + } + } while (n_a--); + if (!gcalc_is_zero(result, result_len)) + result[0]|= GCALC_SIGN(a[0] ^ b[0]); +#ifdef GCALC_CHECK_WITH_FLOAT + GCALC_DBUG_ASSERT(de_check(gcalc_get_double(a, a_len) * + gcalc_get_double(b, b_len), + gcalc_get_double(result, result_len))); +#endif /*GCALC_CHECK_WITH_FLOAT*/ +} + + +inline void gcalc_mul_coord1(Gcalc_coord1 result, + const Gcalc_coord1 a, const Gcalc_coord1 b) +{ + return gcalc_mul_coord(result, GCALC_COORD_BASE2, + a, GCALC_COORD_BASE, b, GCALC_COORD_BASE); +} + + +void gcalc_add_coord(Gcalc_internal_coord *result, int result_len, + const Gcalc_internal_coord *a, + const Gcalc_internal_coord *b) +{ + if (GCALC_SIGN(a[0]) == GCALC_SIGN(b[0])) + do_add(result, result_len, a, b); + else + { + int cmp_res= do_cmp(a, b, result_len); + if (cmp_res == 0) + gcalc_set_zero(result, result_len); + else if (cmp_res > 0) + do_sub(result, result_len, a, b); + else + do_sub(result, result_len, b, a); + } +#ifdef GCALC_CHECK_WITH_FLOAT + GCALC_DBUG_ASSERT(de_check(gcalc_get_double(a, result_len) + + gcalc_get_double(b, result_len), + gcalc_get_double(result, result_len))); +#endif /*GCALC_CHECK_WITH_FLOAT*/ +} + + +void gcalc_sub_coord(Gcalc_internal_coord *result, int result_len, + const Gcalc_internal_coord *a, + const Gcalc_internal_coord *b) +{ + if (GCALC_SIGN(a[0] ^ b[0])) + do_add(result, result_len, a, b); + else + { + int cmp_res= do_cmp(a, b, result_len); + if (cmp_res == 0) + gcalc_set_zero(result, result_len); + else if (cmp_res > 0) + do_sub(result, result_len, a, b); + else + { + do_sub(result, result_len, b, a); + result[0]^= GCALC_COORD_MINUS; + } + } +#ifdef GCALC_CHECK_WITH_FLOAT + GCALC_DBUG_ASSERT(de_check(gcalc_get_double(a, result_len) - + gcalc_get_double(b, result_len), + gcalc_get_double(result, result_len))); +#endif /*GCALC_CHECK_WITH_FLOAT*/ +} + + +inline void gcalc_sub_coord1(Gcalc_coord1 result, + const Gcalc_coord1 a, const Gcalc_coord1 b) +{ + return gcalc_sub_coord(result, GCALC_COORD_BASE, a, b); +} + + +int gcalc_cmp_coord(const Gcalc_internal_coord *a, + const Gcalc_internal_coord *b, int len) +{ + int n_digit= 0; + int result= 0; + + do + { + if (a[n_digit] == b[n_digit]) + { + n_digit++; + continue; + } + if (a[n_digit] > b[n_digit]) + result= GCALC_SIGN(a[0]) ? -1 : 1; + else + result= GCALC_SIGN(b[0]) ? 1 : -1; + break; + + } while (n_digit < len); + +#ifdef GCALC_CHECK_WITH_FLOAT + if (result == 0) + GCALC_DBUG_ASSERT(de_check(gcalc_get_double(a, len), + gcalc_get_double(b, len))); + else if (result == 1) + GCALC_DBUG_ASSERT(de_check(gcalc_get_double(a, len), + gcalc_get_double(b, len)) || + gcalc_get_double(a, len) > gcalc_get_double(b, len)); + else + GCALC_DBUG_ASSERT(de_check(gcalc_get_double(a, len), + gcalc_get_double(b, len)) || + gcalc_get_double(a, len) < gcalc_get_double(b, len)); +#endif /*GCALC_CHECK_WITH_FLOAT*/ + return result; +} + + +#define gcalc_cmp_coord1(a, b) gcalc_cmp_coord(a, b, GCALC_COORD_BASE) + +int gcalc_set_double(Gcalc_internal_coord *c, double d, double ext) +{ + int sign; + double ds= d * ext; + if ((sign= ds < 0)) + ds= -ds; + c[0]= (gcalc_digit_t) (ds / (double) GCALC_DIG_BASE); + c[1]= (gcalc_digit_t) (ds - ((double) c[0]) * (double) GCALC_DIG_BASE); + if (c[1] >= GCALC_DIG_BASE) + { + c[1]= 0; + c[0]++; + } + if (sign && (c[0] | c[1])) + c[0]|= GCALC_COORD_MINUS; +#ifdef GCALC_CHECK_WITH_FLOAT + GCALC_DBUG_ASSERT(de_check(d, gcalc_get_double(c, 2))); +#endif /*GCALC_CHECK_WITH_FLOAT*/ + return 0; +} + + +typedef gcalc_digit_t Gcalc_coord4[GCALC_COORD_BASE*4]; +typedef gcalc_digit_t Gcalc_coord5[GCALC_COORD_BASE*5]; + + +void Gcalc_scan_iterator::intersection_info::do_calc_t() +{ + Gcalc_coord1 a2_a1x, a2_a1y; + Gcalc_coord2 x1y2, x2y1; + + gcalc_sub_coord1(a2_a1x, edge_b->pi->node.shape.ix, edge_a->pi->node.shape.ix); + gcalc_sub_coord1(a2_a1y, edge_b->pi->node.shape.iy, edge_a->pi->node.shape.iy); + + GCALC_DBUG_ASSERT(!gcalc_is_zero(edge_a->dy, GCALC_COORD_BASE) || + !gcalc_is_zero(edge_b->dy, GCALC_COORD_BASE)); + + gcalc_mul_coord1(x1y2, edge_a->dx, edge_b->dy); + gcalc_mul_coord1(x2y1, edge_a->dy, edge_b->dx); + gcalc_sub_coord(t_b, GCALC_COORD_BASE2, x1y2, x2y1); + + + gcalc_mul_coord1(x1y2, a2_a1x, edge_b->dy); + gcalc_mul_coord1(x2y1, a2_a1y, edge_b->dx); + gcalc_sub_coord(t_a, GCALC_COORD_BASE2, x1y2, x2y1); + t_calculated= 1; +} + + +void Gcalc_scan_iterator::intersection_info::do_calc_y() +{ + GCALC_DBUG_ASSERT(t_calculated); + + Gcalc_coord3 a_tb, b_ta; + + gcalc_mul_coord(a_tb, GCALC_COORD_BASE3, + t_b, GCALC_COORD_BASE2, edge_a->pi->node.shape.iy, GCALC_COORD_BASE); + gcalc_mul_coord(b_ta, GCALC_COORD_BASE3, + t_a, GCALC_COORD_BASE2, edge_a->dy, GCALC_COORD_BASE); + + gcalc_add_coord(y_exp, GCALC_COORD_BASE3, a_tb, b_ta); + y_calculated= 1; +} + + +void Gcalc_scan_iterator::intersection_info::do_calc_x() +{ + GCALC_DBUG_ASSERT(t_calculated); + + Gcalc_coord3 a_tb, b_ta; + + gcalc_mul_coord(a_tb, GCALC_COORD_BASE3, + t_b, GCALC_COORD_BASE2, edge_a->pi->node.shape.ix, GCALC_COORD_BASE); + gcalc_mul_coord(b_ta, GCALC_COORD_BASE3, + t_a, GCALC_COORD_BASE2, edge_a->dx, GCALC_COORD_BASE); + + gcalc_add_coord(x_exp, GCALC_COORD_BASE3, a_tb, b_ta); + x_calculated= 1; +} + + +static int cmp_node_isc(const Gcalc_heap::Info *node, + const Gcalc_heap::Info *isc) +{ + GCALC_DBUG_ASSERT(node->type == Gcalc_heap::nt_shape_node); + Gcalc_scan_iterator::intersection_info *inf= i_data(isc); + Gcalc_coord3 exp; + int result; + + inf->calc_t(); + inf->calc_y_exp(); + + gcalc_mul_coord(exp, GCALC_COORD_BASE3, + inf->t_b, GCALC_COORD_BASE2, node->node.shape.iy, GCALC_COORD_BASE); + + result= gcalc_cmp_coord(exp, inf->y_exp, GCALC_COORD_BASE3); +#ifdef GCALC_CHECK_WITH_FLOAT + long double int_x, int_y; + isc->calc_xy_ld(&int_x, &int_y); + if (result < 0) + { + if (!de_check(int_y, node->node.shape.y) && node->node.shape.y > int_y) + GCALC_DBUG_PRINT(("floatcheck cmp_nod_iscy %g < %LG", node->node.shape.y, int_y)); + } + else if (result > 0) + { + if (!de_check(int_y, node->node.shape.y) && node->node.shape.y < int_y) + GCALC_DBUG_PRINT(("floatcheck cmp_nod_iscy %g > %LG", node->node.shape.y, int_y)); + } + else + { + if (!de_check(int_y, node->node.shape.y)) + GCALC_DBUG_PRINT(("floatcheck cmp_nod_iscy %g == %LG", node->node.shape.y, int_y)); + } +#endif /*GCALC_CHECK_WITH_FLOAT*/ + if (result) + goto exit; + + + inf->calc_x_exp(); + gcalc_mul_coord(exp, GCALC_COORD_BASE3, + inf->t_b, GCALC_COORD_BASE2, node->node.shape.ix, GCALC_COORD_BASE); + + result= gcalc_cmp_coord(exp, inf->x_exp, GCALC_COORD_BASE3); +#ifdef GCALC_CHECK_WITH_FLOAT + if (result < 0) + { + if (!de_check(int_x, node->node.shape.x) && node->node.shape.x > int_x) + GCALC_DBUG_PRINT(("floatcheck cmp_nod_iscx failed %g < %LG", + node->node.shape.x, int_x)); + } + else if (result > 0) + { + if (!de_check(int_x, node->node.shape.x) && node->node.shape.x < int_x) + GCALC_DBUG_PRINT(("floatcheck cmp_nod_iscx failed %g > %LG", + node->node.shape.x, int_x)); + } + else + { + if (!de_check(int_x, node->node.shape.x)) + GCALC_DBUG_PRINT(("floatcheck cmp_nod_iscx failed %g == %LG", + node->node.shape.x, int_x)); + } +#endif /*GCALC_CHECK_WITH_FLOAT*/ +exit: + return result; +} + + +static int cmp_intersections(const Gcalc_heap::Info *i1, + const Gcalc_heap::Info *i2) +{ + Gcalc_scan_iterator::intersection_info *inf1= i_data(i1); + Gcalc_scan_iterator::intersection_info *inf2= i_data(i2); + Gcalc_coord5 exp_a, exp_b; + int result; + + inf1->calc_t(); + inf2->calc_t(); + + inf1->calc_y_exp(); + inf2->calc_y_exp(); + + gcalc_mul_coord(exp_a, GCALC_COORD_BASE5, + inf1->y_exp, GCALC_COORD_BASE3, inf2->t_b, GCALC_COORD_BASE2); + gcalc_mul_coord(exp_b, GCALC_COORD_BASE5, + inf2->y_exp, GCALC_COORD_BASE3, inf1->t_b, GCALC_COORD_BASE2); + + result= gcalc_cmp_coord(exp_a, exp_b, GCALC_COORD_BASE5); +#ifdef GCALC_CHECK_WITH_FLOAT + long double x1, y1, x2, y2; + i1->calc_xy_ld(&x1, &y1); + i2->calc_xy_ld(&x2, &y2); + + if (result < 0) + { + if (!de_check(y1, y2) && y2 > y1) + GCALC_DBUG_PRINT(("floatcheck cmp_intersections_y failed %LG < %LG", + y2, y1)); + } + else if (result > 0) + { + if (!de_check(y1, y2) && y2 < y1) + GCALC_DBUG_PRINT(("floatcheck cmp_intersections_y failed %LG > %LG", + y2, y1)); + } + else + { + if (!de_check(y1, y2)) + GCALC_DBUG_PRINT(("floatcheck cmp_intersections_y failed %LG == %LG", + y2, y1)); + } +#endif /*GCALC_CHECK_WITH_FLOAT*/ + + if (result != 0) + return result; + + + inf1->calc_x_exp(); + inf2->calc_x_exp(); + gcalc_mul_coord(exp_a, GCALC_COORD_BASE5, + inf1->x_exp, GCALC_COORD_BASE3, inf2->t_b, GCALC_COORD_BASE2); + gcalc_mul_coord(exp_b, GCALC_COORD_BASE5, + inf2->x_exp, GCALC_COORD_BASE3, inf1->t_b, GCALC_COORD_BASE2); + + result= gcalc_cmp_coord(exp_a, exp_b, GCALC_COORD_BASE5); +#ifdef GCALC_CHECK_WITH_FLOAT + if (result < 0) + { + if (!de_check(x1, x2) && x2 > x1) + GCALC_DBUG_PRINT(("floatcheck cmp_intersectionsx failed %LG < %LG", + x2, x1)); + } + else if (result > 0) + { + if (!de_check(x1, x2) && x2 < x1) + GCALC_DBUG_PRINT(("floatcheck cmp_intersectionsx failed %LG > %LG", + x2, x1)); + } + else + { + if (!de_check(x1, x2)) + GCALC_DBUG_PRINT(("floatcheck cmp_intersectionsx failed %LG == %LG", + x2, x1)); + } +#endif /*GCALC_CHECK_WITH_FLOAT*/ + return result; +} +/* Internal coordinates implementation end */ + + +#define GCALC_SCALE_1 1e18 + +static double find_scale(double extent) +{ + double scale= 1e-2; + while (scale < extent) + scale*= (double ) 10; + return GCALC_SCALE_1 / scale / 10; +} + + +void Gcalc_heap::set_extent(double xmin, double xmax, double ymin, double ymax) +{ + xmin= fabs(xmin); + xmax= fabs(xmax); + ymin= fabs(ymin); + ymax= fabs(ymax); + + if (xmax < xmin) + xmax= xmin; + if (ymax < ymin) + ymax= ymin; + + coord_extent= xmax > ymax ? xmax : ymax; + coord_extent= find_scale(coord_extent); +#ifdef GCALC_CHECK_WITH_FLOAT + gcalc_coord_extent= &coord_extent; +#endif /*GCALC_CHECK_WITH_FLOAT*/ +} + + +void Gcalc_heap::free_point_info(Gcalc_heap::Info *i, + Gcalc_dyn_list::Item **i_hook) +{ + if (m_hook == &i->next) + m_hook= i_hook; + *i_hook= i->next; + free_item(i); + m_n_points--; +} + + +Gcalc_heap::Info *Gcalc_heap::new_point_info(double x, double y, + gcalc_shape_info shape) +{ + Info *result= (Info *)new_item(); + if (!result) + return NULL; + *m_hook= result; + m_hook= &result->next; + result->node.shape.x= x; + result->node.shape.y= y; + result->node.shape.shape= shape; + result->node.shape.top_node= 1; + result->type= nt_shape_node; + gcalc_set_double(result->node.shape.ix, x, coord_extent); + gcalc_set_double(result->node.shape.iy, y, coord_extent); + + m_n_points++; + return result; +} + + +static Gcalc_heap::Info *new_intersection( + Gcalc_heap *heap, Gcalc_scan_iterator::intersection_info *ii) +{ + Gcalc_heap::Info *isc= (Gcalc_heap::Info *)heap->new_item(); + if (!isc) + return 0; + isc->type= Gcalc_heap::nt_intersection; + isc->node.intersection.p1= ii->edge_a->pi; + isc->node.intersection.p2= ii->edge_a->next_pi; + isc->node.intersection.p3= ii->edge_b->pi; + isc->node.intersection.p4= ii->edge_b->next_pi; + isc->node.intersection.data= ii; + return isc; +} + + +static Gcalc_heap::Info *new_eq_point( + Gcalc_heap *heap, const Gcalc_heap::Info *p, + Gcalc_scan_iterator::point *edge) +{ + Gcalc_heap::Info *eqp= (Gcalc_heap::Info *)heap->new_item(); + if (!eqp) + return 0; + eqp->type= Gcalc_heap::nt_eq_node; + eqp->node.eq.node= p; + eqp->node.eq.data= edge; + return eqp; +} + + +void Gcalc_heap::Info::calc_xy(double *x, double *y) const +{ + double b0_x= node.intersection.p2->node.shape.x - node.intersection.p1->node.shape.x; + double b0_y= node.intersection.p2->node.shape.y - node.intersection.p1->node.shape.y; + double b1_x= node.intersection.p4->node.shape.x - node.intersection.p3->node.shape.x; + double b1_y= node.intersection.p4->node.shape.y - node.intersection.p3->node.shape.y; + double b0xb1= b0_x * b1_y - b0_y * b1_x; + 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; + + t/= b0xb1; + + *x= node.intersection.p1->node.shape.x + b0_x * t; + *y= node.intersection.p1->node.shape.y + b0_y * t; +} + + +#ifdef GCALC_CHECK_WITH_FLOAT +void Gcalc_heap::Info::calc_xy_ld(long double *x, long double *y) const +{ + long double b0_x= ((long double) p2->node.shape.x) - p1->node.shape.x; + long double b0_y= ((long double) p2->node.shape.y) - p1->node.shape.y; + long double b1_x= ((long double) p4->node.shape.x) - p3->node.shape.x; + long double b1_y= ((long double) p4->node.shape.y) - p3->node.shape.y; + long double b0xb1= b0_x * b1_y - b0_y * b1_x; + long double ax= ((long double) p3->node.shape.x) - p1->node.shape.x; + long double ay= ((long double) p3->node.shape.y) - p1->node.shape.y; + long double t_a= ax * b1_y - ay * b1_x; + long double hx= (b0xb1 * (long double) p1->node.shape.x + b0_x * t_a); + long double hy= (b0xb1 * (long double) p1->node.shape.y + b0_y * t_a); + + if (fabs(b0xb1) < 1e-15) + { + *x= p1->node.shape.x; + *y= p1->node.shape.y; + return; + } + + *x= hx/b0xb1; + *y= hy/b0xb1; +} +#endif /*GCALC_CHECK_WITH_FLOAT*/ + + +static int cmp_point_info(const Gcalc_heap::Info *i0, + const Gcalc_heap::Info *i1) +{ + int cmp_y= gcalc_cmp_coord1(i0->node.shape.iy, i1->node.shape.iy); + if (cmp_y) + return cmp_y; + return gcalc_cmp_coord1(i0->node.shape.ix, i1->node.shape.ix); +} + + +static inline void trim_node(Gcalc_heap::Info *node, Gcalc_heap::Info *prev_node) +{ + if (!node) + return; + node->node.shape.top_node= 0; + GCALC_DBUG_ASSERT((node->node.shape.left == prev_node) || (node->node.shape.right == prev_node)); + if (node->node.shape.left == prev_node) + node->node.shape.left= node->node.shape.right; + node->node.shape.right= NULL; + GCALC_DBUG_ASSERT(cmp_point_info(node, prev_node)); +} + + +static int compare_point_info(const void *e0, const void *e1) +{ + const Gcalc_heap::Info *i0= (const Gcalc_heap::Info *)e0; + const Gcalc_heap::Info *i1= (const Gcalc_heap::Info *)e1; + return cmp_point_info(i0, i1) > 0; +} + + +void Gcalc_heap::prepare_operation() +{ + Info *cur; + GCALC_DBUG_ASSERT(m_hook); + *m_hook= NULL; + m_hook= NULL; /* just to check it's not called twice */ + m_first= sort_list(compare_point_info, m_first, m_n_points); + + /* TODO - move this to the 'normal_scan' loop */ + for (cur= get_first(); cur; cur= cur->get_next()) + { + trim_node(cur->node.shape.left, cur); + trim_node(cur->node.shape.right, cur); + } +} + + +void Gcalc_heap::reset() +{ + if (m_n_points) + { + if (m_hook) + *m_hook= NULL; + free_list(m_first); + m_n_points= 0; + } + m_hook= &m_first; +} + + +int Gcalc_shape_transporter::int_single_point(gcalc_shape_info Info, + double x, double y) +{ + Gcalc_heap::Info *point= m_heap->new_point_info(x, y, Info); + if (!point) + return 1; + point->node.shape.left= point->node.shape.right= 0; + return 0; +} + + +int Gcalc_shape_transporter::int_add_point(gcalc_shape_info Info, + double x, double y) +{ + Gcalc_heap::Info *point; + Gcalc_dyn_list::Item **hook; + + hook= m_heap->get_cur_hook(); + + if (!(point= m_heap->new_point_info(x, y, Info))) + return 1; + if (m_first) + { + if (cmp_point_info(m_prev, point) == 0) + { + /* Coinciding points, do nothing */ + m_heap->free_point_info(point, hook); + return 0; + } + GCALC_DBUG_ASSERT(!m_prev || m_prev->node.shape.x != x || m_prev->node.shape.y != y); + m_prev->node.shape.left= point; + point->node.shape.right= m_prev; + } + else + m_first= point; + m_prev= point; + m_prev_hook= hook; + return 0; +} + + +void Gcalc_shape_transporter::int_complete() +{ + GCALC_DBUG_ASSERT(m_shape_started == 1 || m_shape_started == 3); + + if (!m_first) + return; + + /* simple point */ + if (m_first == m_prev) + { + m_first->node.shape.right= m_first->node.shape.left= NULL; + return; + } + + /* line */ + if (m_shape_started == 1) + { + m_first->node.shape.right= NULL; + m_prev->node.shape.left= m_prev->node.shape.right; + m_prev->node.shape.right= NULL; + return; + } + + /* polygon */ + if (cmp_point_info(m_first, m_prev) == 0) + { + /* Coinciding points, remove the last one from the list */ + m_prev->node.shape.right->node.shape.left= m_first; + m_first->node.shape.right= m_prev->node.shape.right; + m_heap->free_point_info(m_prev, m_prev_hook); + } + else + { + GCALC_DBUG_ASSERT(m_prev->node.shape.x != m_first->node.shape.x || m_prev->node.shape.y != m_first->node.shape.y); + m_first->node.shape.right= m_prev; + m_prev->node.shape.left= m_first; + } +} + + +inline void calc_dx_dy(Gcalc_scan_iterator::point *p) +{ + gcalc_sub_coord1(p->dx, p->next_pi->node.shape.ix, p->pi->node.shape.ix); + gcalc_sub_coord1(p->dy, p->next_pi->node.shape.iy, p->pi->node.shape.iy); + if (GCALC_SIGN(p->dx[0])) + { + p->l_border= &p->next_pi->node.shape.ix; + p->r_border= &p->pi->node.shape.ix; + } + else + { + p->r_border= &p->next_pi->node.shape.ix; + p->l_border= &p->pi->node.shape.ix; + } +} + + +Gcalc_scan_iterator::Gcalc_scan_iterator(size_t blk_size) : + Gcalc_dyn_list(blk_size, sizeof(point) > sizeof(intersection_info) ? + sizeof(point) : + sizeof(intersection_info)) +{ + state.slice= NULL; + m_bottom_points= NULL; + m_bottom_hook= &m_bottom_points; +} + + +void Gcalc_scan_iterator::init(Gcalc_heap *points) +{ + GCALC_DBUG_ASSERT(points->ready()); + GCALC_DBUG_ASSERT(!state.slice); + + if (!(m_cur_pi= points->get_first())) + return; + m_heap= points; + state.event_position_hook= &state.slice; + state.event_end= NULL; +#ifndef GCALC_DBUG_OFF + m_cur_thread= 0; +#endif /*GCALC_DBUG_OFF*/ + GCALC_SET_TERMINATED(killed, 0); +} + +void Gcalc_scan_iterator::reset() +{ + state.slice= NULL; + m_bottom_points= NULL; + m_bottom_hook= &m_bottom_points; + Gcalc_dyn_list::reset(); +} + + +int Gcalc_scan_iterator::point::cmp_dx_dy(const Gcalc_coord1 dx_a, + const Gcalc_coord1 dy_a, + const Gcalc_coord1 dx_b, + const Gcalc_coord1 dy_b) +{ + Gcalc_coord2 dx_a_dy_b; + Gcalc_coord2 dy_a_dx_b; + gcalc_mul_coord1(dx_a_dy_b, dx_a, dy_b); + gcalc_mul_coord1(dy_a_dx_b, dy_a, dx_b); + + return gcalc_cmp_coord(dx_a_dy_b, dy_a_dx_b, GCALC_COORD_BASE2); +} + + +int Gcalc_scan_iterator::point::cmp_dx_dy(const Gcalc_heap::Info *p1, + const Gcalc_heap::Info *p2, + const Gcalc_heap::Info *p3, + const Gcalc_heap::Info *p4) +{ + Gcalc_coord1 dx_a, dy_a, dx_b, dy_b; + gcalc_sub_coord1(dx_a, p2->node.shape.ix, p1->node.shape.ix); + gcalc_sub_coord1(dy_a, p2->node.shape.iy, p1->node.shape.iy); + gcalc_sub_coord1(dx_b, p4->node.shape.ix, p3->node.shape.ix); + gcalc_sub_coord1(dy_b, p4->node.shape.iy, p3->node.shape.iy); + return cmp_dx_dy(dx_a, dy_a, dx_b, dy_b); +} + + +int Gcalc_scan_iterator::point::cmp_dx_dy(const point *p) const +{ + GCALC_DBUG_ASSERT(!is_bottom()); + return cmp_dx_dy(dx, dy, p->dx, p->dy); +} + + +#ifdef GCALC_CHECK_WITH_FLOAT +void Gcalc_scan_iterator::point::calc_x(long double *x, long double y, + long double ix) const +{ + long double ddy= gcalc_get_double(dy, GCALC_COORD_BASE); + if (fabsl(ddy) < (long double) 1e-20) + { + *x= ix; + } + else + *x= (ddy * (long double) pi->node.shape.x + gcalc_get_double(dx, GCALC_COORD_BASE) * + (y - pi->node.shape.y)) / ddy; +} +#endif /*GCALC_CHECK_WITH_FLOAT*/ + + +static int compare_events(const void *e0, const void *e1) +{ + const Gcalc_scan_iterator::point *p0= (const Gcalc_scan_iterator::point *)e0; + const Gcalc_scan_iterator::point *p1= (const Gcalc_scan_iterator::point *)e1; + return p0->cmp_dx_dy(p1) > 0; +} + + +int Gcalc_scan_iterator::arrange_event(int do_sorting, int n_intersections) +{ + int ev_counter; + point *sp; + point **sp_hook; + + ev_counter= 0; + + *m_bottom_hook= NULL; + for (sp= m_bottom_points; sp; sp= sp->get_next()) + sp->ev_next= sp->get_next(); + + for (sp= state.slice, sp_hook= &state.slice; + sp; sp_hook= sp->next_ptr(), sp= sp->get_next()) + { + if (sp->event) + { + state.event_position_hook= sp_hook; + break; + } + } + + for (sp= *(sp_hook= state.event_position_hook); + sp && sp->event; sp_hook= sp->next_ptr(), sp= sp->get_next()) + { + ev_counter++; + if (sp->get_next() && sp->get_next()->event) + sp->ev_next= sp->get_next(); + else + sp->ev_next= m_bottom_points; + } + +#ifndef GCALC_DBUG_OFF + { + point *cur_p= sp; + for (; cur_p; cur_p= cur_p->get_next()) + GCALC_DBUG_ASSERT(!cur_p->event); + } +#endif /*GCALC_DBUG_OFF*/ + + state.event_end= sp; + + if (ev_counter == 2 && n_intersections == 1) + { + /* If we had only intersection, just swap the two points. */ + sp= *state.event_position_hook; + *state.event_position_hook= sp->get_next(); + sp->next= (*state.event_position_hook)->next; + (*state.event_position_hook)->next= sp; + + /* The list of the events should be restored. */ + (*state.event_position_hook)->ev_next= sp; + sp->ev_next= m_bottom_points; + } + else if (ev_counter == 2 && get_events()->event == scev_two_threads) + { + /* Do nothing. */ + } + else if (ev_counter > 1 && do_sorting) + { + point *cur_p; + *sp_hook= NULL; + sp= (point *) sort_list(compare_events, *state.event_position_hook, + ev_counter); + /* Find last item in the list, it's changed after the sorting. */ + for (cur_p= sp->get_next(); cur_p->get_next(); + cur_p= cur_p->get_next()) + {} + cur_p->next= state.event_end; + *state.event_position_hook= sp; + /* The list of the events should be restored. */ + for (; sp && sp->event; sp= sp->get_next()) + { + if (sp->get_next() && sp->get_next()->event) + sp->ev_next= sp->get_next(); + else + sp->ev_next= m_bottom_points; + } + } + +#ifndef GCALC_DBUG_OFF + { + const event_point *ev= get_events(); + for (; ev && ev->get_next(); ev= ev->get_next()) + { + if (ev->is_bottom() || ev->get_next()->is_bottom()) + break; + GCALC_DBUG_ASSERT(ev->cmp_dx_dy(ev->get_next()) <= 0); + } + } +#endif /*GCALC_DBUG_OFF*/ + return 0; +} + + +int Gcalc_heap::Info::equal_pi(const Info *pi) const +{ + if (type == nt_intersection) + return node.intersection.equal; + if (pi->type == nt_eq_node) + return 1; + if (type == nt_eq_node || pi->type == nt_intersection) + return 0; + return cmp_point_info(this, pi) == 0; +} + +int Gcalc_scan_iterator::step() +{ + int result= 0; + int do_sorting= 0; + int n_intersections= 0; + point *sp; + GCALC_DBUG_ENTER("Gcalc_scan_iterator::step"); + GCALC_DBUG_ASSERT(more_points()); + + if (GCALC_TERMINATED(killed)) + GCALC_DBUG_RETURN(0xFFFF); + + /* Clear the old event marks. */ + if (m_bottom_points) + { + free_list((Gcalc_dyn_list::Item **) &m_bottom_points, + (Gcalc_dyn_list::Item **) m_bottom_hook); + m_bottom_points= NULL; + m_bottom_hook= &m_bottom_points; + } + for (sp= *state.event_position_hook; + sp != state.event_end; sp= sp->get_next()) + sp->event= scev_none; + +//#ifndef GCALC_DBUG_OFF + state.event_position_hook= NULL; + state.pi= NULL; +//#endif /*GCALC_DBUG_OFF*/ + + do + { +#ifndef GCALC_DBUG_OFF + if (m_cur_pi->type == Gcalc_heap::nt_intersection && + m_cur_pi->get_next()->type == Gcalc_heap::nt_intersection && + m_cur_pi->node.intersection.equal) + GCALC_DBUG_ASSERT(cmp_intersections(m_cur_pi, m_cur_pi->get_next()) == 0); +#endif /*GCALC_DBUG_OFF*/ + GCALC_DBUG_CHECK_COUNTER(); + GCALC_DBUG_PRINT_SLICE("step:", state.slice); + GCALC_DBUG_PRINT_PI(m_cur_pi); + if (m_cur_pi->type == Gcalc_heap::nt_shape_node) + { + if (m_cur_pi->is_top()) + { + result= insert_top_node(); + if (!m_cur_pi->is_bottom()) + do_sorting++; + } + else if (m_cur_pi->is_bottom()) + remove_bottom_node(); + else + { + do_sorting++; + result= node_scan(); + } + if (result) + GCALC_DBUG_RETURN(result); + state.pi= m_cur_pi; + } + else if (m_cur_pi->type == Gcalc_heap::nt_eq_node) + { + do_sorting++; + eq_scan(); + } + else + { + /* nt_intersection */ + do_sorting++; + n_intersections++; + intersection_scan(); + if (!state.pi || state.pi->type == Gcalc_heap::nt_intersection) + state.pi= m_cur_pi; + } + + m_cur_pi= m_cur_pi->get_next(); + } while (m_cur_pi && state.pi->equal_pi(m_cur_pi)); + + GCALC_DBUG_RETURN(arrange_event(do_sorting, n_intersections)); +} + + +static int node_on_right(const Gcalc_heap::Info *node, + const Gcalc_heap::Info *edge_a, const Gcalc_heap::Info *edge_b) +{ + Gcalc_coord1 a_x, a_y; + Gcalc_coord1 b_x, b_y; + Gcalc_coord2 ax_by, ay_bx; + int result; + + gcalc_sub_coord1(a_x, node->node.shape.ix, edge_a->node.shape.ix); + gcalc_sub_coord1(a_y, node->node.shape.iy, edge_a->node.shape.iy); + gcalc_sub_coord1(b_x, edge_b->node.shape.ix, edge_a->node.shape.ix); + gcalc_sub_coord1(b_y, edge_b->node.shape.iy, edge_a->node.shape.iy); + gcalc_mul_coord1(ax_by, a_x, b_y); + gcalc_mul_coord1(ay_bx, a_y, b_x); + result= gcalc_cmp_coord(ax_by, ay_bx, GCALC_COORD_BASE2); +#ifdef GCALC_CHECK_WITH_FLOAT + { + long double dx= gcalc_get_double(edge_b->node.shape.ix, GCALC_COORD_BASE) - + gcalc_get_double(edge_a->node.shape.ix, GCALC_COORD_BASE); + long double dy= gcalc_get_double(edge_b->node.shape.iy, GCALC_COORD_BASE) - + gcalc_get_double(edge_a->node.shape.iy, GCALC_COORD_BASE); + long double ax= gcalc_get_double(node->node.shape.ix, GCALC_COORD_BASE) - + gcalc_get_double(edge_a->node.shape.ix, GCALC_COORD_BASE); + long double ay= gcalc_get_double(node->node.shape.iy, GCALC_COORD_BASE) - + gcalc_get_double(edge_a->node.shape.iy, GCALC_COORD_BASE); + long double d= ax * dy - ay * dx; + if (result == 0) + GCALC_DBUG_ASSERT(de_check(d, 0.0)); + else if (result < 0) + GCALC_DBUG_ASSERT(de_check(d, 0.0) || d < 0); + else + GCALC_DBUG_ASSERT(de_check(d, 0.0) || d > 0); + } +#endif /*GCALC_CHECK_WITH_FLOAT*/ + return result; +} + + +static int cmp_tops(const Gcalc_heap::Info *top_node, + const Gcalc_heap::Info *edge_a, const Gcalc_heap::Info *edge_b) +{ + int cmp_res_a, cmp_res_b; + + cmp_res_a= gcalc_cmp_coord1(edge_a->node.shape.ix, top_node->node.shape.ix); + cmp_res_b= gcalc_cmp_coord1(edge_b->node.shape.ix, top_node->node.shape.ix); + + if (cmp_res_a <= 0 && cmp_res_b > 0) + return -1; + if (cmp_res_b <= 0 && cmp_res_a > 0) + return 1; + if (cmp_res_a == 0 && cmp_res_b == 0) + return 0; + + return node_on_right(edge_a, top_node, edge_b); +} + + +int Gcalc_scan_iterator::insert_top_node() +{ + point *sp= state.slice; + point **prev_hook= &state.slice; + point *sp1= NULL; + point *sp0= new_slice_point(); + int cmp_res; + + GCALC_DBUG_ENTER("Gcalc_scan_iterator::insert_top_node"); + if (!sp0) + GCALC_DBUG_RETURN(1); + sp0->pi= m_cur_pi; + sp0->next_pi= m_cur_pi->node.shape.left; +#ifndef GCALC_DBUG_OFF + sp0->thread= m_cur_thread++; +#endif /*GCALC_DBUG_OFF*/ + if (m_cur_pi->node.shape.left) + { + calc_dx_dy(sp0); + if (m_cur_pi->node.shape.right) + { + if (!(sp1= new_slice_point())) + GCALC_DBUG_RETURN(1); + sp1->event= sp0->event= scev_two_threads; + sp1->pi= m_cur_pi; + sp1->next_pi= m_cur_pi->node.shape.right; +#ifndef GCALC_DBUG_OFF + sp1->thread= m_cur_thread++; +#endif /*GCALC_DBUG_OFF*/ + calc_dx_dy(sp1); + /* We have two threads so should decide which one will be first */ + cmp_res= cmp_tops(m_cur_pi, m_cur_pi->node.shape.left, m_cur_pi->node.shape.right); + if (cmp_res > 0) + { + point *tmp= sp0; + sp0= sp1; + sp1= tmp; + } + else if (cmp_res == 0) + { + /* Exactly same direction of the edges. */ + cmp_res= gcalc_cmp_coord1(m_cur_pi->node.shape.left->node.shape.iy, m_cur_pi->node.shape.right->node.shape.iy); + if (cmp_res != 0) + { + if (cmp_res < 0) + { + if (add_eq_node(sp0->next_pi, sp1)) + GCALC_DBUG_RETURN(1); + } + else + { + if (add_eq_node(sp1->next_pi, sp0)) + GCALC_DBUG_RETURN(1); + } + } + else + { + cmp_res= gcalc_cmp_coord1(m_cur_pi->node.shape.left->node.shape.ix, m_cur_pi->node.shape.right->node.shape.ix); + if (cmp_res != 0) + { + if (cmp_res < 0) + { + if (add_eq_node(sp0->next_pi, sp1)) + GCALC_DBUG_RETURN(1); + } + else + { + if (add_eq_node(sp1->next_pi, sp0)) + GCALC_DBUG_RETURN(1); + } + } + } + } + } + else + sp0->event= scev_thread; + } + else + sp0->event= scev_single_point; + + + /* Check if we already have an event - then we'll place the node there */ + for (; sp && !sp->event; prev_hook= sp->next_ptr(), sp=sp->get_next()) + {} + if (!sp) + { + sp= state.slice; + prev_hook= &state.slice; + /* We need to find the place to insert. */ + for (; sp; prev_hook= sp->next_ptr(), sp=sp->get_next()) + { + if (sp->event || gcalc_cmp_coord1(*sp->r_border, m_cur_pi->node.shape.ix) < 0) + continue; + cmp_res= node_on_right(m_cur_pi, sp->pi, sp->next_pi); + if (cmp_res == 0) + { + /* The top node lies on the edge. */ + /* Nodes of that edge will be handled in other places. */ + sp->event= scev_intersection; + } + else if (cmp_res < 0) + break; + } + } + + if (sp0->event == scev_single_point) + { + /* Add single point to the bottom list. */ + *m_bottom_hook= sp0; + m_bottom_hook= sp0->next_ptr(); + state.event_position_hook= prev_hook; + } + else + { + *prev_hook= sp0; + sp0->next= sp; + if (add_events_for_node(sp0)) + GCALC_DBUG_RETURN(1); + + if (sp0->event == scev_two_threads) + { + *prev_hook= sp1; + sp1->next= sp; + if (add_events_for_node(sp1)) + GCALC_DBUG_RETURN(1); + + sp0->next= sp1; + *prev_hook= sp0; + } + } + + GCALC_DBUG_RETURN(0); +} + + +void Gcalc_scan_iterator::remove_bottom_node() +{ + point *sp= state.slice; + point **sp_hook= &state.slice; + point *first_bottom_point= NULL; + + GCALC_DBUG_ENTER("Gcalc_scan_iterator::remove_bottom_node"); + for (; sp; sp= sp->get_next()) + { + if (sp->next_pi == m_cur_pi) + { + *sp_hook= sp->get_next(); + sp->pi= m_cur_pi; + sp->next_pi= NULL; + if (first_bottom_point) + { + first_bottom_point->event= sp->event= scev_two_ends; + break; + } + first_bottom_point= sp; + sp->event= scev_end; + state.event_position_hook= sp_hook; + } + else + sp_hook= sp->next_ptr(); + } + GCALC_DBUG_ASSERT(first_bottom_point); + *m_bottom_hook= first_bottom_point; + m_bottom_hook= first_bottom_point->next_ptr(); + if (sp) + { + *m_bottom_hook= sp; + m_bottom_hook= sp->next_ptr(); + } + + GCALC_DBUG_VOID_RETURN; +} + + +int Gcalc_scan_iterator::add_events_for_node(point *sp_node) +{ + point *sp= state.slice; + int cur_pi_r, sp_pi_r; + + GCALC_DBUG_ENTER("Gcalc_scan_iterator::add_events_for_node"); + + /* Scan to the event point. */ + for (; sp != sp_node; sp= sp->get_next()) + { + GCALC_DBUG_ASSERT(!sp->is_bottom()); + GCALC_DBUG_PRINT(("left cut_edge %d", sp->thread)); + if (sp->next_pi == sp_node->next_pi || + gcalc_cmp_coord1(*sp->r_border, *sp_node->l_border) < 0) + continue; + sp_pi_r= node_on_right(sp->next_pi, sp_node->pi, sp_node->next_pi); + if (sp_pi_r < 0) + continue; + cur_pi_r= node_on_right(sp_node->next_pi, sp->pi, sp->next_pi); + if (cur_pi_r > 0) + continue; + if (cur_pi_r == 0 && sp_pi_r == 0) + { + int cmp_res= cmp_point_info(sp->next_pi, sp_node->next_pi); + if (cmp_res > 0) + { + if (add_eq_node(sp_node->next_pi, sp)) + GCALC_DBUG_RETURN(1); + } + else if (cmp_res < 0) + { + if (add_eq_node(sp->next_pi, sp_node)) + GCALC_DBUG_RETURN(1); + } + continue; + } + + if (cur_pi_r == 0) + { + if (add_eq_node(sp_node->next_pi, sp)) + GCALC_DBUG_RETURN(1); + continue; + } + else if (sp_pi_r == 0) + { + if (add_eq_node(sp->next_pi, sp_node)) + GCALC_DBUG_RETURN(1); + continue; + } + + if (sp->event) + { +#ifndef GCALC_DBUG_OFF + cur_pi_r= node_on_right(sp_node->pi, sp->pi, sp->next_pi); + GCALC_DBUG_ASSERT(cur_pi_r == 0); +#endif /*GCALC_DBUG_OFF*/ + continue; + } + cur_pi_r= node_on_right(sp_node->pi, sp->pi, sp->next_pi); + GCALC_DBUG_ASSERT(cur_pi_r >= 0); + //GCALC_DBUG_ASSERT(cur_pi_r > 0); /* Is it ever violated? */ + if (cur_pi_r > 0 && add_intersection(sp, sp_node, m_cur_pi)) + GCALC_DBUG_RETURN(1); + } + + /* Scan to the end of the slice */ + sp= sp->get_next(); + + for (; sp; sp= sp->get_next()) + { + GCALC_DBUG_ASSERT(!sp->is_bottom()); + GCALC_DBUG_PRINT(("right cut_edge %d", sp->thread)); + if (sp->next_pi == sp_node->next_pi || + gcalc_cmp_coord1(*sp_node->r_border, *sp->l_border) < 0) + continue; + sp_pi_r= node_on_right(sp->next_pi, sp_node->pi, sp_node->next_pi); + if (sp_pi_r > 0) + continue; + cur_pi_r= node_on_right(sp_node->next_pi, sp->pi, sp->next_pi); + if (cur_pi_r < 0) + continue; + if (cur_pi_r == 0 && sp_pi_r == 0) + { + int cmp_res= cmp_point_info(sp->next_pi, sp_node->next_pi); + if (cmp_res > 0) + { + if (add_eq_node(sp_node->next_pi, sp)) + GCALC_DBUG_RETURN(1); + } + else if (cmp_res < 0) + { + if (add_eq_node(sp->next_pi, sp_node)) + GCALC_DBUG_RETURN(1); + } + continue; + } + if (cur_pi_r == 0) + { + if (add_eq_node(sp_node->next_pi, sp)) + GCALC_DBUG_RETURN(1); + continue; + } + else if (sp_pi_r == 0) + { + if (add_eq_node(sp->next_pi, sp_node)) + GCALC_DBUG_RETURN(1); + continue; + } + + if (sp->event) + { +#ifndef GCALC_DBUG_OFF + cur_pi_r= node_on_right(sp_node->pi, sp->pi, sp->next_pi); + GCALC_DBUG_ASSERT(cur_pi_r == 0); +#endif /*GCALC_DBUG_OFF*/ + continue; + } + cur_pi_r= node_on_right(sp_node->pi, sp->pi, sp->next_pi); + GCALC_DBUG_ASSERT(cur_pi_r <= 0); + //GCALC_DBUG_ASSERT(cur_pi_r < 0); /* Is it ever violated? */ + if (cur_pi_r < 0 && add_intersection(sp_node, sp, m_cur_pi)) + GCALC_DBUG_RETURN(1); + } + + GCALC_DBUG_RETURN(0); +} + + +int Gcalc_scan_iterator::node_scan() +{ + point *sp= state.slice; + Gcalc_heap::Info *cur_pi= m_cur_pi; + + GCALC_DBUG_ENTER("Gcalc_scan_iterator::node_scan"); + + /* Scan to the event point. */ + /* Can be avoided if we add link to the sp to the Info. */ + for (; sp->next_pi != cur_pi; sp= sp->get_next()) + {} + + GCALC_DBUG_PRINT(("node for %d", sp->thread)); + /* Handle the point itself. */ + sp->pi= cur_pi; + sp->next_pi= cur_pi->node.shape.left; + sp->event= scev_point; + calc_dx_dy(sp); + + GCALC_DBUG_RETURN(add_events_for_node(sp)); +} + + +void Gcalc_scan_iterator::eq_scan() +{ + point *sp= eq_sp(m_cur_pi); + GCALC_DBUG_ENTER("Gcalc_scan_iterator::eq_scan"); + +#ifndef GCALC_DBUG_OFF + { + point *cur_p= state.slice; + for (; cur_p && cur_p != sp; cur_p= cur_p->get_next()) + {} + GCALC_DBUG_ASSERT(cur_p); + } +#endif /*GCALC_DBUG_OFF*/ + if (!sp->event) + { + sp->event= scev_intersection; + sp->ev_pi= m_cur_pi; + } + + GCALC_DBUG_VOID_RETURN; +} + + +void Gcalc_scan_iterator::intersection_scan() +{ + intersection_info *ii= i_data(m_cur_pi); + GCALC_DBUG_ENTER("Gcalc_scan_iterator::intersection_scan"); + +#ifndef GCALC_DBUG_OFF + { + point *sp= state.slice; + for (; sp && sp != ii->edge_a; sp= sp->get_next()) + {} + GCALC_DBUG_ASSERT(sp); + for (; sp && sp != ii->edge_b; sp= sp->get_next()) + {} + GCALC_DBUG_ASSERT(sp); + } +#endif /*GCALC_DBUG_OFF*/ + + ii->edge_a->event= ii->edge_b->event= scev_intersection; + ii->edge_a->ev_pi= ii->edge_b->ev_pi= m_cur_pi; + free_item(ii); + m_cur_pi->node.intersection.data= NULL; + + GCALC_DBUG_VOID_RETURN; +} + + +int Gcalc_scan_iterator::add_intersection(point *sp_a, point *sp_b, + Gcalc_heap::Info *pi_from) +{ + Gcalc_heap::Info *ii; + intersection_info *i_calc; + int cmp_res; + int skip_next= 0; + + GCALC_DBUG_ENTER("Gcalc_scan_iterator::add_intersection"); + if (!(i_calc= new_intersection_info(sp_a, sp_b)) || + !(ii= new_intersection(m_heap, i_calc))) + GCALC_DBUG_RETURN(1); + + ii->node.intersection.equal= 0; + + for (; + pi_from->get_next() != sp_a->next_pi && + pi_from->get_next() != sp_b->next_pi; + pi_from= pi_from->get_next()) + { + Gcalc_heap::Info *cur= pi_from->get_next(); + if (skip_next) + { + if (cur->type == Gcalc_heap::nt_intersection) + skip_next= cur->node.intersection.equal; + else + skip_next= 0; + continue; + } + if (cur->type == Gcalc_heap::nt_intersection) + { + cmp_res= cmp_intersections(cur, ii); + skip_next= cur->node.intersection.equal; + } + else if (cur->type == Gcalc_heap::nt_eq_node) + continue; + else + cmp_res= cmp_node_isc(cur, ii); + if (cmp_res == 0) + { + ii->node.intersection.equal= 1; + break; + } + else if (cmp_res > 0) + break; + } + + /* Intersection inserted before the equal point. */ + ii->next= pi_from->get_next(); + pi_from->next= ii; + + GCALC_DBUG_RETURN(0); +} + + +int Gcalc_scan_iterator::add_eq_node(Gcalc_heap::Info *node, point *sp) +{ + Gcalc_heap::Info *en; + + GCALC_DBUG_ENTER("Gcalc_scan_iterator::add_intersection"); + en= new_eq_point(m_heap, node, sp); + if (!en) + GCALC_DBUG_RETURN(1); + + /* eq_node inserted after the equal point. */ + en->next= node->get_next(); + node->next= en; + + GCALC_DBUG_RETURN(0); +} + + +void calc_t(Gcalc_coord2 t_a, Gcalc_coord2 t_b, + Gcalc_coord1 dxa, Gcalc_coord1 dxb, + const Gcalc_heap::Info *p1, const Gcalc_heap::Info *p2, + const Gcalc_heap::Info *p3, const Gcalc_heap::Info *p4) +{ + Gcalc_coord1 a2_a1x, a2_a1y; + Gcalc_coord2 x1y2, x2y1; + Gcalc_coord1 dya, dyb; + + gcalc_sub_coord1(a2_a1x, p3->node.shape.ix, p1->node.shape.ix); + gcalc_sub_coord1(a2_a1y, p3->node.shape.iy, p1->node.shape.iy); + + gcalc_sub_coord1(dxa, p2->node.shape.ix, p1->node.shape.ix); + gcalc_sub_coord1(dya, p2->node.shape.iy, p1->node.shape.iy); + gcalc_sub_coord1(dxb, p4->node.shape.ix, p3->node.shape.ix); + gcalc_sub_coord1(dyb, p4->node.shape.iy, p3->node.shape.iy); + + gcalc_mul_coord1(x1y2, dxa, dyb); + gcalc_mul_coord1(x2y1, dya, dxb); + gcalc_sub_coord(t_b, GCALC_COORD_BASE2, x1y2, x2y1); + + + gcalc_mul_coord1(x1y2, a2_a1x, dyb); + gcalc_mul_coord1(x2y1, a2_a1y, dxb); + gcalc_sub_coord(t_a, GCALC_COORD_BASE2, x1y2, x2y1); +} + + +double Gcalc_scan_iterator::get_y() const +{ + if (state.pi->type == Gcalc_heap::nt_intersection) + { + Gcalc_coord1 dxa, dya; + Gcalc_coord2 t_a, t_b; + Gcalc_coord3 a_tb, b_ta, y_exp; + calc_t(t_a, t_b, dxa, dya, + state.pi->node.intersection.p1, state.pi->node.intersection.p2, state.pi->node.intersection.p3, state.pi->node.intersection.p4); + + + gcalc_mul_coord(a_tb, GCALC_COORD_BASE3, + t_b, GCALC_COORD_BASE2, state.pi->node.intersection.p1->node.shape.iy, GCALC_COORD_BASE); + gcalc_mul_coord(b_ta, GCALC_COORD_BASE3, + t_a, GCALC_COORD_BASE2, dya, GCALC_COORD_BASE); + + gcalc_add_coord(y_exp, GCALC_COORD_BASE3, a_tb, b_ta); + + return (get_pure_double(y_exp, GCALC_COORD_BASE3) / + get_pure_double(t_b, GCALC_COORD_BASE2)) / m_heap->coord_extent; + } + else + return state.pi->node.shape.y; +} + + +double Gcalc_scan_iterator::get_event_x() const +{ + if (state.pi->type == Gcalc_heap::nt_intersection) + { + Gcalc_coord1 dxa, dya; + Gcalc_coord2 t_a, t_b; + Gcalc_coord3 a_tb, b_ta, x_exp; + calc_t(t_a, t_b, dxa, dya, + state.pi->node.intersection.p1, state.pi->node.intersection.p2, state.pi->node.intersection.p3, state.pi->node.intersection.p4); + + + gcalc_mul_coord(a_tb, GCALC_COORD_BASE3, + t_b, GCALC_COORD_BASE2, state.pi->node.intersection.p1->node.shape.ix, GCALC_COORD_BASE); + gcalc_mul_coord(b_ta, GCALC_COORD_BASE3, + t_a, GCALC_COORD_BASE2, dxa, GCALC_COORD_BASE); + + gcalc_add_coord(x_exp, GCALC_COORD_BASE3, a_tb, b_ta); + + return (get_pure_double(x_exp, GCALC_COORD_BASE3) / + get_pure_double(t_b, GCALC_COORD_BASE2)) / m_heap->coord_extent; + } + else + return state.pi->node.shape.x; +} + +double Gcalc_scan_iterator::get_h() const +{ + double cur_y= get_y(); + double next_y; + if (state.pi->type == Gcalc_heap::nt_intersection) + { + double x; + state.pi->calc_xy(&x, &next_y); + } + else + next_y= state.pi->next ? state.pi->get_next()->node.shape.y : 0.0; + return next_y - cur_y; +} + + +double Gcalc_scan_iterator::get_sp_x(const point *sp) const +{ + double dy; + if (sp->event & (scev_end | scev_two_ends | scev_point)) + return sp->pi->node.shape.x; + dy= sp->next_pi->node.shape.y - sp->pi->node.shape.y; + if (fabs(dy) < 1e-12) + return sp->pi->node.shape.x; + return sp->pi->node.shape.x + (sp->next_pi->node.shape.x - sp->pi->node.shape.x) * dy; +} + + +double Gcalc_scan_iterator::get_pure_double(const Gcalc_internal_coord *d, + int d_len) +{ + int n= 1; + long double res= (long double) FIRST_DIGIT(d[0]); + do + { + res*= (long double) GCALC_DIG_BASE; + res+= (long double) d[n]; + } while(++n < d_len); + + if (GCALC_SIGN(d[0])) + res*= -1.0; + return res; +} + + +#endif /* HAVE_SPATIAL */ |