! Copyright (c) 1993-1996 Unicomp, Inc.
!
! Developed at Unicomp, Inc.
!
! Permission to use, copy, modify, and distribute this
! software is freely granted, provided that this notice 
! is preserved.

! The module named BIG_INTEGERS defines a new data type
! BIG_INTEGER.  This data type represents nonnegative integers
! up to 10 ** n-1, where n is the parameter (named constant)
! NR_OF_DECIMAL_DIGITS.  This value may be changed, but the
! module must then be recompiled (it is not dynamic).  The
! module also defines the following operations, where b
! represents a big_integer, i represents an ordinary
! Fortran integer, and c represents a Fortran character
! string.
! 
! big_int (i), big_int (b), big_int (c)
! b = i, b = c
! b + i, i + b, b + b
! b - i, b - b
! b * i, i * b, b * b
! b / i
! b ** i
! b ? i, i ? b, b ? b, where ? is ==. /=, <=, <, >=, >
! huge (b)
! modulo (b, i)
! digits (b)
! radix (b)
! range (b)
! sqrt (b)
! prime (b)
! factorial (b)
! call print_big (b)
! call print_factors (b)

module big_integer_module

   public :: big_int, assignment (=), &
             operator (+), operator (-), &
             operator (*), operator (/), operator (**), &
             operator (<), operator (>), &
             operator (==), operator (/=), &
             operator (<=), operator (>=), &
             huge, digits, range, radix, &
             modulo, sqrt, &
             print_big, print_factors, &
             factorial, prime

   private ::  big_int_int, &
               big_int_b, &
               big_int_char, &
               big_gets_int, &
               big_gets_char, &
               big_plus_big, &
               big_plus_int, &
               int_plus_big, &
               big_minus_big, &
               big_minus_int, &
               big_times_big, &
               big_times_int, &
               int_times_big, &
               big_div_big, &
               big_div_int, &
               big_power_int
   private ::  big_eq_big, &
               big_eq_int, &
               int_eq_big, &
               big_ne_big, &
               big_ne_int, &
               int_ne_big, &
               big_le_big, &
               big_le_int, &
               int_le_big, &
               big_lt_big, &
               big_lt_int, &
               int_lt_big, &
               big_ge_big, &
               big_ge_int, &
               int_ge_big, &
               big_gt_big, &
               big_gt_int, &
               int_gt_big, &
               huge_big, &
               modulo_big_int, &
               modulo_big_big, &
               digits_big, &
               radix_big, &
               range_big, &
               sqrt_big, &
               print_big_integer, &
               prime_big, &
               factorial_big, &
               big_base_to_power, &
               print_factors_big

   intrinsic huge, range, digits, radix
   intrinsic modulo, sqrt

   integer, parameter, private :: &
         nr_of_decimal_digits = 50, &
         d = (range (0) - 1) / 2, &
         base = 10 ** d, &
         nr_of_digits = nr_of_decimal_digits / d + 1

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

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

   interface big_int
      module procedure big_int_int, &
                       big_int_char, &
                       big_int_b
   end interface

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

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

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

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

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

   interface operator (**)
      module procedure big_power_int
   end interface

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

   interface operator (/=)
      module procedure big_ne_int, &
                       int_ne_big, &
                       big_ne_big
   end interface

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

   interface operator (>=)
      module procedure big_ge_int, &
                       int_ge_big, &
                       big_ge_big
   end interface

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

   interface operator (>)
      module procedure big_gt_int, &
                       int_gt_big, &
                       big_gt_big
   end interface

   interface huge
      module procedure huge_big
   end interface

   interface modulo
      module procedure modulo_big_int, &
                       modulo_big_big
   end interface

   interface digits
      module procedure digits_big
   end interface

   interface radix
      module procedure radix_big
   end interface

   interface range
      module procedure range_big
   end interface

   interface sqrt
      module procedure sqrt_big
   end interface

   interface print_big
      module procedure print_big_integer
   end interface

   interface prime
      module procedure prime_big
   end interface

   interface factorial
      module procedure factorial_big
   end interface

   interface print_factors
      module procedure print_factors_big
   end interface

contains

function big_int_int (i) result (b_result)

   integer, intent (in) :: i
   type (big_integer) :: b_result

   if (i < 0) then
      b_result = huge (b_result)
   else
      b_result % digit (0) = modulo (i, base)
      b_result % digit (1) = i / base
      b_result % digit (2:) = 0
   end if

end function big_int_int

function big_int_char (c)  result (b_result)

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

   b_result % digit = 0
   do n = 1, len (c)
      temp_digit = index ("0123456789", c (n:n)) - 1
      if (temp_digit < 0) then
         b_result = huge (b_result)
      else
         b_result = b_result * 10 + temp_digit
      end if
   end do

end function big_int_char

function big_int_b (b)  result (b_result)

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

   b_result = b

end function big_int_b

subroutine big_gets_int (b, i)

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

   b = big_int_int (i)

end subroutine big_gets_int

subroutine big_gets_char (b, c)

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

   b = big_int_char (c)

end subroutine big_gets_char

function big_plus_int (b, i) result (b_result)

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

   b_result = b + big_int (i)

end function big_plus_int

function int_plus_big (i, b) result (b_result)

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

   b_result = b + big_int (i)

end function int_plus_big

function big_plus_big (x, y)  result (b_result)

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

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

   if (b_result % digit (nr_of_digits) /= 0 .or.  &
         carry /= 0) then
      b_result = huge (b_result)
   end if

end function big_plus_big

function big_minus_int (b, i)  result (b_result)

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

   b_result = b - big_int (i)

end function big_minus_int

function big_minus_big (x, y)  result (b_result)

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

   temp_big = x
   do n = 0, nr_of_digits - 1
      b_result % digit (n) = temp_big % digit (n) - y % digit (n)
      if (b_result % digit (n) < 0) then
         b_result % digit (n) = b_result % 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
      b_result % digit = 0
   else
      b_result % digit (nr_of_digits) = 0
   end if

end function big_minus_big

function big_times_int (b, i)  result (b_result)

   type (big_integer), intent (in) :: b
   integer, intent (in) :: i
   type (big_integer) :: b_result
   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
      b_result % digit (n) = modulo (temp_digit, base)
      carry = temp_digit / base
   end do

   if (b_result % digit (nr_of_digits) /= 0 .or.  &
         carry /= 0) then
      b_result = huge (b_result)
      return
   end if

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

      if (b_result % digit (nr_of_digits) /= 0 .or.  &
            carry /= 0) then
         b_result = huge (b_result)
      end if
   end if

end function big_times_int

function int_times_big (i, b)  result (b_result)

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

   b_result = b * i

end function int_times_big

function big_times_big (x, y)  result (b_result)

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

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

end function big_times_big

function big_div_int (b, i)  result (b_result)

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

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

end function big_div_int

function big_div_big (x,y)  result (b_result)

   type (big_integer), intent (in) :: x, y
   type (big_integer) :: b_result
   type (big_integer) :: q
   integer :: high, low, mid, b, d, m, n, k

   do n = nr_of_digits -1, 0, -1
      if (y % digit (n) /= 0) then
         exit
      end if
   end do
   if (n == -1) then  ! denominator = 0
      b_result = huge (b_result)
      return
   else if (n <= 1) then
      b_result = x / (base * y % digit (1) + y % digit (0))
      return
   end if

   do m = nr_of_digits -1, 0, -1
      if (x % digit (m) /= 0) then
         exit
      end if
   end do
   if (m < n) then
      b_result = 0
      return
   end if

   d = y % digit (n)
   q = x
   b_result % digit (m - n + 1 :) = 0

   do k = m - n, 0, -1
      b = base * q % digit (m + 1) + q % digit (m)
      high = (b + 1) / d + 1
      low = b / (d + 1)
      do
         if (low >= high - 1) then
            exit 
         end if
         mid = (high + low) / 2
         if (mid * y <= x) then
            low = mid
         else
            high = mid
         end if
      end do
      b_result % digit (k) = low
      q = q - low * y * big_base_to_power (k)
      m = m - 1
   end do

end function big_div_big

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

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

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

end function big_power_int

function big_eq_int (b, i)  result (l_result)

   type (big_integer), intent (in) :: b
   integer, intent (in) :: i
   logical :: l_result

   l_result = b == big_int (i)

end function big_eq_int

function int_eq_big (i, b)  result (l_result)

   integer, intent (in) :: i
   type (big_integer), intent (in) :: b
   logical :: l_result

   l_result = b == big_int (i)

end function int_eq_big

function big_eq_big (x, y)  result (l_result)

   type (big_integer), intent (in) :: x, y
   logical :: l_result

   l_result = all (x % digit (0 : nr_of_digits - 1) ==   &
                   y % digit (0 : nr_of_digits - 1))

end function big_eq_big

function big_ne_int (b, i)  result (l_result)

   type (big_integer), intent (in) :: b
   integer, intent (in) :: i
   logical :: l_result

   l_result = b /= big_int (i)

end function big_ne_int

function int_ne_big (i, b)  result (l_result)

   integer, intent (in) :: i
   type (big_integer), intent (in) :: b
   logical :: l_result

   l_result = b /= big_int (i)

end function int_ne_big

function big_ne_big (x, y)  result (l_result)

   type (big_integer), intent (in) :: x, y
   logical :: l_result

   l_result = any (x % digit (0 : nr_of_digits - 1) /=   &
                   y % digit (0 : nr_of_digits - 1))

end function big_ne_big

function big_le_int (b, i)  result (l_result)

   type (big_integer), intent (in) :: b
   integer, intent (in) :: i
   logical :: l_result

   l_result = b <= big_int (i)
 
end function big_le_int

function int_le_big (i, b)  result (l_result)

   type (big_integer), intent (in) :: b
   integer, intent (in) :: i
   logical :: l_result

   l_result = b <= big_int (i)
 
end function int_le_big

function big_le_big (x, y)  result (l_result)

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

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

end function big_le_big

function big_lt_int (b, i)  result (l_result)

   type (big_integer), intent (in) :: b
   integer, intent (in) :: i
   logical :: l_result

   l_result = b < big_int (i)
 
end function big_lt_int

function int_lt_big (i, b)  result (l_result)

   type (big_integer), intent (in) :: b
   integer, intent (in) :: i
   logical :: l_result

   l_result = big_int (i) < b
 
end function int_lt_big

function big_lt_big (x, y)  result (l_result)

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

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

end function big_lt_big

function big_ge_int (b, i)  result (l_result)

   type (big_integer), intent (in) :: b
   integer, intent (in) :: i
   logical :: l_result

   l_result = b >= big_int (i)
 
end function big_ge_int

function int_ge_big (i, b)  result (l_result)

   type (big_integer), intent (in) :: b
   integer, intent (in) :: i
   logical :: l_result

   l_result = b >= big_int (i)
 
end function int_ge_big

function big_ge_big (x, y)  result (l_result)

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

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

end function big_ge_big

function big_gt_int (b, i)  result (l_result)

   type (big_integer), intent (in) :: b
   integer, intent (in) :: i
   logical :: l_result

   l_result = b > big_int (i)
 
end function big_gt_int

function int_gt_big (i, b)  result (l_result)

   type (big_integer), intent (in) :: b
   integer, intent (in) :: i
   logical :: l_result

   l_result = big_int (i) > b
 
end function int_gt_big

function big_gt_big (x, y)  result (l_result)

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

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

end function big_gt_big

function huge_big (b)  result (b_result)

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

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

end function huge_big

function modulo_big_int (b, i)  result (b_result)

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

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

   b_result = remainder

end function modulo_big_int

function modulo_big_big (a, p)  result (b_result)

   type (big_integer), intent (in) :: a, p
   type (big_integer) :: b_result

!  b_result = a - a / p * p
   type (big_integer) :: q
   integer :: high, low, mid, b, d, m, n, k

   do n = nr_of_digits -1, 0, -1
      if (p % digit (n) /= 0) then
         exit
      end if
   end do
   if (n == -1) then  ! denominator = 0
      b_result = huge (b_result)
      return
   else if (n <= 1) then
      b_result = modulo (a, (base * p % digit (1) + p % digit (0)))
      return
   end if

   do m = nr_of_digits -1, 0, -1
      if (a % digit (m) /= 0) then
         exit
      end if
   end do
   if (m < n) then
      b_result = p
      return
   end if

   d = p % digit (n)
   q = a
   b_result % digit (m - n + 1 :) = 0

   do k = m - n, 0, -1
      b = base * q % digit (m + 1) + q % digit (m)
      high = (b + 1) / d + 1
      low = b / (d + 1)
      do
         if (low >= high - 1) then
            exit 
         end if
         mid = (high + low) / 2
         if (mid * p <= a) then
            low = mid
         else
            high = mid
         end if
      end do
      q = q - low * p * big_base_to_power (k)
      m = m - 1
   end do

   b_result = q

end function modulo_big_big

function digits_big (b)  result (digits_result)

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

   digits_result = nr_of_digits

end function digits_big

function radix_big (b)  result (radix_result)

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

   radix_result = base

end function radix_big

function range_big (b)  result (range_result)

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

   range_result = d * nr_of_digits - 1

end function range_big

function sqrt_big (x)  result (b_result)

   type (big_integer), intent (in) :: x
   type (big_integer) :: b_result
   type (big_integer) :: old_sqrt_big, new_sqrt_big
   integer :: n, i

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

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

      do
         new_sqrt_big = (old_sqrt_big + x / 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
      b_result = old_sqrt_big
   end if

end function sqrt_big

function big_base_to_power (n)  result (b_result)

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

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

end function big_base_to_power

function factorial_big (x)  result (b_result)

   integer, intent (in) :: x
   type (big_integer) :: b_result
   type (big_integer) :: h
   integer :: m

   h = huge (h)
   b_result = 1
   m = 1

   do
      if (m > x) then
         exit
      end if
      if (b_result == h) then
         return
      end if
      b_result = b_result * m
      m = m + 1
   end do

end function factorial_big

function prime_big (x)  result (l_result)

   type (big_integer), intent (in) :: x
   logical :: l_result
   type (big_integer) :: s, div

   if (x <= 1) then
      l_result = .false.
   else if (x == 2) then
      l_result = .true.
   else if (modulo (x, 2) == 0) then
         l_result = .false.
   else
      l_result = .true.
      s = sqrt (x)
      div = 3
      do
         if (div > s) then
            exit
         else if (modulo (x, div) == 0) then
            l_result = .false.
            exit
         end if
         div = div + 2
      end do
   end if

end function prime_big

subroutine print_big_integer (x)

   type (big_integer), intent (in) :: x
   integer :: n, p
   character (len = 100) :: char_d, char_x, big_format

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

   if (n == -1) then
      write (unit=*, fmt="(a)", advance = "no") "0"
   else
      write (unit=char_x, fmt=*) x % digit (n)
      p = verify (char_x, " ")
      char_d = char_x (p :)
      write (unit=*, fmt="(a)", advance = "no") trim (char_d)
      write (unit=char_d, fmt=*) d
      p = verify (char_d, " ")
      char_d = char_d (p :)
      write (unit=big_format, fmt=*) nr_of_digits
      p = verify (big_format, " ")
      big_format = big_format (p :)
      big_format = "(" // trim (big_format) // "i" //  &
                    trim (char_d) // "." // trim (char_d) // ")"
      write (unit=*, fmt=big_format, advance = "no") x % digit (n-1 : 0 : -1)
   end if

end subroutine print_big_integer

subroutine print_factors_big (x)

   type (big_integer), intent (in) :: x
   type (big_integer) :: temp_x, f, sqrt_x

   temp_x = x
   f = 2
   sqrt_x = sqrt (x)
   do
     if (modulo (temp_x, f) == 0) then
        write (unit=*, fmt="(a)", advance="no") "2 "
        temp_x = temp_x / 2
     else
        exit
     end if
   end do
   f = 3
   do
      if (temp_x == 1) then
         exit
      else if (sqrt_x < f) then
         call print_big (temp_x)
         write (unit=*, fmt="(a)", advance="no") " "
         exit
      end if
      do
         if (modulo (temp_x, f) == 0) then
            call print_big (f)
            write (unit=*, fmt="(a)", advance="no") " "
            temp_x = temp_x / f
            sqrt_x = sqrt (temp_x)
         else
            exit
         end if
      end do
      f = f + 2
   end do
   print *

end subroutine print_factors_big

end module big_integer_module

program test_big_module

use big_integer_module

type (big_integer) :: b1

b1 = "1234567654321"
call print_big (b1)
print *
call print_factors (b1)

end program test_big_module
