/*******************************************************************************
+
+  LEDA 3.5
+
+  _rational.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/basic.h>
#include <LEDA/rational.h>
#include <math.h>
#include <ctype.h>


rational& rational::normalize()
{ 
  // divide num and den by their gcd(num,den)
  // precondition: den  > 0

  if (num == den) 
    num = den = 1; 
  else
    if (-num == den) 
      { num = -1; 
        den = 1; 
       }
    else
      { integer ggt = gcd(num, den);
        if (ggt != 1) 
        { num /= ggt;
          den /= ggt;
         }
       }

  return *this;
}


rational& rational::simplify(const integer& a)
{ 
  // divide num and den by a
 
  integer r;
  den = den.div(a,r);

  if (r !=0)
    error_handler(1,"rational::simplify: argument does not divide denominator");

  num = num.div(a,r);

  if (r !=0)
    error_handler(1,"rational::simplify: argument does not divide numerator");

  return *this;
}



// constructors

rational::rational(int n, int d)
{ 
  if (d == 0) error_handler(1,"rational: Zero denominator!");

  if (d < 0) { n = -n; d = -d; }  //wegen reference counting notwendig

  num = integer(n); 
  den = integer(d);
 }


rational::rational(integer n, integer d)
{ 
  if (d.PTR->sign == 0) error_handler(1,"rational: Zero denominator!");
 
  if (d.PTR->sign == -1) 
    { num = -n; den = -d; }
  else
    { num =  n; den =  d; }
}


rational::rational(double x)
  { num = 0; 
    den = 1;

    if (x != 0)
    { int neg = (x < 0);
      if (neg) x = -x;

      const unsigned shift = 15;   // a safe shift per step
      const double width = 32768;  // = 2^shift
      const int maxiter = 20;      // ought not be necessary, but just in case,
                                   // max 300 bits of precision
      int expt;
      double mantissa = frexp(x, &expt);
      long exponent = expt;
      double intpart;
      int k = 0;
      while (mantissa != 0.0 && k++ < maxiter)
      { mantissa *= width; // shift double mantissa
        mantissa = modf(mantissa, &intpart);
        num <<= shift;
        num += (long)intpart;
        exponent -= shift;
      }
      if (exponent > 0)
        num <<= (unsigned)exponent;
      else if (exponent < 0)
        den <<= (unsigned)(-exponent);
      if (neg)
        num = -num;
    }
    normalize();
  }




// operators

rational& rational::operator+= (const rational& r)
{ 
  if (den == r.den)
    num += r.num;
  else
    { num = num*r.den + r.num*den;
      den *= r.den;
     }
  return *this;
}

rational& rational::operator-= (const rational& r)
{ if (den == r.den)
    num -= r.num;
  else
    { num = num*r.den - r.num*den;
      den *= r.den;
     }
  return *this;
}

rational& rational::operator*= (const rational& r)
{ num *= r.num;
  den *= r.den;
  return *this;
}

rational& rational::operator/= (const rational& r)
{ 
  if (r.num.PTR->sign == 0) error_handler(1,"division by 0");

  if (r.num.PTR->sign == -1)
    { num *= -r.den;
      den *= -r.num;
     }
   else
    { num *= r.den;
      den *= r.num;
     }
  return *this;
}


rational& rational::operator= (const rational& r)
{ if (this == &r) return *this;
  num = r.num;
  den = r.den;
  return *this;
}



void rational::invert()
{ 
  if (num.PTR->sign == 0) error_handler(1,"Zero denominator!");

  integer tmp = num;
  if (num.PTR->sign == 1)
    { num = den;
      den = tmp;
     }
  else
    { num = -den;
      den = -tmp;
     }
}


rational rational::inverse()
{ 
  if (num.PTR->sign == 0) error_handler(1,"Zero denominator!");
  
  if (num.PTR->sign == 1)
     return rational(den,num);
  else
     return rational(-den,-num);
}




int rational::cmp(const rational& x, const rational& y)
{ int xsign = sign(x.num);
  int ysign = sign(y.num);

  if (xsign == 0) return -ysign;
  if (ysign == 0) return  xsign;

  if (xsign == ysign) 
    return compare(x.num * y.den, y.num * x.den);
  else 
    return compare(xsign,ysign);
}


// friend functions

rational pow(const rational& r, int i)
{ 
  //r.normalize();

  rational mul(r);
  rational result(1,1);

  if (i < 0) 
  { mul = mul.inverse();
    i = -i;
   }

  while (i--)
  { result.num *= mul.num;
    result.den *= mul.den;
   }

  return result;
 }


rational pow(const rational& r, integer I)
{ 
  //r.normalize();

  rational mul(r);
  rational result(1,1);

  if (I < 0)
  { mul = mul.inverse();
    I=-I;
   }

  while (I-- != 0)
  { result.num *= mul.num;
    result.den *= mul.den;
   }

  return result;
 }


double rational::to_double() const
{ 
  if (num == 0) return 0;
  
  integer p = 1;

  if (den > abs(num)) p = den/num;
  
  integer q = ((num*p) << 53)/den;
  
  return (q.to_double()/p.to_double()) * ldexp(1.0,-53);
 } 


/*
 double rational::to_double() const
 { return div(bigfloat(num),bigfloat(den),53).to_double();
 }
*/



// floor, ceil, round (use schooldiv!)

integer floor(const rational& r)
{ integer x = r.num/r.den;
  if (sign(r.num) == -1 && r.num%r.den != 0) x--;
  return x;
 }


integer ceil(const rational& r)
{ integer x = r.num/r.den; 
  if (sign(r.num) == -1 && r.num%r.den != 0) x++;
  return x;
 }


integer round(const rational& r)
{ integer rem;
  integer quot = (r.num).div(r.den, rem);

  rem <<= 1;
  if (sign(rem) < 0) 
      { if (-rem >= r.den) quot--; }
  else 
      { if ( rem >= r.den) quot++; } 

  return quot;
}



ostream& operator<< (ostream& s, const rational& r)
{  s << r.num << "/" << r.den; 
   return s; 
}

istream& operator>> (istream& in, rational& r)
{ 
   // Format: "r.num / r.den"

   char c;

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

   integer rx;
   in >> rx;

   do in.get(c); while (isspace(c));
   if (c != '/')   error_handler(1,"rational input: syntax error."); 

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

   integer ry;
   in >> ry;

   r = rational(rx,ry);

   return in;
}




  rational& rational::operator++ ()       { num += den; return *this; }
  rational& rational::operator-- ()       { num -= den; return *this; }
  integer   rational::numerator()   const { return num; }
  integer   rational::denominator() const { return den; }
  void      rational::negate()            { num = - num; }


  rational operator+(const rational& x, const rational& y)
  { rational z = x; return z += y; }

  rational operator-(const rational& x, const rational& y)
  { rational z = x; return z -= y; }

  rational operator*(const rational& x, const rational& y)
  { rational z = x; return z *= y; }

  rational operator/(const rational& x, const rational& y)
  { rational z = x; return z /= y; }

  rational operator-(const rational& x)
  { return rational(-x.num,x.den); }

  int sign(const rational& r)
  { return sign(r.num); }

  rational abs(const rational& r)
  { if (sign(r.num) > -1) return r; else return -r; }

  rational sqr(const rational& r)
  { return rational(r.num*r.num, r.den*r.den); }

  integer trunc(const rational& r)
  { return (r.num / r.den); }



// comparison operators

  bool operator==(const rational& x, const rational& y)
  { return x.num * y.den == x.den * y.num; }

  bool operator==(const rational& x, int y)
  { return x.den * y == x.num; }

  bool operator==(int x, const rational& y)
  { return y.den * x == y.num; }

  bool operator==(const rational& x, const integer& y)
  { return x.den * y == x.num; }

  bool operator==(const integer& x, const rational& y)
  { return y.den * x == y.num; }

  bool operator!=(const rational& x, const rational& y)
  { return x.num * y.den != x.den * y.num; }

  bool operator!=(const rational& x, int y)
  { return x.den * y != x.num; }

  bool operator!=(int x, const rational& y)
  { return y.den * x != y.num; }

  bool operator!=(const rational& x, const integer& y)
  { return x.den * y != x.num; }

  bool operator!=(const integer& x, const rational& y)
  { return y.den * x != y.num; }
    
  bool operator<(const rational& x, const rational& y)
  { return rational::cmp(x,y) < 0; }

  bool operator<=(const rational& x, const rational& y)
  { return rational::cmp(x,y) <= 0; }

  bool operator>(const rational& x, const rational& y)
  { return rational::cmp(x,y) > 0; }

  bool operator>=(const rational& x, const rational& y)
  { return rational::cmp(x,y) >= 0; }



rational small_rational_between(const rational& X, const rational& Y)
{  
  // This procedure is similar to Figures 6, 7 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)

  if (Y < X) 
    error_handler(1, "small_rational_between: precondition violated");

  rational floor_y = floor(Y);

  if (floor_y >= X) return floor_y;
 
  integer p0 = 0;
  integer q0 = 1;
  integer p1 = 1;
  integer q1 = 1;
  integer r;

  rational x = X-floor_y;
  rational y = Y-floor_y;
  rational result;

  while(1)
  {  // invariant (I):    0 <= p0/q0 < x <= y < p1/q1 <= 1

    r = floor((rational(p1)-rational(q1)*x)/(rational(q0)*x-rational(p0)));
    p1 = r*p0+p1; 
    q1 = r*q0+q1;
    if ( rational(p1,q1) <= y)
    { result = rational(p1+q1*floor(Y),q1); 
      break;
    }

    // Invariant (I) from above holds again 

    r = floor((rational(q0)*y-rational(p0))/(rational(p1)-rational(q1)*y));
    p0 = r*p1+p0; 
    q0 = r*q1+q0;
    if (rational(p0,q0) >= x)
    { result = rational(p0+q0*floor(Y),q0);
      break;
    }
  } 
  
  if ((result > Y) || (result < X))
         error_handler(1, "small_rational_between: internal error.");

  return result;

}   



rational small_rational_near(const rational& X, rational epsilon)
{  
  // returns a small rational between X-epsilon and X+epsilon
  // precondition: epsilon >= 0

  if (sign(epsilon) < 0)
     error_handler(1,
        "small_rational_near: precondition violated");

  return small_rational_between(X-epsilon,X+epsilon);
}

