Example F Program--Big Integers
! 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
! F integer, and c represents an F 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