/*******************************************************************************
+
+  LEDA 3.5
+
+  _rat_circle.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.
+ 
*******************************************************************************/
#include <LEDA/rat_circle.h>
#include <LEDA/rat_line.h>
#include <LEDA/circle.h>
#include <math.h>


//------------------------------------------------------------------------------
// rat_circles
//
// S. N"aher (1996)
//------------------------------------------------------------------------------


// 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

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

/*
Kurt:: The constructor checks whether the points p1, p2, p3 are admissable. 
If all three of them are equal it defines the center.
*/

rat_circle_rep::rat_circle_rep(const rat_point& p1, 
                               const rat_point& p2, 
                               const rat_point& p3) : 
a(p1), b(p2), c(p3), cp(0), first_side_of(true)
{ orient = orientation(p1,p2,p3);
  if (orient == 0 && (p1 == p2 || p1 == p3 || p2 == p3) )
  { if (p1 == p2 && p1 == p3)
       cp = new rat_point(p1);
    else
       error_handler(1,"rat_circle::rat_circle:non-admissable triple");
  }
} 
                              



rat_circle::rat_circle(const rat_point& a,const rat_point& b,const rat_point& c)
{ PTR = new rat_circle_rep(a,b,c); }


rat_circle::rat_circle(const rat_point& a, const rat_point& b0)
{ // center a and b0 on circle
  rat_point b1 = b0.rotate90(a);
  rat_point b2 = b1.rotate90(a);
  PTR = new rat_circle_rep(b0,b1,b2);
  ptr()->cp = new rat_point(a);
 }


rat_circle::rat_circle(const rat_point& a) 
{ PTR = new rat_circle_rep(a,a,a);
  ptr()->cp = new rat_point(a);
}

rat_circle::rat_circle() 
{ rat_point p(0,0,1);
  PTR = new rat_circle_rep(p,p,p);
  ptr()->cp = new rat_point(p);
 }


circle rat_circle::to_circle() const
{ point p = ptr()->a.to_point();
  point q = ptr()->b.to_point();
  point r = ptr()->c.to_point();
  return circle(p,q,r);
 }


rat_point rat_circle::center() const
{ 
  if (ptr()->cp) return *(ptr()->cp);

  if (is_trivial()) 
  { ptr()->cp = new rat_point(ptr()->a);
    return *(ptr()->cp);
   }

   if (ptr()->orient == 0)
     error_handler(1,"rat_circle::center(): points are collinear.");

  // The formula is analogous to the PhD dissertation of
  // C. Burnikel (page 30 above).

  rat_point p = ptr()->a;
  rat_point q = ptr()->b;
  rat_point r = ptr()->c;

 
  integer xp = p.X(), yp = p.Y(), zp = p.W();
  integer xq = q.X(), yq = q.Y(), zq = q.W();
  integer xr = r.X(), yr = r.Y(), zr = r.W();
 
  integer cx  =  zq*zr*(xp*xp+yp*yp)*(yq*zr-yr*zq)
               + zr*zp*(xq*xq+yq*yq)*(yr*zp-yp*zr)
               + zp*zq*(xr*xr+yr*yr)*(yp*zq-yq*zp);
 
  integer cy  = -zq*zr*(xp*xp+yp*yp)*(xq*zr-xr*zq)
                -zr*zp*(xq*xq+yq*yq)*(xr*zp-xp*zr)
                -zp*zq*(xr*xr+yr*yr)*(xp*zq-xq*zp);

  integer cw  =  2*zp*zq*zr*(zp*(xq*yr-xr*yq) 
               + zq*(xr*yp-xp*yr) + zr*(xp*yq-xq*yp));

  ptr()->cp = new rat_point(cx,cy,cw);


/*
  rat_line l1 = p_bisector(ptr()->a,ptr()->b);
  rat_line l2 = p_bisector(ptr()->b,ptr()->c);
  rat_point m;
  l1.intersection(l2,m);
  ptr()->cp = new rat_point(m);
*/


  return *(ptr()->cp);
}


/*
Kurt: outside muss auch wahr liefern, wenn orient == 0, 
alle drei Punkte gleich und p ungleich diesem Punkt.
*/

bool rat_circle::outside(const rat_point& p) const 
{ if ( ptr()->orient != 0 )
    return (ptr()->orient * side_of(p)) < 0;
  else
    return (!is_trivial() && p != ptr()->a); 
 }


bool rat_circle::inside(const rat_point& p) const 
{ return (ptr()->orient * side_of(p)) > 0; }


bool rat_circle::contains(const rat_point& p) const 
{ if (is_trivial()) 
     return p == ptr()->a;
  else
     return side_of(p) == 0;
}


rat_circle rat_circle::translate(integer dx, integer dy, integer dw) const
{ rat_point a1 = ptr()->a.translate(dx,dy,dw);
  rat_point b1 = ptr()->b.translate(dx,dy,dw);
  rat_point c1 = ptr()->c.translate(dx,dy,dw);
  return rat_circle(a1,b1,c1);
 }

rat_circle rat_circle::translate(const rational& dx, const rational& dy) const
{ rat_point a1 = ptr()->a.translate(dx,dy);
  rat_point b1 = ptr()->b.translate(dx,dy);
  rat_point c1 = ptr()->c.translate(dx,dy);
  return rat_circle(a1,b1,c1);
 }

rat_circle rat_circle::translate(const rat_vector& v) const
{ rat_point a1 = ptr()->a.translate(v);
  rat_point b1 = ptr()->b.translate(v);
  rat_point c1 = ptr()->c.translate(v);
  return rat_circle(a1,b1,c1);
 }



rat_circle  rat_circle::rotate90(const rat_point& o)  const
{ rat_point a1 = ptr()->a.rotate90(o);
  rat_point b1 = ptr()->b.rotate90(o);
  rat_point c1 = ptr()->c.rotate90(o);
  return rat_circle(a1,b1,c1);
}


rat_circle  rat_circle::reflect(const rat_point& p, const rat_point& q)  const
{ rat_point a1 = ptr()->a.reflect(p,q);
  rat_point b1 = ptr()->b.reflect(p,q);
  rat_point c1 = ptr()->c.reflect(p,q);
  return rat_circle(a1,b1,c1);
}


rat_circle  rat_circle::reflect(const rat_point& p)  const
{ rat_point a1 = ptr()->a.reflect(p);
  rat_point b1 = ptr()->b.reflect(p);
  rat_point c1 = ptr()->c.reflect(p);
  return rat_circle(a1,b1,c1);
}


bool rat_circle::operator==(const rat_circle& c)  const
{ if (!contains(c.ptr()->a)) return false;
  if (!contains(c.ptr()->b)) return false;
  if (!contains(c.ptr()->c)) return false;
  return (is_trivial() || orientation() == c.orientation());
}


bool equal_as_sets(const rat_circle& c1, const rat_circle& c2)
{ if (!c1.contains(c2.point1())) return false;
  if (!c1.contains(c2.point2())) return false;
  if (!c1.contains(c2.point3())) return false;
  return true;
}


ostream& operator<<(ostream& out, const rat_circle& c) 
{ out << c.ptr()->a << " " << c.ptr()->b << " " << c.ptr()->c;
  return out;
 } 

istream& operator>>(istream& in,  rat_circle& c) 
{ rat_point p,q,r;
  in >> p >> q >> r;
  c = rat_circle(p,q,r);
  return in;
}



int rat_circle::side_of(const rat_point& p) const
{ 
  integer AX = ptr()->a.X();
  integer AY = ptr()->a.Y();
  integer AW = ptr()->a.W();

  if (ptr()->first_side_of)
  {
    integer BX = ptr()->b.X();
    integer BY = ptr()->b.Y();
    integer BW = ptr()->b.W();
  
    integer CX = ptr()->c.X();
    integer CY = ptr()->c.Y();
    integer CW = ptr()->c.W();

    integer b1 = BX*AW - AX*BW;
    integer b2 = BY*AW - AY*BW;
    integer bx = b1 * AW * BW;
    integer by = b2 * AW * BW;
    integer bw = b1*b1 + b2*b2;

    integer c1 = CX*AW - AX*CW;
    integer c2 = CY*AW - AY*CW;
    integer cx = c1 * AW * CW;
    integer cy = c2 * AW * CW;
    integer cw = c1*c1 + c2*c2;

    ptr()->D1 = by*cx - bx*cy;
    ptr()->D2 = cy*bw - by*cw;
    ptr()->D3 = bx*cw - cx*bw;

    ptr()->first_side_of = false;
  }

  // floating point filter :

  double pxaw = p.XD() * ptr()->a.WD();   
  double pyaw = p.YD() * ptr()->a.WD();
  double axpw = ptr()->a.XD() * p.WD(); 
  double aypw = ptr()->a.YD() * p.WD(); 

  double pax = pxaw - axpw;
  double pay = pyaw - aypw;
  double paw = p.WD()*(ptr()->a.WD());

  double Pxi = pax * paw;
  double Pyi = pay * paw;
  double Pwi = pax*pax + pay*pay;

  double D1 = (ptr()->D1).todouble();
  double D2 = (ptr()->D2).todouble();
  double D3 = (ptr()->D3).todouble();

  double E = Pwi*D1 + Pxi*D2 + Pyi*D3;
        
  /* ----------------------------------------------------------------
   * ERROR BOUND:
   * ----------------------------------------------------------------
   *
   * mes(E)   = mes(D1*Pwi + D2*Pxi + D3*Pyi)
   *          = 2*( mes(D1*Pwi) + 2*mes(D2*Pxi) + 2*mes(D3*Pyi) )
   *          = 2*( mes(D1)*mes(Pwi) + 2*mes(D2)*mes(Pxi) + 2*mes(D3)*mes(Pyi) )
   *        
   * mes(Pyi) = mes(pay) * mes(paw) 
   *          = mes(pyaw - aypw)* mes(paw)
   *          = 2*(mes(pyaw) + mes(aypw))* mes(paw)
   *          = 2*(fabs(pyaw) + fabs(aypw))*paw
   *            ( because a.WD() and p.WD() are positive )
   *
   * mes(Pxi) = 2*(fabs(pxaw) + fabs(axpw))*paw
   *
   * mes(Pwi) = 2*(mes(pax)*mes(pax) + mes(pay)*mes(pay))
   *          = 2*( 2*(fabs(pxaw) + fabs(axpw)) * 
   *                2*(fabs(pxaw) + fabs(axpw)) +
   *                2*(fabs(pyaw) + fabs(aypw)) * 
   *                2*(fabs(pyaw) + fabs(aypw)))
   *          = 8*((fabs(pxaw)+fabs(axpw))*(fabs(pxaw)+fabs(axpw)) +
   *               (fabs(pyaw)+fabs(aypw))*(fabs(pyaw)+fabs(aypw)))
   *
   * mes(E) = 2*( fabs(D1) *
   *              8*((fabs(pxaw)+fabs(axpw))*(fabs(pxaw)+fabs(axpw)) +
   *                 (fabs(pyaw)+fabs(aypw))*(fabs(byaw)+fabs(aybw)))
   *              + 
   *              fabs(D2)*4*(fabs(pxaw) + fabs(axpw))*paw
   *              +
   *              fabs(D3)*4*(fabs(pyaw) + fabs(aypw))*paw )
   *
   * ind(E) = ( ind(D1*Pwi) + ind(D2*Pxi + D3*Pyi) +1 ) /2
   *        = ( ind(D1*Pwi) + (ind(D2*Pxi) + ind(D3*Pyi) + 1)/2 +1) /2
   *        = ( ind(D1) + ind(Pwi) + 0.5 + (ind(D2) + ind(Pxi) + 0.5 + 
   *                                  ind(D3) +ind(Pyi) + 0.5 + 1)/2 +1) /2
   *
   * ind(Pyi) = ind(pay)+ind(paw)+0.5
   *          = ind(pxaw-axpw)+ind(baw)+0.5
   *          = (1.5 + 1.5 +1)/2 + 1.5 +0.5
   *          = 4 = ind(Pxi)
   *
   * ind(Pwi) = ind(pax*pax + pay*pay)
   *          = (ind(pax) + ind(pax) +0.5 + ind(pay) + ind(pay)+0.5+1)/2
   *          = (1.5 + 1.5 +1)/2 + (1.5 + 1.5 +1)/2 + 1
   *          = 5 
   *
   * ind(E) = (0.5  + 5 + 0.5 + (0.5 + 4 + 0.5 + 
   *                             0.5 + 4 + 0.5 + 1)/2 +1) /2
   *        = (6.0 + 5.5 + 1)/2
   *        = 25/4
   *
   * eps(E) = ind(E) * mes(E) * eps0
   *        = 50*( fabs(D1) *
   *              2*((fabs(pxaw)+fabs(axpw))*(fabs(pxaw)+fabs(axpw)) +
   *                 (fabs(pyaw)+fabs(aypw))*(fabs(pyaw)+fabs(aypw)))
   *              + 
   *              fabs(D2)*(fabs(pxaw) + fabs(axpw))*paw
   *              +
   *              fabs(D3)*(fabs(pyaw) + fabs(aypw))*paw ) * eps0
   *             
   * -----------------------------------------------------------------*/
  
  FABS(D1); 
  FABS(D2); 
  FABS(D3);
  FABS(pxaw); 
  FABS(axpw); 
  FABS(pyaw); 
  FABS(aypw);

  double eps = 50*( D1*2*((pxaw+axpw)*(pxaw+axpw) +
                          (pyaw+aypw)*(pyaw+aypw)) + 
                    D2*(pxaw+axpw)*paw +D3*(pyaw+aypw)*paw ) * circle_eps0;
  
  if (E   >  eps)  return  1;
  if (E   < -eps)  return -1;
  if (eps <  1)    return  0;

  integer p1 = p.X()*AW - AX*p.W();
  integer p2 = p.Y()*AW - AY*p.W();

  integer px = p1 * AW * p.W();
  integer py = p2 * AW * p.W();
  integer pw = p1*p1 + p2*p2;

  return sign(pw*ptr()->D1 + px*ptr()->D2  + py*ptr()->D3);
}


// We implement the procedure rat_point_on_circle that computes 
// for a circle C given by rational points p, q, r and an angle
// alpha a rational point s on C such that the angle of s 
// at C differs from alpha by at most a given delta


// The following procedure computes a rat_point phi_alpha
// on the circle C whose angle at C satisfies
//     |angle(phi_alpha)-alpha| < epsilon
// alpha and the corresponding error bound epsilon are in radians

// It is similar to the method shown in the paper by Canny et al:
// "A Rational Rotation Method for Robust Geometric Algorithms"
// (Proc. of the 8th ACM Symposium on Computational Geometry, 
//  pages 251-260, 1992)

// preconditions: 0 <= alpha <= 2*pi; epsilon >= 2^{-48}


static
void compute_phi_C( const integer& m,  const integer& n,
                    const integer& cx, const integer& cy, const integer& cw,
                    integer& phi_x, integer& phi_y, integer& phi_w)
{
  // C is a circle with center c_trans(cx,cy,cw) passing 
  // through the origin.
  // Every rational point on circle C has the form 
  //        phi(n/m) = (phi_x/phi_w,phi_y/phi_w)
  // for integeral l, m where (phi_x, phi_y, phi_w) are as below

  phi_x = (n*n-m*m) * cx + 2*m*n * cy;
  phi_y = (m*m-n*n) * cy + 2*m*n * cx;
  phi_w = (m*m+n*n) * cw;
}




// For given alpha we determine a value t = t(alpha) such that 
//       angle(phi(t)) = alpha.
// The solution is
//       t = (cy + r*sin(alpha)) / (cx + r*cos(alpha))
// where
//       r = sqrt(cx^2+cy^2)

static double compute_t(double alpha, const rat_point& center)
{
  // preconditions: 0 <= alpha <= pi/2, 
  //                0 <= angle(center) <= pi/4

  double si = sin(alpha), co = cos(alpha);
  double sinalpha, cosalpha;
  if (co < si) {  sinalpha = co; cosalpha = si; }
  else         {  sinalpha = si; cosalpha = co; }

  double cx = center.X().todouble();
  double cy = center.Y().todouble();
  
  // compute result

  integer icx = center.X(), icy = center.Y();
  integer rsq = icx*icx+icy*icy;
  double r = sqrt(rsq.todouble());
  double t = (cy + r*sinalpha) / (cx + r*cosalpha); 

  return(t);
}




rat_point rat_circle::point_on_circle(double alpha, double epsilon) const
{
  // precondition: epsilon >= 2^{-48}

  rat_point phi_alpha;
  rat_point c = center();
  rat_point c_trans;
  
  {
    rat_point p = point1();

    integer pfactor = c.W()/p.W();

    integer dx = c.X()-pfactor*p.X();
    integer dy = c.Y()-pfactor*p.Y();
   
    if (abs(dx) >= abs(dy))
       c_trans = rat_point(abs(dx),abs(dy),c.W());
    else
       c_trans = rat_point(abs(dy),abs(dx),c.W()); 
  }

  integer ctx = c_trans.X();
  integer cty = c_trans.Y();
  integer ctw = c_trans.W();

  double alpha_normal, alpha_trans;
  {
    int sector_2_pi    = (int) floor(alpha/(4*asin(1.0)));
    alpha_normal = alpha - sector_2_pi * 4* asin(1.0);
    int sector_pi_half = (int) floor(alpha_normal/asin(1.0));
    alpha_trans  = alpha_normal - sector_pi_half * asin(1.0);
  }

  double sinalpha = sin(alpha_normal);
  double cosalpha = cos(alpha_normal);

  double t = compute_t(alpha_trans,c_trans);
  double delta;
  {
    delta = epsilon - 26*t*ldexp(1.0,-53);
    delta = ldexp(delta,-1);

    if (delta <= 0) 
      error_handler(1,"rat_direction called with epsilon_degree too small");
  }


  rational rt_low = t - delta, rt_high = t + delta;
  if (rt_high > 1) rt_high = 1;
  if (rt_low  < 0) rt_low  = 0;
  rational rt = small_rational_between(rt_low, rt_high);
  integer m = rt.numerator(), n = rt.denominator();

  integer phi_x, phi_y, phi_w;
  compute_phi_C(m,n,ctx,cty,ctw,phi_x,phi_y,phi_w);
  
  {
    integer xresult, yresult, wresult;

    if (fabs(cosalpha) >= fabs(sinalpha))
    {
      xresult = phi_x;
      yresult = phi_y;
    }
    else
    {
      xresult = phi_y;
      yresult = phi_x;
    }

    wresult = phi_w;

    if (cosalpha < 0) xresult = -xresult;
    if (sinalpha < 0) yresult = -yresult;

    phi_alpha = rat_point(
      xresult+(m*m+n*n)*c.X(),yresult+(m*m+n*n)*c.Y(),wresult);
  }


  
  {
    if (!contains(phi_alpha))
      error_handler(1,"rat_point_on_circle: result not on circle");

    rat_point phi_C = phi_alpha.translate(-c.X(),-c.Y(),c.W());

    double x = phi_C.xcoordD(), y = phi_C.ycoordD();
    double radius = sqrt(x*x+y*y);

    if (fabs(x/radius - cosalpha) > epsilon)
      error_handler(1,"rat_point_on_circle: phi_x/phi_w not correct");
    if (fabs(y/radius - sinalpha) > epsilon)
      error_handler(1,"rat_point_on_circle: phi_y/phi_w not correct");
  }


  return phi_alpha;
}



