big_integers.f90

用来求取大的整数

module big_integers_module

!  Copyright (c) 1993-2000 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 following operations are implemented.
!  b represents a big integer, c a character string,
!  and i an ordinary integer.

!  integer, parameter :: nr_of_decimal_digits

!  big (i)
!  big (c)
!  int (b)
!  int (c)
!  char (b)
!  char (i)

!  b = i
!  b = c
!  i = b
!  i = c
!  c = b
!  c = i

!  b ? i, i ? b, and b ? b, where ? is
!    +, -, *, /, 
!    <, <=, >, >=, ==, or /=

!  b ** i

!  modulo (b, i)  [result is integer]
!  modulo (i, b)  [result is big_integer]
!  modulo (b, b)  [result is big_integer]

!  huge (b)
!  sqrt (b)
!  prime (b) [result is logical]

!  call print_big (b)
!  call random_number (b, low, high)
!  call print_factors (b)

!  Many operations of the form b ? i, where i < base,
!  are implemented to be efficient as special cases.

   implicit none

   public :: big
   public :: int
   public :: char
   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 :: random_number
   public :: operator (==)
   public :: operator (/=)
   public :: operator (<=)
   public :: operator (<)
   public :: operator (>=)
   public :: operator (>)

   private :: big_gets_int, &
              big_int, &
              big_gets_char, &
              big_char, &
              int_gets_big, &
              int_big, &
              int_gets_char, &
              int_char, &
              char_gets_big, &
              char_big, &
              char_gets_int, &
              char_int

   private :: big_plus_int, &
              int_plus_big, &
              big_plus_big, &
              big_minus_int, &
              int_minus_big, &
              big_minus_big, &
              big_times_int, &
              int_times_big, &
              big_times_big, &
              big_div_int, &
              int_div_big, &
              big_div_big, &
              big_power_int, &
              modulo_big_int, &
              modulo_int_big, &
              modulo_big_big

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

   private :: huge_big, &
              big_base_to_power, &
              print_big_base, &
              sqrt_big, &
              msd, &
              random_number_big

   intrinsic char
   intrinsic int
   intrinsic modulo
   intrinsic huge
   intrinsic sqrt
   intrinsic random_number

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

   integer, parameter, public :: 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

   interface big
      module procedure big_char, &
                       big_int
   end interface

   interface int
      module procedure int_char, &
                       int_big
   end interface

   interface char
      module procedure char_big, &
                       char_int
   end interface

   interface assignment (=)
      module procedure big_gets_int, &
                       big_gets_char, &
                       int_gets_big, &
                       int_gets_char, &
                       char_gets_big, &
                       char_gets_int
   end interface

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

   interface operator (-)
      module procedure big_minus_int, &
                       int_minus_big, &
                       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, &
                       int_div_big, &
                       big_div_big
   end interface

   interface operator (**)
      module procedure big_power_int
   end interface

   interface modulo
      module procedure modulo_big_int, &
                       modulo_int_big, &
                       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

   interface random_number
      module procedure random_number_big
   end interface

contains

function big_char (c) result (b)

   character (len=*), intent (in) :: c
   type (big_integer) :: b
   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
         b = huge (b)
      end if
      b = b * 10 + temp_digit
   end do

end function big_char

subroutine big_gets_char (b, c)

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

   b = big_char (trim (c))

end subroutine big_gets_char

function big_int (i) result (b)

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

   if (i < 0) then
      b = huge (b)
   end if

   b % digit = 0
   temp_i = i
   do n = 0, nr_of_digits - 1
      if (temp_i == 0) then
         return
      end if
      b % digit (n) = modulo (temp_i, base)
      temp_i = temp_i / base
   end do

   if (temp_i /= 0) then
      b = huge (b)
   end if

end function big_int

subroutine big_gets_int (b, i)

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

   b = big (i)

end subroutine big_gets_int

function int_char (c) result (i)

   character (len=*), intent (in) :: c
   integer :: i

   i = int (big (trim (c)))

end function int_char

subroutine int_gets_char (i, c)

   integer, intent (out) :: i
   character (len=*), intent (in) :: c

   i = int_char (trim (c))

end subroutine int_gets_char

function int_big (b) result (i)

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

   if (msd (b) > 1) then
      i = huge (i)
   else
      i = base * b % digit (1) + b % digit (0)
   end if

end function int_big

subroutine int_gets_big (i, b)

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

   i = int (b)

end subroutine int_gets_big

function char_big (b) result (c)

   type (big_integer), intent (in) :: b
   character (len=nr_of_decimal_digits+9) :: c 
   type (big_integer) :: temp_big
   integer :: n, remainder
   character (len = 10), parameter :: digit_chars = "0123456789"

   temp_big = b
   c = repeat (" ", len(c)-1) // "0"
   do n = len (c), 1, -1
      if (temp_big == 0) then
         exit
      end if
      remainder = modulo (temp_big, 10) + 1
      temp_big = temp_big / 10
      c (n:n) = digit_chars (remainder:remainder)
   end do

   c = adjustl (c)

end function char_big

subroutine char_gets_big (c, b)

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

   c = char (b)

end subroutine char_gets_big

function char_int (i) result (c)

   integer, intent (in) :: i
   character (len=nr_of_decimal_digits+9) :: c

   c = big (i)

end function char_int

subroutine char_gets_int (c, i)

   integer, intent (in) :: i
   character (len=*), intent (out) :: c

   c = big (i)

end subroutine char_gets_int

function msd (x) result (msd_result)

! Find most significan digit of x

   type (big_integer), intent (in) :: x
   integer :: msd_result
   integer :: n

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

   msd_result = -1

end function msd

function big_plus_int (b, i) result (bi)

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

   if (i < base) then
      carry = i
      do n = 0, nr_of_digits - 1
         summ = b % digit (n) + carry
         bi % digit (n) = modulo (summ, base)
         carry = summ / base
         if (carry == 0) then
            bi % digit (n+1:) = b % digit (n+1:)
            return
         end if
      end do
      if (n==nr_of_digits) then
         bi = huge (bi)
      end if
   else
      bi = b + big (i)
   end if

end function big_plus_int

function int_plus_big (i, b) result (bi)

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

   bi = b + i

end function int_plus_big

function big_plus_big (x, y) result (bb)

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

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

   bb % digit (m+1) = carry
   bb % digit (m+2:nr_of_digits) = 0
   if (bb % digit (nr_of_digits) /= 0) then
      bb = huge (bb)
   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
   integer :: n, borrow, diff, msdb

   bi % digit = 0
   msdb = msd (b)
   if (msdb<1 .and. b % digit (0) < i) then
      return
   end if
 
   if (i < base) then
      borrow = i
      do n = 0, nr_of_digits - 1
         diff = b % digit (n) - borrow
         bi % digit (n) = modulo (diff, base)
         borrow = (base - diff) / base
         if (borrow == 0) then
            bi % digit (n+1:msdb) = b % digit (n+1:msdb)
            return
         end if
      end do
   else
      bi = b - big (i)
   end if

end function big_minus_int

function int_minus_big (i, b) result (ib)

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

   ib = big (i) - b

end function int_minus_big

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 :: ib, prod, carry

   if (i < base) then
      bi % digit = 0
      carry = 0
      do ib = 0, msd (b)
         prod = b % digit (ib) * i + carry
         bi % digit (ib) = modulo (prod, base)
         carry = prod / base
      end do
      if (ib==nr_of_digits .and. carry /= 0) then
         bi = huge (bi)
      else
         bi % digit (ib) = carry
      end if
   else
      bi = b * big (i)
   end if

end function big_times_int

function int_times_big (i, b) result (bi)

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

   bi = b * i

end function int_times_big

function big_times_big (x, y) result (bb)

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

   integer :: ix, iy, ib, carry, prod

   bb % digit = 0

   do ix = 0, msd (x)
      carry = 0
      ib = ix
      do iy = 0, msd (y)
         prod = x % digit (ix) * y % digit (iy) + bb % digit (ib) + carry
         carry = prod / base
         bb % digit (ib) = modulo (prod, base)
         if (ib == nr_of_digits) then
            bb = huge (bb)
            return
         end if
         ib = ib + 1
      end do
      bb % digit (ib) = bb % digit (ib) + carry