CoCalc Public Files2020-04-29-MoMath share.sagews
Author: Paul Zeitz
Views : 19
#Now identify the primes #ulam spiral # find n, given coords, first normalized # so that (0,0)-->1, (1,0)-->2, (1,1)-->3, etc. def pointToNum(x,y): if (x==0) & (y==0): return 1 #first determine quadrant based on y= x, y=-x if (y>-x) & (y<=x): q=3 if (y>x) & (y>=-x): q=2 if (y<x) & (y<=-x): q=0 if (y<-x) & (y>=x): q=1 #next, determine which concentric square level we are at a = max(abs(x),abs(y)) #list the coordinates of each corner point corners = [[a,-a],[-a,-a],[-a,a],[a,a]] #find the value of each quadrants corner point. #start with southeast corner (#0) se se=(2*a+1)^2 cornerVals=[se, se-2*a, se-4*a, se-6*a] #now find taxicab distance from given point to corner point dist = abs(x-corners[q][0]) + abs(y-corners[q][1]) #finally, subtract this from corner value return cornerVals[q]-dist # Now we make a matrix of size 2a+1 where the origin (0,0) is now (a,a). # Hence the point (x,y) in matrix coords (with (0,0) at lower left) # will be translated to (x-a, y-a) for pointToNum computation a= 30 # Make a value array ulamArrayVals = [[pointToNum(x-a, y-a) for x in range(2*a+1)] for y in range(2*a+1)] #Also, make an array with values measuring "how prime" using sum of divisor; # if it is prime the value will be zero ulamArrayDivisors = [[sigma(pointToNum(x-a, y-a),0)-2 for x in range(2*a+1)] for y in range(2*a+1)] #create array of abundant numbers ulamArrayAbundant = [[float(sigma(pointToNum(x-a, y-a))/pointToNum(x-a, y-a)) >2 for x in range(2*a+1)] for y in range(2*a+1)] #abundancy quotient ulamArrayAbundantQ = [[float(sigma(pointToNum(x-a, y-a))/pointToNum(x-a, y-a)) for x in range(2*a+1)] for y in range(2*a+1)] def plotFromArray(array,colormap): return matrix_plot(array,origin='lower',frame=False,cmap=colormap) def ulamPrimeMake(size,name,colormap): a = size #Now identify the primes ulamArrayPrimes = [[pointToNum(x-a, y-a) not in Primes() for x in range(2*a+1)] for y in range(2*a+1)] plot = plotFromArray(ulamArrayPrimes,colormap) plotName=name+'.pdf' plot.save(filename=plotName) def ulamAbundantMake(size,name,colormap): a = size #Now identify the primes ulamArray = [[float(sigma(pointToNum(x-a, y-a))/pointToNum(x-a, y-a))>2 for x in range(2*a+1)] for y in range(2*a+1)] plot = plotFromArray(ulamArray,colormap) plotName=name+'.pdf' plot.save(filename=plotName) def ulamPrimeMake(size,name,colormap): a = size #Now identify the primes ulamArray = [[pointToNum(x-a, y-a) not in Primes() for x in range(2*a+1)] for y in range(2*a+1)] plot = plotFromArray(ulamArray,colormap) plotName=name+'.pdf' plot.save(filename=plotName) def ulamNotPrimeMake(size,name,colormap): a = size #Now identify the primes ulamArray = [[pointToNum(x-a, y-a) in Primes() for x in range(2*a+1)] for y in range(2*a+1)] plot = plotFromArray(ulamArray,colormap) plotName=name+'.pdf' plot.save(filename=plotName) def relPrimeMake(size,name,colormap): relPrimeArray=[[gcd(x,y)!=1 for x in [1..size]] for y in [1..size]] plot=matrix_plot(relPrimeArray,origin='lower',frame=False,cmap=colormap) plotName=name+'.pdf' plot.save(filename=plotName,frame=False) def notRelPrimeMake(size,name,colormap): relPrimeArray=[[gcd(x,y)==1 for x in [1..size]] for y in [1..size]] plot=matrix_plot(relPrimeArray,origin='lower',frame=False,cmap=colormap) plotName=name+'.pdf' plot.save(filename=plotName,frame=False)
ulamPrimeMake(100,'prime100','winter')
def ulamRandomMake(size,name,colormap): a = size #Now identify the primes ulamArrayPrimes = [[ (random()>1/float(log(2+pointToNum(x-a, y-a)))) for x in range(2*a+1)] for y in range(2*a+1)] plot = plotFromArray(ulamArrayPrimes,colormap) plotName=name+'.pdf' plot.save(filename=plotName)
for k in [1..40]: print k, k^2-k+41, (k^2-k+41).is_prime()
1 41 True 2 43 True 3 47 True 4 53 True 5 61 True 6 71 True 7 83 True 8 97 True 9 113 True 10 131 True 11 151 True 12 173 True 13 197 True 14 223 True 15 251 True 16 281 True 17 313 True 18 347 True 19 383 True 20 421 True 21 461 True 22 503 True 23 547 True 24 593 True 25 641 True 26 691 True 27 743 True 28 797 True 29 853 True 30 911 True 31 971 True 32 1033 True 33 1097 True 34 1163 True 35 1231 True 36 1301 True 37 1373 True 38 1447 True 39 1523 True 40 1601 True
n=1000000 sum([(k^2-k+41).is_prime() for k in [1..n]])
a=2 n=10000 powers = [a^k for k in range(n)] lastDigits = [powers[i].digits()[-1] for i in range(n)] [lastDigits.count(j) for j in [1..9]]
[3010, 1761, 1249, 970, 791, 670, 579, 512, 458]
mylist=[7,7,7] mlen=len(mylist) matchlist=[] for n in [10..999]: testlist=powers[n].digits() tlen=len(testlist) result=false for i in [0..tlen-mlen-1]: if testlist[i:i+mlen]==mylist: result=true matchlist.append(n) break matchlist
[24, 40, 75, 152, 166, 179, 181, 191, 194, 199, 214, 230, 235, 260, 296, 304, 311, 317, 323, 326, 332, 342, 345, 363, 370, 374, 390, 417, 424, 426, 443, 455, 468, 471, 474, 475, 483, 489, 490, 505, 512, 523, 524, 536, 540, 559, 567, 581, 584, 585, 588, 593, 608, 609, 611, 618, 620, 626, 632, 643, 645, 649, 650, 651, 653, 660, 663, 669, 670, 678, 689, 692, 693, 694, 695, 696, 702, 704, 706, 710, 715, 717, 718, 721, 722, 743, 746, 756, 759, 763, 770, 775, 777, 781, 782, 787, 791, 799, 801, 807, 816, 818, 819, 824, 828, 830, 837, 851, 863, 867, 880, 883, 902, 906, 909, 922, 925, 927, 930, 931, 944, 954, 958, 967, 969, 972, 973, 974, 975, 978, 981, 984, 985, 987, 991, 992, 997]