subroutine aptpfft (ifact, npmax, iprime, ipow, np, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTPFFT
c
c call aptpfft (ifact, npmax, iprime, ipow, np, nerr)
c
c Version: aptpfft Updated 2003 September 26 14:00.
c aptpfft Originated 2003 September 26 14:00.
c
c Authors: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: For the factorial function ifact!, find all of the prime factors
c iprime(n), and their powers ipow(n), n = 1, np <= npmax, such
c that if np is not limited by npmax,
c ifact! = product {iprime(n)^ipow(n)}, n = 1, np,
c with iprime(np) <= ifact.
c
c Flag nerr indicates any input error or limited results.
c
c Input: ifact, npmax
c
c Output: iprime, np, nerr
c
c Calls: (none)
c
c Glossary:
c
c ifact Input The argument of the factorial function ifact!, for
c which all prime factors and their powers are to be
c found, up to a limit of npmax prime factors.
c Unless limited by npmax, the largest prime factor
c will be the largest one less than or equat to ifact.
c If np is not limited by npmax,
c ifact! = product {iprime(n)^ipow(n)}, n = 1, np,
c with iprime(np) <= ifact.
c
c ipow(n) Output The power of prime factor iprime(n).
c Size at least npmax.
c
c iprime(n) Output The n'th prime number. Size at least npmax.
c None will be returned if ifact < 1 (nerr = 1) or
c npmax < 1 (nerr = 2).
c
c The density of prime numbers from 1 to iprime(n) is
c d = n/iprime(n). For large iprime(n),
c d ~ dlim = 1/ln(iprime(n)), and
c n ~ nlim = iprime(n)/ln (iprime(n)), and
c d/dlim = n/nlim = n*ln(iprime(n))/iprime(n).
c
c d = dlim = d/dlim =
c n iprime n/iprime 1/ln(iprime) n/nlim
c ----- ------ -------- ------------ --------
c 1 1 1.000 infinite 0.000
c 2 2 1.000 1.443 0.683
c 3 3 1.000 0.910 1.099
c 4 5 0.800 0.621 1.288
c 5 7 0.714 0.514 1.390
c 10 23 0.435 0.319 1.363
c 30 109 0.275 0.213 1.291
c 100 523 0.191 0.160 1.197
c 300 1979 0.152 0.132 1.150
c 1000 7907 0.126 0.111 1.135
c 3000 27437 0.109 0.098 1.117
c 10000 104723 0.095 0.087 1.104
c 30000 350351 0.086 0.078 1.093
c
c nerr Output Indicates an input error, if not zero:
c -1 if npmax prime numbers and their powers are found
c without reaching the highest prime number not
c exceeding ifact.
c 1 if ifact is less than 1.
c 2 if npmax is less than 1.
c
c npmax Input The maximum number of sequential prime numbers to find,
c from 1 up to but not exceeding ifact.
c For large ifact, the necessary npmax is at least
c nlim = iprime/ln (iprime). See the table under
c "iprime(n)" above for factors n/nlim to estimate a
c sufficient value of npmax.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c.... Declarations for arguments.
integer ifact ! Argument of factorial function.
integer ipow(1) ! An array of up to npmax powers.
integer iprime(1) ! An array of up to npmax prime numbers.
integer np ! The actual number of prime numbers.
integer npmax ! Maximum number of prime numbers allowed.
integer nerr ! Indicates error, if not zero.
c.... Declarations for local variables.
integer int ! A number to be tested.
integer i ! Index of preceding prime number.
integer iadd ! Integer added to ipow.
integer n ! Index of current prime number.
integer nterm ! Term in addition to ipow.
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptpfft finding prime factors and powers of'
cbug & ' factorial',i13 /
cbug $ 'but no more than',i13,' prime factors.' )
cbug write ( 3, 9901) ifact, npmax
cbugc***DEBUG ends.
c.... Initialize.
nerr = 0
np = 0
if (npmax .gt. 0) then
do 100 n = 1, npmax
iprime(n) = -999999999
ipow(n) = -999999999
100 continue
endif
c.... Test for input errors.
if (ifact .lt. 1) then
nerr = 1
go to 210
endif
if (npmax .lt. 1) then
nerr = 2
go to 210
endif
c.... Set the first prime factor and its power to 1.
np = 1
iprime(1) = 1
ipow(1) = 1
if (ifact .eq. 1) go to 210
if (npmax .eq. 1) go to 210
c.... Set the second prime factor to 2, and find its power.
np = 2
iprime(2) = 2
ipow(2) = 0
nterm = 0
105 nterm = nterm + 1
iadd = ifact / iprime(2)**nterm
ipow(2) = ipow(2) + iadd
if (iadd .gt. 0) go to 105
if (ifact .eq. 2) go to 210
if (npmax .eq. 2) go to 210
c.... Loop over higher integers to find prime numbers and powers.
int = 2
110 int = int + 1
c.... See if maximum prime number has been reached.
if (int .gt. ifact) go to 210
c.... See if the current test integer has any prime factors among the prime
c.... numbers less than the square root of the integer.
do 120 i = 2, np
if (iprime(i)**2 .gt. int) go to 130 ! Integer int is prime.
if (mod (int, iprime(i)) .eq. 0) go to 110 ! Integer int is not prime.
120 continue
c.... The current test integer is prime.
130 np = np + 1
iprime(np) = int
c.... Find the power for this prime.
ipow(np) = 0
nterm = 0
140 nterm = nterm + 1
iadd = ifact / iprime(np)**nterm
ipow(np) = ipow(np) + iadd
if (iadd .gt. 0) go to 140
c.... See if npmax primes have been found.
if (np .lt. npmax) go to 110
c.... Reached npmax without reaching largest prime not exceeding ifact.
nerr = -1
210 continue
cbugc***DEBUG begins.
cbug 9902 format (/ 'aptpfft found',i12,' primes and powers. nerr =',
cbug & i3,'.' //
cbug & ' np iprime ipow' / )
cbug 9903 format (3i13)
cbug write ( 3, 9902) np, nerr
cbug if (nerr .le. 0) then
cbug write ( 3, 9903) (i, iprime(i), ipow(i), i = 1, np)
cbug endif
cbugc***DEBUG ends.
return
c.... End of subroutine aptpfft. (+1 line.)
end
UCRL-WEB-209832