subroutine aptpfac (inum, nfactm, ifact, ipow, nfact,
& itot, inumr, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTPFAC
c
c call aptpfac (inum, nfactm, ifact, ipow, nfact, itot, inumr, nerr)
c
c Version: aptpfac Updated 2003 September 4 15:00.
c aptpfac Originated 2003 August 29 13:40.
c
c Authors: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: For the integer inum, find any prime factors ifact(n) > 1,
c and their exponents ipow(n), for n = 1, nfact, with nfact
c limited to nfactm. If not all prime factors of inum are
c returned, because nfactm is insufficient, then inumr, the
c remaining unfactored part of inum, is returned with a value
c other than 1, with nerr = 3.
c
c Also find itot, Euler's totient function for modulo inum
c arithmetic, the number of integers in the set (1, ..., inum - 1)
c that are relatively prime to inum (have no common factor other
c than 1).
c
c The method is slow, but avoids storing a large number of prime
c numbers. Each successive integer, 2, 3, 5, 7, 9, ..., is
c tested as a factor multiple times, and the result of dividing
c inum by each factor used to limit the search to integers no
c larger than the square root of that result.
c
c Flag nerr indicates any input error or limited results.
c
c Input: inum, nfactm.
c
c Output: ifact, ipow, nfact, itot, inumr, nerr.
c
c Calls: (none)
c
c Glossary:
c
c ifact Output Prime factor ifact(n) is the n'th smallest prime
c factor of inum, with an exponent of ipow(n).
c If inum = 1, then ifact(1) = ipow(1) = 1, itot = 0.
c Size at least nfactm. Initially set to -999999999.
c
c inum Input Any positive integer, for which all prime factors and
c their exponents are to be found, up to a maximum of
c nfactm prime factors.
c
c inumr Output Any remaining unfactored part of inum, when nfact
c is prevented from exceeding nfactm. Normally 1.
c
c ipow Output Exponent ipow(n) is the largest power of ifact(n)
c that divides inum. Initially set to -999999999.
c Size at least nfactm.
c
c itot Output Euler's totient function for modulo inum arithmetic.
c The number of integers in the set
c (1, ..., inum - 1) that are relatively prime to
c inum (have no common factors other then 1).
c Initially set to -999999999.
c
c nerr Output If not zero, indicates an error.
c 1 if inum is not 1 or more.
c 2 if nfactm is not positive.
c 3 if more than nfactm prime factors exist.
c
c nfact Output The number of prime factors found for integer inum.
c If nfactm are found, and more exist, nerr = 3.
c
c nfactm Input The maximum allowed size of arrays ifact and ipow.
c On a machine with 64-bit integers, 18 is enough.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c.... Declarations for arguments.
integer ifact(1) ! A prime factor of inum.
integer inum ! Any positive integer.
integer inumr ! Result of inum divided by all factors.
integer ipow(1) ! Largest power of ifact(n) dividing inum.
integer itot ! Euler's totient function for mod inum.
integer nerr ! Indicates an error, if not zero.
integer nfact ! How many prime factors of inum found.
integer nfactm ! The maximum allowed value of nfact.
c.... Declarations for local variables.
integer itest ! Integer to test as factor.
integer itestm ! Maximum integer to test as factor.
integer n ! Index in arrays ifact, ipow.
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptpfac finding prime factors of inum =',i20 /
cbug & ' but no more than',i6,' prime factors.' )
cbug write ( 3, 9901) inum, nfactm
cbugc***DEBUG ends.
c.... Initialize.
nerr = 0
nfact = 0
inumr = inum
if (nfactm .gt. 0) then
do n = 1, nfactm
ifact(n) = -999999999
ipow(n) = -999999999
enddo
endif
itot = -999999999
c.... Test for input errors.
if (inum .lt. 1) then
nerr = 1
go to 210
endif
if (nfactm .le. 0) then
nerr = 2
go to 210
endif
if (inum .eq. 1) then
nfact = 1
ifact(1) = 1
ipow(1) = 1
itot = 0
go to 210
endif
itot = inum
c.... Find maximum possible factor.
itestm = sqrt (float (inumr))
cbugcbugc***DEBUG begins.
cbugcbug 9700 format ('DEBUG 0. nfact, itestm =',i3,21x, i18)
cbugcbug write ( 3, 9700) nfact, itestm
cbugcbugc***DEBUG ends.
c.... Test 2.
itest = 2
if (itest .gt. itestm) then ! 2 is too high.
nfact = nfact + 1
cbugcbugc***DEBUG begins.
cbugcbug 9801 format ('DEBUG A. nfact, itestm =',i3,21x, i18)
cbugcbug write ( 3, 9801) nfact, itestm
cbugcbugc***DEBUG ends.
ifact(nfact) = inumr
ipow(nfact) = 1
itot = (itot / inumr) * (inumr - 1)
inumr = 1
go to 210
endif ! 2 is too high.
if (mod (inumr, itest) .eq. 0) then ! Residual inumr is even.
if (nfact .ge. nfactm) then
nerr = 3
go to 210
endif
nfact = nfact + 1
cbugcbugc***DEBUG begins.
cbugcbug 9802 format ('DEBUG B. nfact, itestm =',i3,21x, i18)
cbugcbug write ( 3, 9802) nfact, itestm
cbugcbugc***DEBUG ends.
ifact(nfact) = itest
ipow(nfact) = 1
itot = (itot / itest) * (itest - 1)
inumr = inumr / itest
if (inumr .eq. 1) go to 210
itestm = sqrt (float (inumr))
cbugcbugc***DEBUG begins.
cbugcbug 9701 format ('DEBUG 1. nfact,ifact,itestm =',i3,i18,'^',i2,i18)
cbugcbug write ( 3, 9701) nfact, ifact(nfact), ipow(nfact), itestm
cbugcbugc***DEBUG ends.
120 if (mod (inumr, itest) .eq. 0) then
ipow(nfact) = ipow(nfact) + 1
inumr = inumr / itest
if (inumr .eq. 1) go to 210
itestm = sqrt (float (inumr))
cbugcbugc***DEBUG begins.
cbugcbug 9702 format ('DEBUG 2. nfact,ifact,itestm =',i3,i18,'^',i2,i18)
cbugcbug write ( 3, 9702) nfact, ifact(nfact), ipow(nfact), itestm
cbugcbugc***DEBUG ends.
go to 120
endif
endif
if (nfact .ge. nfactm) then
nerr = 3
go to 210
endif
c.... Test each higher odd integer.
itest = 3
130 if (itest .gt. itestm) then
nfact = nfact + 1
cbugcbugc***DEBUG begins.
cbugcbug 9803 format ('DEBUG C. nfact, itestm =',i3,21x, i18)
cbugcbug write ( 3, 9803) nfact, itestm
cbugcbugc***DEBUG ends.
ifact(nfact) = inumr
ipow(nfact) = 1
itot = (itot / inumr) * (inumr - 1)
inumr = 1
go to 210
endif
if (mod (inumr, itest) .eq. 0) then ! Test next odd integer.
nfact = nfact + 1
cbugcbugc***DEBUG begins.
cbugcbug 9804 format ('DEBUG C. nfact, itestm =',i3,21x, i18)
cbugcbug write ( 3, 9804) nfact, itestm
cbugcbugc***DEBUG ends.
ifact(nfact) = itest
ipow(nfact) = 1
itot = (itot / itest) * (itest - 1)
inumr = inumr / itest
if (inumr .eq. 1) go to 210
itestm = sqrt (float (inumr))
cbugcbugc***DEBUG begins.
cbugcbug 9703 format ('DEBUG 3. nfact,ifact,itestm =',i3,i18,'^',i2,i18)
cbugcbug write ( 3, 9703) nfact, ifact(nfact), ipow(nfact), itestm
cbugcbugc***DEBUG ends.
140 if (mod (inumr, itest) .eq. 0) then
ipow(nfact) = ipow(nfact) + 1
inumr = inumr / itest
if (inumr .eq. 1) go to 210
itestm = sqrt (float (inumr))
cbugcbugc***DEBUG begins.
cbugcbug 9704 format ('DEBUG 4. nfact,ifact,itestm =',i3,i18,'^',i2,i18)
cbugcbug write ( 3, 9704) nfact, ifact(nfact), ipow(nfact), itestm
cbugcbugc***DEBUG ends.
go to 140
endif
endif ! Tested next odd integer.
if (nfact .ge. nfactm) then
nerr = 3
go to 210
endif
itest = itest + 2
go to 130
210 continue
cbugc***DEBUG begins.
cbug 9902 format (/ 'aptpfac found',i12,' factors. nerr =',i3,'.' )
cbug 9903 format (' n =',i3,' ifact(n) =',i22,'^',i2)
cbug 9904 format (' Totient function = ',i20)
cbug 9905 format (' Unfactored residue =',i20)
cbug write ( 3, 9902) nfact, nerr
cbug if (nfact .gt. 0) then
cbug write ( 3, 9903) (n, ifact(n), ipow(n), n = 1, nfact)
cbug endif
cbug write ( 3, 9904) itot
cbug if (nerr .eq. 3) then
cbug write ( 3, 9905) inumr
cbug endif
cbugc***DEBUG ends.
return
c.... End of subroutine aptpfac. (+1 line.)
end
UCRL-WEB-209832