E.1.1 Prime_Factors_Math_Utils Procedure

The main documentation of the Prime_Factors_Math_Utils Procedure contains additional explanation of this code listing.

  subroutine Prime_Factors_Math_Utils (Number, NFactors, Factors, Verbose)

    ! Input variables.

    type(integer), intent(in) :: Number            ! The number to be factored.
    type(logical), intent(in), optional :: Verbose ! Verbosity.

    ! Output variables.

    ! The number of prime factors found.
    type(integer), intent(out) :: NFactors 
    ! The vector of prime factors.
    type(integer), dimension(32), intent(out) :: Factors

    ! Internal variables.

    type(logical) :: A_Verbose        ! Actual verbosity.
    type(integer) :: i                ! Loop counter.
    type(integer) :: Remaining_Factor ! The current number to be factored.

    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    ! Verify requirements.

    VERIFY(Valid_State(Number),5)             ! Number is valid.

    ! Set verbosity toggle.
    
    if (PRESENT(Verbose)) then
      A_Verbose = Verbose
    else
      A_Verbose = .false.
    end if

    ! Initializations.

    Factors = 0
    NFactors = 0

    ! Set Remaining_Factor based on sign of Number, pull out -1 
    ! factor if negative, treat special cases of 0 and 1.

    select case (Number)
    case (:-1)
      Remaining_Factor = -Number
      NFactors = 1
      Factors(1) = -1
    case (2:)
      Remaining_Factor = Number
    case (1)
      NFactors = 1
      Factors(1) = 1
      go to 1  ! Output if requested and return.
    case (0)
      NFactors = 1
      Factors(1) = 0
      go to 1  ! Output if requested and return.
    end select

    ! Pull out all factors of 2 and 3. Note that the problem size 
    ! is reduced with every factor that is found.

    call Extract_Factor (Remaining_Factor, Factors, NFactors, 2)
    call Extract_Factor (Remaining_Factor, Factors, NFactors, 3)

    ! Next, notice that the set of integers may be written:
    !
    !   <Integers> = {6i, 6i+1, 6i+2, 6i+3, 6i+4, 6i+5}
    !
    ! Since all factors of 2 and 3 have already been removed from 
    ! Remaining_Factor, we only need to check numbers that are
    ! not divisible by 2 and 3. This rules out 6i, 6i+2, 6i+3, and
    ! 6i+4. Therefore, we only need to check numbers in the set {6i+1,
    ! 6i+5}, or equivalently, the set {6i-1, 6i+1}. The next prime
    ! number to check after 3 is 5, so we start with i=1.
    !
    ! The problem size is again reduced every time a factor is found. 
    ! When the current i value becomes greater than the square root of 
    ! the current problem size, we are finished. Note that the upper 
    ! limit on this do-loop is unimportant -- the exit statements are 
    ! always executed first.

    do i = 6, ABS(Number), 6
      call Extract_Factor (Remaining_Factor, Factors, NFactors, i-1)
      if (i - 1 > INT(sqrt(REAL(Remaining_Factor)))) exit
      call Extract_Factor (Remaining_Factor, Factors, NFactors, i+1)
      if (i + 1 > INT(sqrt(REAL(Remaining_Factor)))) exit
    end do

    ! Include the Remaining_Factor if necessary.

    if (Remaining_Factor /= 1) then
      NFactors = NFactors + 1
      Factors(NFactors) = Remaining_Factor
    end if

    ! Output.

1   if (A_Verbose) then 
      write (6,*)
      write (6,100) 'Number:        ', Number
      write (6,100) 'Prime Factors: ', Factors(1:NFactors)
100   format (a,6i11,/,(15x,6i11))
    end if
  
    ! Verify guarantees.

    VERIFY(NFactors<=32,5)         ! Factors array is big enough.
    ! Product of factors is correct.
    VERIFY(Number==PRODUCT(Factors(1:NFactors)),5)  

    return
  end subroutine Prime_Factors_Math_Utils

  ! Auxiliary procedure for Prime_Factors_Math_Utils.

  subroutine Extract_Factor_Math_Utils (Remaining_Factor, Factors, NFactors, &
                                        Divisor)

    implicit none

    ! Variables described in Prime_Factors_Math_Utils.

    integer, intent(inout) :: Factors(32), NFactors, Remaining_Factor

    ! Internal variables.

    integer, intent(in) :: Divisor  ! Divisor to be checked.
    integer :: Remainder   ! The remainder when Remaining_Factor is 
                           ! divided by Divisor.

    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    ! Pull out all multiples of Divisor.

    do
      Remainder = Remaining_Factor - Divisor*(Remaining_Factor/Divisor)
      if (Remainder == 0) then 
        NFactors = NFactors + 1
        Factors(NFactors) = Divisor
        Remaining_Factor = Remaining_Factor/Divisor
        cycle
      end if
      exit
    end do
    return

  end subroutine Extract_Factor_Math_Utils