subroutine aptpfac (inum, nfactm, ifact, ipow, nfact, & itot, inumr, nerr) ccbeg. cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c SUBROUTINE APTPFAC c c call aptpfac (inum, nfactm, ifact, ipow, nfact, itot, inumr, nerr) c c Version: aptpfac Updated 2003 September 4 15:00. c aptpfac Originated 2003 August 29 13:40. c c Authors: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123. c c c Purpose: For the integer inum, find any prime factors ifact(n) > 1, c and their exponents ipow(n), for n = 1, nfact, with nfact c limited to nfactm. If not all prime factors of inum are c returned, because nfactm is insufficient, then inumr, the c remaining unfactored part of inum, is returned with a value c other than 1, with nerr = 3. c c Also find itot, Euler's totient function for modulo inum c arithmetic, the number of integers in the set (1, ..., inum - 1) c that are relatively prime to inum (have no common factor other c than 1). 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 inum 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: inum, nfactm. c c Output: ifact, ipow, nfact, itot, inumr, nerr. c c Calls: (none) c c Glossary: c c ifact Output Prime factor ifact(n) is the n'th smallest prime c factor of inum, with an exponent of ipow(n). c If inum = 1, then ifact(1) = ipow(1) = 1, itot = 0. c Size at least nfactm. Initially set to -999999999. c c inum Input Any positive integer, for which all prime factors and c their exponents are to be found, up to a maximum of c nfactm prime factors. c c inumr Output Any remaining unfactored part of inum, when nfact c is prevented from exceeding nfactm. Normally 1. c c ipow Output Exponent ipow(n) is the largest power of ifact(n) c that divides inum. Initially set to -999999999. c Size at least nfactm. c c itot Output Euler's totient function for modulo inum arithmetic. c The number of integers in the set c (1, ..., inum - 1) that are relatively prime to c inum (have no common factors other then 1). c Initially set to -999999999. c c nerr Output If not zero, indicates an error. c 1 if inum 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 integer inum. c If nfactm are found, and more exist, nerr = 3. c c nfactm Input The maximum allowed size of arrays ifact and ipow. c On a machine with 64-bit integers, 18 is enough. c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ccend. c.... Declarations for arguments. integer ifact(1) ! A prime factor of inum. integer inum ! Any positive integer. integer inumr ! Result of inum divided by all factors. integer ipow(1) ! Largest power of ifact(n) dividing inum. integer itot ! Euler's totient function for mod inum. integer nerr ! Indicates an error, if not zero. integer nfact ! How many prime factors of inum found. integer nfactm ! The maximum allowed value of nfact. c.... Declarations for local variables. integer itest ! Integer to test as factor. integer itestm ! Maximum integer to test as factor. integer n ! Index in arrays ifact, ipow. cbugc***DEBUG begins. cbug 9901 format (/ 'aptpfac finding prime factors of inum =',i20 / cbug & ' but no more than',i6,' prime factors.' ) cbug write ( 3, 9901) inum, nfactm cbugc***DEBUG ends. c.... Initialize. nerr = 0 nfact = 0 inumr = inum if (nfactm .gt. 0) then do n = 1, nfactm ifact(n) = -999999999 ipow(n) = -999999999 enddo endif itot = -999999999 c.... Test for input errors. if (inum .lt. 1) then nerr = 1 go to 210 endif if (nfactm .le. 0) then nerr = 2 go to 210 endif if (inum .eq. 1) then nfact = 1 ifact(1) = 1 ipow(1) = 1 itot = 0 go to 210 endif itot = inum c.... Find maximum possible factor. itestm = sqrt (float (inumr)) cbugcbugc***DEBUG begins. cbugcbug 9700 format ('DEBUG 0. nfact, itestm =',i3,21x, i18) cbugcbug write ( 3, 9700) nfact, itestm cbugcbugc***DEBUG ends. c.... Test 2. itest = 2 if (itest .gt. itestm) then ! 2 is too high. nfact = nfact + 1 cbugcbugc***DEBUG begins. cbugcbug 9801 format ('DEBUG A. nfact, itestm =',i3,21x, i18) cbugcbug write ( 3, 9801) nfact, itestm cbugcbugc***DEBUG ends. ifact(nfact) = inumr ipow(nfact) = 1 itot = (itot / inumr) * (inumr - 1) inumr = 1 go to 210 endif ! 2 is too high. if (mod (inumr, itest) .eq. 0) then ! Residual inumr is even. if (nfact .ge. nfactm) then nerr = 3 go to 210 endif nfact = nfact + 1 cbugcbugc***DEBUG begins. cbugcbug 9802 format ('DEBUG B. nfact, itestm =',i3,21x, i18) cbugcbug write ( 3, 9802) nfact, itestm cbugcbugc***DEBUG ends. ifact(nfact) = itest ipow(nfact) = 1 itot = (itot / itest) * (itest - 1) inumr = inumr / itest if (inumr .eq. 1) go to 210 itestm = sqrt (float (inumr)) cbugcbugc***DEBUG begins. cbugcbug 9701 format ('DEBUG 1. nfact,ifact,itestm =',i3,i18,'^',i2,i18) cbugcbug write ( 3, 9701) nfact, ifact(nfact), ipow(nfact), itestm cbugcbugc***DEBUG ends. 120 if (mod (inumr, itest) .eq. 0) then ipow(nfact) = ipow(nfact) + 1 inumr = inumr / itest if (inumr .eq. 1) go to 210 itestm = sqrt (float (inumr)) cbugcbugc***DEBUG begins. cbugcbug 9702 format ('DEBUG 2. nfact,ifact,itestm =',i3,i18,'^',i2,i18) cbugcbug write ( 3, 9702) nfact, ifact(nfact), ipow(nfact), itestm cbugcbugc***DEBUG ends. go to 120 endif endif if (nfact .ge. nfactm) then nerr = 3 go to 210 endif c.... Test each higher odd integer. itest = 3 130 if (itest .gt. itestm) then nfact = nfact + 1 cbugcbugc***DEBUG begins. cbugcbug 9803 format ('DEBUG C. nfact, itestm =',i3,21x, i18) cbugcbug write ( 3, 9803) nfact, itestm cbugcbugc***DEBUG ends. ifact(nfact) = inumr ipow(nfact) = 1 itot = (itot / inumr) * (inumr - 1) inumr = 1 go to 210 endif if (mod (inumr, itest) .eq. 0) then ! Test next odd integer. nfact = nfact + 1 cbugcbugc***DEBUG begins. cbugcbug 9804 format ('DEBUG C. nfact, itestm =',i3,21x, i18) cbugcbug write ( 3, 9804) nfact, itestm cbugcbugc***DEBUG ends. ifact(nfact) = itest ipow(nfact) = 1 itot = (itot / itest) * (itest - 1) inumr = inumr / itest if (inumr .eq. 1) go to 210 itestm = sqrt (float (inumr)) cbugcbugc***DEBUG begins. cbugcbug 9703 format ('DEBUG 3. nfact,ifact,itestm =',i3,i18,'^',i2,i18) cbugcbug write ( 3, 9703) nfact, ifact(nfact), ipow(nfact), itestm cbugcbugc***DEBUG ends. 140 if (mod (inumr, itest) .eq. 0) then ipow(nfact) = ipow(nfact) + 1 inumr = inumr / itest if (inumr .eq. 1) go to 210 itestm = sqrt (float (inumr)) cbugcbugc***DEBUG begins. cbugcbug 9704 format ('DEBUG 4. nfact,ifact,itestm =',i3,i18,'^',i2,i18) cbugcbug write ( 3, 9704) nfact, ifact(nfact), ipow(nfact), itestm cbugcbugc***DEBUG ends. go to 140 endif endif ! Tested next odd integer. if (nfact .ge. nfactm) then nerr = 3 go to 210 endif itest = itest + 2 go to 130 210 continue cbugc***DEBUG begins. cbug 9902 format (/ 'aptpfac found',i12,' factors. nerr =',i3,'.' ) cbug 9903 format (' n =',i3,' ifact(n) =',i22,'^',i2) cbug 9904 format (' Totient function = ',i20) cbug 9905 format (' Unfactored residue =',i20) cbug write ( 3, 9902) nfact, nerr cbug if (nfact .gt. 0) then cbug write ( 3, 9903) (n, ifact(n), ipow(n), n = 1, nfact) cbug endif cbug write ( 3, 9904) itot cbug if (nerr .eq. 3) then cbug write ( 3, 9905) inumr cbug endif cbugc***DEBUG ends. return c.... End of subroutine aptpfac. (+1 line.) end UCRL-WEB-209832