! dzielenie, odejmowanie, minus, **, sign, sin, ..., ==, !=

module ARITHMETIC

  implicit none

  Type Float
     Real :: f
  end type Float

  Integer :: s, m_low, m_high, baza_do_s
  Logical :: chopped ! true = chopping, false = rounding
  Real    :: u_error

  interface assignment (=)
     module procedure a_eq_a, a_eq_r, r_eq_a, i_eq_a, a_eq_i
  end interface assignment (=)

  interface operator (+)
     module procedure add_aa,add_ra,add_ar,add_ai,add_ia
  end interface assignment (+)

  interface operator (*)
     module procedure mul_aa,mul_ra,mul_ar,mul_ai,mul_ia
  end interface assignment (*)

  interface operator (/)
     module procedure div_aa,div_ia,div_ai
  end interface assignment (/)


CONTAINS
  
  ! (=)

  subroutine a_eq_a(z,x)
    Type(Float), intent(out) :: z
    Type(Float), intent(in) :: x
    z%f = x%f
    return
  end subroutine a_eq_a
  subroutine a_eq_r(z,x)
    Type(Float), intent(out) :: z
    Real, intent(in) :: x
    z%f = fl(x)
    return
  end subroutine a_eq_r
  subroutine r_eq_a(z,x)
    Real, intent(out) :: z
    Type(Float), intent(in) :: x
    z = x%f
    return
  end subroutine r_eq_a
  subroutine a_eq_i(z,x)
    Type(Float), intent(out) :: z
    Integer, intent(in) :: x
    z%f = fl(real(x))
    return
  end subroutine a_eq_i
  subroutine i_eq_a(z,x)
    Integer, intent(out) :: z
    Type(Float), intent(in) :: x
    z = int(x%f)
    return
  end subroutine i_eq_a


  ! (+) wynik Float

  function add_aa(x,y) Result(z)
    Type(Float), intent(in) :: x
    Type(Float), intent(in) :: y
    Type(Float) :: z
    z = fl(x%f + y%f)
    return
  end function add_aa

  function add_ra(x,y) Result(z)
    Real, intent(in) :: x
    Type(Float), intent(in) :: y
    Type(Float) :: z
    z = fl(fl(x) + y%f)
    return
  end function add_ra

  function add_ar(x,y) Result(z)
    Real, intent(in) :: y
    Type(Float), intent(in) :: x
    Type(Float) :: z
    z = fl(x%f + fl(y))
    return
  end function add_ar

  function add_ai(x,i) Result(z)
    Integer, intent(in) :: i
    Type(Float), intent(in) :: x
    Type(Float) :: z
    z = fl(x%f + fl(real(i)))
    return
  end function add_ai

  function add_ia(i,y) Result(z)
    Integer, intent(in) :: i
    Type(Float), intent(in) :: y
    Type(Float) :: z
    z = fl(fl(real(i)) + y%f)
    return
  end function add_ia


  ! (*) wynik Float

  function mul_aa(x,y) Result(z)
    Type(Float), intent(in) :: x
    Type(Float), intent(in) :: y
    Real :: z
    z = fl(x%f * y%f)
    return
  end function mul_aa

  function mul_ra(x,y) Result(z)
    Real, intent(in) :: x
    Type(Float), intent(in) :: y
    Real :: z
    z = fl(fl(x) * y%f)
    return
  end function mul_ra

  function mul_ar(x,y) Result(z)
    Real, intent(in) :: y
    Type(Float), intent(in) :: x
    Real :: z
    z = fl(x%f * fl(y))
    return
  end function mul_ar

  function mul_ai(x,y) Result(z)
    Integer, intent(in) :: y
    Type(Float), intent(in) :: x
    Real :: z
    z = fl(x%f * fl(real(y)))
    return
  end function mul_ai

  function mul_ia(x,y) Result(z)
    Integer, intent(in) :: x
    Type(Float), intent(in) :: y
    Real :: z
    z = fl(y%f * fl(real(x)))
    return
  end function mul_ia


  ! (/) Wynik = float

  function div_aa(x,y) Result(z)
    Type(Float), intent(in) :: x,y
    Real :: z
    if (y%f==0) then
       stop "Error w div_aa: Dzielenie przez 0!"
    end if
    z = fl(x%f/y%f)
    return
  end function div_aa

  function div_ia(x,y) Result(z)
    Integer, intent(in) :: x
    Type(Float), intent(in) :: y
    Real :: z
    if (y%f==0) then
       stop "Error w div_ia: Dzielenie przez 0!"
    end if
    z = fl(fl(real(x))/y%f)
    return
  end function div_ia

  function div_ai(x,y) Result(z)
    Integer, intent(in) :: y
    Type(Float), intent(in) :: x
    Real :: z
    if (y==0) then
       stop "Error w div_ai: Dzielenie przez 0!"
    end if
    z = fl(x%f/fl(real(y)))
    return
  end function div_ai


  subroutine Arytmetyka(mantysa, mdolne, mgorne, chopping)
    implicit none
    Integer, intent(in) :: mantysa, mdolne, mgorne
    Logical, intent(in) :: chopping
    s         = mantysa
    m_low     = mdolne
    m_high    = mgorne
    baza_do_s = 10**s
    chopped   = chopping
    u_error   = 1./baza_do_s
    if (.not. chopping) then 
       u_error = 0.5*u_error 
    end if
    return
  end subroutine Arytmetyka

  function fl(x)
    implicit none
    Real :: fl
    Real, intent(in) :: x
    Integer :: norm, sigx
    Real :: absx
    absx = abs(x)
    sigx = sign(1.,x)
    norm = 10**(floor(log10(absx)) + 1)
    if (chopped) then
       fl = floor( (absx/norm)*baza_do_s )
    else
       fl = floor( ((absx/norm+u_error))*baza_do_s )
    end if
    fl = sigx*real(norm)/baza_do_s * fl
    return
  end function fl

end module ARITHMETIC

program ary

  use ARITHMETIC

  Integer :: mantysa, mdolne, mgorne
  Type(Float) :: xf, yf
  Real :: x, y
  Logical :: chopping

  mantysa = 4
  mdolne  = 0
  mgorne  = 0

  chopping = .true.
  call arytmetyka(mantysa, mdolne, mgorne, chopping)
  x  = 3.25676289
  xf = x
  write(*,*) u_error, 3*xf, fl(xf*3), xf/xf

  chopping = .false.
  call arytmetyka(mantysa, mdolne, mgorne, chopping)
  x  = -x
  xf = x
  write(*,*) u_error, 3*xf, fl(xf*3), xf/xf

end program ary

