# how many numbers less than n are prime?
prime_pi(2017)

n=1000
stairs=[prime_pi(k) for k in [1..n]]
list_plot(stairs, plotjoined=true)

next_prime(2017)

n=2000

n.is_square()

nsq=0
npr=0
for n in [1..1000]:
    if n.is_square():
        nsq += 1
    if n.is_prime():
        npr +=1

npr, nsq    

prime_pi(1000)

#the algorithm that never stops outputting primes.  
#we do it for just 10 steps, and it gives us 10 distinct primes
steps = 10
atoms =[2] #the first prime
for k in range(steps):
    #take product of all atoms, add 1; call it q
    pr=1
    for a in atoms:
        pr *= a
    q = pr +1
    if q.is_prime():
        # if q is a prime, append it to our list
        atoms.append(q)
    else: #otherwise, q factors, and append its least prime factor to our list
        atoms.append(min(q.prime_factors()))
#Now, print out the list of primes we found
atoms

#the song!
n=8675309
n.is_prime()

n=15 #size of grid
nticks =15 #number of ticks to mark on axes
d=ceil(n/nticks) #spacing between ticks

#create a matrix of 0,1 values depending on whether x is relatively prime to y
#recall that "!=" is the python symbol for not equals.  We could use "==" for logical equals but I like the colors better this way.
m = [[gcd(x,y) != 1 for x in [1..n]] for y in [1..n]]

#now we draw the matrix.  The 'origin' option places the lower values on the lower left (default is upper left)
#the 'ticks' and 'tick-formatter' stuff is designed to get the labeling right.  If you leave it off, you get Python
#defaults starting from zero
#Also note the "cmap" command which picks a nice color scheme.  Sage has a few default color schemes; look them up.
#You can load a huge collection of color schemes as well with the command 'import matplotlib.cm' and if you
#want to see the list of color maps, type 'matplotlib.cm.datad.keys()'
#finally, the option 'figsize' allows you to get a larger figure, duh.  good for large n. If you omit it you get a default size which is fine.
matrix_plot(m,origin='lower',ticks=[range(0,n,d),range(0,n,d)],tick_formatter=[range(1,n+1,d),range(1,n+1,d)],cmap='winter', figsize=8)

n=29 #size of grid


#create a matrix of 0,1 values depending on whether x is relatively prime to y
m = [[gcd(x,y)==1 for x in [1..n]]for y in [1..n]]
#here is a matrix plot with no tickmarks.  the way you do this is by declaring ticks that are outside the actual range
#(i.e., greater than n).  It is stupid, but it works.  It took me a long time to discover this trick.
matrix_plot(m,origin='lower',ticks=[range(2*n,3*n),range(2*n,3*n)])

n=10 #size of grid



m = [[gcd(x,y) for x in [1..n]]for y in [1..n]]
#here is a matrix plot with no tickmarks.  the way you do this is by declaring ticks that are outside the actual range
#(i.e., greater than n).  It is stupid, but it works.  It took me a long time to discover this trick.
matrix_plot(m,origin='lower',ticks=[range(2*n,3*n),range(2*n,3*n)], cmap='rainbow',colorbar=true)

################
# GCD calculations
gcd(30,200)

xgcd(30,200)

#implement Euclidean Algorithm and count number of steps
def EALength(a,b):
    #count number of steps in Euclidean Algorithm
    #if a<b then switch em
    if a < b:
        temp = a
        a = b
        b = temp
        
    r=b
    steps =0
    while r>0:
        temp = b
        r = a% b
        a = temp
        b = r
        steps += 1
    return steps

n=90 #size of grid


#create a matrix of values for EALength
m = [[EALength(x,y) for x in [1..n]]for y in [1..n]]
#here is a matrix plot with no tickmarks.  the way you do this is by declaring ticks that are outside the actual range
#(i.e., greater than n).  It is stupid, but it works.  It took me a long time to discover this trick.
matrix_plot(m,origin='lower',ticks=[range(2*n,3*n),range(2*n,3*n)],cmap='winter',colorbar=true)

sigma(5,1)

sigma(10,0)

Mod(5, 12)

a= mod(3,7)

a

a^100

3^100

1/a

b=mod(-1,17)

sqrt(b)

c=mod(3,12)

1/c

a=mod(2, 645)

a^644

factor(644)

a^14

a^7

a^28

euler_phi(645)

a^336

factor(336)

xgcd(27720,13)

13*6397

mod(83161,27720)

-27720*3

mod(-83160,27720*13)

mod(277200,13)

crt(1,2,3,4)

a=crt(0,-1,2,3)

is_prime(219)

factor(219)

crt(32,-5,30,7)

210*11+2

for k in [888..888+19]:
    print k, factor(k)

length = 35
for start in range(10000):
    prod=1
    for k in [start..start+length-1]:
        prod *= 1-is_prime(k)
    if prod==1:
        print start

crt(11, 19, 17,50)    

f(x)=x/(2-3*x+x^2)
f.taylor(x,0,18)

u = (-1+sqrt(5))/2
v = (-1-sqrt(5))/2
k=11
float((u^k-v^k)/sqrt(5))

solve(1-x-x^2-x^3==0,x)

mod(24^37,77)

mod(73^13,77)