subroutine aptpfab (idiga, ndiga, nfactm, idigb, idigc,
& nfact, idigd, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTPFAB
c
c call aptpfab (idiga, ndiga, nfactm, idigb, idigd,
c nfact, idigd, nerr)
c
c Version: aptpfab Updated 2006 April 4, 17:30.
c aptpfab Originated 2006 April 4, 17:30.
c
c Authors: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: For the big integer idiga(i), i = 1, ndiga, find any prime
c factors idigb(n) > 1, and their exponents idigd(n), for
c n = 1, nfact, with nfact limited to nfactm.
c If not all prime factors of idiga are returned, because nfactm
c is insufficient, then idigd, the remaining unfactored part of
c idiga, is returned with a value other than 1, with nerr = 3.
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 idiga 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: idiga, ndiga, nfactm.
c
c Output: idigb, idigd, nfact, idigd, nerr.
c
c Calls: (none)
c
c Glossary:
c
c idiga Input Any big integer, for which all prime factors and
c their exponents are to be found, up to a maximum of
c nfactm prime factors. Size at least ndiga.
c
c idigb Output Prime factor idigb(n) is the n'th smallest prime
c factor of idiga, with an exponent of idigc(n).
c If idiga = 1, then idigb(1) = idigc(1) = 0.
c Size at least nfactm. Initially set to -999999999.
c
c idigb Input Any trial factor.
c
c idigc Output Exponent idigc(n) is the largest power of
c idigb(n) that divides idiga.
c Size at least nfactm. Initially set to -999999999.
c
c idigd Output Any remaining unfactored part of idiga, when nfact
c is prevented from exceeding nfactm. Normally 1.
c
c ndiga Input The number of digits in idiga.
c
c nerr Output If not zero, indicates an error.
c 1 if idiga 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 big integer
c idiga. If more than nfactm exist, nerr = 3.
c
c nfactm Input The maximum allowed number of prime factors.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c.... Declarations for arguments.
integer idiga(1) ! A big integer to be factored.
integer idigb(1) ! A prime factor of idiga.
integer idigc(1) ! Power of idigb(n) dividing idiga.
integer idigd(1) ! Result of idiga divided by all factors.
integer nerr ! Indicates an error, if not zero.
integer nfact ! Number of prime factors of idiga found.
integer nfactm ! The maximum allowed value of nfact.
c.... Declarations for local variables.
integer idigb(ndiga) ! Any big integer.
integer idigc(ndiga) ! Any big integer.
integer idigd(ndiga) ! Any big integer.
integer idige(ndiga) ! Any big integer.
integer idigf(ndiga) ! Any big integer.
integer idigg(ndiga) ! Any big integer.
integer idigh(ndiga) ! Any big integer.
integer idigi(ndiga) ! Any big integer.
integer idigj(ndiga) ! Any big integer.
integer idigk(ndiga) ! Any big integer.
integer idigl(ndiga) ! Any big integer.
integer idigm(ndiga) ! Any big integer.
integer idign(ndiga) ! Any big integer.
integer idigo(ndiga) ! Any big integer.
integer idigp(ndiga) ! Any big integer.
integer idigq(ndiga) ! Any big integer.
integer idigr(ndiga) ! Any big integer.
integer idigs(ndiga) ! Any big integer.
integer idigt(ndiga) ! Any big integer.
integer idigu(ndiga) ! Any big integer.
integer idigv(ndiga) ! Any big integer.
integer idigw(ndiga) ! Any big integer.
integer idigx(ndiga) ! Any big integer.
integer idigy(ndiga) ! Any big integer.
integer idigz(ndiga) ! Any big integer.
integer n ! Index in arrays idigb, idigc.
c***DEBUG begins.
9901 format ('' / 'aptpfab finding prime factors of big integer',
& ' with length',i6,' digits,' /
& ' But no more than',i6,' prime factors.' )
9902 format ('big ',11x,i6,' = ',50i1)
write ( 3, 9901) ndiga, nfactm
do nin = 1, ndiga, 50
nlast = min (nin + 49, ndiga)
write ( 3, 9902) nin, (idiga(n), n = nin, nlast)
enddo
c***DEBUG ends.
c.... Initialize.
nerr = 0
nfact = 0
idigd = idiga
if (nfactm .gt. 0) then
do n = 1, nfactm
idigb(n) = -999999999
idigc(n) = -999999999
enddo
endif
c.... Test for input errors.
if (idiga .lt. 1) then
nerr = 1
go to 210
endif
if (nfactm .le. 0) then
nerr = 2
go to 210
endif
if (idiga .eq. 1) then
nfact = 1
idigb(1) = 1
idigc(1) = 1
go to 210
endif
c.... Find maximum possible factor.
idigd = sqrt (idigd)
call aptbrtn (2, 10, idiga, ndiga, idmax, iemax,
& idigs, idigt, idigu, idigv,
& idigw, idigx, idigy, idigz,
& ndigm, idigd, ndigc, isign, nerr)
c***DEBUG begins.
9700 format ('DEBUG 0. nfact, idigd =',i3,21x, i18)
write ( 3, 9700) nfact, idigd
c***DEBUG ends.
c.... Test 2.
idigb = 2
if (idigb .gt. idigd) then ! 2 is too high.
nfact = nfact + 1
c***DEBUG begins.
9801 format ('DEBUG A. nfact, idigd =',i3,21x, i18)
write ( 3, 9801) nfact, idigd
c***DEBUG ends.
idigb(nfact) = idigd
idigc(nfact) = 1
idigd = 1
go to 210
endif ! 2 is too high.
if (mod (idigd, idigb) .eq. 0) then ! Residual idigd is even.
if (nfact .ge. nfactm) then
nerr = 3
go to 210
endif
nfact = nfact + 1
c***DEBUG begins.
9802 format ('DEBUG B. nfact, idigd =',i3,21x, i18)
write ( 3, 9802) nfact, idigd
c***DEBUG ends.
idigb(nfact) = idigb
idigc(nfact) = 1
idigd = idigd / idigb
if (idigd .eq. 1) go to 210
idigd = sqrt (float (idigd))
c***DEBUG begins.
9701 format ('DEBUG 1. nfact,idigb,idigc =',i3,i18,'^',i2,i18)
write ( 3, 9701) nfact, idigb(nfact), idigd(nfact), idigd
c***DEBUG ends.
120 if (mod (idigd, idigb) .eq. 0) then
idigc(nfact) = idigc(nfact) + 1
idigd = idigd / idigb
if (idigd .eq. 1) go to 210
idigd = sqrt (float (idigd))
c***DEBUG begins.
9702 format ('DEBUG 2. nfact,idigb,idigd =',i3,i18,'^',i2,i18)
write ( 3, 9702) nfact, idigb(nfact), idigc(nfact), idigc
c***DEBUG ends.
go to 120
endif
endif
if (nfact .ge. nfactm) then
nerr = 3
go to 210
endif
c.... Test each higher odd integer.
idigb = 3
130 if (idigb .gt. idigd) then
nfact = nfact + 1
c***DEBUG begins.
9803 format ('DEBUG C. nfact, idigd =',i3,21x, i18)
write ( 3, 9803) nfact, idigd
c***DEBUG ends.
idigb(nfact) = idigd
idigc(nfact) = 1
idigd = 1
go to 210
endif
if (mod (idigd, idigb) .eq. 0) then ! Test next odd integer.
nfact = nfact + 1
c***DEBUG begins.
9804 format ('DEBUG C. nfact, idigd =',i3,21x, i18)
write ( 3, 9804) nfact, idigd
c***DEBUG ends.
idigb(nfact) = idigb
idigc(nfact) = 1
idigd = idigd / idigb
if (idigd .eq. 1) go to 210
idigd = sqrt (float (idigd))
c***DEBUG begins.
9703 format ('DEBUG 3. nfact,idigb,idigd =',i3,i18,'^',i2,i18)
write ( 3, 9703) nfact, idigb(nfact), idigc(nfact), idigc
c***DEBUG ends.
140 if (mod (idigd, idigb) .eq. 0) then
idigd(nfact) = idigd(nfact) + 1
idigc = idigc / idigb
if (idigd .eq. 1) go to 210
idigd = sqrt (float (idigd))
c***DEBUG begins.
9704 format ('DEBUG 4. nfact,idigb,idigd =',i3,i18,'^',i2,i18)
write ( 3, 9704) nfact, idigb(nfact), idigc(nfact), idigc
c***DEBUG ends.
go to 140
endif
endif ! Tested next odd integer.
if (nfact .ge. nfactm) then
nerr = 3
go to 210
endif
call aptbadd (10, idigb, ndigb, idigb, ndigb,
& ndigm, idigb, ndigb, nerr)
go to 130
210 continue
c***DEBUG begins.
9902 format (/ 'aptpfab found',i12,' factors. nerr =',i3,'.' )
9903 format (' n =',i3,' idigb(n) =',i22,'^',i2)
9905 format (' Unfactored residue =',i20)
write ( 3, 9902) nfact, nerr
if (nfact .gt. 0) then
write ( 3, 9903) (n, idigb(n), idigc(n), n = 1, nfact)
endif
if (nerr .eq. 3) then
write ( 3, 9905) idigd
endif
c***DEBUG ends.
return
c.... End of subroutine aptpfab. (+1 line.)
end
UCRL-WEB-209832