module big_integers_module
   implicit none
   public :: print_big
   public :: prime
   public :: assignment (=)
   public :: operator (+)
   public :: operator (-)
   public :: operator (*)
   public :: operator (/)
   public :: operator (**)
   public :: modulob
   public :: huge
   public :: sqrtb
   public :: operator (==)
   public :: operator (<=)
   public :: operator (<)
   private ::  big_gets_int, &
               big_gets_char, &
               f_gets_big, &
               big_gets_f, &
               big_plus_int, &
               big_plus_big, &
               f_plus_f, &
               big_minus_int, &
               big_minus_big, &
               f_minus_f, &
               big_times_int, &
               big_times_big, &
               f_times_int, &
               f_times_f, &
               big_div_int, &
               big_div_big, &
               big_power_int, &
               modulo_big_int, &
               modulo_big_big, &
               big_eq_int, &
               big_eq_big, &
               big_le_int, &
               big_le_big, &
               f_le_f, &
               big_lt_int, &
               big_lt_big, &
               huge_big, &
               print_f, &
               f_base_to_power, &
               big_base_to_power, &
               f_reciprocal, &
               print_factors_big, &
               print_big_base, &
               sqrt_big

!!!   intrinsic modulo
   intrinsic huge
!!!   intrinsic sqrt

   integer, parameter, private :: &
         d = (digits (0) - 1) / 2, &
         r = radix (0), &
         base = r ** d, &
         nr_of_digits = 20, &
         f_nr_of_digits = 2 * nr_of_digits, &
         nr_of_decimal_digits = 100

   ! base of number system is r ** d,
   ! so that each "digit" is 0 to r**d - 1

   type, public :: big_integer
      private
      integer, dimension (0 : nr_of_digits) &
            :: digit
   end type big_integer

   type, public :: big_fraction
      private
      integer, dimension (-f_nr_of_digits : nr_of_digits) &
            :: digit
   end type big_fraction

   interface assignment (=)
      module procedure big_gets_int, &
                       big_gets_char, &
                       f_gets_big, &
                       big_gets_f
   end interface

   interface operator (+)
      module procedure big_plus_int, &
                       big_plus_big, &
                       f_plus_f
   end interface

   interface operator (-)
      module procedure big_minus_int, &
                       big_minus_big, &
                       f_minus_f
   end interface

   interface operator (*)
      module procedure big_times_int, &
                       big_times_big, &
                       f_times_int, &
                       f_times_f
   end interface

   interface operator (/)
      module procedure big_div_int, &
                       big_div_big
   end interface

   interface operator (**)
      module procedure big_power_int
   end interface

   interface modulob
      module procedure modulo_big_int, &
                       modulo_big_big
   end interface

   interface operator (==)
      module procedure big_eq_int, &
                       big_eq_big
   end interface

   interface operator (<=)
      module procedure big_le_int, &
                       big_le_big, &
                       f_le_f
   end interface

   interface operator (<)
      module procedure big_lt_int, &
                       big_lt_big
   end interface

   interface huge
      module procedure huge_big
   end interface

   interface sqrtb
      module procedure sqrt_big
   end interface

   type (big_fraction), private :: eps

contains

subroutine big_gets_int (b, i)

   type (big_integer), intent (out) :: b
   integer, intent (in) :: i

   if (i < 0) then
      print *, "attempt to assign negative number"
      stop
   end if
   b % digit (0) = modulo (i, base)
   b % digit (1) = i / base
   b % digit (2:) = 0

end subroutine big_gets_int

subroutine big_gets_char (b, c)

   type (big_integer), intent (out) :: b
   character (len=*), intent (in) :: c
   integer :: temp_digit, n

   b % digit = 0
   do n = 1, len (c)
      temp_digit = index ("0123456789", c (n:n)) - 1
      if (temp_digit < 0) then
         print *, "attempt to assign nonnumeric string"
         stop
      end if
      b = b * 10 + temp_digit
   end do

end subroutine big_gets_char

subroutine f_gets_big (f, b)

   type (big_fraction), intent (out) :: f
   type (big_integer), intent (in) :: b

   f % digit (-f_nr_of_digits : -1) = 0
   f % digit (0 : nr_of_digits) = b % digit

end subroutine f_gets_big

subroutine big_gets_f (b, f)

   type (big_integer), intent (out) :: b
   type (big_fraction), intent (in) :: f

   b % digit = f % digit (0 : nr_of_digits)

end subroutine big_gets_f

function big_base_to_power (n)  result (b)

   integer, intent (in) :: n
   type (big_integer) :: b

   if (n < 0) then
     b = 0
   else if (n >= nr_of_digits) then
      b = huge (b)
   else
      b % digit = 0
      b % digit (n) = 1
   end if

end function big_base_to_power

function f_base_to_power (n)  result (f)

   integer, intent (in) :: n
   type (big_fraction) :: f
   integer :: d

   d = n
   if (d < -f_nr_of_digits) then
      d = -f_nr_of_digits
   else if (d >= 2) then
      print *, "overflow in f_base_to_power"
      stop
   end if
   f % digit = 0
   f % digit (d) = 1

end function f_base_to_power

function big_plus_int (b, i) result (bi)

   type (big_integer), intent (in) :: b
   integer, intent (in) :: i
   type (big_integer) :: bi
   type (big_integer) :: temp_big

   temp_big = i
   bi = b + temp_big

end function big_plus_int

function big_plus_big (x, y) result (bb)

   type (big_integer), intent (in) :: x, y
   type (big_integer) :: bb
   integer :: carry, temp_digit, n

   carry = 0
   do n = 0, nr_of_digits
      temp_digit = &
         x % digit (n) + y % digit (n) + carry
      bb % digit (n) = modulo (temp_digit, base)
      carry = temp_digit / base
   end do

   if (bb % digit (nr_of_digits) /= 0 .or.  &
         carry /= 0) then
      print *, "overflow adding two big integers"
      stop
   end if

end function big_plus_big

function f_plus_f (x, y) result (ff)

   type (big_fraction), intent (in) :: x, y
   type (big_fraction) :: ff
   integer :: carry, temp_digit, n

   carry = 0
   do n = -f_nr_of_digits, nr_of_digits
      temp_digit = &
         x % digit (n) + y % digit (n) + carry
      ff % digit (n) = modulo (temp_digit, base)
      carry = temp_digit / base
   end do

   if (ff % digit (nr_of_digits) /= 0 .or.  &
         carry /= 0) then
      print *, "overflow adding two fractions"
      stop
   end if

end function f_plus_f

function big_minus_int (b, i) result (bi)

   type (big_integer), intent (in) :: b
   integer, intent (in) :: i
   type (big_integer) :: bi
   type (big_integer) :: temp_big

   temp_big = i
   bi = b - temp_big

end function big_minus_int

function big_minus_big (x, y) result (bb)

   type (big_integer), intent (in) :: x, y
   type (big_integer) :: bb
   type (big_integer) :: temp_big
   integer :: n

   temp_big = x
   do n = 0, nr_of_digits - 1
      bb % digit (n) = temp_big % digit (n) - y % digit (n)
      if (bb % digit (n) < 0) then
         bb % digit (n) = bb % digit (n) + base
         temp_big % digit (n + 1) = temp_big % digit (n + 1) - 1
      end if
   end do

   if (temp_big % digit (nr_of_digits) < 0) then
      bb % digit = 0
   else
      bb % digit (nr_of_digits) = 0
   end if

end function big_minus_big

function f_minus_f (x, y) result (ff)

   type (big_fraction), intent (in) :: x, y
   type (big_fraction) :: ff
   type (big_fraction) :: temp_f
   integer :: n

   temp_f = x
   do n = -f_nr_of_digits, nr_of_digits - 1
      ff % digit (n) = temp_f % digit (n) - y % digit (n)
      if (ff % digit (n) < 0) then
         ff % digit (n) = ff % digit (n) + base
         temp_f % digit (n + 1) = temp_f % digit (n + 1) - 1
      end if
   end do

   if (temp_f % digit (nr_of_digits) < 0) then
      ff % digit = 0
   else
      ff % digit (nr_of_digits) = 0
   end if

end function f_minus_f

function big_times_int (b, i) result (bi)

   type (big_integer), intent (in) :: b
   integer, intent (in) :: i
   type (big_integer) :: bi
   integer :: i0, i1, carry, n, temp_digit

   i0 = modulo (i, base)
   carry = 0
   do n = 0, nr_of_digits
      temp_digit = b % digit (n) * i0 + carry
      bi % digit (n) = modulo (temp_digit, base)
      carry = temp_digit / base
   end do

   if (bi % digit (nr_of_digits) /= 0 .or.  &
         carry /= 0) then
      print *, "overflow multiplying by an integer"
      stop
   end if

   if (i1 >= base) then
      i1 = i / base
      carry = 0
      do n = 1, nr_of_digits
         temp_digit = b % digit (n-1) * i1 + bi % digit (n) + carry
         bi % digit (n) = modulo (temp_digit, base)
         carry = temp_digit / base
      end do

      if (bi % digit (nr_of_digits) /= 0 .or.  &
            carry /= 0) then
         print *, "overflow multiplying by an integer"
         stop
      end if
   end if

end function big_times_int

function big_times_big (x, y) result (bb)

   type (big_integer), intent (in) :: x, y
   type (big_integer) :: bb
   type (big_integer) :: temp_big
   integer :: n

   bb % digit = 0
   do n = nr_of_digits, 0, -1
      if (x % digit (n) /= 0) then
         temp_big = y * (x % digit (n))
         temp_big % digit = eoshift (temp_big % digit, -n)
         bb = bb + temp_big
      end if
   end do

end function big_times_big

function f_times_int (f, i) result (fi)

   type (big_fraction), intent (in) :: f
   integer, intent (in) :: i
   type (big_fraction) :: fi
   integer :: i0, i1, carry, n, temp_digit

   i0 = modulo (i, base)
   carry = 0
   do n = -f_nr_of_digits, nr_of_digits
      temp_digit = f % digit (n) * i0 + carry
      fi % digit (n) = modulo (temp_digit, base)
      carry = temp_digit / base
   end do

   if (fi % digit (nr_of_digits) /= 0 .or.  &
         carry /= 0) then
      print *, "overflow multiplying by an integer"
      stop
   end if

   if (i1 >= base) then
      i1 = i / base
      carry = 0
      do n = -f_nr_of_digits + 1, nr_of_digits
         temp_digit = f % digit (n-1) * i1 + fi % digit (n) + carry
         fi % digit (n) = modulo (temp_digit, base)
         carry = temp_digit / base
      end do

      if (fi % digit (nr_of_digits) /= 0 .or.  &
            carry /= 0) then
         print *, "overflow multiplying by an integer"
         stop
      end if
   end if

end function f_times_int

function f_times_f (x, y) result (ff)

   type (big_fraction), intent (in) :: x, y
   type (big_fraction) :: ff
   type (big_fraction) :: temp_f
   integer :: n

   ff % digit = 0
   do n = nr_of_digits, -f_nr_of_digits, -1
      if (x % digit (n) > 0) then
         temp_f = y * (x % digit (n))
         temp_f % digit = eoshift (temp_f % digit, -n)
         ff = ff + temp_f
      end if
   end do

end function f_times_f

function big_div_int (b, i) result (bi)

   type (big_integer), intent (in) :: b
   integer, intent (in) :: i
   type (big_integer) :: bi
   integer :: n, temp_int, remainder

   remainder = 0
   do n = nr_of_digits, 0, -1
      temp_int = base * remainder + b % digit (n)
      bi % digit (n) = temp_int / i
      remainder = modulo (temp_int, i)
   end do

end function big_div_int

function big_div_big (x, y) result (bb)

   type (big_integer), intent (in) :: x, y
   type (big_integer) :: bb
   type (big_integer) :: temp_b
   type (big_fraction) :: temp_f
   integer :: m, n

   if (x < y) then
      bb = 0
   else if (y == 0) then
      bb = huge (y)
   else
       temp_f = x
       eps % digit = 0
       do m = nr_of_digits, 0, -1
          if (x % digit (m) > 0) then
             exit
          end if
       end do
       do n = nr_of_digits, 0, -1
          if (y % digit (n) > 0) then
             exit
          end if
       end do
       eps % digit (-n - max(m,n) - 1) = 1
       bb = temp_f * f_reciprocal (y)
       temp_b = bb + 1
       if (y * temp_b <= x) then
          bb = temp_b
       end if
   end if

end function big_div_big

function modulo_big_int (b, i) result (bi)

   type (big_integer), intent (in) :: b
   integer, intent (in) :: i
   integer :: bi
   integer :: n, remainder

   remainder = 0
   do n = nr_of_digits, 0, -1
      remainder = modulo (base * remainder + b % digit (n), i)
   end do

   bi = remainder

end function modulo_big_int

function modulo_big_big (x, y) result (bb)

   type (big_integer), intent (in) :: x, y
   type (big_integer) :: bb

   bb = x - x / y * y

end function modulo_big_big

function big_eq_int (b, i) result (bi)

   type (big_integer), intent (in) :: b
   integer, intent (in) :: i
   logical :: bi
   integer :: i0, i1

   i0 = modulo (i, base)
   i1 = i / base
   
   bi = b % digit (0) == i0 .and. &
        b % digit (1) == i1 .and. &
        all (b % digit (2 : nr_of_digits) == 0)

end function big_eq_int

function big_eq_big (x, y) result (bb)

   type (big_integer), intent (in) :: x, y
   logical :: bb
   integer :: n

   bb = .true.
   do n = nr_of_digits, 0, -1
      if (x % digit (n) /= y % digit (n)) then
         bb = .false.
         exit
      end if
   end do

end function big_eq_big

function big_le_int (b, i) result (bi)

   type (big_integer), intent (in) :: b
   integer, intent (in) :: i
   logical :: bi
   integer :: i0, i1

   i0 = modulo (i, base)
   i1 = i / base
   
   if (any (b % digit (2 : nr_of_digits) /= 0)) then
      bi = .false.
   else if (b % digit (1) > i1) then
      bi = .false.
   else if (b % digit (1) == i1 .and.  &
            b % digit (0) <= i0) then
      bi = .true.
   else
      bi = .false.
   end if
 
end function big_le_int

function big_le_big (x, y) result (bb)

   type (big_integer), intent (in) :: x, y
   logical :: bb
   integer :: n

   bb = .true.
   do n = nr_of_digits, 0, -1
      if (x % digit (n) /= y % digit (n)) then
         bb = (x % digit (n) < y % digit (n))
         exit
      end if
   end do

end function big_le_big

function f_le_f (x, y) result (ff)

   type (big_fraction), intent (in) :: x, y
   logical :: ff
   integer :: n

   ff = .true.
   do n = nr_of_digits, -f_nr_of_digits, -1
      if (x % digit (n) /= y % digit (n)) then
         ff = (x % digit (n) < y % digit (n))
         exit
      end if
   end do

end function f_le_f

function big_lt_int (b, i) result (bi)

   type (big_integer), intent (in) :: b
   integer, intent (in) :: i
   logical :: bi
   integer :: i0, i1

   i0 = modulo (i, base)
   i1 = i / base
   
   if (any (b % digit (2 : nr_of_digits) /= 0)) then
      bi = .false.
   else if (b % digit (1) > i1) then
      bi = .false.
   else if (b % digit (1) == i1 .and.  &
            b % digit (0) < i0) then
      bi = .true.
   else
      bi = .false.
   end if
 
end function big_lt_int

function big_lt_big (x, y) result (bb)

   type (big_integer), intent (in) :: x, y
   logical :: bb
   integer :: n

   bb = .false.
   do n = nr_of_digits, 0, -1
      if (x % digit (n) /= y % digit (n)) then
         bb = (x % digit (n) < y % digit (n))
         exit
      end if
   end do

end function big_lt_big

function huge_big (b) result (hb)

   type (big_integer), intent (in) :: b
   type (big_integer) :: hb

   hb % digit = base - 1
   hb % digit (nr_of_digits) = 0

end function huge_big

function sqrt_big (b) result (sb)

   type (big_integer), intent (in) :: b
   type (big_integer) :: sb
   type (big_integer) :: old_sqrt_big, new_sqrt_big
   integer :: i, n

   n = -1
   do i = nr_of_digits, 0, -1
      if (b % digit (i) /= 0) then
         n = i
         exit
      end if
   end do

   if (n == -1) then
      sb = 0
   else if (n == 0) then
      sb = int (sqrt (real (b % digit (0))))
   else
      old_sqrt_big = 0
      if (modulo (n, 2) == 0) then
         old_sqrt_big % digit (n / 2) = int (sqrt (real (b % digit (n))))
      else
         old_sqrt_big % digit ((n - 1) / 2) =  &
               int (sqrt (real (base * b % digit (n) + b % digit (n-1))))
      end if

      do
         new_sqrt_big = (old_sqrt_big + b / old_sqrt_big) / 2
         if (new_sqrt_big == old_sqrt_big .or.  &
             new_sqrt_big == old_sqrt_big + 1 .or.  &
             new_sqrt_big == 0) then
            exit
         else
            old_sqrt_big = new_sqrt_big
         end if
      end do
      sb = old_sqrt_big
   end if

end function sqrt_big

recursive function big_power_int (b, i)  &
      result (big_power_int_result)

   type (big_integer), intent (in) :: b
   integer, intent (in) :: i
   type (big_integer) :: big_power_int_result
   type (big_integer) :: temp_big

   if (i <= 0) then
      big_power_int_result = 1
   else
      temp_big = big_power_int (b, i / 2)
      if (modulo (i, 2) == 0) then
         big_power_int_result = temp_big * temp_big
      else
         big_power_int_result = temp_big * temp_big * b
      end if
   end if

end function big_power_int

function prime (b) result (pb)

   type (big_integer), intent (in) :: b
   logical :: pb
   type (big_integer) :: s, div

   if (b <= 1) then
      pb = .false.
   else if (b == 2) then
      pb = .true.
   else if (modulob (b, 2) == 0) then
         pb = .false.
   else
      pb = .true.
      s = sqrtb (b)
      div = 3
      do
         if (.not. s <= div) then
            exit
         end if
         if (modulob (b, div) == 0) then
            pb = .false.
            exit
         end if
         div = div + 2
      end do
   end if

end function prime

subroutine print_big (b)

   type (big_integer), intent (in) :: b
   type (big_integer) :: temp_big
   integer :: n, remainder
   character (len=1), dimension (nr_of_decimal_digits) :: &
         decimal_digits 
   character (len = 10), parameter :: digit_chars = "0123456789"

   temp_big = b
   decimal_digits = " "
   do n = 1, nr_of_decimal_digits
      if (all (temp_big % digit == 0)) then
         exit
      end if
      remainder = modulob (temp_big, 10) + 1
      temp_big = temp_big / 10
      decimal_digits (n) = digit_chars (remainder:remainder)
   end do

   write (unit = *, fmt = "(a)", advance = "no") " "
   if (n == 1) then
      write (unit = *, fmt = "(a)", advance = "no") "0"
   else
      do n = n - 1, 1, -1
         write (unit = *, fmt = "(a)", advance = "no") decimal_digits (n)
      end do
   end if

end subroutine print_big

subroutine print_f (f)

   type (big_fraction), intent (in) :: f
   type (big_fraction) :: temp_f
   type (big_integer) :: temp_b
   integer :: n, digit0
   character (len = 10), parameter :: digit_chars = "0123456789"

   temp_b = f
   call print_big (temp_b)
   write (unit = *, fmt = "(a)", advance = "no")  "."
   temp_f = f
   temp_f % digit (1 :) = 0
   do n = 1, nr_of_decimal_digits
      temp_f % digit (0) = 0
      temp_f = temp_f * 10
      digit0 = temp_f % digit (0) + 1
      write (unit = *, fmt = "(a)", advance = "no") &
            digit_chars (digit0 : digit0)
   end do

end subroutine print_f

subroutine print_big_base (b)

   type (big_integer), intent (in) :: b
   integer :: n, k

   print *, "base: ", base
   do n = nr_of_digits, 1, -1
      if (b % digit (n) /= 0) then
         exit
      end if
   end do
   print "(99i9)", (b % digit (k), k = n, 0, -1)

end subroutine print_big_base

function f_reciprocal (b) result (fr)

   type (big_integer), intent (in) :: b
   type (big_fraction) :: fr
   type (big_fraction) :: two_f, b_f, old_x, new_x
   type (big_integer) :: two_big
   integer, parameter :: base_2 = base * base
   integer :: n, b_approx

   two_big = 2
   two_f = two_big
   b_f = b
   do n = nr_of_digits - 1, 0, -1
      if (b % digit (n) /= 0) then
         exit
      end if
   end do
   if (n == -1) then
      print *, "attempt to divide by zero."
      stop
   else
      old_x % digit = 0
      b_approx = base * b_f % digit(n) + b_f % digit (n-1)
      old_x % digit (-n-1) = base_2 / b_approx
      old_x % digit (-n-2) =  base * modulo (base_2, b_approx) / b_approx
      do
         new_x = old_x * (two_f - old_x * b_f)
         if (new_x - old_x <= eps) then
            exit
         end if
         old_x = new_x
      end do
   end if

   fr = new_x

end function f_reciprocal

subroutine print_factors_big (b)

   type (big_integer), intent (in) :: b
   type (big_integer) :: temp_b, f, sqrt_b

   temp_b = b
   f = 2
   sqrt_b = sqrtb (b)
   do
     if (modulob (temp_b, f) == 0) then
        write (unit = *, fmt = "(a)", advance="no") "2 "
        temp_b = temp_b / 2
     else
        exit
     end if
   end do
   f = 3
   do
      if (temp_b == 1) then
         exit
      end if
      if (sqrt_b < f) then
         call print_big (temp_b)
         write (unit = *, fmt = "(a)", advance="no") " "
         exit
      end if
      do
         if (modulob (temp_b, f) == 0) then
            call print_big (f)
            write (unit = *, fmt = "(a)", advance="no") " "
            temp_b = temp_b / f
            sqrt_b = sqrtb (temp_b)
         else
            exit
         end if
      end do
      f = f + 2
   end do
   print *

end subroutine print_factors_big

end module big_integers_module

!program test_big_5
!
!   use big_integers_module
!   type (big_integer) :: a, b, c
!
!   a = "1"
!   b = "9999999999999999999"
!   c = "9999999999999999999"
!   call print_big (a + b * c)
!   print*
!
!end program test_big_5