
/* rational.q: rational numbers
   $Id: rational.q,v 1.15 2007/09/27 10:09:06 agraef Exp $ */

/* This is based on Rob Hubbard's "Q+Q" rational number module (v. 160).
   Copyright (C) 2005-2006 by Rob Hubbard */

/* This file is part of the Q programming system.

   The Q programming system is free software; you can redistribute it and/or
   modify it under the terms of the GNU General Public License as published by
   the Free Software Foundation; either version 2, or (at your option)
   any later version.

   The Q programming system 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 the
   GNU General Public License for more details.

   You should have received a copy of the GNU General Public License
   along with this program; if not, write to the Free Software
   Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */

include clib, complex, stdlib;

/* The Rational type. The standard arithmetic operators, abs, sgn, round,
   trunc, int, frac, and the pow function from clib are overloaded on this
   type. Float conversion is also supported, so that all Real functions
   defined in prelude.q work, too. Access functions for numerator and
   denominator and an exact division operator (%) are also provided. The
   latter is declared as a virtual constructor of the type which can be used
   in pattern-matching. All operations work consistently on both Rational and
   Int arguments, assuming that an integer represents a rational value with
   denominator 1.

   The exact division division operator (%) returns a rational or complex
   rational for any combination of integer, rational and complex integer/
   rational arguments. For other Num arguments it works just like (/). This
   can be used to construct rational numbers, and is also employed as the
   default view of rational numbers in the format `N%D', where N and D are the
   (normalised) numerator and denominator. */

public type Rational : Real = virtual (%) P Q @ (/) | private const rat N D;

@-0x80000000
view (rat N:Int D:Int)		= '(N % D);
@0

/* Construction: The rational function constructs a Rational from either an
   Int or a pair of Int values (numerator,denominator). */

public rational X;

/* Deconstruction: The num and den functions return the normalised numerator
   and denominator of a rational number, respectively (i.e., numerator and
   denominator are coprime, and the sign is always stored in the numerator).
   The num_den function returns both numerator and denominator as a pair of
   Int values. */

public num Q, den Q, num_den Q;

/* Implementation. *********************************************************/

private mkrat N D;

rational (N:Int,D:Int)		= mkrat N D if D<>0;
rational N:Int			= rat N 1;

num (rat N:Int _:Int)		= N;
num N:Int			= N;

den (rat _:Int D:Int)		= D;
den _:Int			= 1;

num_den (rat N:Int D:Int)	= (N,D);
num_den N:Int			= (N,1);

float (rat N:Int D:Int)		= N/D if D<>0;

/* Construct a normalised Rational value. Optimize the cases with zero
   numerators and unit numerators or denominators. Also handle the case of a
   zero denominator (this shouldn't happen during normal operation) by just
   returning the number as is. */

mkrat N:Int 0			= rat N 0;
mkrat 0 D:Int			= rat 0 1;
mkrat N:Int 1			= rat N 1;
mkrat 1 D:Int			= rat (sgn D) (abs D);
mkrat N:Int (-1)		= rat (-N) 1;
mkrat (-1) D:Int		= rat (-sgn D) (abs D);
mkrat N:Int D:Int		= rat N D
  where G = gcd N D, S = sgn D*sgn G, G = S*G, N = N div G, D = D div G;

/* The following definitions take priority over the stuff in prelude.q. */

/* Basic rounding and conversion functions. Note that floor and ceil are taken
   care of in math.q. */

trunc (rat N:Int D:Int)		= N div D if D<>0;
round (rat N:Int D:Int)		= -round (rat (-N) D) if N<0;
				= N div 2+1 if D=2;
				= (2*N+D) div (2*D) if D<>0 otherwise;
int Q:Rational			= float IP where IP:Int = trunc Q;
frac Q:Rational			= Q-IP where IP:Float = int Q;

/* Absolute value and sign. Note that the sign is always kept in the numerator
   of the normalised representation. */

abs (rat N:Int D:Int)		= rat (abs N) D;
sgn (rat N:Int _:Int)		= sgn N;

/* Extend the pow function to negative integer exponents and rational values
   in the first argument. */

pow (rat N:Int D:Int) I:Int	= rat 1 1 if (I=0) and then (N<>0);
				= mkrat (pow N I) (pow D I) if I>0;
				= mkrat (pow D (-I)) (pow N (-I))
				    if (I<0) and then (N<>0);
pow N:Int I:Int			= rat 1 1 if (I=0) and then (N<>0);
				= mkrat 1 (pow N (-I))
				    if (I<0) and then (N<>0);

/* Exact division operator. */

rat N1:Int D1:Int % rat N2:Int D2:Int
				= mkrat (N1*D2) (D1*N2) if N2<>0;
rat N1:Int D1:Int % N2:Int	= mkrat N1 (D1*N2) if N2<>0;
N1:Int % rat N2:Int D2:Int	= mkrat (N1*D2) N2 if N2<>0;
P:Int % Q:Int			= mkrat P Q if Q<>0;

Z1:Complex % Z2:Complex		= ((X1*X2+Y1*Y2)%(X2*X2+Y2*Y2)) +:
				  ((Y1*X2-X1*Y2)%(X2*X2+Y2*Y2))
				    where X1:Real+:Y1:Real = Z1,
				      X2:Real+:Y2:Real = Z2;
X1:Real % Z2:Complex		= (X1+:0) % Z2;
Z1:Complex % X2:Real		= Z1 % (X2+:0);

P:Num % Q:Num			= R where R:Num = P/Q otherwise;

/* Standard arithmetic operators. Note that (/) and (^) are not defined here,
   instead we rely on the default implementations of the prelude which always
   give inexact results. Use (%) and `pow' with integer exponents if you want
   exact results. */

-rat N:Int D:Int			= rat (-N) D;

rat N1:Int D1:Int + rat N2:Int D2:Int	= mkrat (N1*D2+N2*D1) (D1*D2);
rat N1:Int D1:Int - rat N2:Int D2:Int	= mkrat (N1*D2-N2*D1) (D1*D2);
rat N1:Int D1:Int * rat N2:Int D2:Int	= mkrat (N1*N2) (D1*D2);

/* Mixed rational/int arithmetic. */

rat N1:Int D1:Int + N2:Int		= mkrat (N1+N2*D1) D1;
rat N1:Int D1:Int - N2:Int		= mkrat (N1-N2*D1) D1;
rat N1:Int D1:Int * N2:Int		= mkrat (N1*N2) D1;

N1:Int + rat N2:Int D2:Int		= mkrat (N1*D2+N2) D2;
N1:Int - rat N2:Int D2:Int		= mkrat (N1*D2-N2) D2;
N1:Int * rat N2:Int D2:Int		= mkrat (N1*N2) D2;

/* Comparisons. */

(rat N1:Int D1:Int = rat N2:Int D2:Int)	= (N1*D2 = N2*D1);
(rat N1:Int D1:Int <> rat N2:Int D2:Int)= (N1*D2 <> N2*D1);
(rat N1:Int D1:Int < rat N2:Int D2:Int)	= (N1*D2 < N2*D1);
(rat N1:Int D1:Int <= rat N2:Int D2:Int)= (N1*D2 <= N2*D1);
(rat N1:Int D1:Int > rat N2:Int D2:Int)	= (N1*D2 > N2*D1);
(rat N1:Int D1:Int >= rat N2:Int D2:Int)= (N1*D2 >= N2*D1);

/* Mixed rational/int comparisons. */

(rat N1:Int D1:Int = N2:Int)		= (N1 = N2*D1);
(rat N1:Int D1:Int <> N2:Int)		= (N1 <> N2*D1);
(rat N1:Int D1:Int < N2:Int)		= (N1 < N2*D1);
(rat N1:Int D1:Int <= N2:Int)		= (N1 <= N2*D1);
(rat N1:Int D1:Int > N2:Int)		= (N1 > N2*D1);
(rat N1:Int D1:Int >= N2:Int)		= (N1 >= N2*D1);

(N1:Int = rat N2:Int D2:Int)		= (N1*D2 = N2);
(N1:Int <> rat N2:Int D2:Int)		= (N1*D2 <> N2);
(N1:Int < rat N2:Int D2:Int)		= (N1*D2 < N2);
(N1:Int <= rat N2:Int D2:Int)		= (N1*D2 <= N2);
(N1:Int > rat N2:Int D2:Int)		= (N1*D2 > N2);
(N1:Int >= rat N2:Int D2:Int)		= (N1*D2 >= N2);
