big_integers.f90

用来求取大的整数

   end do
   
end function big_times_big

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_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)
   else if (i < base) then
      bi % digit = 0
      remainder = 0
      do n = msd(b), 0, -1
         temp_int = base * remainder + b % digit (n)
         bi % digit (n) = temp_int / i
         remainder = modulo (temp_int, i)
      end do
   else
      bi = b / big (i)
   end if

end function big_div_int

function int_div_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_div_big

function big_div_big (x, y) result (bb)

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

   type (big_integer) :: tx, ty

   integer :: msdx, msdy, ix, iy
   integer :: v1, v2, u0, u1, u2
   integer :: dd, bi, car, bar, prd

   if (y == 0) then
      bb = huge (bb)
      return
   end if

   msdx = msd(x)
   msdy = msd(y)

   if (msdy == 0) then
      bb = x / y % digit (0)
      return
   end if

   bb % digit = 0

   if (msdy < msdy) then
      return
   end if

   tx = x
   ty = y

   car = 0
   bar = 0
   prd = 0
   dd = base / (ty % digit (msdy) + 1)
   if (dd /= 1) then
      do ix = 0, msdx
         tx % digit (ix) = tx % digit (ix) * dd + car
         car = tx % digit (ix) / base
         tx % digit (ix) = tx % digit (ix) - base * car
      end do
      tx % digit (msdx+1) = car
      car = 0
      do iy = 0, msdy
         ty % digit (iy) = ty % digit (iy) * dd + car
         car = ty % digit (iy) / base
         ty % digit (iy) = ty % digit (iy) - base * car
      end do
   end if

   msdx = msdx + 1

   v1 = ty % digit (msdy)
   v2 = ty % digit (msdy-1)
   bb % digit = 0

   do msdx = msdx, msdy + 1, -1

      u0 = tx % digit (msdx)
      u1 = tx % digit (msdx-1)
      u2 = tx % digit (msdx-2)

      if (u0 == v1) then
         bi = base - 1
      else 
         bi = (u0*base + u1) / v1
      end if

      do
         if (v2*bi <= (u0*base + u1 - bi*v1) * base + u2) then
            exit
         end if
         bi = bi - 1
      end do

      if (bi > 0) then
         car = 0
         bar = 0
         ix = msdx - msdy - 1
         do iy = 0, msdy
            prd = bi * ty % digit (iy) + car
            car = prd / base
            prd = prd - base * car
            tx % digit (ix) = tx % digit (ix) - (prd + bar)
            if (tx % digit (ix) < 0) then
               bar = 1
               tx % digit (ix) = tx % digit (ix) + base
            else
               bar = 0
            end if
            ix = ix + 1
         end do
         if (tx % digit (msdx) < car + bar) then
            car = 0
            bi = bi -1
            ix = msdx - msdy - 1
            do iy = 0, msdy
               tx % digit (ix) = tx % digit (ix) + ty % digit (iy) + car
               if (tx % digit (ix) > base) then
                  car = 1
                  tx % digit (ix) = tx % digit (ix) - base
               else
                  car = 0
               end if
               ix = ix + 1
            end do
         end if
      end if
      tx % digit (msdx) = 0
      bb % digit (1:nr_of_digits) = bb % digit (0:nr_of_digits-1)
      bb % digit (0) = bi
   end do

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

   if (i == 0) then
      bi = huge (bi)
   else if (i < base) then
      bi = 0
      do n = msd(b), 0, -1
         bi = modulo (base * bi + b % digit (n), i)
      end do
   else
      bi = modulo (b, big (i))
   end if

end function modulo_big_int

function modulo_int_big (ii, b) result (ib)

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

   ib = modulo (big (ii), b)

end function modulo_int_big

function modulo_big_big (x, y) result (bb)

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

   bb = x - x / y * y

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

   bi = int (b) == i

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

   bi = int (b) == i

end function int_eq_big

function big_eq_big (x, y) result (bb)

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

   bb = all (x % digit == y % digit)

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

   bi = int (b) /= i

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

   bi = int (b) /= i

end function int_ne_big

function big_ne_big (x, y) result (bb)

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

   bb = any (x % digit /= y % digit)

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

   bi = int (b) <= i
 
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

   bi = i <= int (b)
 
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

   bi = int (b) > i
 
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

   bi = i > int (b)
 
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

   bi = int (b) < i

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

   bi = i < int (b)
 
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

   bi = int (b) >= i
 
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

   bi = i >= int (b)
 
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 (0) = b % digit (0) ! to avoid diagnostic
   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

   write (unit = *, fmt = "(a)", advance = "no") trim (char (b))

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

!  Use ordinary integer fi as a test factor.

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

subroutine random_number_big (r, low, high)
 
!  Generate by linear congruence x' = ax + c mod m
!  where m is huge (b) + 1 = base ** nr_of_digits

   type (big_integer), intent (out) :: r
   type (big_integer), intent (in) :: low, high
   integer :: n, i, carry, prod, summ
   type (big_integer), save :: x = big_integer ( (/ (1, i=0,nr_of_digits-1), 0 /) )
   type (big_integer), parameter :: h = big_integer ( (/ (base-1, i=0,nr_of_digits-1), 0 /) )
   integer, parameter :: a = 16907, c = 8191

!  Multiply by a
   carry = 0
   do n = 0, nr_of_digits - 1
      prod = x % digit (n) * a + carry
      x % digit (n) = modulo (prod, base)
      carry = prod / base
   end do

!  Add c
   carry = c
   do n = 0, nr_of_digits - 1
      summ = x % digit (n) + carry
      x % digit (n) = modulo (summ, base)
      carry = summ / base
   end do

   r = x / (h / (high -low + 1)) + low

end subroutine random_number_big

end module big_integers_module

program test_big

   use big_integers_module
   implicit none

   type (big_integer) :: a

   a = "12345678987654321"
   call print_factors (a)

end program test_big