#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()

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

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