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