! RationalArithmetic.f90: rational arithmetic data type and operations (i e, class)

    MODULE RationalArithmetic
        IMPLICIT NONE
        PRIVATE ! the default for all names in this module, except those specified below
        PUBLIC  :: rational, OPERATOR(+), OPERATOR(-), OPERATOR(*), OPERATOR(/), OPERATOR(**), &
          OPERATOR(.NE.), OPERATOR(.EQ.), OPERATOR(.LT.), OPERATOR(.LE.), OPERATOR(.GT.), &
          ASSIGNMENT(=), RationalString, MATMULRat, Reduce


        TYPE rational
!			PRIVATE
            INTEGER :: num,denom
        END TYPE rational

        INTERFACE OPERATOR(+)
            MODULE PROCEDURE AddRational, AddRationalInt
        END INTERFACE

        INTERFACE OPERATOR(-)
            MODULE PROCEDURE SubtractRational
        END INTERFACE

        INTERFACE OPERATOR(*)
            MODULE PROCEDURE MultRational
        END INTERFACE

        INTERFACE OPERATOR(**)
            MODULE PROCEDURE RealToRationalPower
            MODULE PROCEDURE RealToRationalPower16
        END INTERFACE

        INTERFACE OPERATOR(/)
            MODULE PROCEDURE DivRational
        END INTERFACE

        INTERFACE OPERATOR(.EQ.)
            MODULE PROCEDURE EQ_Rational
            MODULE PROCEDURE EQ_RationalInt
        END INTERFACE

        INTERFACE OPERATOR(.NE.)
            MODULE PROCEDURE NE_Rational
            MODULE PROCEDURE NE_RationalInt
        END INTERFACE

        INTERFACE OPERATOR(.LT.)
            MODULE PROCEDURE LT_Rational
            MODULE PROCEDURE LT_RationalInt
        END INTERFACE

        INTERFACE OPERATOR(.LE.)
            MODULE PROCEDURE LE_Rational
            MODULE PROCEDURE LE_RationalInt
        END INTERFACE

        INTERFACE OPERATOR(.GT.)
            MODULE PROCEDURE GT_Rational
            MODULE PROCEDURE GT_RationalInt
        END INTERFACE

        INTERFACE ASSIGNMENT(=)
            MODULE PROCEDURE IntToRational
            MODULE PROCEDURE RationalToRational
            MODULE PROCEDURE RationalToReal
            MODULE PROCEDURE RationalToReal16
            MODULE PROCEDURE RealToRational
        END INTERFACE

        CONTAINS

        ! rat=int
        ELEMENTAL SUBROUTINE IntToRational(rat,int)
!!DEC$ ATTRIBUTES STDCALL :: IntToRational
!!DEC$ ATTRIBUTES REFERENCE :: rat,int
        INTEGER, INTENT(IN) :: int
        TYPE(rational), INTENT(OUT) :: rat
            rat%num=int
            rat%denom=1
        END SUBROUTINE IntToRational

        ! rat=rat
        ELEMENTAL SUBROUTINE RationalToRational(ratL,ratR)
!!DEC$ ATTRIBUTES STDCALL :: RationalToRational
!!DEC$ ATTRIBUTES REFERENCE :: ratL,ratR
        TYPE(rational), INTENT(IN) :: ratR
        TYPE(rational), INTENT(OUT) :: ratL
            ratL%num=ratR%num
            ratL%denom=ratR%denom
        END SUBROUTINE RationalToRational

        ! real=rat
        ELEMENTAL SUBROUTINE RationalToReal(realV,rat)
!!DEC$ ATTRIBUTES STDCALL :: RationalToReal
!!DEC$ ATTRIBUTES REFERENCE :: realV,rat
        TYPE(rational), INTENT(IN) :: rat
        REAL(8), INTENT(OUT) :: realV
            realV=REAL(rat%num,8)/REAL(rat%denom,8)
        END SUBROUTINE RationalToReal

        ! real(16)=rat
        ELEMENTAL SUBROUTINE RationalToReal16(realV,rat)
        TYPE(rational), INTENT(IN) :: rat
        REAL(16), INTENT(OUT) :: realV
            realV=REAL(rat%num,16)/REAL(rat%denom,16)
        END SUBROUTINE RationalToReal16

        ! rat=real: make a real into the closest rational
        ELEMENTAL SUBROUTINE RealToRational(rat,realV)
!!DEC$ ATTRIBUTES STDCALL :: RealToRational
!!DEC$ ATTRIBUTES REFERENCE :: rat,realV
        REAL(8), INTENT(IN) :: realV
        TYPE(rational), INTENT(OUT) :: rat
        REAL(8) :: realVA,fract,x
        INTEGER :: a,b,n,d
        LOGICAL :: neg
        neg=(realV<0)
        realVA=ABS(realV)

        ! Is it already an integer?
        N=FLOOR(realVA); fract=realVA-N
        IF(fract==0.0) THEN
            rat%num=N; rat%denom=1
        ELSE
            ! Try small integers on fractional part-allow slight tolerance
            rat%denom=0 ! "found" flag
    Loop_a: DO a=1,10
                DO b=a+1,10
                    x=DBLE(a)/DBLE(b)
                    IF(ABS(x-fract)<1d-16) THEN
                        rat%num=N*b+a; rat%denom=b
                        EXIT Loop_a
                    END IF
                END DO
            END DO Loop_a

            IF(rat%denom.EQ.0) THEN
                ! Reduce decimal fraction to lowest terms
                d=DIGITS(a) ! significant bits in an integer
                x=SCALE(realVA,d-3)
                rat%num=INT(ANINT(x)); rat%denom=RADIX(a)**(d-3)
                rat=Reduce(rat)
            END IF
        END IF

        IF(neg) rat%num=-rat%num
        END SUBROUTINE RealToRational

        ! rat+int
        TYPE(rational) ELEMENTAL FUNCTION AddRationalInt(rat, int)
!!DEC$ ATTRIBUTES STDCALL :: AddRationalInt
!!DEC$ ATTRIBUTES REFERENCE :: rat,int
        TYPE(rational), INTENT(IN) :: rat
        INTEGER, INTENT(IN) :: int
            AddRationalInt%num = rat%num + int*rat%denom ! a/b + c = a/b + c/1 = a + bc / b.
            AddRationalInt%denom = rat%denom
            ! If a/b is reduced, a+bc will be also, since b divides bc but not a
            !AddRationalInt=Reduce(AddRationalInt)
    END FUNCTION AddRationalInt

    ! rat1+rat2
        TYPE(rational) ELEMENTAL FUNCTION AddRational(rat1,rat2)
!!DEC$ ATTRIBUTES STDCALL :: AddRational
!!DEC$ ATTRIBUTES REFERENCE :: rat1,rat2
        TYPE(rational), INTENT(IN) :: rat1,rat2
            AddRational%num=rat1%num*rat2%denom+rat2%num*rat1%denom
            AddRational%denom=rat1%denom*rat2%denom
            AddRational=Reduce(AddRational)
        END FUNCTION AddRational

        ! rat1-rat2
        TYPE(rational) ELEMENTAL FUNCTION SubtractRational(rat1,rat2)
!!DEC$ ATTRIBUTES STDCALL :: SubtractRational
!!DEC$ ATTRIBUTES REFERENCE :: rat1,rat2
        TYPE(rational), INTENT(IN) :: rat1,rat2
            SubtractRational%num=rat1%num*rat2%denom-rat2%num*rat1%denom
            SubtractRational%denom=rat1%denom*rat2%denom
            SubtractRational=Reduce(SubtractRational)
        END FUNCTION SubtractRational

        ! rat1*rat2
        TYPE(rational) ELEMENTAL FUNCTION MultRational(rat1,rat2)
!!DEC$ ATTRIBUTES STDCALL :: MultRational
!!DEC$ ATTRIBUTES REFERENCE :: rat1,rat2
        TYPE(rational), INTENT(IN) :: rat1,rat2
            MultRational%num=rat1%num*rat2%num
            MultRational%denom=rat1%denom*rat2%denom
            MultRational=Reduce(MultRational)
        END FUNCTION MultRational

        ! real**rat
        REAL(8) ELEMENTAL FUNCTION RealToRationalPower(realBase,ratPower)
!!DEC$ ATTRIBUTES STDCALL :: RealToRationalPower
!!DEC$ ATTRIBUTES REFERENCE :: realBase,ratPower
        REAL(8), INTENT(IN) :: realBase
        TYPE(rational), INTENT(IN) :: ratPower
        TYPE(rational) :: ratR
        ratR=Reduce(ratPower)
        SELECT CASE(ratR%num)
            CASE(0)
                RealToRationalPower=1D0
            CASE(1)
                RealToRationalPower=realBase
            CASE(-1)
                RealToRationalPower=1D0/realBase
            CASE DEFAULT
                RealToRationalPower=realBase**ratR%num
        END SELECT
        SELECT CASE(ratR%denom) ! in reduced form, this is >0
            CASE(1)
                ! does nothing
            CASE(2)
                RealToRationalPower=SQRT(RealToRationalPower)
            CASE DEFAULT
                RealToRationalPower=RealToRationalPower**(1D0/ratR%denom)
        END SELECT
        END FUNCTION RealToRationalPower

        ! real16**rat
        REAL(16) ELEMENTAL FUNCTION RealToRationalPower16(realBase,ratPower)
        REAL(16), INTENT(IN) :: realBase
        TYPE(rational), INTENT(IN) :: ratPower
        TYPE(rational) :: ratR
        !REAL(16) :: x
        ratR=Reduce(ratPower)
        SELECT CASE(ratR%num)
            CASE(0)
                RealToRationalPower16=1.0_16 !Q0
            CASE(1)
                RealToRationalPower16=realBase
            CASE(-1)
                !RealToRationalPower16=1Q0/realBase
                RealToRationalPower16=1.0_16/realBase
            CASE DEFAULT
                RealToRationalPower16=realBase**ratR%num
        END SELECT
        SELECT CASE(ratR%denom) ! in reduced form, this is >0
            CASE(1)
                ! does nothing
            CASE(2)
                RealToRationalPower16=SQRT(RealToRationalPower16)
            CASE DEFAULT
                !RealToRationalPower16=RealToRationalPower16**(1Q0/ratR%denom)
                RealToRationalPower16=RealToRationalPower16**(1.0_16/ratR%denom)
        END SELECT
        END FUNCTION RealToRationalPower16

        ! rat1/rat2
        TYPE(rational) ELEMENTAL FUNCTION DivRational(rat1,rat2)
!!DEC$ ATTRIBUTES STDCALL :: DivRational
!!DEC$ ATTRIBUTES REFERENCE :: rat1,rat2
        TYPE(rational), INTENT(IN) :: rat1,rat2
            DivRational%num=rat1%num*rat2%denom
            DivRational%denom=rat2%num*rat1%denom
            DivRational=Reduce(DivRational)
        END FUNCTION DivRational

        ! rat1.EQ.rat2
        LOGICAL ELEMENTAL FUNCTION EQ_Rational(rat1,rat2)
!!DEC$ ATTRIBUTES STDCALL :: EQ_Rational
!!DEC$ ATTRIBUTES REFERENCE :: rat1,rat2
        TYPE(rational), INTENT(IN) :: rat1,rat2
        TYPE(rational) :: rat1R,rat2R
        rat1R=Reduce(rat1); rat2R=Reduce(rat2) ! put them into canonical form
        EQ_Rational=(rat1R%num.EQ.rat2R%num) .AND. (rat1R%denom.EQ.rat2R%denom)
        END FUNCTION EQ_Rational

        ! rat1.EQ.int
        LOGICAL ELEMENTAL FUNCTION EQ_RationalInt(rat1,int)
!!DEC$ ATTRIBUTES STDCALL :: EQ_RationalInt
!!DEC$ ATTRIBUTES REFERENCE :: rat1,int
        TYPE(rational), INTENT(IN) :: rat1
        INTEGER, INTENT(IN) :: int
        TYPE(rational) :: rat1R
        rat1R=Reduce(rat1) ! put it into canonical form
        EQ_RationalInt=(rat1R%num.EQ.int) .AND. (rat1R%denom.EQ.1)
        END FUNCTION EQ_RationalInt

        ! rat1.NE.rat2
        LOGICAL ELEMENTAL FUNCTION NE_Rational(rat1,rat2)
!!DEC$ ATTRIBUTES STDCALL :: NE_Rational
!!DEC$ ATTRIBUTES REFERENCE :: rat1,rat2
        TYPE(rational), INTENT(IN) :: rat1,rat2
        TYPE(rational) :: rat1R,rat2R
        rat1R=Reduce(rat1); rat2R=Reduce(rat2) ! put them into canonical form
        NE_Rational=(rat1R%num.NE.rat2R%num) .OR. (rat1R%denom.NE.rat2R%denom)
        END FUNCTION NE_Rational

        ! rat1.NE.int
        LOGICAL ELEMENTAL FUNCTION NE_RationalInt(rat1,int)
!!DEC$ ATTRIBUTES STDCALL :: NE_RationalInt
!!DEC$ ATTRIBUTES REFERENCE :: rat1,int
        TYPE(rational), INTENT(IN) :: rat1
        INTEGER, INTENT(IN) :: int
        TYPE(rational) :: rat1R
        rat1R=Reduce(rat1) ! put it into canonical form
        NE_RationalInt=(rat1R%num.NE.int) .OR. (rat1R%denom.NE.1)
        END FUNCTION NE_RationalInt

        ! rat1.LT.rat2
        LOGICAL ELEMENTAL FUNCTION LT_Rational(rat1,rat2)
!!DEC$ ATTRIBUTES STDCALL :: LT_Rational
!!DEC$ ATTRIBUTES REFERENCE :: rat1,rat2
        TYPE(rational), INTENT(IN) :: rat1,rat2
        LT_Rational=(rat1-rat2)<0
        END FUNCTION LT_Rational

        ! rat1.LT.int
        LOGICAL ELEMENTAL FUNCTION LT_RationalInt(rat1,int)
!!DEC$ ATTRIBUTES STDCALL :: LT_RationalInt
!!DEC$ ATTRIBUTES REFERENCE :: rat1,int
        TYPE(rational), INTENT(IN) :: rat1
        INTEGER, INTENT(IN) :: int
        TYPE(rational) :: rat2,rat3
        rat2%num=int; rat2%denom=1
        rat3=rat1-rat2
        LT_RationalInt=rat3%num<0
        END FUNCTION LT_RationalInt

        ! rat1.LE.rat2
        LOGICAL ELEMENTAL FUNCTION LE_Rational(rat1,rat2)
!!DEC$ ATTRIBUTES STDCALL :: LE_Rational
!!DEC$ ATTRIBUTES REFERENCE :: rat1,rat2
        TYPE(rational), INTENT(IN) :: rat1,rat2
        LE_Rational=(rat1-rat2)<=0
        END FUNCTION LE_Rational

        ! rat1.LE.int
        LOGICAL ELEMENTAL FUNCTION LE_RationalInt(rat1,int)
!!DEC$ ATTRIBUTES STDCALL :: LE_RationalInt
!!DEC$ ATTRIBUTES REFERENCE :: rat1,int
        TYPE(rational), INTENT(IN) :: rat1
        INTEGER, INTENT(IN) :: int
        TYPE(rational) :: rat2,rat3
        rat2%num=int; rat2%denom=1
        rat3=rat1-rat2
        LE_RationalInt=rat3%num<=0
        END FUNCTION LE_RationalInt

        ! rat1.GT.rat2
        LOGICAL ELEMENTAL FUNCTION GT_Rational(rat1,rat2)
!!DEC$ ATTRIBUTES STDCALL :: GT_Rational
!!DEC$ ATTRIBUTES REFERENCE :: rat1,rat2
        TYPE(rational), INTENT(IN) :: rat1,rat2
        GT_Rational=(rat1-rat2)>0
        END FUNCTION GT_Rational

        ! rat1.GT.int
        LOGICAL ELEMENTAL FUNCTION GT_RationalInt(rat1,int)
!!DEC$ ATTRIBUTES STDCALL :: GT_RationalInt
!!DEC$ ATTRIBUTES REFERENCE :: rat1,int
        TYPE(rational), INTENT(IN) :: rat1
        INTEGER, INTENT(IN) :: int
        TYPE(rational) :: rat2,rat3
        rat2%num=int; rat2%denom=1
        rat3=rat1-rat2
        GT_RationalInt=rat3%num>0
        END FUNCTION GT_RationalInt

        ! Matrix multiplication
        FUNCTION MATMULRat(rat1,rat2)
!!DEC$ ATTRIBUTES STDCALL :: MATMULRat
!!DEC$ ATTRIBUTES REFERENCE :: rat1,rat2
        TYPE(rational), DIMENSION(:,:), INTENT(IN) :: rat1,rat2
        TYPE(rational), DIMENSION(SIZE(rat1,1),SIZE(rat2,2)) :: MATMULRat ! the function is array-valued
        TYPE(rational) :: sum
        INTEGER :: i,j,k,N
        N=SIZE(rat1,2)
        IF(N.NE.SIZE(rat2,1)) THEN
            MATMULRat=0 ! non-conforming matrices cannot be multiplied
            RETURN
        END IF
        DO i=1,SIZE(rat1,1)
            DO j=1,SIZE(rat2,2)
                sum=0
                DO k=1,N
                    sum=sum+rat1(i,k)*rat2(k,j)
                END DO
                MATMULRat(i,j)=Reduce(sum)
            END DO
        END DO
        END FUNCTION MATMULRat

        ! Reduce to lowest terms and put into canonical form
        TYPE(rational) ELEMENTAL FUNCTION Reduce(rat)
!!DEC$ ATTRIBUTES STDCALL :: Reduce
!!DEC$ ATTRIBUTES REFERENCE :: rat
        TYPE(rational), INTENT(IN) :: rat
        INTEGER :: a,b,t
        IF(rat%num==0) THEN
            Reduce%num=0; Reduce%denom=1
        ELSE
            ! Find LCD
            a=ABS(rat%num); b=ABS(rat%denom)
            IF(a<b) THEN
                t=a; a=b; b=t
            END IF
            DO
                IF(b<=0) EXIT
                a=MOD(a,b)
                IF(a<b) THEN
                    t=a; a=b; b=t
                END IF
            END DO
            ! a=LCD (a>0)
            Reduce%num=rat%num/a; Reduce%denom=rat%denom/a
            IF(Reduce%denom<0) THEN
                Reduce%denom=-Reduce%denom ! make denominator +
                Reduce%num=-Reduce%num
            END IF
        END IF
        END FUNCTION Reduce

        ! Display rational as a string
        CHARACTER(LEN=25) ELEMENTAL FUNCTION RationalString(rat)
!!DEC$ ATTRIBUTES STDCALL :: RationalString
!!DEC$ ATTRIBUTES REFERENCE :: rat
        TYPE(rational), INTENT(IN) :: rat
        TYPE(rational) :: ratR
        ratR=Reduce(rat) ! get canonical form
        IF(ratR%denom==1) THEN
            WRITE(RationalString,"(I0)") ratR%num
        ELSE
            WRITE(RationalString,"(I0,'/',I0)") ratR%num,ratR%denom
        END IF
        END FUNCTION RationalString

    END MODULE RationalArithmetic