2020-04-29-MoMath share.sagews

#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]