------------------------->  GNU Sather - sourcefile  <-------------------------
-- Copyright (C) 1994 by International Computer Science Institute            --
-- This file is part of the GNU Sather library. It is free software; you may --
-- redistribute  and/or modify it under the terms of the GNU Library General --
-- Public  License (LGPL)  as published  by the  Free  Software  Foundation; --
-- either version 2 of the license, or (at your option) any later version.   --
-- This  library  is distributed  in the  hope that it will  be  useful, but --
-- WITHOUT ANY WARRANTY without even the implied warranty of MERCHANTABILITY --
-- or FITNESS FOR A PARTICULAR PURPOSE. See Doc/LGPL for more details.       --
-- The license text is also available from:  Free Software Foundation, Inc., --
-- 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA                     --
-------------->  Please email comments to <bug-sather@gnu.org>  <--------------

-- Implementation of rational numbers.
----------------------------------------------------------------------
immutable class RAT < $IS_LT{RAT}, $STR is
   -- Rational numbers. A rational number is represented by by two
   -- integers u and v (numerator and denominator) and always
   -- normalized, i.e. u.gcd(v) = 1 and v > 0.
   
   include COMPARABLE;
   
   readonly attr u, v: INTI; -- numerator, denominator

   -------------------------------------------------- object creation

   create (u: INTI): SAME is
      r: SAME; return r.u(u).v(#INTI(1))
   end;

   create (u: INT): SAME is
      return create(#INTI(u))
   end;

   create (u, v: INTI): SAME is
      assert ~v.is_zero;
      if u.is_zero then v := #INTI(1)
      else g ::= u.gcd(v);
         if v.is_neg then g := -g end;
         u := u/g; v := v/g
      end;
      assert v.is_pos;
      r: SAME; return r.u(u).v(v)
   end;

   create (u, v: INT): SAME is
      return create(#INTI(u), #INTI(v))
   end;

-------------------------------------------------- binary arithmetics

   plus (y: SAME): SAME is return #SAME(u*y.v + v*y.u, v*y.v) end;
   minus (y: SAME): SAME is return #SAME(u*y.v - v*y.u, v*y.v) end;
   times (y: SAME): SAME is return #SAME(u*y.u, v*y.v) end;
   div (y: SAME): SAME is assert ~y.u.is_zero; return #SAME(u*y.v, v*y.u) end;
   mod (y: SAME): SAME is return self - #SAME((self/y).floor)*y end;

   pow (i: INT): SAME is
   -- Returns self raised to the power i. Returns 1 for i < 0.
   --
      x ::= self; z ::= #SAME(1);
      loop while!(i > 0);
         -- z * x^i = self ^ i0
         if i.is_odd then z := z*x end;
         x := x.square; i := i/2
      end;
      return z
   end;

-------------------------------------------------------------------
   
   hash:INT is return u.hash.bxor(v.hash) end;

-------------------------------------------------- binary relations

   is_eq (y: SAME): BOOL is return (u = y.u) and (v = y.v) end;
   is_lt (y: SAME): BOOL is return u*y.v < v*y.u end;

-------------------------------------------------- unary predicates

   is_pos: BOOL is return u.is_pos end;
   is_neg: BOOL is return u.is_neg end;
   is_zero: BOOL is return u.is_zero end;
   is_int: BOOL is return v = #INTI(1) end;

-------------------------------------------------- unary functions

   abs: SAME is r: SAME; return r.u(u.abs).v(v) end;
   negate: SAME is r: SAME; return r.u(-u).v(v) end;
   sign: INT is return u.sign end;
   square: SAME is r: SAME; return r.u(u.square).v(v.square) end;
   cube: SAME is r: SAME; return r.u(u.cube).v(v.cube) end;
   floor: INTI is return u/v end;
   ceiling: INTI is raise "RAT::ceiling not implemented" end;

-------------------------------------------------- miscellaneous

   float (n, b: INT): TUP{BOOL, INTI, INT, BOOL} pre (n > 0) and (b > 1) is
   -- Returns the floating point number x = s * m * b^e that is
   -- closest to self. x is returned in form of a 4-tupel (s, m, e, a).
   -- n specifies the length of the mantissa m in no. of digits to the
   -- base b. The mantissa is either 0 (self = 0) or normalized, i.e.
   -- it holds: b^(n-1) <= mantissa < b^n. The sign s indicates a
   -- negative number. The flag a indicates that x is accurate
   -- (i.e. x = self). Rounding to the nearest is used; ties are
   -- broken by rounding to even (IEEE rounding).
   --
   -- The implementation is essentialy AlgorithmM described by
   -- W.D. Clinger in "How to Read Floating Point Numbers
   -- Accurately", PLDI 1990 proceedings, p. 92-101.
   --
      if u.is_zero then return #(false, #INTI(0), 0, true)
      else
         beta ::= #INTI(b);
         beta_n1 ::= beta^(n-1);
         beta_n ::= beta_n1 * beta;
         s ::= self.u.is_neg;
         u ::= self.u.abs;
         v ::= self.v; assert v.is_pos;
         -- u/v = |self|
         -- approximation for e
         e ::= u.log2 - v.log2;
         if b /= 2 then e := (e.flt / b.flt.log2).round.int end;
         e := e-n;
         if e < 0 then u := u * beta^(-e)
         else v := v * beta^e
         end;
         -- normalize
         -- (usually 0 or 1 iteration steps required)
         m: INTI;
         loop
            -- u/v * b^e = |self|
            m := u/v;
            if m < beta_n1 then u := u*beta; e := e-1
            elsif m >= beta_n then v := v*beta; e := e+1
            else break!
            end
         end;
         -- (u/v * b^e = |self|) & (b^(n-1) <= u/v < b^n)
         -- conversion to float & rounding to nearest (IEEE)
         r ::= u%v; v_r ::= v-r;
         if (r > v_r) or ((r = v_r) and m.is_odd) then -- next float
            m := m + #INTI(1);
            if m = beta_n then m := beta_n1; e := e+1 end
         end;
         return #(s, m, e, r.is_zero)
      end
   end;

-------------------------------------------------- output

   str: STR is
   -- Returns a string representation of self in the form
   -- u/v or u (if self.is_int).
   --
      if is_int then return u.str
      else return u.str + '/' + v.str
      end
   end;

   str (n: INT): STR pre n > 0 is
   -- Returns a string representation of self in decimal floating
   -- point form [sign] digit '.' {digit} ['e' [sign] digit {digit}] where
   -- n specifies the number of digits of the mantissa.
   --
	x ::= float(n, 10); s ::= "";
	if x.t1 then s := s + '-' end;
	m::=x.t2.str;
	digits::=m.length;
	trailing_zeros:INT;
	loop
	    trailing_zeros:=0.upto!(digits-1);
	    while!(m[digits-1-trailing_zeros]='0');
	end;
	s := s + m[0] + '.';
	after_decimal::=digits-1-trailing_zeros;
	if after_decimal>0 then
	    s := s + m.substring(1,after_decimal); 
	end;
	e ::= x.t3 + digits-1;
	if e /= 0 then s := s + 'e' + e end;
	return s
   end;

end
