Back to F Example Codes Page
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?