/*******************************************************************************
+
+  LEDA 3.5
+
+  _rat_point.c
+
+  This file is part of the LEDA research version (LEDA-R) that can be 
+  used free of charge in academic research and teaching. Any commercial
+  use of this software requires a license which is distributed by the
+  LEDA Software GmbH, Postfach 151101, 66041 Saarbruecken, FRG
+  (fax +49 681 31104).
+
+  Copyright (c) 1991-1997  by  Max-Planck-Institut fuer Informatik
+  Im Stadtwald, 66123 Saarbruecken, Germany     
+  All rights reserved.
+ 
*******************************************************************************/
//------------------------------------------------------------------------------
// rat_points :  points with rational (homogeneous) coordinates
//
// by S. Naeher (1995)
//------------------------------------------------------------------------------


#include <math.h>
#include <ctype.h>
#include <LEDA/rat_point.h>

// the fabs function is used very often therefore we provide
// a fast version for sparc machines (should work on all big-endian 
// architectures) that simply clears the sign bit to zero  
//
// FABS(x) clears the sign bit of (double) floating point number x

#if defined(sparc)
#define FABS(x) (*(unsigned long*)&x) &= 0x7FFFFFFF
#else
#define FABS(x) x=fabs(x)
#endif


leda_mutex rat_point_rep::mutex_id_counter;

unsigned long rat_point_rep::id_counter = 0;

// constructors

   rat_point_rep::rat_point_rep() : x(0), y(0), w(1)
   { xd = 0;
     yd = 0;
     wd = 1;
     mutex_id_counter.lock();
     id = id_counter++;
     mutex_id_counter.unlock();
    }

   rat_point_rep::rat_point_rep(const rational& a, const rational& b) 
       : x(a.numerator()*b.denominator()),
         y(b.numerator()*a.denominator()),
         w(a.denominator()*b.denominator())
   { xd = x.todouble(); 
     yd = y.todouble(); 
     wd = w.todouble(); 
     mutex_id_counter.lock();
     id = id_counter++;
     mutex_id_counter.unlock();
    }

   rat_point_rep::rat_point_rep(integer a, integer b) : x(a), y(b), w(1)
   { xd = a.todouble(); 
     yd = b.todouble(); 
     wd = 1;
     mutex_id_counter.lock();
     id = id_counter++;
     mutex_id_counter.unlock();
    }

   rat_point_rep::rat_point_rep(integer a, integer b, integer c) 
                                 : x(a), y(b), w(c)
   { if (c < 0)
     { x = -a;
       y = -b;
       w = -c;
      }
     xd = x.todouble(); 
     yd = y.todouble(); 
     wd = w.todouble();
     mutex_id_counter.lock();
     id = id_counter++;
     mutex_id_counter.unlock();
    }



// static members used for statistics
// (MOVED TO _pglobal.c)
//

/*
int rat_point::orient_count= 0;
int rat_point::exact_orient_count = 0;

int rat_point::cmp_count= 0;
int rat_point::exact_cmp_count = 0;

int rat_point::soc_count= 0;
int rat_point::exact_soc_count = 0;

// use_filter flag controls whether the floating point filter is used 

int rat_point::use_filter = 1;

*/




point rat_point::to_point() const
{ return point(xcoordD(),ycoordD()); }


rat_point rat_point::rotate90(const rat_point& p) const
{ integer x0 = p.X();
  integer y0 = p.Y();
  integer w0 = p.W();
  integer x1 = X();
  integer y1 = Y();
  integer w1 = W();
  integer x = (x0 + y0) * w1 - y1 * w0;
  integer y = (y0 - x0) * w1 + x1 * w0;
  return rat_point(x,y,w0*w1);
 }


rat_point rat_point::rotate90() const
{ return rat_point(-Y(),X(),W()); }



rat_point rat_point::reflect(const rat_point& q) const
{ // reflect point across point q

  integer x1 = X();
  integer y1 = Y();
  integer w1 = W();
  integer x2 = 2*q.X();
  integer y2 = 2*q.Y();
  integer w2 = q.W();

  return rat_point(x2*w1-x1*w2 , y2*w2-y1*w1, w1*w2);
}




rat_point rat_point::reflect(const rat_point& p, const rat_point& q) const
{  // reflect point across line through p and q

  integer dx = p.X();
  integer dy = p.Y();
  integer dw = p.W();

  rat_point p0 = translate(-dx,-dy,dw);
  rat_point q0 = q.translate(-dx,-dy,dw);

  integer x1 = p0.X();
  integer y1 = p0.Y();
  integer w1 = p0.W();

  integer x2 = q0.X();
  integer y2 = q0.Y();

  integer cosfi = (x1*x2 + y1*y2);
  integer sinfi = (x1*y2 - x2*y1); 

  integer cos2 = (cosfi*cosfi - sinfi*sinfi);
  integer sin2 = 2*cosfi*sinfi;
  integer W    = w1 * (x1*x1 + y1*y1) * (x2*x2 + y2*y2);

  rat_point r0(x1*cos2-y1*sin2,x1*sin2+y1*cos2,W);

  return r0.translate(dx,dy,dw);
}
  

rat_point rat_point::translate(integer dx, integer dy, integer dw) const
{ integer x1 = X();
  integer y1 = Y();
  integer w1 = W();
  integer x = dx * w1 + x1 * dw;
  integer y = dy * w1 + y1 * dw;
  return rat_point(x,y,dw*w1);
 }

rat_point rat_point::translate(const rational& dx, const rational& dy) const
{ integer x = dx.numerator() * dy.denominator();
  integer y = dy.numerator() * dx.denominator();
  integer w = dx.denominator() * dy.denominator();
  return translate(x,y,w);
}


rational rat_point::sqr_dist(const rat_point& p) const
{ 
  integer x1 = X();
  integer y1 = Y();
  integer w1 = W();
  integer x2 = p.X();
  integer y2 = p.Y();
  integer w2 = p.W();

  integer dx = x1*w2 - x2*w1;
  integer dy = y1*w2 - y2*w1;

  return rational(dx*dx+dy*dy,w1*w1*w2*w2);
}

rational rat_point::xdist(const rat_point& p) const
{ integer x1  = X();
  integer w1  = W();
  integer x2  = p.X();
  integer w2  = p.W();
  return rational(abs(x1*w2-x2*w1), w1*w2);
}

rational rat_point::ydist(const rat_point& p) const
{ integer y1  = Y();
  integer w1  = W();
  integer y2  = p.Y();
  integer w2  = p.W();
  return rational(abs(y1*w2-y2*w1), w1*w2);
}



// basic stream I/O operations

ostream& operator<<(ostream& out, const rat_point& p)
{ out << "(" << p.X() << "," << p.Y() << "," << p.W() << ")";
  return out;
 } 

istream& operator>>(istream& in, rat_point& p) 
{ // syntax: {(} x {,} y {,} w {)}   

  integer x,y,w;
  char c;

  do in.get(c); while (in && isspace(c));

  if (!in) return in;

  if (c != '(') in.putback(c);

  in >> x;

  do in.get(c); while (isspace(c));
  if (c != ',') in.putback(c);

  in >> y; 

  do in.get(c); while (isspace(c));
  if (c != ',') in.putback(c);

  in >> w; 

  do in.get(c); while (c == ' ');
  if (c != ')') in.putback(c);

  p = rat_point(x,y,w ); 
  return in; 

 } 



// machine precision

const double point_eps0 = ldexp(1.0,-53);


// cmp - primitive

int rat_point::cmp_xy(const rat_point& a, const rat_point& b)
{ 
  rat_point::cmp_count++;

  if (use_filter != 0)
  { 
    double axbw = a.XD()*b.WD();
    double bxaw = b.XD()*a.WD();
    double E    = axbw - bxaw;       // floating point result

    //----------------------------------------------------------------
    // ERROR BOUND:
    //----------------------------------------------------------------
    //
    // mes(E) = mes(ax*bw - bx*aw)
    //        = 2*(mes(ax)*mes(bw) + mes(bx)*mes(aw))
    //        = 2*(fabs(ax)*fabs(bw) + fabs(bx)*fabs(aw)) 
    //        = 2*(fabs(axbw) + fabs(bxaw)) 
    //
    // ind(E) = ind(ax*bw - bx*aw)
    //        = (ind(ax*bw) + ind(bx*aw) + 1)/2
    //        = ((ind(ax)+ind(bw)+0.5) + (ind(bx)+ind(aw)+0.5) + 1)/2
    //        = (1.5 + 1.5 + 1)/2
    //        = 2
    //
    // eps(E) = ind(E) * mes(E) * eps0 
    //        = 4 * (fabs(axbw) + fabs(bxaw)) * eps0
    //----------------------------------------------------------------


    FABS(axbw);
    FABS(bxaw);
    double eps = 4 * (axbw+bxaw) * point_eps0;

    if (E > +eps) return +1;
    if (E < -eps) return -1;
   
    if (eps < 1)  // compare y-coordinates
    { double aybw = a.YD()*b.WD();
      double byaw = b.YD()*a.WD();
      double E    = aybw - byaw;
      FABS(aybw);
      FABS(byaw);
      double eps = 4 * (aybw+byaw) * point_eps0;
      if (E > +eps) return +1;
      if (E < -eps) return -1;
      if (eps < 1)  return  0; 
     }
   }
  
    //use big integer arithmetic

    rat_point::exact_cmp_count++;

    integer axbw = a.X()*b.W();
    integer bxaw = b.X()*a.W();
    if (axbw > bxaw) return  1;
    if (axbw < bxaw) return -1;

    integer aybw = a.Y()*b.W();
    integer byaw = b.Y()*a.W();
    if (aybw > byaw) return  1;
    if (aybw < byaw) return -1;
    return 0;
}


int rat_point::cmp_yx(const rat_point& a, const rat_point& b)
{ 
  rat_point::cmp_count++;

  if (use_filter != 0)
  { 
    double aybw = a.YD()*b.WD();
    double byaw = b.YD()*a.WD();
    double E    = aybw - byaw;       // floating point result

    FABS(aybw);
    FABS(byaw);
    double eps = 4 * (aybw+byaw) * point_eps0;

    if (E > +eps) return +1;
    if (E < -eps) return -1;
   
    if (eps < 1)  // compare x-coordinates
    { double axbw = a.XD()*b.WD();
      double bxaw = b.XD()*a.WD();
      double E    = axbw - bxaw;
      FABS(axbw);
      FABS(bxaw);
      double eps = 4 * (axbw+bxaw) * point_eps0;
      if (E > +eps) return +1;
      if (E < -eps) return -1;
      if (eps < 1)  return  0; 
     }
   }
  
    //use big integer arithmetic

    rat_point::exact_cmp_count++;

    integer aybw = a.Y()*b.W();
    integer byaw = b.Y()*a.W();
    if (aybw > byaw) return  1;
    if (aybw < byaw) return -1;

    integer axbw = a.X()*b.W();
    integer bxaw = b.X()*a.W();
    if (axbw > bxaw) return  1;
    if (axbw < bxaw) return -1;
    return 0;
}






rational area(const rat_point& a, const rat_point& b, const rat_point& c)
{ return ((a.xcoord()-b.xcoord()) * (a.ycoord()-c.ycoord()) -
          (a.ycoord()-b.ycoord()) * (a.xcoord()-c.xcoord()))/2; }


int orientation(const rat_point& a, 
                const rat_point& b, 
                const rat_point& c)
{  
    rat_point::orient_count++;

    if (rat_point::use_filter != 0)
    { 
      double ax =  a.ptr()->xd; 
      double bx =  b.ptr()->xd; 
      double cx =  c.ptr()->xd; 
      double ay =  a.ptr()->yd;
      double by =  b.ptr()->yd;
      double cy =  c.ptr()->yd;
      double aw =  a.ptr()->wd;
      double bw =  b.ptr()->wd;
      double cw =  c.ptr()->wd;
 
      double aybw = ay*bw;
      double byaw = by*aw;
      double axcw = ax*cw;
      double cxaw = cx*aw;
      double axbw = ax*bw;
      double bxaw = bx*aw;
      double aycw = ay*cw;
      double cyaw = cy*aw;
 
      double E = (axbw-bxaw) * (aycw-cyaw) - (aybw-byaw) * (axcw-cxaw);
 
  //---------------------------------------------------------------------------
  // ERROR BOUNDS
  //---------------------------------------------------------------------------
  // mes(E) = 2*(mes(aybw-byaw)*mes(axcw-cxaw) + mes(axbw-bxaw)*mes(aycw-cyaw))
  //        = 2*(4*(fabs(aybw)+fabs(byaw)) * (fabs(axcw)+fabs(cxaw)) +
  //             4*(fabs(axbw)+fabs(bxaw)) * (fabs(aycw)+fabs(cyaw)))
  //        = 8*((fabs(aybw)+fabs(byaw)) * (fabs(axcw)+fabs(cxaw)) + 
  //             (fabs(axbw)+fabs(bxaw)) * (fabs(aycw)+fabs(cyaw)))
  //
  // ind(E) = ((ind(aybw-byaw) + ind(axcw-cxaw) +0.5) +
  //           (ind(axbw-bxaw) + ind(aycw-cyaw) +0.5) + 1 ) / 2
  //        = (4.5 + 4.5 + 1) / 2  =  5
  //
  // eps(E) = ind(E) * mes(E) * eps0
  //        = 40 * ((fabs(aybw)+fabs(byaw))*(fabs(axcw)-fabs(cxaw)) +
  //               (fabs(axbw)-fabs(bxaw))*(fabs(aycw)-fabs(cyaw))) * eps0;
  //---------------------------------------------------------------------------
 
      FABS(aybw);
      FABS(byaw);
      FABS(axcw);
      FABS(cxaw);
      FABS(axbw);
      FABS(bxaw);
      FABS(aycw);
      FABS(cyaw);
 
      double eps = 40*((aybw+byaw)*(axcw+cxaw)+(axbw+bxaw)*(aycw+cyaw))*point_eps0;
   
      if (E > eps)  return  1;
      if (E < -eps) return -1;
      if (eps < 1)  return  0;
     }
 
     // big integer arithmetic

     rat_point::exact_orient_count++;

     integer AX = a.X(); integer AY = a.Y(); integer AW = a.W();
     integer BX = b.X(); integer BY = b.Y(); integer BW = b.W();
     integer CX = c.X(); integer CY = c.Y(); integer CW = c.W();
     integer D = (AX*BW-BX*AW) * (AY*CW-CY*AW) - (AY*BW-BY*AW) * (AX*CW-CX*AW);
     return sign(D);
}


int side_of_circle(const rat_point& a, 
                   const rat_point& b, 
                   const rat_point& c, 
                   const rat_point& d)
{ 
  rat_point::soc_count++;

  if (rat_point::use_filter != 0)
    {
      double bxaw = b.XD() * a.WD();   
      double byaw = b.YD() * a.WD();
      double axbw = a.XD() * b.WD(); 
      double aybw = a.YD() * b.WD();    

      double bax = bxaw - axbw;
      double bay = byaw - aybw;   
      double baw = b.WD()*a.WD();

      double cxaw = c.XD() * a.WD();   
      double cyaw = c.YD() * a.WD();
      double axcw = a.XD() * c.WD(); 
      double aycw = a.YD() * c.WD(); 
     
      double cax = cxaw - axcw;
      double cay = cyaw - aycw;
      double caw = c.WD()*a.WD();

      double dxaw = d.XD() * a.WD();   
      double dyaw = d.YD() * a.WD();
      double axdw = a.XD() * d.WD(); 
      double aydw = a.YD() * d.WD(); 
     
      double dax = dxaw - axdw;
      double day = dyaw - aydw;
      double daw = d.WD()*a.WD();
      
      double Bxi = bax * baw;
      double Byi = bay * baw;
      double Bwi = bax*bax + bay*bay;

      double Cxi = cax * caw;
      double Cyi = cay * caw;
      double Cwi = cax*cax + cay*cay;

      double Dxi = dax * daw;
      double Dyi = day * daw;
      double Dwi = dax*dax + day*day;

      double E = (Byi*Cxi - Bxi*Cyi)*Dwi + 
        (Bwi*Cyi - Byi*Cwi)*Dxi + (Bxi*Cwi - Bwi*Cxi)*Dyi;
      
      /* ----------------------------------------------------------------
       * ERROR BOUND:
       * ----------------------------------------------------------------
       *
       * mes(E) = mes((Byi*Cxi - Bxi*Cyi)*Dwi + (Bwi*Cyi - Byi*Cwi)*Dxi + 
       *              (Bxi*Cwi - Bwi*Cxi)*Dyi)
       *        = 2*( mes((Byi*Cxi - Bxi*Cyi)*Dwi) + 
       *              2*mes((Bwi*Cyi - Byi*Cwi)*Dxi) +
       *              2*mes((Bxi*Cwi - Bwi*Cxi)*Dyi) )
       *        = 4*( (mes(Byi)*mes(Cxi) + mes(Bxi)*mes(Cyi))*mes(Dwi) + 
       *              2*(mes(Bwi)*mes(Cyi) + mes(Byi)*mes(Cwi))*mes(Dxi) +
       *              2*(mes(Bxi)*mes(Cwi) + mes(Bwi)*mes(Cxi))*mes(Dyi) )
       *
       * mes(Byi) = mes(bay) * mes(baw) 
       *          = mes(byaw - aybw)* mes(baw)
       *          = 2*(mes(byaw) + mes(aybw))* mes(Bw)*mes(Aw)
       *          = 2*(fabs(byaw) + fabs(aybw))*baw
       *            ( because a.WD() and b.WD() are positive )
       *
       * mes(Bxi) = 2*(fabs(bxaw) + fabs(axbw))*baw
       *
       * mes(Bwi) = 2*(mes(bax)*mes(bax) + mes(bay)*mes(bay))
       *          = 2*( 2*(fabs(bxaw) + fabs(axbw)) * 
       *                2*(fabs(bxaw) + fabs(axbw)) +
       *                2*(fabs(byaw) + fabs(aybw)) * 
       *                2*(fabs(byaw) + fabs(aybw)))
       *          = 8*((fabs(bxaw)+fabs(axbw))*(fabs(bxaw)+fabs(axbw)) +
       *               (fabs(byaw)+fabs(aybw))*(fabs(byaw)+fabs(aybw)))
       * 
       * mes(Cyi) = 2*(fabs(cyaw) + fabs(aycw))*caw
       * 
       * mes(Cxi) = 2*(fabs(cxaw) + fabs(axcw))*caw
       *
       * mes(Cwi) = 8*((fabs(cxaw)+fabs(axcw))*(fabs(cxaw)+fabs(axcw)) +
       *               (fabs(cyaw)+fabs(aycw))*(fabs(cyaw)+fabs(aycw)))
       *
       * mes(Dyi) = 2*(fabs(dyaw) + fabs(aydw))*daw
       * 
       * mes(Dxi) = 2*(fabs(dxaw) + fabs(axdw))*daw
       *
       * mes(Dwi) = 8*((fabs(dxaw)+fabs(axdw))*(fabs(dxaw)+fabs(axdw)) +
       *               (fabs(dyaw)+fabs(aydw))*(fabs(dyaw)+fabs(aydw)))
       *
       *
       * mes(E) = 4*(  ( 2*(fabs(byaw) + fabs(aybw))*baw *
       *                 2*(fabs(cxaw) + fabs(axcw))*caw + 
       *                 2*(fabs(bxaw) + fabs(axbw))*baw *
       *                 2*(fabs(cyaw) + fabs(aycw))*caw ) * 
       *               8*((fabs(dxaw)+fabs(axdw))*(fabs(dxaw)+fabs(axdw)) +
       *               (fabs(dyaw)+fabs(aydw))*(fabs(dyaw)+fabs(aydw)))
       *             + 
       *               ( 8*((fabs(bxaw)+fabs(axbw))*(fabs(bxaw)+fabs(axbw))+
       *                 (fabs(byaw)+fabs(aybw))*(fabs(byaw)+fabs(aybw))) *
       *                 2*(fabs(cyaw) + fabs(aycw))*caw + 
       *                 2*(fabs(byaw) + fabs(aybw))*baw *
       *                 8*((fabs(cxaw)+fabs(axcw))*(fabs(cxaw)+fabs(axcw))+
       *                 (fabs(cyaw)+fabs(aycw))*(fabs(cyaw)+fabs(aycw))) )*
       *               4*(fabs(dxaw) + fabs(axdw))*daw 
       *             +
       *               ( 2*(fabs(bxaw) + fabs(axbw))*baw *
       *                 8*((fabs(cxaw)+fabs(axcw))*(fabs(cxaw)+fabs(axcw))+
       *                 (fabs(cyaw)+fabs(aycw))*(fabs(cyaw)+fabs(aycw))) + 
       *                 8*((fabs(bxaw)+fabs(axbw))*(fabs(bxaw)+fabs(axbw))+
       *                 (fabs(byaw)+fabs(aybw))*(fabs(byaw)+fabs(aybw))) * 
       *                 2*(fabs(cxaw) + fabs(axcw))*caw ) *
       *               4*(fabs(dyaw) + fabs(aydw))*daw )
       *
       *
       * ind(E) = ( ind((Byi*Cxi-Bxi*Cyi)*Dwi) +
       *            ind((Bwi*Cyi-Byi*Cwi)*Dxi + (Bxi*Cwi-Bwi*Cxi)*Dyi) +
       *            1 ) /2
       *        = ( ind((Byi*Cxi-Bxi*Cyi)*Dwi) +
       *            (ind((Bwi*Cyi-Byi*Cwi)*Dxi)+
       *             ind((Bxi*Cwi-Bwi*Cxi)*Dyi)+1)/2
       *            + 1 ) / 2
       *        = ( (ind(Byi)+ind(Cxi)+0.5+ind(Bxi)+ind(Cyi)+0.5+1)/2 
       *            + ind(Dwi)+0.5+
       *            ((ind(Bwi)+ind(Cyi)+0.5+ind(Byi)+ind(Cwi)+0.5+1)/2
       *             +ind(Dxi)+0.5+(ind(Bxi)+ind(Cwi)+0.5+ind(Bwi)+ind(Cxi)
       *             +0.5+1)/2+ind(Dyi)+0.5+1)/2 +1)/2
       *
       * ind(Byi) = ind(bay)+ind(baw)+0.5
       *          = ind(bxaw-axbw)+ind(baw)+0.5
       *          = (1.5 + 1.5 +1)/2 + 1.5 +0.5
       *          = 4
       *          = ind(Bxi) = ind(Cxi) = ind(Cyi) = ind(Dxi) = ind(Dyi)
       *
       * ind(Dwi) = ind(dax*dax + day*day)
       *          = (ind(dax) + ind(dax) +0.5 + ind(day) + ind(day)+0.5+1)/2
       *          = ind(dax) + ind(day) + 1
       *          = (1.5 + 1.5 +1)/2 + (1.5 + 1.5 +1)/2 + 1
       *          = 5 
       *          = ind(Bwi) = ind(Cwi)
       *
       * ind(E) = ( 4 + 4 + 0.5 + 4 + 4 + 0.5 + 1 )/2 + 5 + 0.5 +
       *            ((5 + 4 + 0.5 + 4 + 5 + 0.5 + 1)/2 + 4 + 0.5 +
       *             (4 + 5 + 0.5 + 5 + 4 + 0.5 + 1)/2 + 4 + 0.5 + 1)/2 +
       *            1 ) /2
       *        = (9 + 5 + 0.5 + (10 + 4 + 0.5 + 10 + 4 + 0.5 + 1)/2 + 1)/2
       *        = (9 + 5 + 0.5 + 15 + 1)/2 = 30.5/2 = 61/4
       *
       * eps(E) = ind(E) * mes(E) * eps0
       *        = 32*61*
       *            (  ( (fabs(byaw) + fabs(aybw))*baw *
       *                 (fabs(cxaw) + fabs(axcw))*caw + 
       *                 (fabs(bxaw) + fabs(axbw))*baw *
       *                 (fabs(cyaw) + fabs(aycw))*caw ) * 
       *               ((fabs(dxaw)+fabs(axdw))*(fabs(dxaw)+fabs(axdw)) +
       *                (fabs(dyaw)+fabs(aydw))*(fabs(dyaw)+fabs(aydw)))
       *             + 
       *               ( ((fabs(bxaw)+fabs(axbw))*(fabs(bxaw)+fabs(axbw))+
       *                  (fabs(byaw)+fabs(aybw))*(fabs(byaw)+fabs(aybw))) *
       *                 (fabs(cyaw) + fabs(aycw))*caw + 
       *                 (fabs(byaw) + fabs(aybw))*baw *
       *                 ((fabs(cxaw)+fabs(axcw))*(fabs(cxaw)+fabs(axcw))+
       *                  (fabs(cyaw)+fabs(aycw))*(fabs(cyaw)+fabs(aycw))))*
       *               (fabs(dxaw) + fabs(axdw))*daw*2
       *             +
       *               ( (fabs(bxaw) + fabs(axbw))*baw *
       *                 ((fabs(cxaw)+fabs(axcw))*(fabs(cxaw)+fabs(axcw))+
       *                 (fabs(cyaw)+fabs(aycw))*(fabs(cyaw)+fabs(aycw))) + 
       *                 ((fabs(bxaw)+fabs(axbw))*(fabs(bxaw)+fabs(axbw))+
       *                  (fabs(byaw)+fabs(aybw))*(fabs(byaw)+fabs(aybw)))* 
       *                 (fabs(cxaw) + fabs(axcw))*caw ) *
       *               (fabs(dyaw) + fabs(aydw))*daw*2 )
       *           * eps0
       *             
       * -----------------------------------------------------------------*/
      
      FABS(bxaw); FABS(axbw); FABS(byaw); FABS(aybw);
      FABS(cxaw); FABS(axcw); FABS(cyaw); FABS(aycw);
      FABS(dxaw); FABS(axdw); FABS(dyaw); FABS(aydw);

      double eps = 1952 *
        (((byaw + aybw)*baw *(cxaw + axcw)*caw + 
          (bxaw + axbw)*baw *(cyaw + aycw)*caw ) * 
         ((dxaw + axdw)*(dxaw + axdw) + (dyaw + aydw)*(dyaw +aydw))
         + 
         (((bxaw + axbw)*(bxaw + axbw) + (byaw + aybw)*(byaw + aybw)) *
          (cyaw +aycw)*caw   +  (byaw + aybw)*baw *
          ((cxaw + axcw)*(cxaw + axcw) + (cyaw + aycw)*(cyaw + aycw))) *
         (dxaw + axdw)*daw*2
         +
         ((bxaw + axbw)*baw *
          ((cxaw + axcw)*(cxaw + axcw) + (cyaw + aycw)*(cyaw + aycw)) + 
          ((bxaw + axbw)*(bxaw + axbw) + (byaw + aybw)*(byaw + aybw)) * 
          (cxaw + axcw)*caw ) *
         (dyaw + aydw)*daw*2 )
        * point_eps0;


      if (E   >  eps)  return  1;
      if (E   < -eps)  return -1;
      if (eps <  1)    return  0;
    }

  rat_point::exact_soc_count++;

  integer AX = a.X();
  integer AY = a.Y();
  integer AW = a.W();

  integer BX = b.X();
  integer BY = b.Y();
  integer BW = b.W();

  integer CX = c.X();
  integer CY = c.Y();
  integer CW = c.W();

  integer DX = d.X();
  integer DY = d.Y();
  integer DW = d.W();

  integer bx,by,bw,cx,cy,cw,dx,dy,dw;

  if (AW==1 && BW==1 && CW==1 && DW==1)
  { bx = BX - AX;
    by = BY - AY;
    bw = bx*bx + by*by;

    cx = CX - AX;
    cy = CY - AY;
    cw = cx*cx + cy*cy;
    
    dx = DX - AX;
    dy = DY - AY;
    dw = dx*dx + dy*dy;
   }
  else
  { integer b1 = BX*AW - AX*BW;
    integer b2 = BY*AW - AY*BW;
    integer c1 = CX*AW - AX*CW;
    integer c2 = CY*AW - AY*CW;
    integer d1 = DX*AW - AX*DW;
    integer d2 = DY*AW - AY*DW;

    bx = b1 * AW * BW;
    by = b2 * AW * BW;
    bw = b1*b1 + b2*b2;

    cx = c1 * AW * CW;
    cy = c2 * AW * CW;
    cw = c1*c1 + c2*c2;

    dx = d1 * AW * DW;
    dy = d2 * AW * DW;
    dw = d1*d1 + d2*d2;
   }

  return sign((by*cx-bx*cy)*dw + (cy*bw-by*cw)*dx  + (bx*cw-cx*bw)*dy);
}


void rat_point::print_statistics()
{
  float cmp_percent = 0;
  float orient_percent = 0;
  float soc_percent = 0;
  
  if (cmp_count > 0) 
      cmp_percent= 100*float(exact_cmp_count)/cmp_count;
  
  if (orient_count > 0) 
      orient_percent= 100*float(exact_orient_count)/orient_count;
  
  if (soc_count > 0) 
      soc_percent= 100*float(exact_soc_count)/soc_count;
  
  cout << endl;
  cout << string("compare:        %6d /%6d   (%.2f %%)",
                  exact_cmp_count,cmp_count,cmp_percent) << endl;
  cout << string("orientation:    %6d /%6d   (%.2f %%)",
                  exact_orient_count,orient_count,orient_percent) << endl;
  cout << string("side of circle: %6d /%6d   (%.2f %%)",
                  exact_soc_count,soc_count,soc_percent) << endl;
  cout << endl;
}
  
