################################### C:\python27\main2.py ############################## import cFract2, time, random tBegin = time.clock() pList = cFract2.primes(pow(2,22)) # Get list of primes < pow(2,22) topPrime,modulus = 183,pow(2,106) - 1 # topPrime,modulus = 29,2979553936 print "Number to be factored: ",modulus i,pSlice = 1,[pList[0]] # The prime 2 is handled separately. while pList[i] <= topPrime: if pow(modulus, (pList[i]-1)/2, pList[i]) <= 1: pSlice.append(pList[i]) i += 1 print "Factor base: primes <= ",topPrime print "Length of factor base is: ",len(pSlice) print i,iEnd,outNums = 0,len(pList),[] while i < iEnd: # Extract small factors quotient = modulus/pList[i] # Use pList, NOT pSlice! if modulus == pList[i]*quotient: outNums.append(pList[i]) modulus = quotient else: i = i+1 wrkStack = [modulus] recNo = [0] while len(wrkStack) > 0: modulus = wrkStack.pop() if cFract2.prCheck(modulus,10)==1: if modulus > 1: outNums.append(modulus) else: factor = cFract2.getFact(modulus,pSlice,recNo,topPrime) wrkStack.append(factor) modulus = modulus/factor wrkStack.append(modulus) outNums.sort() tEnd = time.clock() print "Elapsed time: ",tEnd-tBegin," sec" print "Record number: ",recNo[0] print "Factors: ",outNums[:] ###################################### C:\python27\cFract2.py ######################### #- Programs: primes, getFlist, newRec, factSqrt, combo, getFact, prCheck # Programmer: Terry Cooper # Last modified: 11/05/2012 import math,random,time,fractions,shelve def primes(n): # Return a list of primes less than or equal to some integer n if n==2: # Make n = 2 a valid request return [2] # Return list of primes elif n<2: # Make n < 2 a valid reques return [] # Return list of primes s = range(3,n+1,2) # This gets the space required for the sieve, and initializes it mroot = n**0.5 # Result is rounded down since n is an integer half = (n+1)/2 - 1 # Length of sieve, since only odd numbers are sieved i, m = 0, 3 # Initialize i and m for "while" loop below while m <= mroot: # "3" is first prime processed; "2" is handled separately if s[i]: # If number in location i is composite, don't bother processing it j = (m*m-3)/2 # Lower composite numbers have been dealt with while j < half: # Set composite locations equal to zero s[j] = 0 # Multiples of m are composite, so set sieve value equal to zero j += m # Next value of loop index i += 1 # Not a loop index, but it needs to be incremented anyway m = 2*i + 3 # Loop index. "Points" to beginning of next work area return [2]+[x for x in s if x] # Return only prime numbers def newRec(dBasis,topPrime,fList,modulus,recNo): i = 1 # We do NOT create a static constant for len(fList), since len(fList) can increase. while i < len(fList): # "combo" never messes up values of i that have already been looked at, so we can actually do this. if fList[i][1] % 2: iSub = fList[i][0] if iSub in dBasis: good = diagnose(dBasis[iSub],modulus,recNo) if good: combine(fList,dBasis,iSub,modulus,recNo) else: del dBasis[iSub] i -= 1 else: good = diagnose(fList,modulus,recNo) if good: dBasis[iSub] = fList break i += 1 if i >= len(fList): # Normal exit from "while" loop if fList[0][1] == -1: if 1 in dBasis: good = diagnose(dBasis[1],modulus,recNo) if good: combine(fList,dBasis,1,modulus,recNo) else: del dBasis[1] good = diagnose(fList,modulus,recNo) if good: dBasis[1] = fList else: good = diagnose(fList,modulus,recNo) if good: dBasis[1] = fList if fList[0][1] == 1: # fList[0][1] may very well have changed since the last test. good = diagnose(fList,modulus,recNo) if good: dBasis[0] = fList return(fList) def factSqrt(fList,modulus,recNo): good = diagnose(fList,modulus,recNo) if not good: abort() iEnd = len(fList) # At this point in time, len(fList) is stable i = fSqrt = 1 # More initialization while i < iEnd: # Get fSqrt (mod modulus) fSqrt *= pow(fList[i][0],fList[i][1]/2,modulus) i += 1 if fSqrt >= fList[0][0]: # Get absolute value arg2 = fSqrt - fList[0][0] # Guarantees arg2 < modulus else: arg2 = fList[0][0] - fSqrt return fractions.gcd(modulus,arg2) # 1 => solution is worthless; won't return modulus since arg2 < modulus def getFact(modulus,pSlice,recNo,topPrime): dBasis = shelve.open('cFract2.tmp','n') i,testVal,dBasis,oldTime = 0,1,{},int(time.clock()) # Initialize g = int(math.sqrt(modulus)) # If modulus = g*g, return g as factor if (g*g == modulus): return g capQnM2, capQnM1, qNm1, rNm2, rNm1, aNm2, aNm1, sign = modulus, 1, g, g, 0, 1, g, -1 while ((testVal == 1) | (testVal == modulus)): while 0 not in dBasis: # Possible solution not yet obtained capQn = capQnM2 + qNm1 * (rNm1 - rNm2) capGn = 2 * g - rNm1 qN = capGn / capQn rN = capGn - qN * capQn aN = (qN * aNm1 + aNm2)%modulus if (recNo[0] == 0): aNsave,qNsave = aN,qN else: if ((aNsave == aN) & (qNsave == qN)): print "Run aborted: continued fraction is repeating, with a cycle of length ",recNo[0] abort() k = capQn # capQn = number from which small factors are extracted fList = [[aNm1,sign]] # aN = independent square root i,j = len(pSlice)-1, 0 while i >= 0: if k%pSlice[i]==0: j += 1 # j = location of appended item fList.append([pSlice[i],1]) k /= pSlice[i] while k%pSlice[i]==0: fList[j][1] += 1 k /= pSlice[i] i -= 1 if k==1: fList = newRec(dBasis,topPrime,fList,modulus,recNo) capQnM2, rNm2, aNm2, sign = capQnM1, rNm1, aNm1, -sign capQnM1, rNm1, aNm1, qNm1 = capQn, rN, aN, qN newTime = int(time.clock()) if (newTime > oldTime): print "recNo,newTime: ", recNo[0],newTime," sec., ",len(dBasis)," records on dBasis" oldTime = newTime recNo[0] += 1 testVal = factSqrt(dBasis[0],modulus,recNo) del dBasis[0] del dBasis return (testVal) def prCheck(tNum,iter): if tNum < 4: return 1 i = 0 while i < iter: k = random.randint(2,tNum - 1) if pow(k,tNum - 1,tNum) != 1: return 0 i += 1 return 1 def diagnose(fList,modulus,recNo): #Checks for error on record i,accum,fLen = 1,1,len(fList) while i < fLen: if fList[i][1] > 32767: print "Integer overflow on record",recNo print "Record contents:",fList[:] print abort() accum *= pow(fList[i][0],fList[i][1]) i += 1 accum %= modulus test = (fList[0][1]*fList[0][0]*fList[0][0])%modulus # fList[0][1] = plus or minus sign if test != accum: print "Bad record:",recNo[0] return False else: return True def combine(gList1,dBasis,iSub,modulus,recNo): # Merge gList2 into gList1 to get product gList2 = dBasis[iSub] gList1[0][0] = (gList1[0][0]*gList2[0][0]) % modulus gList1[0][1] *= gList2[0][1] i1,i2,len1,len2 = 1,1,len(gList1),len(gList2) while i2 < len2: if i1 >= len1: gList1.append(gList2[i2]) i2 += 1 elif gList1[i1][0] < gList2[i2][0]: gList1.append(gList2[i2]) i2 += 1 elif gList1[i1][0] == gList2[i2][0]: if (recNo[0]==150377) & (i2==2): print "gList1 (before): ",gList1 print "dBasis[163] (before):",dBasis[163] gList1[i1][1] += gList2[i2][1] if (recNo[0]==150377) & (i2==2): print "gList1 (after): ",gList1 print "dBasis[163] (after): ",dBasis[163] i2 += 1 else: i1 += 1 if gList1[0][0] > gList1[1][0]: # The sane situation gList1.sort(reverse = True) else: # Insane, but we check for it anyway gList3 = sorted(gList[1:],reverse=True) del gList1[1:] gList1.extend(gList3) 