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