OS/2 Shareware BBS: 10 Tools/10-Tools.zip/ledar34.zip/leda-r-3_4_tar/LEDA-3.4/src/numbers/_rational.c

C/C++ Source or Header | 1996-09-03 | 10KB | 504 lines

/*******************************************************************************
+
+  LEDA 3.4
+
+  _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 [phone]).
+
+  Copyright (c) 1991-1996  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,"/ expected"); }

   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");

  integer floor_y = floor(Y);
  if (floor_y >= X) return floor_y;

  integer p0 = 0, q0 = 1, p1 = 1, q1 = 1, r;
  rational x = X-floor_y, y = Y-floor_y;
  rational z, result;

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

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

   // Invariant (I) from above holds again

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

  if ((result > Y) || (result < X))
         error_handler(1,
           "small_rational_between: result wrong");
  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);
}