Example F Program--Big Integers

module big_integers_module

   public :: print_big
   public :: print_factors
   public :: prime
   public :: assignment (=)
   public :: operator (+)
   public :: operator (-)
   public :: operator (*)
   public :: operator (/)
   public :: operator (**)
   public :: modulo
   public :: huge
   public :: sqrt
   public :: operator (==)
   public :: operator (/=)
   public :: operator (<=)
   public :: operator (<)
   public :: operator (>=)
   public :: operator (>)
   private ::  big_gets_int, &
               big_gets_char, &
               big_plus_int, &
               big_plus_big, &
               big_minus_int, &
               big_minus_big, &
               big_times_int, &
               big_times_big, &
               big_div_int, &
               big_div_big, &
               big_power_int, &
               modulo_big_int, &
               modulo_big_big, &
               big_eq_int, &
               int_eq_big, &
               big_eq_big, &
               big_ne_int, &
               int_ne_big, &
               big_ne_big, &
               big_le_int, &
               int_le_big, &
               big_le_big, &
               big_ge_int, &
               int_ge_big, &
               big_ge_big, &
               big_lt_int, &
               int_lt_big, &
               big_lt_big, &
               big_gt_int, &
               int_gt_big, &
               big_gt_big, &
               huge_big, &
               big_base_to_power, &
               print_big_base, &
               sqrt_big

   intrinsic modulo
   intrinsic huge
   intrinsic sqrt

   ! This indicates the maximum number of decimal digits
   ! that a big integer may contain.

   integer, parameter, private :: nr_of_decimal_digits = 50

   ! If the radix (returned by "radix(0)" of the integers on
   ! your system is not 2 change the following constant to
   ! the logarithm in the base 10 of the radix: log10(radix)

   real, parameter, private :: log_base_10_of_radix = 0.30103

   integer, parameter, private :: &
         d = digits (0) / 2, &
         r = radix (0), &
         base = r ** d, &
         nr_of_digits = nr_of_decimal_digits / (log_base_10_of_radix * d) + 1

   ! The base of the 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

   integer, parameter, private :: f_nr_of_digits = 3 * nr_of_digits

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

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

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

   interface operator (*)
      module procedure big_times_int, &
                       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 modulo
      module procedure modulo_big_int, &
                       modulo_big_big
   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 sqrt
      module procedure sqrt_big
   end interface

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

   if (len (c) > nr_of_decimal_digits) then
      b = huge (b)
      return
   end if
   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

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 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 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 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 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

   if (i == 0) then
      bi = huge (bi)
      return
   end if
   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) :: high, low
   integer :: m, n

   if (x < y) then
      bb = 0
   else if (y == 0) then
      bb = huge (y)
   else if (y == 1) then
      bb = x
   else
      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
      high % digit = 0
      if (m - n + 1 >= nr_of_digits) then
         high = huge (high)
      else
         high % digit (m - n + 1) = 1
      end if
      low % digit = 0
      low % digit (max(0, m - n - 1)) = 1
      do
         bb = (high + low) / 2
         if (x < bb * y) then
            high = bb
         else if ((bb + 1) * y <= x) then
            low = bb
         else
            exit
         end if
      end do
   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

   if (i==0) then
      bi = huge (bi)
      return
   end if
   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

   if (y == 0) then
      bb = huge (bb)
   else
      bb = x - x / y * y
   end if

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 int_eq_big (i, b) 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 int_eq_big

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_ne_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)

   bi = .not. bi

end function big_ne_int

function int_ne_big (i, b) 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)

   bi = .not. bi

end function int_ne_big

function big_ne_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

   bb = .not. bb

end function big_ne_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 int_le_big (i, b) 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 int_le_big

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 big_gt_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

   bi = .not. bi
 
end function big_gt_int

function int_gt_big (i, b) 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

   bi = .not. bi
 
end function int_gt_big

function big_gt_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

   bb = .not. bb

end function big_gt_big

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 int_lt_big (i, b) 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 int_lt_big

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 big_ge_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

   bi = .not. bi
 
end function big_ge_int

function int_ge_big (i, b) 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

   bi = .not. bi
 
end function int_ge_big

function big_ge_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

   bb = .not. bb

end function big_ge_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 (modulo (b, 2) == 0) then
         pb = .false.
   else
      pb = .true.
      s = sqrt (b)
      div = 3
      do
         if (.not. s <= div) then
            exit
         end if
         if (modulo (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 = modulo (temp_big, 10) + 1
      temp_big = temp_big / 10
      decimal_digits (n) = digit_chars (remainder:remainder)
   end do

   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_big_base (b)

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

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

end subroutine print_big_base

subroutine print_factors (b)

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

   temp_b = b
   f = 2
   sqrt_b = sqrt (b)
   do
     if (modulo (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 (modulo (temp_b, f) == 0) then
            call print_big (f)
            write (unit = *, fmt = "(a)", advance="no") " "
            temp_b = temp_b / f
            sqrt_b = sqrt (temp_b)
         else
            exit
         end if
      end do
      f = f + 2
   end do

end subroutine print_factors

end module big_integers_module

program test_big

   use big_integers_module

   type (big_integer) :: b1

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

end program test_big

Back to F Example Codes Page

Back to F Homepage


Questions or Comments?