home *** CD-ROM | disk | FTP | other *** search
/ Black Box 4 / BlackBox.cdr / progc / djlsr106.arj / RATIONAL.CC < prev    next >
C/C++ Source or Header  |  1992-03-30  |  8KB  |  401 lines
  1. /* 
  2. Copyright (C) 1988 Free Software Foundation
  3.     written by Doug Lea (dl@rocky.oswego.edu)
  4.  
  5. This file is part of the GNU C++ Library.  This library is free
  6. software; you can redistribute it and/or modify it under the terms of
  7. the GNU Library General Public License as published by the Free
  8. Software Foundation; either version 2 of the License, or (at your
  9. option) any later version.  This library is distributed in the hope
  10. that it will be useful, but WITHOUT ANY WARRANTY; without even the
  11. implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  12. PURPOSE.  See the GNU Library General Public License for more details.
  13. You should have received a copy of the GNU Library General Public
  14. License along with this library; if not, write to the Free Software
  15. Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
  16. */
  17.  
  18. #ifdef __GNUG__
  19. #pragma implementation "Rational.h"
  20. #endif
  21. #include <Rational.h>
  22. #include <std.h>
  23. #include <math.h>
  24. #include <values.h>
  25. #include <builtin.h>
  26.  
  27.  
  28. void Rational::error(const char* msg) const
  29. {
  30.   (*lib_error_handler)("Rational", msg);
  31. }
  32.  
  33. static const Integer _Int_One(1);
  34.  
  35. void Rational::normalize()
  36. {
  37.   int s = sign(den);
  38.   if (s == 0)
  39.     error("Zero denominator.");
  40.   else if (s < 0)
  41.   {
  42.     den.negate();
  43.     num.negate();
  44.   }
  45.  
  46.   Integer g = gcd(num, den);
  47.   if (ucompare(g, _Int_One) != 0)
  48.   {
  49.     num /= g;
  50.     den /= g;
  51.   }
  52. }
  53.  
  54. void      add(const Rational& x, const Rational& y, Rational& r)
  55. {
  56.   if (&r != &x && &r != &y)
  57.   {
  58.     mul(x.num, y.den, r.num);
  59.     mul(x.den, y.num, r.den);
  60.     add(r.num, r.den, r.num);
  61.     mul(x.den, y.den, r.den);
  62.   }
  63.   else
  64.   {
  65.     Integer tmp;
  66.     mul(x.den, y.num, tmp);
  67.     mul(x.num, y.den, r.num);
  68.     add(r.num, tmp, r.num);
  69.     mul(x.den, y.den, r.den);
  70.   }
  71.   r.normalize();
  72. }
  73.  
  74. void      sub(const Rational& x, const Rational& y, Rational& r)
  75. {
  76.   if (&r != &x && &r != &y)
  77.   {
  78.     mul(x.num, y.den, r.num);
  79.     mul(x.den, y.num, r.den);
  80.     sub(r.num, r.den, r.num);
  81.     mul(x.den, y.den, r.den);
  82.   }
  83.   else
  84.   {
  85.     Integer tmp;
  86.     mul(x.den, y.num, tmp);
  87.     mul(x.num, y.den, r.num);
  88.     sub(r.num, tmp, r.num);
  89.     mul(x.den, y.den, r.den);
  90.   }
  91.   r.normalize();
  92. }
  93.  
  94. void      mul(const Rational& x, const Rational& y, Rational& r)
  95. {
  96.   mul(x.num, y.num, r.num);
  97.   mul(x.den, y.den, r.den);
  98.   r.normalize();
  99. }
  100.  
  101. void      div(const Rational& x, const Rational& y, Rational& r)
  102. {
  103.   if (&r != &x && &r != &y)
  104.   {
  105.     mul(x.num, y.den, r.num);
  106.     mul(x.den, y.num, r.den);
  107.   }
  108.   else
  109.   {
  110.     Integer tmp;
  111.     mul(x.num, y.den, tmp);
  112.     mul(y.num, x.den, r.den);
  113.     r.num = tmp;
  114.   }
  115.   r.normalize();
  116. }
  117.  
  118.  
  119.  
  120.  
  121. void Rational::invert()
  122. {
  123.   Integer tmp = num;  
  124.   num = den;  
  125.   den = tmp;  
  126.   int s = sign(den);
  127.   if (s == 0)
  128.     error("Zero denominator.");
  129.   else if (s < 0)
  130.   {
  131.     den.negate();
  132.     num.negate();
  133.   }
  134. }
  135.  
  136. int compare(const Rational& x, const Rational& y)
  137. {
  138.   int xsgn = sign(x.num);
  139.   int ysgn = sign(y.num);
  140.   int d = xsgn - ysgn;
  141.   if (d == 0 && xsgn != 0) d = compare(x.num * y.den, x.den * y.num);
  142.   return d;
  143. }
  144.  
  145. Rational::Rational(double x)
  146. {
  147.   num = 0;
  148.   den = 1;
  149.   if (x != 0.0)
  150.   {
  151.     int neg = x < 0;
  152.     if (neg)
  153.       x = -x;
  154.  
  155.     const long shift = 15;         // a safe shift per step
  156.     const double width = 32768.0;  // = 2^shift
  157.     const int maxiter = 20;        // ought not be necessary, but just in case,
  158.                                    // max 300 bits of precision
  159.     int expt;
  160.     double mantissa = frexp(x, &expt);
  161.     long exponent = expt;
  162.     double intpart;
  163.     int k = 0;
  164.     while (mantissa != 0.0 && k++ < maxiter)
  165.     {
  166.       mantissa *= width;
  167.       mantissa = modf(mantissa, &intpart);
  168.       num <<= shift;
  169.       num += (long)intpart;
  170.       exponent -= shift;
  171.     }
  172.     if (exponent > 0)
  173.       num <<= exponent;
  174.     else if (exponent < 0)
  175.       den <<= -exponent;
  176.     if (neg)
  177.       num.negate();
  178.   }
  179.   normalize();
  180. }
  181.  
  182.  
  183. Integer trunc(const Rational& x)
  184. {
  185.   return x.num / x.den ;
  186. }
  187.  
  188.  
  189. Rational pow(const Rational& x, const Integer& y)
  190. {
  191.   long yy = long(y);
  192.   return pow(x, yy);
  193. }               
  194.  
  195. #if defined(__GNUG__) && !defined(NO_NRV)
  196.  
  197. Rational operator - (const Rational& x) return r(x)
  198. {
  199.   r.negate();
  200. }
  201.  
  202. Rational abs(const Rational& x) return r(x)
  203. {
  204.   if (sign(r.num) < 0) r.negate();
  205. }
  206.  
  207.  
  208. Rational sqr(const Rational& x) return r
  209. {
  210.   mul(x.num, x.num, r.num);
  211.   mul(x.den, x.den, r.den);
  212.   r.normalize();
  213. }
  214.  
  215. Integer floor(const Rational& x) return q
  216. {
  217.   Integer r;
  218.   divide(x.num, x.den, q, r);
  219.   if (sign(x.num) < 0 && sign(r) != 0) q--;
  220. }
  221.  
  222. Integer ceil(const Rational& x) return q
  223. {
  224.   Integer  r;
  225.   divide(x.num, x.den, q, r);
  226.   if (sign(x.num) >= 0 && sign(r) != 0) q++;
  227. }
  228.  
  229. Integer round(const Rational& x) return q
  230. {
  231.   Integer r;
  232.   divide(x.num, x.den, q, r);
  233.   r <<= 1;
  234.   if (ucompare(r, x.den) >= 0)
  235.   {
  236.     if (sign(x.num) >= 0)
  237.       q++;
  238.     else
  239.       q--;
  240.   }
  241. }
  242.  
  243. // power: no need to normalize since num & den already relatively prime
  244.  
  245. Rational pow(const Rational& x, long y) return r
  246. {
  247.   if (y >= 0)
  248.   {
  249.     pow(x.num, y, r.num);
  250.     pow(x.den, y, r.den);
  251.   }
  252.   else
  253.   {
  254.     y = -y;
  255.     pow(x.num, y, r.den);
  256.     pow(x.den, y, r.num);
  257.     if (sign(r.den) < 0)
  258.     {
  259.       r.num.negate();
  260.       r.den.negate();
  261.     }
  262.   }
  263. }
  264.  
  265. #else
  266.  
  267. Rational operator - (const Rational& x) 
  268. {
  269.   Rational r(x); r.negate(); return r;
  270. }
  271.  
  272. Rational abs(const Rational& x) 
  273. {
  274.   Rational r(x);
  275.   if (sign(r.num) < 0) r.negate();
  276.   return r;
  277. }
  278.  
  279.  
  280. Rational sqr(const Rational& x)
  281. {
  282.   Rational r;
  283.   mul(x.num, x.num, r.num);
  284.   mul(x.den, x.den, r.den);
  285.   r.normalize();
  286.   return r;
  287. }
  288.  
  289. Integer floor(const Rational& x)
  290. {
  291.   Integer q;
  292.   Integer r;
  293.   divide(x.num, x.den, q, r);
  294.   if (sign(x.num) < 0 && sign(r) != 0) q--;
  295.   return q;
  296. }
  297.  
  298. Integer ceil(const Rational& x)
  299. {
  300.   Integer q;
  301.   Integer  r;
  302.   divide(x.num, x.den, q, r);
  303.   if (sign(x.num) >= 0 && sign(r) != 0) q++;
  304.   return q;
  305. }
  306.  
  307. Integer round(const Rational& x) 
  308. {
  309.   Integer q;
  310.   Integer r;
  311.   divide(x.num, x.den, q, r);
  312.   r <<= 1;
  313.   if (ucompare(r, x.den) >= 0)
  314.   {
  315.     if (sign(x.num) >= 0)
  316.       q++;
  317.     else
  318.       q--;
  319.   }
  320.   return q;
  321. }
  322.  
  323. Rational pow(const Rational& x, long y)
  324. {
  325.   Rational r;
  326.   if (y >= 0)
  327.   {
  328.     pow(x.num, y, r.num);
  329.     pow(x.den, y, r.den);
  330.   }
  331.   else
  332.   {
  333.     y = -y;
  334.     pow(x.num, y, r.den);
  335.     pow(x.den, y, r.num);
  336.     if (sign(r.den) < 0)
  337.     {
  338.       r.num.negate();
  339.       r.den.negate();
  340.     }
  341.   }
  342.   return r;
  343. }
  344.  
  345. #endif
  346.  
  347. ostream& operator << (ostream& s, const Rational& y)
  348. {
  349.   if (y.denominator() == 1)
  350.     s << y.numerator();
  351.   else
  352.   {
  353.     s << y.numerator();
  354.     s << "/";
  355.     s << y.denominator();
  356.   }
  357.   return s;
  358. }
  359.  
  360. istream& operator >> (istream& s, Rational& y)
  361. {
  362. #ifdef _OLD_STREAMS
  363.   if (!s.good())
  364.   {
  365.     return s;
  366.   }
  367. #else
  368.   if (!s.ipfx(0))
  369.   {
  370.     s.set(ios::failbit); // Redundant if using GNU iostreams.
  371.     return s;
  372.   }
  373. #endif
  374.   Integer n = 0;
  375.   Integer d = 1;
  376.   if (s >> n)
  377.   {
  378.     char ch = 0;
  379.     s.get(ch);
  380.     if (ch == '/')
  381.     {
  382.       s >> d;
  383.     }
  384.     else
  385.     {
  386.       s.putback(ch);
  387.     }
  388.   }
  389.   y = Rational(n, d);
  390.   return s;
  391. }
  392.  
  393. int Rational::OK() const
  394. {
  395.   int v = num.OK() && den.OK(); // have valid num and denom
  396.   v &= sign(den) > 0;           // denominator positive;
  397.   v &=  ucompare(gcd(num, den), _Int_One) == 0; // relatively prime
  398.   if (!v) error("invariant failure");
  399.   return v;
  400. }
  401.