## 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
```

Michael L. Hall