!  Rational arithmetic module
!  Defines basic integer type, then
!  rational operations, =,+,-,*,/
!  Ralph Frisbie
!  rfrisbie@vcnet.com
!  28 July 1996
! Further mods and logicals added
!  20 October 1996
!

module ratiomod

public ::  assignment (=), &
             operator (+), &
             operator (-), &
             operator (*), &
             operator (/), &
             operator (**), &
             operator (==), &
             operator (/=), &
             operator (<=), &
             operator (>=), &
             operator (>), &
             operator (<), &
             GCD,  reduce, &
             ratin, ratout

private ::   rat_equals_int, rat_equals_rat, &
             array_equals_rat,rat_equals_array, &
             rat_plus_rat, rat_plus_int,int_plus_rat, &
             rat_minus_rat,rat_minus_int,int_minus_rat, &
             rat_times_rat,rat_times_int,int_times_rat, &
             rat_div_rat,rat_div_int,int_div_rat, &
             rat_to_int, rat_eq_rat,rat_ne_rat,&
             rat_le_rat,rat_ge_rat,rat_gt_rat,rat_lt_rat

integer, parameter, public  ::  rat_kind = selected_int_kind(9)

type, public  :: rational
      integer (kind=rat_kind)  :: numerator
      integer (kind=rat_kind)  :: denominator
end type rational

type (rational), public, parameter   :: &
                       one =  rational (1_rat_kind,1_rat_kind )
type (rational), public, parameter   :: &
                        zero = rational( 0_rat_kind,1_rat_kind )


interface assignment (=)
     module procedure rat_equals_int, &
                      rat_equals_rat, &
                      array_equals_rat, &
                      rat_equals_array
end interface

interface operator (+)
     module procedure   rat_plus_rat, &
                        rat_plus_int, &
                        int_plus_rat
end interface

interface operator (-)
      module procedure  rat_minus_rat, &
                        rat_minus_int, &
                        int_minus_rat
end interface

interface operator (*)
     module procedure   rat_times_rat, &
                        rat_times_int, &
                        int_times_rat
end interface

interface operator (/)
     module procedure   rat_div_rat, &
                        rat_div_int, &
                        int_div_rat
end interface

interface operator (**)
     module procedure  rat_to_int
end interface

interface operator (==)
     module procedure  rat_eq_rat
end interface

interface operator (/=)
     module procedure  rat_ne_rat
end interface

interface operator (>=)
     module procedure  rat_ge_rat
end interface

interface operator (<=)
     module procedure  rat_le_rat
end interface

interface operator (<)
     module procedure  rat_lt_rat
end interface

interface operator (>)
     module procedure  rat_gt_rat
end interface


contains

  subroutine rat_equals_int(r,i)
     type (rational),intent (out) :: r
     integer (kind=rat_kind),intent(in)  :: i
     r = rational( i,1_rat_kind)
  end subroutine rat_equals_int


  subroutine rat_equals_rat (r1,r2)
    type (rational), intent(out)   :: r1
    type (rational), intent(in)   :: r2
      r1%numerator = r2%numerator
      r1%denominator = r2%denominator
  end subroutine rat_equals_rat

  subroutine array_equals_rat (ar,r)
    integer (kind=rat_kind), intent(out),dimension(2)         :: ar
    type (rational),intent(in)   :: r
    ar = (/ r%numerator,r%denominator /)
  end subroutine array_equals_rat


  subroutine rat_equals_array(r,ar)
    type (rational),intent(out)  :: r
    integer (kind=rat_kind), dimension(2), intent(in)  :: ar
      r = rational( ar(1),ar(2))
  end subroutine rat_equals_array

  function  rat_plus_rat (r1,r2) result (summ)
    type (rational),intent(in)    :: r1,r2
    type (rational)           :: summ
    integer (kind=rat_kind)       :: num, denom
    num = r1%numerator*r2%denominator + r1%denominator*r2%numerator
    denom = r1%denominator*r2%denominator
    summ = rational( num,denom)
  end function rat_plus_rat


  function rat_plus_int(r,i) result(summ)
   type (rational),intent(in)   :: r
   integer (kind=rat_kind),intent(in) :: i
   type (rational)         :: summ
   summ = rational( r%numerator + i,r%denominator )
  end function rat_plus_int

  function int_plus_rat(i,r) result(summ)
   type (rational), intent(in)    :: r
   integer (kind=rat_kind),intent(in)  :: i
   type (rational)       :: summ
    summ = rational(r%numerator + i*r%denominator, r%denominator)
   end function int_plus_rat

  function  rat_minus_rat (r1,r2) result (diff)
   type (rational),intent(in)    :: r1,r2
   type (rational)              :: diff
     diff = rational(r1%numerator*r2%denominator -r2%numerator*r1%denominator, &
                      r1%denominator * r2%denominator)
   end function rat_minus_rat

  function rat_minus_int(r,i) result(diff)
   type (rational),intent(in)  :: r
   integer (kind=rat_kind),intent(in)  :: i
   type (rational)                :: diff
    diff = r + (-i)
  end function rat_minus_int

  function int_minus_rat(i,r) result(diff)
   type (rational),intent(in)  :: r
   integer (kind=rat_kind),intent(in)  :: i
   type (rational)                   :: diff
     diff = rational(i*r%denominator - r%numerator, r%denominator)
  end function int_minus_rat

  function  rat_times_rat(r1,r2) result(prod)
   type (rational), intent(in)  :: r1,r2
   type (rational)                  :: prod
   prod = rational(r1%numerator*r2%numerator, r1%denominator*r2%denominator)
  end function rat_times_rat

  function rat_times_int(r,i) result(prod)
   type (rational), intent(in) :: r
   integer (kind=rat_kind), intent(in)  :: i
   type (rational)                 :: prod
   prod = rational(i*r%numerator,r%denominator)
  end function rat_times_int

  function int_times_rat(i,r) result(prod)
   type (rational), intent(in) :: r
   integer (kind=rat_kind), intent(in) :: i
   type (rational)                 :: prod
   prod = r*i
  end function int_times_rat

  function  rat_div_rat(r1,r2) result(quot)
   type (rational), intent(in) :: r1,r2
   type (rational)                 :: quot
    quot = r1*rational(r2%denominator,r2%numerator)
  end function rat_div_rat

  function rat_div_int(r,i) result(quot)
   type(rational), intent(in)  :: r
   integer (kind=rat_kind), intent(in) :: i
   type (rational)                :: quot
    quot = rational(r%numerator,i*r%denominator)
  end function rat_div_int

  function int_div_rat(i,r) result(quot)
   type(rational), intent(in)  :: r
   integer (kind=rat_kind), intent(in) :: i
   type (rational)                 :: quot
    quot = rational(i*r%denominator,r%numerator)
  end function int_div_rat

  function rat_to_int(r,i) result(pow)
    type(rational), intent(in)  :: r
    integer, intent(in) :: i
    type(rational)                  :: pow
    integer                         :: absi
    absi = abs(i)
    if( i >= 0) then
      pow%numerator = r%numerator**absi
      pow%denominator = r%denominator**absi
    else
      pow%denominator = r%numerator**absi
      pow%numerator = r%denominator**absi 
    end if
  end function rat_to_int

  function rat_eq_rat(r1,r2) result(equ)
    type(rational),intent(in)  :: r1,r2
    logical                    :: equ
    type (rational)            :: x1,x2
     x1 = reduce(r1)
     x2 = reduce(r2)
     equ = x1%numerator == x2%numerator .and. x1%denominator == x2%denominator
   end function rat_eq_rat

   function rat_ne_rat(r1,r2) result(nequ)
     type (rational), intent (in)  :: r1,r2
     logical                       :: nequ
     nequ = .not. rat_eq_rat(r1,r2)
    end function rat_ne_rat

   function rat_le_rat(r1,r2) result(ler)
     type (rational),intent(in)   :: r1,r2
     logical                       :: ler
     type (rational)              :: x1,x2
       x1 = reduce(r1)
       x2 = reduce(r2)
       ler = x1%numerator*x2%denominator <= x2%numerator * x1%denominator
    end function rat_le_rat

   function rat_lt_rat(r1,r2) result(ler)
     type (rational),intent(in)   :: r1,r2
     logical                       :: ler
     type (rational)              :: x1,x2
       x1 = reduce(r1)
       x2 = reduce(r2)
       ler = x1%numerator * x2%denominator < x2%numerator * x1%denominator
    end function rat_lt_rat

   function rat_gt_rat(r1,r2) result(ler)
     type (rational),intent(in)   :: r1,r2
     logical                       :: ler
     type (rational)              :: x1,x2
       x1 = reduce(r1)
       x2 = reduce(r2)
       ler = x1%numerator * x2%denominator > x2%numerator * x1%denominator
    end function rat_gt_rat

   function rat_ge_rat(r1,r2) result(ler)
     type (rational),intent(in)   :: r1,r2
     logical                       :: ler
     type (rational)              :: x1,x2
       x1 = reduce(r1)
       x2 = reduce(r2)
       ler = x1%numerator * x2%denominator >= x2%numerator * x1%denominator
    end function rat_ge_rat

  recursive function GCD(i1,i2) result(ans)
   integer (kind=rat_kind) ,intent(in)    :: i1,i2
   integer (kind=rat_kind)        :: ans
   integer (kind=rat_kind)        :: work

     work = modulo(i1,i2)
   if(work == 0) then
      ans = i2
   else
      ans = GCD(i2,work)
   end if
  end function GCD

  function reduce(r) result(red)
   type (rational), intent(in)   :: r
   type (rational)               :: red
   integer (kind=rat_kind)       :: mult
     mult = GCD(abs(r%numerator),abs(r%denominator))
    red = rational(r%numerator/mult,r%denominator/mult)
    if(red%denominator < 0) then
       red%numerator = -red%numerator
    end if
  end function reduce

  subroutine ratin(r)
   type (rational), intent(out)  :: r
   character (len=80)           :: rin
   integer                      :: slash_loc
   integer (kind=rat_kind)      :: numer,denom

   read "(a)",rin
   slash_loc = index(rin,"/")
   read(unit=rin(1:slash_loc-1),fmt=*) numer
   if(slash_loc > 1) then
     read(unit=rin(slash_loc+1:),fmt=*) denom
   else
     denom = 1_rat_kind
   end if
   r = rational(numer,denom)
  end subroutine ratin

  subroutine ratout(r,rout)
   type(rational), intent(in)   :: r
   character (len=*), intent(in out) :: rout
   integer                      :: numchar
   character (len=36)           :: cnum,cden
   character (len=73)           :: temp

   write(unit=cnum,fmt=*) r%numerator
   write(unit=cden,fmt=*) r%denominator
   temp = trim(adjustl(cnum))//"/"//trim(adjustl(cden))
   numchar = len(rout)
   if(len(trim(temp)) > numchar) then
      rout = repeat("*",numchar)
   else
      rout = temp
      rout = adjustr(rout)
   end if
  end subroutine ratout


end module ratiomod


program testrat
use ratiomod

type (rational)   :: a,b,c,diff,prod,quot

integer (kind=rat_kind)  :: i,j
integer                  :: k
character (len=79)       :: line,head
k = 1
head = &
"    a     b      a+b   (a+b)      a-b    (a-b)        a*b    (a*b)      a/b "
do i = 1,4
   do j = 1,10
     if(modulo(k,20) == 0) then
       read *
       print "(a)",head
     end if
     a = rational( i,j )
     b = rational( 2*i - 1,2*j+5 )
     c = a + b
     diff = a - b
     prod = a*b

     quot =a/b
     line = " "
     call ratout(a,line(1:5))
!     line(6:6) = "|"
     call ratout(b,line(7:11))
!     line(12:12) = "|"
     call ratout(c,line(13:20))
!     line(21:21) = "|"
     call ratout(reduce(c),line(22:28))
     call ratout(diff,line(30:38))
     call ratout(reduce(diff),line(40:47))
     call ratout(prod,line(49:58))
     call ratout(reduce(prod),line(59:67))
     call ratout(quot,line(68:76))
     print "(a)", line
     k = k + 1

   end do
end do
   a = rational(2,3)
   b = a**3
   quot = a**(-3)
   line = " "
   call ratout(a,line(1:6))
   call ratout(b,line(8:12))
   call ratout(quot,line(14:20))
   call ratout(reduce(b*quot),line(22:26))
   line(27:)= " for 2/3,(2/3)**3, **-3, and product"
   print "(a)",line
   print *, (a==b),(a<b),(a>b)
   do j = 1,4
     a = rational(1,j)
     b = rational(j,4)
     print *,j,a<b,a<=b,a==b,a>=b, a>b, a/=b
   end do


end program testrat
