from __future__ import division from time import clock from gmpy import * from pysymbolicext import * import psyco psyco.full() __all__=["carmichael","carmichael1","carmichael2","carmichael3","iscarmichael","clambda"] #Carmichael numbers are odd composites C with at least 3 prime #factors p(n) such that for each p(n) of C, p(n)-1 divides C-1 def carmichael1(a,b,limit=1000): """ carmichael(a,b) For a,b (odd primes) returns array of c such that a*b*c is a Carmichael number """ if a==b or a==2 or b==2: return [] if not is_prime(a) or not is_prime(b): return [] ab=a*b l=lcm(a-1,b-1) m1=ab%l m2=invert(m1,l) abl=ab*l c=m2 while c<=max(a,b): c+=l xm1=ab*c-1 solutions=[] while c<=limit: if (xm1)%(c-1)==0 and is_prime(c): solutions.append(c) c+=l xm1+=abl return solutions def carmichael2(a,b,limit=1000): """ carmichael(a,b) For a,b (odd primes) returns array of c such that a*b*c is a Carmichael number """ if a==b or a==2 or b==2: return [] if not is_prime(a) or not is_prime(b): return [] am1=a-1 bm1=b-1 am1bm1=am1*bm1 g,s,t=gcdext(am1bm1,-(am1+b)) c=((t*(am1+bm1)//g))%bm1+1 while c<=max(a,b): c+=bm1 xm1=(a)*(b)*(c)-1 inc=a*b*bm1 solutions=[] while c<=limit: if (xm1)%(c-1)==0 and is_prime(c): solutions.append(c) c+=bm1 xm1+=inc return solutions #A different approach,but seems slower on many inputs def carmichael3(*args): """ carmichael(*args) (last arg is search limit) For p1...pn (odd distinct primes) returns array of c such that p1*p2...pn*c is a Carmichael number. """ limit=args[-1] args=args[:-1] if 2 in args: return [] for arg in args: if not is_prime(arg): return [] product=reduce(lambda a,b : a*b, args) l=reduce(lambda a,b : lcm(a,b), (x-1 for x in args)) m1=product%l m2=invert(m1,l) productl=product*l c=m2 while c<=max(args): c+=l xm1=product*c-1 solutions=[] while c<=limit: if (xm1)%(c-1)==0 and is_prime(c): solutions.append(c) c+=l xm1+=productl return solutions def carmichael(*args): """ carmichael(*args) (last arg is search limit) For p1...pn (odd distinct primes) returns array of c such that p1*p2...pn*c is a Carmichael number. Redirection function - calls carmichael 1 or 2 depending on number of arguments. """ if len(args)==3: return carmichael1(*args) else: return carmichael3(*args) def iscarmichael(*args): """ iscarmichael(*args). args are prime factors. Naively tests if the product of the given factors is a Carmichael number """ for arg in args: if not is_prime(arg): return False xm1=reduce(lambda a,b : a*b, args)-1 for arg in args: if xm1%(arg-1)!=0: return False return True def cgenerate(limit=1000000,a=6,b=12,c=18): """ cgenerate(a,b,c,limit) Return a list of all Carmichael numbers of the form (ak+1)*(bk+1)*(ck+1) up to limit """ print "Using search limit",limit result=[] k=1 while True: a1=a*k+1 b1=b*k+1 c1=c*k+1 n=a1*b1*c1 if n>limit: break if is_prime(a1) and is_prime(b1) and is_prime(c1): result.append((a1,b1,c1,n)) #result.append(n) k+=1 return result def clambda(n): """ clambda(n) Returns Carmichael's lambda function for positive integer n. Relies on factoring n """ smallvalues=[1,1,2,2,4,2,6,2,6,4,10,2,12,6,4,4,16,6,18,4,6,10,22,2,20,12,18,\ 6,28,4,30,8,10,16,12,6,36,18,12,4,40,6,42,10,12,22,46,4,42,20,16,12,52,18,\ 20,6,18,28,58,4,60,30,6,16,12,10,66,16,22,12,70,6,72,36,20,18,30,12,78,4,54,\ 40,82,6,16,42,28,10,88,12,12,22,30,46,36,8,96,42,30,20] if n<=100: return smallvalues[n-1] factors=factor(n) l1=[] for p,e in factors: if p==2 and e>2: l1.append(2**(e-2)) else: l1.append((p-1)*p**(e-1)) return reduce(lambda a,b : lcm(a,b), l1) numbers=cgenerate(100000000000000000000,2,4,14) for x in numbers: print x[3], if iscarmichael(x[0],x[1],x[2]): print "is a Carmichael number" else: print #----------------------------------TEST CODE----------------------------------# if __name__=="__main__": def naive(a,b,limit): """ naive(a,b,limit). Naive implementation for benchmarking comparison. """ if a==b or a==2 or b==2: return [] if not is_prime(a) or not is_prime(b): return [] c=b+2 xm1=a*b*c-1 am1=a-1 bm1=b-1 inc=2*a*b solutions=[] while c<limit: if is_prime(c): if xm1%am1==0 and xm1%bm1==0 and xm1%(c-1)==0: solutions.append(c) c+=2 xm1+=inc return solutions #--------------------------------------------------------------------------- #a=727 #b=1453 #limit=500000 #trials=10 #start=clock() #for _ in xrange(trials): # result=naive(a,b,limit) #end=clock() #print "Naive in",end-start,"seconds" #print result #start=clock() #for _ in xrange(trials): # result=carmichael1(a,b,limit) #end=clock() #print "Carmichael1 in",end-start,"seconds" #print result #start=clock() #for _ in xrange(trials): # result=carmichael2(a,b,limit) #end=clock() #print "Carmichael2 in",end-start,"seconds" #print result #start=clock() #for _ in xrange(trials): # result=carmichael3(a,b,limit) #end=clock() #print "Carmichael3 in",end-start,"seconds" #print result #--------------------------------------------------------------------------- #899331780481 = 239 * 409 * 1021 * 9011 #a=239 #b=409 #c=1021 #limit=10000 #start=clock() #for _ in xrange(trials): # result=carmichael(a,b,c,limit) #end=clock() #print "Carmichael (redirection) in",end-start,"seconds" #print result #Naive in 2.19843932679 seconds #[2179, 3631, 10891, 70423, 352111] #Carmichael in 0.0035508004509 seconds #[mpz(2179), mpz(3631), mpz(10891), mpz(70423), mpz(352111)] #Carmichael2 in 0.00424956244439 seconds #[mpz(2179), mpz(3631), mpz(10891), mpz(70423), mpz(352111)] #Carmichael3 in 0.00353201314692 seconds #[mpz(2179), mpz(3631), mpz(10891), mpz(70423), mpz(352111)] #--------------------------------------------------------------------------- #from carmichaels import cdic #k=cdic.keys() #k.sort() #testlist=[] #for c in k: # if len(cdic[c])==3: # testlist.append((cdic[c][0],cdic[c][1])) #start=clock() #limit=1000000 #result=[] #for test in testlist: # a=test[0] # b=test[1] # result.append((a,b,carmichael1(a,b,limit))) #end=clock() #print result #print "Carmichael1: All triples < 10^12 ( 1000 numbers ) in",end-start,"seconds" #result=[] #for test in testlist: #a=test[0] #b=test[1] #result.append((a,b,naive(a,b,limit))) #end=clock() #print result #print "Naive: All triples < 10^12 ( 1000 numbers ) in",end-start,"seconds" #Carmichael1: All triples < 10^12 in 0.764858592363 seconds #Naive: All triples < 10^12 in 373.157365461 seconds #--------------------------------------------------------------------------- #Print maximum number of factors found in cdic #from carmichaels import cdic #print reduce(lambda a,b: max(a,b), (len(cdic[c]) for c in cdic))