subroutine aptprim (ipmax, npmax, nd, ip, fp, np, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTPRIM
c
c call aptprim (ipmax, npmax, nd, ip, fp, np, nerr)
c
c Version: aptprim Updated 2003 June 20 15:00.
c aptprim Originated 2003 June 10 16:00.
c
c Authors: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: Find the prime numbers ip(n), n = 1, np <= npmax, such that
c ip(np) <= ipmax, and if power nd > 0, find the prime number
c functions fp(n) = Product {1 - 1/ip(i)**nd}, i = 2, n,
c for n = 1, np.
c
c For nd = 1, fp(n) is the approximate density of prime numbers
c among the integers from ip(n) + 1 to ip(n+1)**2 - 1, or less
c accurately, from 1 to ip(n)**2.
c
c For nd > 1, fp(n) is the approximate density of irreducible
c proportions (no common factors) between all possible
c combinations of nd integers (k, l, m, ...) from (1, 1, 1, ...)
c to (K, L, M, ...), with ip(n) no larger then the largest
c possible common prime factor of k, l, m, ...
c
c Flag nerr indicates any input error or limited results.
c
c Input: ipmax, npmax, nd
c
c Output: ip, fp, np, nerr
c
c Calls: (none)
c
c Glossary:
c
c fp(n) Output If nd > 0, function fp(1) = 1, and fp(n) is the
c product of all factors (1 - 1/ip(i)**nd), i = 2, n,
c for n = 2, np <= npmax.
c Size at least npmax, if nd > 0. Otherwise, zero.
c
c For nd = 1 and very large n, fp(n) ~ 0.5/ln(ip(n)).
c For nd = 2 and n > 10000, fp(n) ~ 0.60792756,
c For nd = 3 and n > 900, fp(n) ~ 0.831907373.
c For nd = 4 and n > 90, fp(n) ~ 0.923938403.
c For nd = 5 and n > 40, fp(n) ~ 0.964387340.
c
c For nd = 1, fp(n) is the approximate density of
c prime numbers among the integers from ip(n) + 1
c to ip(n+1)**2 - 1, or less accurately, from 1 to
c ip(n)**2.
c
c For nd > 1, fp(n) is the approximate density of
c irreducible proportions (no common factors) between
c all possible combinations of nd integers
c (k, l, m, ...) from (1, 1, 1, ...) to (K, L, M, ...),
c with ip(n) no larger then the largest possible
c common prime factor of k, l, m, ...
c
c If R(k, l, m, ...) is a point in nd-dimensional
c space, with R(1, 1, 1, ...) defined as the origin,
c and any increment in one of the indices k, l, m, ...
c represents a spatial displacement unique to that
c index, then fp(n) is the approximate fraction of all
c points R for which a line from point R to the origin
c does not pass through another point nearer the
c origin. Each such point may define a unique
c direction from the origin through multiple points
c farther from the origin.
c
c ip(n) Output The n'th prime number. Size at least npmax.
c None will be returned if nd < 1 (nerr = 1), or
c ipmax < 1 (nerr = 2) or npmax < 1 (nerr = 3).
c
c The density of prime numbers from 1 to ip(n) is
c d = n/ip(n). For large ip(n),
c d ~ dlim = 1/ln(ip(n)), and
c n ~ nlim = ip(n)/ln (ip(n)), and
c d/dlim = n/nlim = n*ln(ip(n))/ip(n).
c
c d = dlim = d/dlim =
c n ip n/ip 1/ln(ip) 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 ipmax Input The limit on the size of the maximum prime number
c ip(n) to find, if possible without n exceeding npmax.
c If ipmax is not prime, the next lower prime number
c will be the actual limit.
c
c nd Input The power of ip to be used in function fp.
c If nd = 0, no function fp will be calculated.
c It may be the dimension of a symmetric array of
c points for which unique directions are to be found.
c
c nerr Output Indicates an input error, if not zero:
c -1 if npmax prime numbers are found without reaching
c the highest prime number not exceeding ipmax.
c 1 if nd is less than 0.
c 2 if ipmax is less than 1.
c 3 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 ipmax.
c For large ipmax, the necessary npmax is at least
c nlim = ip/ln (ip). See the table under "ip(n)"
c above for factors n/nlim to estimate a sufficient
c value of npmax.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c.... Declarations for arguments.
integer ip(1) ! An array of up to npmax prime numbers.
integer ipmax ! Upper limit on size of ip.
integer np ! The actual number of prime numbers.
integer nd ! Power of ip in function fp.
integer npmax ! Maximum number of prime numbers to find.
integer nerr ! Indicates error, if not zero.
real fp(1) ! Product (1 - 1/ip(i)**nd)), i = 2, n.
c.... Declarations for local variables.
integer int ! A number to be tested.
integer i ! Index of preceding prime number.
integer n ! Index of current prime number.
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptprim finding prime numbers from 1 to',i13 /
cbug $ 'but no more than',i13,' prime numbers, and the products' /
cbug & 'fp(n) = (1-1/ip(i)**k), i = 2 to n, for nd =',i3)
cbug write ( 3, 9901) ipmax, npmax, nd
cbugc***DEBUG ends.
c.... Initialize.
nerr = 0
np = 0
if (npmax .gt. 0) then
do 100 n = 1, npmax
ip(n) = -999999999
100 continue
if (nd .gt. 0) then
do 105 n = 1, npmax
fp(n) = -1.e99
105 continue
endif
endif
c.... Test for input errors.
if (nd .lt. 0) then
nerr = 2
go to 210
endif
if (ipmax .lt. 1) then
nerr = 2
go to 210
endif
if (npmax .lt. 1) then
nerr = 3
go to 210
endif
c.... Find the first prime number.
np = 1
ip(1) = 1
if (ipmax .eq. 1) go to 150
if (npmax .eq. 1) go to 150
c.... Find the second prime number.
np = 2
ip(2) = 2
if (ipmax .eq. 2) go to 150
if (npmax .eq. 2) go to 150
c.... Loop over higher integers to find prime numbers and functions.
int = 2
110 int = int + 1
c.... See if maximum prime number has been reached.
if (int .gt. ipmax) go to 150
c.... Search for factors among prime numbers less than square root of
c.... current integer.
do 120 i = 2, np
if (ip(i)**2 .gt. int) go to 130 ! Integer int is prime.
if (mod (int, ip(i)) .eq. 0) go to 110 ! Integer int is not prime.
120 continue
c.... Found a prime.
130 np = np + 1
ip(np) = int
c.... See if npmax primes have been found.
if (np .lt. npmax) go to 110
c.... Reached npmax without reaching largest prime not exceeding ipmax.
nerr = -1
c.... Find the prime number functions.
150 if (nd .eq. 0) go to 210
fp(1) = 1.0
do 160 n = 2, np
fip = ip(n)
fp(n) = fp(n-1) * (1.0 - 1.0 / fip**nd)
160 continue
210 continue
cbugc***DEBUG begins.
cbug 9902 format (/ 'aptprim found',i12,' primes. nerr =',i3,'.' //
cbug & 'np, ip=' / )
cbug 9903 format (3(i6,i7,f13.9))
cbug write ( 3, 9902) np, nerr
cbug if (nerr .le. 0) then
cbug nn = 1
cbug 220 nmax = nn + 2
cbug if (nmax .gt. np) nmax = np
cbug write ( 3, 9903) (i, ip(i), fp(i), i = nn, nmax)
cbug nn = nn + 3
cbug if (nn .le. np) go to 220
cbug endif
cbugc***DEBUG ends.
return
c.... End of subroutine aptprim. (+1 line.)
end
UCRL-WEB-209832