subroutine aptcrts (p, r, np, pp, ppop, t, x, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTCRTS
c
c call aptcrts (p, r, np, pp, ppop, t, x, nerr)
c
c Version: aptcrts Updated 2003 July 29 12:00.
c aptcrts Originated 2003 July 29 12:00.
c
c Authors: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: Find the smallest integer x which has the non-negative integer
c remainders r(1), r(2), r(3), ..., r(np), when divided by the
c integers p(1), p(2), p(3), ..., p(np). No two p values may have
c common factors. Best if all are prime.
c
c Flag nerr indicates any input error.
c
c Input: p, r, np
c
c Output: x, nerr
c
c Calls: (none)
c
c Glossary:
c
c nerr Output Indicates an input error or incomplete result,
c if not zero:
c 1 if np is less than 2.
c 2 if any of the p values are less than 2.
c 3 if any of the r values are negative.
c 4 if any of the r values exceed p - 1.
c 5 if any of the p values have a common factor.
c 6 if the product of the p values exceeds 10^18.
c
c np Input The number of p values. Must be positive.
c
c p(n) Input An integer divisor of x that produces a non-negative
c integer remainder of r(n), n = 1, np. Size => np.
c None may have common factors.
c
c pp Output The product of all of the p(n) values, n = 1, np.
c Must not exceed the maximum integer value for
c the machine. Limited here to < 10^18.
c
c ppop(n) Output The integer ratio pp / p(n). Size => np.
c
c r(n) Input The non-negative integer remainder when dividing
c x by p(n), n = 1, np. Size => np.
c
c t(n) Output Integer satisfying mod (t(n) * ppop(n), p(n)) = 1,
c n = 1, np. Size => np.
c
c x Output The smallest integer that has a remainder of r(n)
c when divided by p(n), n = 1, np.
c For the same divisors and any other set of
c remainders r(n), n = 1, np,
c x = mod {sum [r(n)*ppop(n)*t(n), n = 1, np], pp}.
c All bigger integer solutions may be obtained by
c repeatedly adding pp.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c.... Declarations for arguments.
integer nerr ! Indicates an error, if not zero.
integer np ! The number of divisor-remainder pairs.
integer p(1) ! An array of np divisors.
integer pp ! Product [p(n), n = 1, np].
integer ppop(1) ! Ratio pp / p(n), n = 1, np.
integer r(1) ! An array of np remainders.
integer t(1) ! Satisfies mod (t(n)*ppop(n), p(n)) = 1.
integer x ! The desired integer.
c.... Declarations for local variables.
integer imax ! Integer <= square root of biggest p(n).
integer n ! Index of a divisor, n = 1, np.
integer nn ! Index of a divisor, n = 1, np - 1.
real fp ! Floating point value of p(n).
real fpp ! Floating point value of pp.
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptcrts finding solution to the Chinese',
cbug & ' Remainder Theorem. np =',i5 )
cbug 9902 format ('Divisor =',i9,' Remainder =',i9)
cbug write ( 3, 9901) np
cbug if (np .ge. 2) then
cbug write ( 3, 9902) (p(n), r(n), n = 1, np)
cbug endif
cbugc***DEBUG ends.
c.... Initialize.
nerr = 0
x = -999999
pp = -999999
c.... Test for input errors.
if (np .lt. 2) then
nerr = 1
go to 210
endif
do 105 n = 1, np
if (p(n) .lt. 2) then
nerr = 2
go to 210
endif
if (r(n) .ge. p(n)) then
nerr = 4
go to 210
endif
105 continue
imax = p(1)
do 115 n = 2, np
if (p(n) .gt. imax) imax = p(n)
115 continue
do 130 n = 2, np
do 125 nn = 1, n - 1
do 120 nnn = 2, imax
if ((mod (p(n), nnn) .eq. 0) .and.
& (mod (p(nn), nnn) .eq. 0)) then
nerr = 5
go to 210
endif
120 continue
125 continue
130 continue
do 140 n = 1, np
if (r(n) .lt. 0) then
nerr = 3
go to 210
endif
140 continue
c.... Find the product pp of all the divisors.
pp = 1
fpp = 1.0
cbugcbugc***DEBUG begins.
cbugcbug 9701 format ('aptcrts DEBUG 1. pp =',i9)
cbugcbug write ( 3, 9701) pp
cbugcbugc***DEBUG ends.
do 150 n = 1, np
fp = p(n)
fpp = fpp * fp
if (fpp .ge. 1.e18) then
nerr = 6
go to 210
endif
pp = pp * p(n)
cbugcbugc***DEBUG begins.
cbugcbug 9702 format ('aptcrts DEBUG 2. pp =',i9)
cbugcbug write ( 3, 9702) pp
cbugcbugc***DEBUG ends.
150 continue
c.... Find the ratios ppop(n) = pp/p(n) for n = 1, np.
do 170 n = 1, np
ppop(n) = pp / p(n)
cbugcbugc***DEBUG begins.
cbugcbug 9703 format ('aptcrts DEBUG 3. n =',i5,' ppop =',i9)
cbugcbug write ( 3, 9703) n, ppop(n)
cbugcbugc***DEBUG ends.
t(n) = 0
160 t(n) = t(n) + 1
if (mod (t(n) * ppop(n), p(n)) .ne. 1) go to 160
cbugcbugc***DEBUG begins.
cbugcbug 9708 format ('aptcrts DEBUG 8. n =',i5,' ppop =',i9,' t =',i9)
cbugcbug write ( 3, 9708) n, ppop(n), t(n)
cbugcbugc***DEBUG ends.
170 continue
c.... Find an x that satisfies requirements.
x = 0
cbugcbugc***DEBUG begins.
cbugcbug 9704 format ('aptcrts DEBUG 4. x =',i9)
cbugcbug write ( 3, 9704) x
cbugcbugc***DEBUG ends.
do 180 n = 1, np
x = x + r(n) * ppop(n) * t(n)
cbugcbugc***DEBUG begins.
cbugcbug 9705 format ('aptcrts DEBUG 5. x =',i9)
cbugcbug write ( 3, 9705) x
cbugcbugc***DEBUG ends.
180 continue
c.... Reduce x to the smallest possible value.
x = mod (x, pp)
cbugcbugc***DEBUG begins.
cbugcbug 9706 format ('aptcrts DEBUG 6. x =',i9)
cbugcbug write ( 3, 9706) x
cbugcbugc***DEBUG ends.
210 continue
cbugc***DEBUG begins.
cbug 9903 format (/ 'aptcrts done. nerr =',i2)
cbug 9904 format ('n =',i5', ppop =',i9,' t =',i9)
cbug 9905 format ('x =',i9)
cbug write ( 3, 9903) nerr
cbug if (np .ge. 2) then
cbug write ( 3, 9904) (n, ppop(n), t(n), n = 1, np)
cbug write ( 3, 9905) x
cbug endif
cbugc***DEBUG ends.
return
c.... End of subroutine aptcrts. (+1 line.)
end
UCRL-WEB-209832