diff options
author | Alexey Botchkov <holyfoot@askmonty.org> | 2011-09-21 00:04:41 +0500 |
---|---|---|
committer | Alexey Botchkov <holyfoot@askmonty.org> | 2011-09-21 00:04:41 +0500 |
commit | 0249413a6ae4a5900b51f434b361b9efa6133ad1 (patch) | |
tree | ce34c73e2d9ad3d7b9d08612c35b65141b55821a /sql/gcalc_slicescan.cc | |
parent | 5a04ac7bf0ee359fb36f4ebc00f0aebc142d36b7 (diff) | |
download | mariadb-git-0249413a6ae4a5900b51f434b361b9efa6133ad1.tar.gz |
several bugs fixed here.
849789 Second assertion `m_poly_borders->next' failed in Gcalc_operation_reducer::count_slice in maria-5.3-gis
849791 Fourth assertion `n > 0 && n < SINUSES_CALCULATED*2+1' in get_n_sincos
849789 Second assertion `m_poly_borders->next' failed in Gcalc_operation_reducer::count_slice in maria-5.3-gis
848901 Assertion `fabs(cur_isc->x-m_cur_intersection->x) + fabs(cur_isc->y-m_cur_intersection->y) < 0.000000000001' failed in Gcalc_scan_iterator::intersection_scan() in maria-5.3-gis
per-file comments:
mysql-test/r/gis-precise.result
test result updated.
mysql-test/r/gis.result
test result updated.
sql/gcalc_slicescan.cc
bugfixes.
sql/gcalc_slicescan.h
bugfixes.
sql/gcalc_tools.cc
bugfixes.
sql/gcalc_tools.h
bugfixes.
sql/item_geofunc.cc
bugfixes.
sql/spatial.cc
bugfixes.
Diffstat (limited to 'sql/gcalc_slicescan.cc')
-rw-r--r-- | sql/gcalc_slicescan.cc | 1326 |
1 files changed, 1189 insertions, 137 deletions
diff --git a/sql/gcalc_slicescan.cc b/sql/gcalc_slicescan.cc index 6f41a5fc371..07120d2e87e 100644 --- a/sql/gcalc_slicescan.cc +++ b/sql/gcalc_slicescan.cc @@ -116,6 +116,450 @@ void Gcalc_dyn_list::reset() } +/* Internal coordinate operations implementations */ + +void Gcalc_internal_coord::set_zero() +{ + int n_res= 0; + do + { + digits[n_res++]= 0; + } while (n_res < n_digits); + sign= 0; +} + + +int Gcalc_internal_coord::is_zero() const +{ + int n_res= 0; + do + { + if (digits[n_res++] != 0) + return 0; + } while (n_res < n_digits); + return 1; +} + + +#ifdef GCALC_CHECK_WITH_FLOAT +long double Gcalc_internal_coord::get_double() const +{ + int n= 0; + long double res= 0.0; + do + { + res*= (long double) DIG_BASE; + res+= digits[n]; + } while(++n < n_digits); + + n= 0; + do + { + if (n & 1) + res/= (long double) C_SCALE; + } while(++n < n_digits); + + if (sign) + res*= -1.0; + return res; +} +#endif /*GCALC_CHECK_WITH_FLOAT*/ + + +static void do_add(Gcalc_internal_coord *result, + const Gcalc_internal_coord *a, + const Gcalc_internal_coord *b) +{ + DBUG_ASSERT(a->n_digits == b->n_digits); + DBUG_ASSERT(a->n_digits == result->n_digits); + int n_digit= a->n_digits-1; + coord_digit_t carry= 0; + + do + { + if ((result->digits[n_digit]= + a->digits[n_digit] + b->digits[n_digit] + carry) >= DIG_BASE) + { + carry= 1; + result->digits[n_digit]-= DIG_BASE; + } + else + carry= 0; + } while (n_digit--); + DBUG_ASSERT(carry == 0); + result->sign= a->sign; +} + + +static void do_sub(Gcalc_internal_coord *result, + const Gcalc_internal_coord *a, + const Gcalc_internal_coord *b) +{ + DBUG_ASSERT(a->n_digits == b->n_digits); + DBUG_ASSERT(a->n_digits == result->n_digits); + int n_digit= a->n_digits-1; + coord_digit_t carry= 0; + + do + { + if ((result->digits[n_digit]= + a->digits[n_digit] - b->digits[n_digit] - carry) < 0) + { + carry= 1; + result->digits[n_digit]+= DIG_BASE; + } + else + carry= 0; + } while (n_digit--); + DBUG_ASSERT(carry == 0); + if (a->sign && result->is_zero()) + result->sign= 0; + else + result->sign= a->sign; +} + + +static int do_cmp(const Gcalc_internal_coord *a, + const Gcalc_internal_coord *b) +{ + DBUG_ASSERT(a->n_digits == b->n_digits); + int n_digit= 0; + + do + { + coord_digit_t d= a->digits[n_digit] - b->digits[n_digit]; + if (d > 0) + return 1; + if (d < 0) + return -1; + n_digit++; + } while (n_digit < a->n_digits); + + return 0; +} + + +#ifdef GCALC_CHECK_WITH_FLOAT +static int de_check(long double a, long double b) +{ + long double d= a - b; + if (d < (long double) 1e-6 && d > (long double) -1e-6) + return 1; + printf("xxx\n"); + return 0; +} + +static int de_check1(long double a, long double b) +{ + long double d= a - b; + if (d < (long double) 1e-6 && d > (long double) -1e-6) + return 1; + return 0; +} +static int de_weak(long double a, long double b) +{ + long double d= a - b; + if (d < (long double) 1 && d > (long double) -1) + return 1; + return 0; +} +#endif /*GCALC_CHECK_WITH_FLOAT*/ + + +void gcalc_mul_coord(Gcalc_internal_coord *result, + const Gcalc_internal_coord *a, + const Gcalc_internal_coord *b) +{ + DBUG_ASSERT(result->n_digits == a->n_digits + b->n_digits); + int n_a, n_b, n_res; + coord_digit_t carry= 0; + + result->set_zero(); + if (a->is_zero() || b->is_zero()) + return; + + n_a= a->n_digits - 1; + do + { + n_b= b->n_digits - 1; + do + { + coord2 mul= ((coord2) a->digits[n_a]) * b->digits[n_b] + carry + + result->digits[n_a + n_b + 1]; + result->digits[n_a + n_b + 1]= mul % DIG_BASE; + carry= mul / DIG_BASE; + } while (n_b--); + if (carry) + { + for (n_res= n_a; (result->digits[n_res]+= carry) >= DIG_BASE; n_res--) + { + result->digits[n_res]-= DIG_BASE; + carry= 1; + } + carry= 0; + } + } while (n_a--); + result->sign= a->sign != b->sign; +#ifdef GCALC_CHECK_WITH_FLOAT + DBUG_ASSERT(de_check(a->get_double() * b->get_double(), + result->get_double())); +#endif /*GCALC_CHECK_WITH_FLOAT*/ +} + + +void gcalc_add_coord(Gcalc_internal_coord *result, + const Gcalc_internal_coord *a, + const Gcalc_internal_coord *b) +{ + if (a->sign == b->sign) + do_add(result, a, b); + else + { + int cmp_res= do_cmp(a, b); + if (cmp_res == 0) + result->set_zero(); + else if (cmp_res > 0) + do_sub(result, a, b); + else + do_sub(result, b, a); + } +#ifdef GCALC_CHECK_WITH_FLOAT + DBUG_ASSERT(de_check(a->get_double() + b->get_double(), + result->get_double())); +#endif /*GCALC_CHECK_WITH_FLOAT*/ +} + + +void gcalc_sub_coord(Gcalc_internal_coord *result, + const Gcalc_internal_coord *a, + const Gcalc_internal_coord *b) +{ + if (a->sign != b->sign) + do_add(result, a, b); + else + { + int cmp_res= do_cmp(a, b); + if (cmp_res == 0) + result->set_zero(); + else if (cmp_res > 0) + do_sub(result, a, b); + else + { + do_sub(result, b, a); + result->sign= 1 - result->sign; + } + } +#ifdef GCALC_CHECK_WITH_FLOAT + DBUG_ASSERT(de_check(a->get_double() - b->get_double(), + result->get_double())); +#endif /*GCALC_CHECK_WITH_FLOAT*/ +} + + +int gcalc_cmp_coord(const Gcalc_internal_coord *a, + const Gcalc_internal_coord *b) +{ + int result; + if (a->sign != b->sign) + return a->sign ? -1 : 1; + result= a->sign ? do_cmp(b, a) : do_cmp(a, b); +#ifdef GCALC_CHECK_WITH_FLOAT + if (result == 0) + DBUG_ASSERT(de_check(a->get_double(), b->get_double())); + else if (result == 1) + DBUG_ASSERT(de_check1(a->get_double(), b->get_double()) || + a->get_double() > b->get_double()); + else + DBUG_ASSERT(de_check1(a->get_double(), b->get_double()) || + a->get_double() < b->get_double()); +#endif /*GCALC_CHECK_WITH_FLOAT*/ + return result; +} + + +int Gcalc_coord1::set_double(double d) +{ + double ds= d * C_SCALE; + init(); + if ((sign= ds < 0)) + ds= -ds; + c[0]= (coord_digit_t) (ds / (double) DIG_BASE); + c[1]= (coord_digit_t) (ds - ((double) c[0]) * DIG_BASE); +#ifdef GCALC_CHECK_WITH_FLOAT + DBUG_ASSERT(de_check(d, get_double())); +#endif /*GCALC_CHECK_WITH_FLOAT*/ + return 0; +} + + +void Gcalc_coord1::copy(const Gcalc_coord1 *from) +{ + c[0]= from->c[0]; + c[1]= from->c[1]; + sign= from->sign; +} + +/* Internal coordinates implementation end */ + + +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; + m_n_points++; + result->x= x; + result->y= y; + result->shape= shape; + result->ix.set_double(x); + result->iy.set_double(y); + return result; +} + + +Gcalc_heap::Intersection_info *Gcalc_heap::new_intersection( + const Info *p1, const Info *p2, + const Info *p3, const Info *p4) +{ + Intersection_info *isc= (Intersection_info *)new_item(); + if (!isc) + return 0; + isc->p1= p1; + isc->p2= p2; + isc->p3= p3; + isc->p4= p4; + *m_intersection_hook= isc; + m_intersection_hook= &isc->next; + return isc; +} + + +class Gcalc_coord3 : public Gcalc_internal_coord +{ + coord_digit_t c[COORD_BASE*3]; + public: + void init() + { + n_digits= COORD_BASE*3; + digits= c; + } +}; + + +static void calc_t(Gcalc_coord2 *t_a, Gcalc_coord2 *t_b, + Gcalc_coord1 *b1x, + Gcalc_coord1 *b1y, + 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_coord1 b2x, b2y; + Gcalc_coord2 x1y2, x2y1; + + a2_a1x.init(); + a2_a1y.init(); + x1y2.init(); + x2y1.init(); + t_a->init(); + t_b->init(); + b1y->init(); + b1x->init(); + b2x.init(); + b2y.init(); + + gcalc_sub_coord(&a2_a1x, &p3->ix, &p1->ix); + gcalc_sub_coord(&a2_a1y, &p3->iy, &p1->iy); + gcalc_sub_coord(b1x, &p2->ix, &p1->ix); + gcalc_sub_coord(b1y, &p2->iy, &p1->iy); + gcalc_sub_coord(&b2x, &p4->ix, &p3->ix); + gcalc_sub_coord(&b2y, &p4->iy, &p3->iy); + + gcalc_mul_coord(&x1y2, b1x, &b2y); + gcalc_mul_coord(&x2y1, &b2x, b1y); + gcalc_sub_coord(t_b, &x1y2, &x2y1); + + + gcalc_mul_coord(&x1y2, &a2_a1x, &b2y); + gcalc_mul_coord(&x2y1, &a2_a1y, &b2x); + gcalc_sub_coord(t_a, &x1y2, &x2y1); +} + + +inline void calc_t(Gcalc_coord2 *t_a, Gcalc_coord2 *t_b, + Gcalc_coord1 *b1x, + Gcalc_coord1 *b1y, + const Gcalc_heap::Intersection_info *isc) +{ + calc_t(t_a, t_b, b1x, b1y, isc->p1, isc->p2, isc->p3, isc->p4); +} + + +void Gcalc_heap::Intersection_info::calc_xy(double *x, double *y) const +{ + double b0_x= p2->x - p1->x; + double b0_y= p2->y - p1->y; + double b1_x= p4->x - p3->x; + double b1_y= p4->y - p3->y; + double b0xb1= b0_x * b1_y - b0_y * b1_x; + double t= (p3->x - p1->x) * b1_y - (p3->y - p1->y) * b1_x; + + t/= b0xb1; + + *x= p1->x + b0_x * t; + *y= p1->y + b0_y * t; +} + + +#ifdef GCALC_CHECK_WITH_FLOAT +void Gcalc_heap::Intersection_info::calc_xy_ld(long double *x, + long double *y) const +{ + long double b0_x= p2->x - p1->x; + long double b0_y= p2->y - p1->y; + long double b1_x= p4->x - p3->x; + long double b1_y= p4->y - p3->y; + long double b0xb1= b0_x * b1_y - b0_y * b1_x; + long double t= (p3->x - p1->x) * b1_y - (p3->y - p1->y) * b1_x; + long double cx, cy; + + t/= b0xb1; + + cx= (long double) p1->x + b0_x * t; + cy= (long double) p1->y + b0_y * t; + + Gcalc_coord2 t_a, t_b; + Gcalc_coord1 yb, xb; + Gcalc_coord3 m1, m2, sum; + + calc_t(&t_a, &t_b, &xb, &yb, this); + if (t_b.is_zero()) + { + *x= p1->x; + *y= p1->y; + return; + } + + m1.init(); + m2.init(); + sum.init(); + gcalc_mul_coord(&m1, &p1->ix, &t_b); + gcalc_mul_coord(&m2, &xb, &t_a); + gcalc_add_coord(&sum, &m1, &m2); + *x= sum.get_double() / t_b.get_double(); + + gcalc_mul_coord(&m1, &p1->iy, &t_b); + gcalc_mul_coord(&m2, &yb, &t_a); + gcalc_add_coord(&sum, &m1, &m2); + *y= sum.get_double() / t_b.get_double(); +} +#endif /*GCALC_CHECK_WITH_FLOAT*/ + + static inline void trim_node(Gcalc_heap::Info *node, Gcalc_heap::Info *prev_node) { if (!node) @@ -127,13 +571,21 @@ static inline void trim_node(Gcalc_heap::Info *node, Gcalc_heap::Info *prev_node } +static int cmp_point_info(const Gcalc_heap::Info *i0, + const Gcalc_heap::Info *i1) +{ + int cmp_y= gcalc_cmp_coord(&i0->iy, &i1->iy); + if (cmp_y) + return cmp_y; + return gcalc_cmp_coord(&i0->ix, &i1->ix); +} + + 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; - if (!coord_eq(i0->y, i1->y)) - return i0->y > i1->y; - return i0->x > i1->x; + return cmp_point_info(i0, i1) > 0; } @@ -155,6 +607,7 @@ void Gcalc_heap::prepare_operation() void Gcalc_heap::reset() { +#ifdef TMP_BLOCK if (!m_hook) { m_hook= &m_first; @@ -164,10 +617,18 @@ void Gcalc_heap::reset() *m_hook= m_free; m_free= m_first; +#endif /*TMP_BLOCK*/ + if (m_n_points) + { + free_list(m_first); + free_list(m_first_intersection, m_intersection_hook); + m_intersection_hook= (Gcalc_dyn_list::Item **) &m_first_intersection; + m_n_points= 0; + } m_hook= &m_first; - m_n_points= 0; } + int Gcalc_shape_transporter::int_single_point(gcalc_shape_info Info, double x, double y) { @@ -229,6 +690,7 @@ void Gcalc_shape_transporter::int_complete() } +#ifdef TMP_BLOCK inline int GET_DX_DY(double *dxdy, const Gcalc_heap::Info *p0, const Gcalc_heap::Info *p1) { @@ -238,6 +700,13 @@ inline int GET_DX_DY(double *dxdy, (*dxdy/= dy)>DBL_MAX || (*dxdy)<-DBL_MAX; } +#endif /*TMP_BLOCK*/ +inline void calc_dx_dy(Gcalc_scan_iterator::point *p) +{ + gcalc_sub_coord(&p->dx, &p->next_pi->ix, &p->pi->ix); + gcalc_sub_coord(&p->dy, &p->next_pi->iy, &p->pi->iy); +} + Gcalc_scan_iterator::Gcalc_scan_iterator(size_t blk_size) : Gcalc_dyn_list(blk_size, @@ -268,14 +737,20 @@ void Gcalc_scan_iterator::init(Gcalc_heap *points) if (!(m_cur_pi= points->get_first())) return; + m_heap= points; m_cur_thread= 0; m_intersections= NULL; + m_cur_intersection= NULL; m_next_is_top_point= true; m_events= NULL; current_state= &state0; next_state= &state1; saved_state= &state_s; +#ifdef TMP_BLOCK next_state->y= m_cur_pi->y; +#endif /*TMP_BLOCK*/ + next_state->intersection_scan= 0; + next_state->pi= m_cur_pi; } void Gcalc_scan_iterator::reset() @@ -286,21 +761,36 @@ void Gcalc_scan_iterator::reset() } -void Gcalc_scan_iterator::point::copy_core(point *from) +void Gcalc_scan_iterator::point::copy_core(const point *from) { +#ifdef TMP_BLOCK dx_dy= from->dx_dy; horiz_dir= from->horiz_dir; +#endif /*TMP_BLOCK*/ pi= from->pi; next_pi= from->next_pi; thread= from->thread; + dx.copy(&from->dx); + dy.copy(&from->dy); #ifdef TO_REMOVE from->next_link= this; #endif /*TO_REMOVE*/ } -int Gcalc_scan_iterator::point::compare_dx_dy(int horiz_dir_a, double dx_dy_a, - int horiz_dir_b, double dx_dy_b) +void Gcalc_scan_iterator::point::copy_all(const point *from) +{ + pi= from->pi; + next_pi= from->next_pi; + thread= from->thread; + dx.copy(&from->dx); + dy.copy(&from->dy); + intersection_link= from->intersection_link; + event= from->event; +} +#ifdef TMP_BLOCK +int Gcalc_scan_iterator::point::cmp_dx_dy(int horiz_dir_a, double dx_dy_a, + int horiz_dir_b, double dx_dy_b) { if (!horiz_dir_a && !horiz_dir_b) { @@ -323,7 +813,169 @@ int Gcalc_scan_iterator::point::cmp_dx_dy(const point *p) const return p->is_bottom() ? 0 : -1; if (p->is_bottom()) return 1; - return compare_dx_dy(horiz_dir, dx_dy, p->horiz_dir, p->dx_dy); + return cmp_dx_dy(horiz_dir, dx_dy, p->horiz_dir, p->dx_dy); +} +#endif /*TMP_BLOCK*/ +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; + dx_a_dy_b.init(); + dy_a_dx_b.init(); + gcalc_mul_coord(&dx_a_dy_b, dx_a, dy_b); + gcalc_mul_coord(&dy_a_dx_b, dy_a, dx_b); + + return gcalc_cmp_coord(&dx_a_dy_b, &dy_a_dx_b); +} + + +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; + dx_a.init(); + dx_b.init(); + dy_a.init(); + dy_b.init(); + gcalc_sub_coord(&dx_a, &p2->ix, &p1->ix); + gcalc_sub_coord(&dy_a, &p2->iy, &p1->iy); + gcalc_sub_coord(&dx_b, &p4->ix, &p3->ix); + gcalc_sub_coord(&dy_b, &p4->iy, &p3->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 +{ + if (is_bottom()) + return p->is_bottom() ? 0 : -1; + if (p->is_bottom()) + return 1; + 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 +{ + if (fabsl(dy.get_double()) < (long double) 1e-15) + { + *x= ix; + } + else + *x= (long double) pi->x + dx.get_double()/dy.get_double() * (y - pi->y); +} +#endif /*GCALC_CHECK_WITH_FLOAT*/ + + +static int cmp_sp_pi(const Gcalc_scan_iterator::point *sp, + const Gcalc_heap::Info *pi) +{ + Gcalc_coord1 dx_pi, dy_pi; + Gcalc_coord2 dx_sp_dy_pi, dy_sp_dx_pi; + + if (!sp->next_pi) + return cmp_point_info(sp->pi, pi); + + dx_pi.init(); + dy_pi.init(); + dx_sp_dy_pi.init(); + dy_sp_dx_pi.init(); + + + gcalc_sub_coord(&dx_pi, &pi->ix, &sp->pi->ix); + gcalc_sub_coord(&dy_pi, &pi->iy, &sp->pi->iy); + gcalc_mul_coord(&dx_sp_dy_pi, &sp->dx, &dy_pi); + gcalc_mul_coord(&dy_sp_dx_pi, &sp->dy, &dx_pi); + + int result= gcalc_cmp_coord(&dx_sp_dy_pi, &dy_sp_dx_pi); +#ifdef GCALC_CHECK_WITH_FLOAT + long double sp_x; + sp->calc_x(&sp_x, pi->y, pi->x); + if (result == 0) + DBUG_ASSERT(de_check1(sp->dy.get_double(), 0.0) || + de_check(sp_x, pi->x)); + if (result < 0) + DBUG_ASSERT(de_check1(sp->dy.get_double(), 0.0) || + de_check1(sp_x, pi->x) || + sp_x < pi->x); + if (result > 0) + DBUG_ASSERT(de_check1(sp->dy.get_double(), 0.0) || + de_check1(sp_x, pi->x) || + sp_x > pi->x); +#endif /*GCALC_CHECK_WITH_FLOAT*/ + return result; +} + + +static void calc_cmp_sp_sp_exp(Gcalc_coord3 *result, + const Gcalc_scan_iterator::point *a, + const Gcalc_scan_iterator::point *b, + const Gcalc_coord1 *ly) +{ + Gcalc_coord2 x_dy, dx_ly, sum; + + x_dy.init(); + dx_ly.init(); + sum.init(); + result->init(); + + gcalc_mul_coord(&x_dy, &a->pi->ix, &a->dy); + gcalc_mul_coord(&dx_ly, &a->dx, ly); + gcalc_add_coord(&sum, &x_dy, &dx_ly); + gcalc_mul_coord(result, &sum, &b->dy); +} + + +static int cmp_sp_sp_cnt(const Gcalc_scan_iterator::point *a, + const Gcalc_scan_iterator::point *b, + const Gcalc_coord1 *y) +{ + Gcalc_coord1 lya, lyb; + Gcalc_coord3 a_exp, b_exp; + + lya.init(); + lyb.init(); + gcalc_sub_coord(&lya, y, &a->pi->iy); + gcalc_sub_coord(&lyb, y, &b->pi->iy); + + calc_cmp_sp_sp_exp(&a_exp, a, b, &lya); + calc_cmp_sp_sp_exp(&b_exp, b, a, &lyb); + + int result= gcalc_cmp_coord(&a_exp, &b_exp); +#ifdef GCALC_CHECK_WITH_FLOAT + long double a_x, b_x; + a->calc_x(&a_x, y->get_double(), 0); + b->calc_x(&b_x, y->get_double(), 0); + if (result == 0) + DBUG_ASSERT(de_check(a_x, b_x)); + if (result < 0) + DBUG_ASSERT(de_check1(a_x, b_x) || a_x < b_x); + if (result > 0) + DBUG_ASSERT(de_check1(a_x, b_x) || a_x > b_x); +#endif /*GCALC_CHECK_WITH_FLOAT*/ + return result; +} + + +static int cmp_sp_sp(const Gcalc_scan_iterator::point *a, + const Gcalc_scan_iterator::point *b, + const Gcalc_heap::Info *pi) +{ + if (a->event == scev_none && b->event == scev_none) + return cmp_sp_sp_cnt(a, b, &pi->iy); + if (a->event == scev_none) + return cmp_sp_pi(a, pi); + if (b->event == scev_none) + return -1 * cmp_sp_pi(b, pi); + + return 0; } @@ -366,7 +1018,7 @@ int Gcalc_scan_iterator::arrange_event() continue; if (!(new_sp= new_slice_point())) return 1; - *new_sp= *sp; + new_sp->copy_all(sp); *ae_hook= new_sp; ae_hook= &new_sp->next; #ifdef TO_REMOVE @@ -420,10 +1072,12 @@ int Gcalc_scan_iterator::insert_top_point() sp0->pi= m_cur_pi; sp0->next_pi= m_cur_pi->left; sp0->thread= m_cur_thread++; +#ifdef TMP_BLOCK sp0->x= coord_to_float(m_cur_pi->x); +#endif /*TMP_BLOCK*/ if (m_cur_pi->left) { - sp0->horiz_dir= GET_DX_DY(&sp0->dx_dy, m_cur_pi, m_cur_pi->left); + calc_dx_dy(sp0); sp0->event= scev_thread; /*Now just to increase the size of m_slice0 to be same*/ @@ -436,15 +1090,13 @@ int Gcalc_scan_iterator::insert_top_point() if (!(sp1= new_slice_point())) return 1; sp1->event= sp0->event= scev_two_threads; -#ifdef TO_REMOVE - sp1->event_pair= sp0; - sp0->event_pair= sp1; -#endif /*TO_REMOVE*/ sp1->pi= m_cur_pi; sp1->next_pi= m_cur_pi->right; sp1->thread= m_cur_thread++; +#ifdef TMP_BLOCK sp1->x= sp0->x; - sp1->horiz_dir= GET_DX_DY(&sp1->dx_dy, m_cur_pi, m_cur_pi->right); +#endif /*TMP_BLOCK*/ + calc_dx_dy(sp1); /* We have two threads so should decide which one will be first */ if (sp0->cmp_dx_dy(sp1)>0) { @@ -463,17 +1115,20 @@ int Gcalc_scan_iterator::insert_top_point() else { sp0->event= scev_single_point; +#ifdef TMP_BLOCK sp0->horiz_dir= 0; sp0->dx_dy= 0.0; +#endif /*TMP_BLOCK*/ } /* We need to find the place to insert. */ - for (; sp && (sp->x < sp0->x); prev_hook= &sp->next, sp=sp->get_next()) + for (; sp && cmp_sp_pi(sp, m_cur_pi) < 0; + prev_hook= &sp->next, sp=sp->get_next()) {} next_state->event_position_hook= prev_hook; - if (sp && coord_eq(sp->x, sp0->x)) + if (sp && cmp_sp_pi(sp, m_cur_pi) == 0) { next_state->event_position= sp; do @@ -482,7 +1137,7 @@ int Gcalc_scan_iterator::insert_top_point() sp->event= scev_intersection; prev_hook= &sp->next; sp= sp->get_next(); - } while (sp && coord_eq(sp->x, sp0->x)); + } while (sp && cmp_sp_pi(sp, m_cur_pi) == 0); } else next_state->event_position= sp0; @@ -503,44 +1158,88 @@ int Gcalc_scan_iterator::insert_top_point() } +#ifndef NO_TESTING +const char *pev(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"; + }; + return "fck"; +} +extern int ca_counter; +#endif /*NO_TESTING*/ int Gcalc_scan_iterator::normal_scan() { point *sp; Gcalc_dyn_list::Item **sp_hook; Gcalc_heap::Info *next_pi; - point *first_bottom_point= NULL; + point *first_bottom_point; if (m_next_is_top_point && insert_top_point()) return 1; for (next_pi= m_cur_pi->get_next(); - next_pi && - coord_eq(m_cur_pi->x, next_pi->x) && - coord_eq(m_cur_pi->y, next_pi->y); + next_pi && cmp_point_info(m_cur_pi, next_pi) == 0; next_pi= next_pi->get_next()) { - DBUG_ASSERT(coord_eq(next_state->event_position->x, next_pi->x)); next_state->clear_event_position(); m_next_is_top_point= true; + first_bottom_point= NULL; for (sp= next_state->slice, sp_hook= (Gcalc_dyn_list::Item **) &next_state->slice; sp; sp_hook= &sp->next, sp= sp->get_next()) { +#ifndef NO_TESTING + // if (ca_counter == 21) + printf("%s%d\t", pev(sp->event), sp->thread); +#endif /*NO_TESTING*/ if (sp->next_pi == next_pi) /* End of the segment */ { +#ifdef TMP_BLOCK sp->x= coord_to_float(next_pi->x); sp->pi= next_pi; +#endif /*TMP_BLOCK*/ + if (cmp_point_info(sp->pi, next_pi)) + sp->pi= next_pi; sp->next_pi= next_pi->left; m_next_is_top_point= false; if (next_pi->is_bottom()) { - if (first_bottom_point) + if (sp->event == scev_thread) + sp->event= scev_single_point; + else if (sp->event == scev_two_threads) + { + if (sp->get_next() && sp->get_next()->pi == sp->pi) + sp->get_next()->event= scev_thread; + else if (sp != next_state->slice) + { + point *fnd_sp; + for (fnd_sp= next_state->slice; fnd_sp->get_next() != sp; + fnd_sp= fnd_sp->get_next()); + DBUG_ASSERT(fnd_sp->pi == sp->pi); + fnd_sp->event= scev_thread; + } + sp->event= scev_single_point; + } + else if (first_bottom_point) { first_bottom_point->event= sp->event= scev_two_ends; -#ifdef TO_REMOVE - first_bottom_point->event_pair= sp; - sp->event_pair= first_bottom_point; -#endif /*TO_REMOVE*/ } else { @@ -550,21 +1249,52 @@ int Gcalc_scan_iterator::normal_scan() } else { - sp->event= scev_point; - sp->horiz_dir= GET_DX_DY(&sp->dx_dy, next_pi, next_pi->left); + if ((sp->event & (scev_point | scev_thread | scev_two_threads)) == 0) + sp->event= scev_point; + calc_dx_dy(sp); } mark_event_position1(sp, sp_hook); } - else if (coord_eq(sp->x, next_pi->x)) + else if (sp->event || (cmp_sp_pi(sp, next_pi) == 0)) { if (!sp->event) sp->event= scev_intersection; mark_event_position1(sp, sp_hook); } } +#ifndef NO_TESTING + //if (ca_counter == 21) + printf("\n"); +#endif /*NO_TESTING*/ m_cur_pi= next_pi; - if (m_next_is_top_point && insert_top_point()) - return 1; + if (m_next_is_top_point) + { + if (insert_top_point()) + return 1; + /* Set proper values to the event position */ + /* TODO: can be done faster */ + next_state->clear_event_position(); + if (next_state->slice->event) + mark_event_position1(next_state->slice, + (Gcalc_dyn_list::Item **) &next_state->slice); +#ifndef NO_TESTING + //if (ca_counter == 21) + printf("*%s%d\t", pev(next_state->slice->event), next_state->slice->thread); +#endif /*NO_TESTING*/ + for (sp= next_state->slice; sp->get_next(); sp= sp->get_next()) + { +#ifndef NO_TESTING + //if (ca_counter == 21) + printf("%s%d\t", pev(sp->get_next()->event), sp->get_next()->thread); +#endif /*NO_TESTING*/ + if (sp->get_next()->event) + mark_event_position1(sp->get_next(), &sp->next); + } +#ifndef NO_TESTING + //if (ca_counter == 21) + printf("\n"); +#endif /*NO_TESTING*/ + } } /* Swap current <-> next */ @@ -580,6 +1310,19 @@ int Gcalc_scan_iterator::normal_scan() point *sp0= current_state->slice; point *sp1= next_state->slice; point *prev_sp1= NULL; +#ifndef NO_TESTING + //if (ca_counter == 21) + { + point *sp= current_state->slice; + printf("After arrange\n"); + for (; sp; sp= sp->get_next()) + printf("%s%d\t", pev(sp->event), sp->thread); + printf("\nEvent\n"); + for (sp= m_events; sp; sp= sp->get_next()) + printf("%s%d\t", pev(sp->event), sp->thread); + printf("\n"); + } +#endif /*NO_TESTING*/ if (!(m_cur_pi= next_pi)) { @@ -592,8 +1335,12 @@ int Gcalc_scan_iterator::normal_scan() return 0; } + next_state->intersection_scan= 0; + next_state->pi= m_cur_pi; Gcalc_heap::Info *cur_pi= m_cur_pi; +#ifdef TMP_BLOCK next_state->y= coord_to_float(cur_pi->y); +#endif /*TMP_BLOCK*/ first_bottom_point= NULL; @@ -603,9 +1350,12 @@ int Gcalc_scan_iterator::normal_scan() for (; sp0; sp0= sp0->get_next()) { + DBUG_ASSERT(!sp0->is_bottom()); if (sp0->next_pi == cur_pi) /* End of the segment */ { +#ifdef TMP_BLOCK sp1->x= coord_to_float(cur_pi->x); +#endif /*TMP_BLOCK*/ sp1->pi= cur_pi; sp1->thread= sp0->thread; sp1->next_pi= cur_pi->left; @@ -625,16 +1375,12 @@ int Gcalc_scan_iterator::normal_scan() else { first_bottom_point->event= sp1->event= scev_two_ends; -#ifdef TO_REMOVE - sp1->event_pair= first_bottom_point; - first_bottom_point->event_pair= sp1; -#endif /*TO_REMOVE*/ } } else { sp1->event= scev_point; - sp1->horiz_dir= GET_DX_DY(&sp1->dx_dy, m_cur_pi, m_cur_pi->left); + calc_dx_dy(sp1); } mark_event_position1(sp1, prev_sp1 ? &prev_sp1->next : @@ -644,10 +1390,12 @@ int Gcalc_scan_iterator::normal_scan() { /* Cut current string with the height of the new point*/ sp1->copy_core(sp0); +#ifdef TMP_BLOCK sp1->x= sp1->horiz_dir ? coord_to_float(cur_pi->x) : (coord_to_float(sp1->pi->x) + (next_state->y-coord_to_float(sp1->pi->y)) * sp1->dx_dy); - if (coord_eq(sp1->x, cur_pi->x)) +#endif /*TMP_BLOCK*/ + if (cmp_sp_pi(sp1, cur_pi) == 0) { mark_event_position1(sp1, prev_sp1 ? &prev_sp1->next : @@ -659,7 +1407,7 @@ int Gcalc_scan_iterator::normal_scan() } intersections_found= intersections_found || - (prev_sp1 && prev_sp1->x > sp1->x); + (prev_sp1 && cmp_sp_sp(prev_sp1, sp1, cur_pi) > 0); prev_sp1= sp1; sp1= sp1->get_next(); @@ -674,6 +1422,16 @@ int Gcalc_scan_iterator::normal_scan() free_list(sp1); } +#ifndef NO_TESTING + //if (ca_counter == 21) + { + point *sp= next_state->slice; + printf("After slice\n"); + for (; sp; sp= sp->get_next()) + printf("%s%d\t", pev(sp->event), sp->thread); + printf("\n"); + } +#endif /*NO_TESTING*/ if (intersections_found) return handle_intersections(); @@ -681,16 +1439,20 @@ int Gcalc_scan_iterator::normal_scan() } -#define INTERSECTION_ZERO 0.000000000001 - +#ifndef NO_TESTING +int isc_counter= 0; +#endif /*NO_TESTING*/ int Gcalc_scan_iterator::add_intersection(int n_row, const point *a, const point *b, Gcalc_dyn_list::Item ***p_hook) { + const point *a0= a->intersection_link; + const point *b0= b->intersection_link; intersection *isc= new_intersection(); if (!isc) return 1; + m_n_intersections++; **p_hook= isc; *p_hook= &isc->next; @@ -698,34 +1460,17 @@ int Gcalc_scan_iterator::add_intersection(int n_row, isc->thread_a= a->thread; isc->thread_b= b->thread; - /* intersection_normal */ - const point *a0= a->intersection_link; - const point *b0= b->intersection_link; - DBUG_ASSERT(!a0->horiz_dir || !b0->horiz_dir); - - if (!a0->horiz_dir && !b0->horiz_dir) + isc->ii= m_heap->new_intersection(a0->pi, a0->next_pi, + b0->pi, b0->next_pi); +#ifndef NO_TESTING + //if (isc_counter == 40) { - double b0_x= a0->next_pi->x - a0->pi->x; - double b0_y= a0->next_pi->y - a0->pi->y; - double b1_x= b0->next_pi->x - b0->pi->x; - double b1_y= b0->next_pi->y - b0->pi->y; - double b1xb0= b1_x * b0_y - b1_y * b0_x; - double t= (a0->pi->x - b0->pi->x) * b0_y - (a0->pi->y - b0->pi->y) * b0_x; - if (fabs(b1xb0) < INTERSECTION_ZERO) - { - isc->y= current_state->y; - isc->x= a0->x; - return 0; - } - - t/= b1xb0; - isc->x= b0->pi->x + b1_x*t; - isc->y= b0->pi->y + b1_y*t; - return 0; + long double ix, iy; + isc->ii->calc_xy_ld(&ix, &iy); + printf("%d\t%d\t%.20LG\t%.20LG\t%d\n", isc->thread_a, isc->thread_b, ix, iy, isc->n_row); } - isc->y= next_state->y; - isc->x= a0->horiz_dir ? b->x : a->x; - return 0; +#endif /*NO_TESTING*/ + return isc->ii == NULL; } @@ -733,6 +1478,10 @@ int Gcalc_scan_iterator::find_intersections() { Gcalc_dyn_list::Item **hook; +#ifndef NO_TESTING + ++isc_counter; + printf("Looking for intersections\n"); +#endif /*NO_TESTING*/ m_n_intersections= 0; { /* Set links between slicepoints */ @@ -743,6 +1492,15 @@ int Gcalc_scan_iterator::find_intersections() DBUG_ASSERT(!sp0->is_bottom()); DBUG_ASSERT(sp0->thread == sp1->thread); sp1->intersection_link= sp0; +#ifndef NO_TESTING + //if (isc_counter == 40) + { + long double spx, spx1; + sp0->calc_x(&spx, current_state->pi->y, current_state->pi->x); + sp1->calc_x(&spx1, m_cur_pi->y, m_cur_pi->x); + printf("%d\t%.20LG\t%.20LG\n", sp0->thread, spx, spx1); + } +#endif /*NO_TESTING*/ } } @@ -761,7 +1519,7 @@ int Gcalc_scan_iterator::find_intersections() point *s1= prev_s1->get_next(); if (!s1) break; - if (prev_s1->x <= s1->x) + if (cmp_sp_sp(prev_s1, s1, m_cur_pi) <= 0) { pprev_s1= (point **) &prev_s1->next; continue; @@ -783,19 +1541,250 @@ int Gcalc_scan_iterator::find_intersections() } -static int compare_intersections(const void *e0, const void *e1) +class Gcalc_coord4 : public Gcalc_internal_coord +{ + coord_digit_t c[COORD_BASE*4]; + public: + void init() + { + n_digits= COORD_BASE*4; + digits= c; + } +}; + +class Gcalc_coord5 : public Gcalc_internal_coord +{ + coord_digit_t c[COORD_BASE*5]; + public: + void init() + { + n_digits= COORD_BASE*5; + digits= c; + } +}; + + +static void calc_isc_exp(Gcalc_coord5 *exp, + const Gcalc_coord2 *bb2, + const Gcalc_coord1 *ya1, + const Gcalc_coord2 *bb1, + const Gcalc_coord1 *yb1, + const Gcalc_coord2 *a21_b1) +{ + Gcalc_coord3 p1, p2, sum; + p1.init(); + p2.init(); + sum.init(); + exp->init(); + + gcalc_mul_coord(&p1, ya1, bb1); + gcalc_mul_coord(&p2, yb1, a21_b1); + gcalc_add_coord(&sum, &p1, &p2); + gcalc_mul_coord(exp, bb2, &sum); +} + + +#ifdef TMP_BLOCK +static void calc_aa1_b(Gcalc_coord2 *res, + const Gcalc_heap::Info *a0, + const Gcalc_heap::Info *a1, + const Gcalc_coord1 *xb, + const Gcalc_coord1 *yb) +{ + Gcalc_coord1 aa1_x, aa1_y; + Gcalc_coord2 p1, p2; + res->init(); + aa1_x.init(); + aa1_y.init(); + p1.init(); + p2.init(); + + gcalc_sub_coord(&aa1_x, &a1->ix, &a0->ix); + gcalc_sub_coord(&aa1_y, &a1->iy, &a0->iy); + gcalc_mul_coord(&p1, &aa1_x, yb); + gcalc_mul_coord(&p2, &aa1_y, xb); + gcalc_sub_coord(res, &p1, &p2); +} + + +static int cmp_intersections_y(const Gcalc_heap::Intersection_info *i1, + const Gcalc_heap::Intersection_info *i2) +{ + Gcalc_coord2 t_a1, t_b1; + Gcalc_coord2 t_a2, t_b2; + Gcalc_coord1 yb1, yb2; + Gcalc_coord1 xb1, xb2; + Gcalc_coord5 exp_a, exp_b; + + calc_t(&t_a1, &t_b1, &xb1, &yb1, i1); + calc_t(&t_a2, &t_b2, &xb2, &yb2, i2); + + calc_isc_exp(&exp_a, &t_b2, &i1->p1->iy, &t_b1, &yb1, &t_a1); + calc_isc_exp(&exp_b, &t_b1, &i2->p1->iy, &t_b2, &yb2, &t_a2); + + int result= gcalc_cmp_coord(&exp_a, &exp_b); +#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) + DBUG_ASSERT(de_check(y1, y2)); + if (result < 0) + DBUG_ASSERT(de_check1(y1, y2) || y1 < y2); + if (result > 0) + DBUG_ASSERT(de_check1(y1, y2) || y1 > y2); +#endif /*GCALC_CHECK_WITH_FLOAT*/ + return result; +} + + +static int compare_intersections(const void *e1, const void *e2) { - Gcalc_scan_iterator::intersection *i0= (Gcalc_scan_iterator::intersection *)e0; Gcalc_scan_iterator::intersection *i1= (Gcalc_scan_iterator::intersection *)e1; + Gcalc_scan_iterator::intersection *i2= (Gcalc_scan_iterator::intersection *)e2; + int result= cmp_intersections_y(i1->ii, i2->ii); + if (result != 0) + return result > 0; + return (i1->n_row > i2->n_row); +} + +#endif /*TMP_BLOCK*/ + +#ifndef NO_TESTING +extern int ca_counter; +#endif /*NO_TESTING*/ +static int cmp_intersections(const Gcalc_heap::Intersection_info *i1, + const Gcalc_heap::Intersection_info *i2) +{ + Gcalc_coord2 t_a1, t_b1; + Gcalc_coord2 t_a2, t_b2; + Gcalc_coord1 yb1, yb2; + Gcalc_coord1 xb1, xb2; + Gcalc_coord5 exp_a, exp_b; + int result; + + calc_t(&t_a1, &t_b1, &xb1, &yb1, i1); + calc_t(&t_a2, &t_b2, &xb2, &yb2, i2); + + calc_isc_exp(&exp_a, &t_b2, &i1->p1->iy, &t_b1, &yb1, &t_a1); + calc_isc_exp(&exp_b, &t_b1, &i2->p1->iy, &t_b2, &yb2, &t_a2); + + result= gcalc_cmp_coord(&exp_a, &exp_b); +#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) + DBUG_ASSERT(de_weak(y1, y2)); +#ifndef NO_TESTING + if (ca_counter == 7) + printf("77"); +#endif /*NO_TESTING*/ + if (result < 0) + DBUG_ASSERT(de_weak(y1, y2) || y1 < y2); + if (result > 0) + DBUG_ASSERT(de_weak(y1, y2) || y1 > y2); +#endif /*GCALC_CHECK_WITH_FLOAT*/ + + if (result != 0) + return result; + + + calc_isc_exp(&exp_a, &t_b2, &i1->p1->ix, &t_b1, &xb1, &t_a1); + calc_isc_exp(&exp_b, &t_b1, &i2->p1->ix, &t_b2, &xb2, &t_a2); + result= gcalc_cmp_coord(&exp_a, &exp_b); +#ifdef GCALC_CHECK_WITH_FLOAT + if (result == 0) + DBUG_ASSERT(de_weak(x1, x2)); + if (result < 0) + DBUG_ASSERT(de_weak(x1, x2) || x1 < x2); + if (result > 0) + DBUG_ASSERT(de_weak(x1, x2) || x1 > x2); +#endif /*GCALC_CHECK_WITH_FLOAT*/ + return result; +} - if (fabs(i0->y - i1->y) > 0.00000000000001) - return i0->y > i1->y; +static int compare_intersections(const void *e1, const void *e2) +{ + Gcalc_scan_iterator::intersection *i1= (Gcalc_scan_iterator::intersection *)e1; + Gcalc_scan_iterator::intersection *i2= (Gcalc_scan_iterator::intersection *)e2; + int result= cmp_intersections(i1->ii, i2->ii); + if (result != 0) + return result > 0; + return (i1->n_row > i2->n_row); +} - if (i0->n_row != i1->n_row) - return i0->n_row > i1->n_row; - if (!coord_eq(i0->y, i1->y)) - return i0->y > i1->y; - return i0->x > i1->x; +inline int intersections_eq(const Gcalc_heap::Intersection_info *i1, + const Gcalc_heap::Intersection_info *i2) +{ + return cmp_intersections(i1, i2) == 0; +} + + +#ifdef TMP_BLOCK +static void calc_isc_sp_exp(Gcalc_coord4 *exp, + const Gcalc_coord2 *bb, + const Gcalc_coord1 *x1, + const Gcalc_coord1 *y1, + const Gcalc_coord1 *x2, + const Gcalc_coord1 *y2) +{ + Gcalc_coord2 p1, p2, sum; + p1.init(); + p2.init(); + sum.init(); + exp->init(); + + gcalc_mul_coord(&p1, x1, y1); + gcalc_mul_coord(&p2, x2, y2); + gcalc_add_coord(&sum, &p1, &p2); + gcalc_mul_coord(exp, bb, &sum); +} +#endif /*TMP_BLOCK*/ + + +static int sp_isc_eq(const Gcalc_scan_iterator::point *sp, + const Gcalc_heap::Intersection_info *isc) +{ + + Gcalc_coord4 exp_a, exp_b; + Gcalc_coord1 xb1, yb1; + Gcalc_coord2 t_a, t_b; + Gcalc_coord2 t_sp_a, t_sp_b; +#ifdef TMP_BLOCK + xa1a0.init(); + ya1a0.init(); + + gcalc_sub_coord(&xa1a0, &isc->p1->ix, &sp->pi->ix); + gcalc_sub_coord(&ya1a0, &isc->p1->iy, &sp->pi->iy); + calc_isc_sp_exp(&exp_a, &t_b, &sp->pi->ix, &ya1a0, &sp->pi->iy, &xa1a0); + calc_isc_sp_exp(&exp_b, &t_a, &sp->pi->iy, &isc->p1->ix, + &sp->pi->ix, &isc->p1->iy); + return gcalc_cmp_coord(&exp_a, &exp_b); +#endif /*TMP_BLOCK*/ + exp_a.init(); + exp_b.init(); + calc_t(&t_a, &t_b, &xb1, &yb1, isc); + calc_t(&t_sp_a, &t_sp_b, &xb1, &yb1, isc->p1, isc->p2, sp->pi, sp->next_pi); + + gcalc_mul_coord(&exp_a, &t_a, &t_sp_b); + gcalc_mul_coord(&exp_b, &t_b, &t_sp_a); + int result= gcalc_cmp_coord(&exp_a, &exp_b); +#ifdef GCALC_CHECK_WITH_FLOAT + long double int_x, int_y, sp_x; + isc->calc_xy_ld(&int_x, &int_y); + sp->calc_x(&sp_x, int_y, int_x); + if (result == 0) + DBUG_ASSERT(de_check1(sp->dy.get_double(), 0.0) || + de_check(sp_x, int_x)); +#ifdef TMP_BLOCK + else + DBUG_ASSERT(!de_check1(sp_x, int_x)); +#endif /*TMP_BLOCK*/ +#endif /*GCALC_CHECK_WITH_FLOAT*/ + return result == 0; } @@ -823,6 +1812,18 @@ int Gcalc_scan_iterator::handle_intersections() /* We need the next slice to be just equal */ next_state->slice= new_slice(saved_state->slice); m_cur_intersection= m_intersections; +#ifndef NO_TESTING + //if (isc_counter == 40) + { + printf("Sorted\n"); + for (intersection *isc= m_intersections; isc; isc= isc->get_next()) + { + long double ix, iy; + isc->ii->calc_xy_ld(&ix, &iy); + printf("%d\t%d\t%.20LG\t%.20LG\t%d\n", isc->thread_a, isc->thread_b, ix, iy, isc->n_row); + } + } +#endif /*NO_TESTING*/ return intersection_scan(); } @@ -832,7 +1833,6 @@ int Gcalc_scan_iterator::intersection_scan() point *sp0, *sp1; Gcalc_dyn_list::Item **hook; intersection *next_intersection= NULL; - int met_equal= 0; if (m_cur_intersection != m_intersections) { @@ -846,31 +1846,72 @@ int Gcalc_scan_iterator::intersection_scan() if (arrange_event()) return 1; +#ifdef TMP_BLOCK if (!m_cur_intersection) { saved_state->event_position_hook= (Gcalc_dyn_list::Item **) &saved_state->slice; -#ifdef TO_REMOVE - for (sp0= current_state->slice, sp1= saved_state->slice; - sp0; - sp0= sp0->get_next(), sp1= sp1->get_next()) + + for (sp1= saved_state->slice; sp1; sp1= sp1->get_next()) { - sp0->next_link= sp1; if (sp1->get_next() == saved_state->event_position) saved_state->event_position_hook= &sp1->next; } -#endif /*TO_REMOVE*/ - for (sp1= saved_state->slice; sp1; sp1= sp1->get_next()) + /* Swap saved <-> next */ { - if (sp1->get_next() == saved_state->event_position) - saved_state->event_position_hook= &sp1->next; + slice_state *tmp= next_state; + next_state= saved_state; + saved_state= tmp; } + free_list(saved_state->slice); + saved_state->slice= NULL; + + free_list(m_intersections); + m_intersections= NULL; + return 0; + } +#endif /*TMP_BLOCK*/ + if (!m_cur_intersection) + { /* Swap saved <-> next */ { slice_state *tmp= next_state; next_state= saved_state; saved_state= tmp; } + Gcalc_dyn_list::Item **n_hook= + (Gcalc_dyn_list::Item **) &next_state->slice; + + next_state->clear_event_position(); + sp0= current_state->slice; + sp1= next_state->slice; + + for (; sp0; + sp0= sp0->get_next(), n_hook= &sp1->next, sp1= sp1->get_next()) + { + if (sp0->thread != sp1->thread) + { + point *fnd_s= sp1->get_next(); + Gcalc_dyn_list::Item **fnd_hook= &sp1->next; + for (; fnd_s && fnd_s->thread != sp0->thread; + fnd_hook= &fnd_s->next, fnd_s= fnd_s->get_next()); + DBUG_ASSERT(fnd_s && fnd_s == *fnd_hook); + /* Now swap items of the next_state->slice */ + *n_hook= fnd_s; + *fnd_hook= fnd_s->next; + fnd_s->next= sp1; + sp1= fnd_s; + } + if (sp1->event) + mark_event_position1(sp1, n_hook); + } +#ifndef DBUG_OFF + sp0= current_state->slice; + sp1= next_state->slice; + for (; sp0; sp0= sp0->get_next(), sp1= sp1->get_next()) + DBUG_ASSERT(sp0->thread == sp1->thread); + DBUG_ASSERT(!sp1); +#endif /*DBUG_OFF*/ free_list(saved_state->slice); saved_state->slice= NULL; @@ -880,12 +1921,12 @@ int Gcalc_scan_iterator::intersection_scan() } } - next_state->y= m_cur_intersection->y; - sp0= current_state->slice; hook= (Gcalc_dyn_list::Item **) &next_state->slice; sp1= next_state->slice; next_state->clear_event_position(); + next_state->intersection_scan= 1; + next_state->isc= m_cur_intersection->ii; for (; sp0; hook= &sp1->next, sp1= sp1->get_next(), sp0= sp0->get_next()) @@ -893,27 +1934,20 @@ int Gcalc_scan_iterator::intersection_scan() if (sp0->thread == m_cur_intersection->thread_a || sp0->thread == m_cur_intersection->thread_b) { -#ifdef REACTIVATE_THIS DBUG_ASSERT(sp0->thread != m_cur_intersection->thread_a || - sp0->get_next()->thread == m_cur_intersection->thread_b); -#endif /*REACTIVATE_THIS*/ + sp0->get_next()->thread == m_cur_intersection->thread_b || + sp_isc_eq(sp0->get_next(), m_cur_intersection->ii)); sp1->copy_core(sp0); - sp1->x= m_cur_intersection->x; sp1->event= scev_intersection; mark_event_position1(sp1, hook); } else { sp1->copy_core(sp0); - sp1->x= sp1->horiz_dir ? - m_cur_intersection->x : - coord_to_float(sp1->pi->x) + - (next_state->y-coord_to_float(sp1->pi->y)) * sp1->dx_dy; - if (coord_eq(sp1->x, m_cur_intersection->x)) + if (sp_isc_eq(sp1, m_cur_intersection->ii)) { sp1->event= scev_intersection; mark_event_position1(sp1, hook); - met_equal= 1; } else sp1->event= scev_none; @@ -926,45 +1960,10 @@ int Gcalc_scan_iterator::intersection_scan() *hook= NULL; } - if (met_equal) - { - /* Remove superfluous intersections. */ - /* Double operations can produce unexact result, so it's needed. */ - for (sp0= next_state->event_position; - sp0 != *next_state->event_end_hook; - sp0= sp0->get_next()) - { - for (sp1= sp0->get_next(); - sp1 != *next_state->event_end_hook; - sp1= sp1->get_next()) - { - intersection *isc= m_cur_intersection; - while (isc->get_next()) - { - intersection *cur_isc= isc->get_next(); - if ((cur_isc->thread_a == sp0->thread && - cur_isc->thread_b == sp1->thread) || - (cur_isc->thread_a == sp1->thread && - cur_isc->thread_b == sp0->thread)) - { - /* The superfluous intersection should be close to the current. */ - DBUG_ASSERT(fabs(cur_isc->x-m_cur_intersection->x) + - fabs(cur_isc->y-m_cur_intersection->y) < - 0.00000000001); - isc->next= isc->next->next; - free_item(cur_isc); - } - else - isc= isc->get_next(); - } - } - } - } - + /* Now check equal intersections */ for (next_intersection= m_cur_intersection->get_next(); next_intersection && - coord_eq(next_intersection->x, m_cur_intersection->x) && - coord_eq(next_intersection->y, m_cur_intersection->y); + intersections_eq(next_intersection->ii, m_cur_intersection->ii); next_intersection= next_intersection->get_next()) { /* Handle equal intersections. We only need to set proper events */ @@ -990,5 +1989,58 @@ int Gcalc_scan_iterator::intersection_scan() return 0; } + +double Gcalc_scan_iterator::get_y() const +{ + if (current_state->intersection_scan) + { + double x, y; + current_state->isc->calc_xy(&x, &y); + return y; + } + else + return current_state->pi->y; +} + + +double Gcalc_scan_iterator::get_event_x() const +{ + if (current_state->intersection_scan) + { + double x, y; + current_state->isc->calc_xy(&x, &y); + return x; + } + else + return current_state->pi->x; +} + +double Gcalc_scan_iterator::get_h() const +{ + double cur_y= get_y(); + double next_y; + if (next_state->intersection_scan) + { + double x; + next_state->isc->calc_xy(&x, &next_y); + } + else + next_y= next_state->pi->y; + 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->x; + dy= sp->next_pi->y - sp->pi->y; + if (fabs(dy) < 1e-12) + return sp->pi->x; + return (sp->next_pi->x - sp->pi->x) * dy; +} + + #endif /* HAVE_SPATIAL */ |