Number Theory 17S Scratch Paper.sagews

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

306

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

next_prime(2017)

2027

n=2000

n.is_square()

False

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

npr, nsq    

(168, 31)

prime_pi(1000)

168

#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

[2, 3, 7, 43, 13, 53, 5, 6221671, 38709183810571, 139, 2801]

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

True

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)

10

xgcd(30,200)

(10, 7, -1)

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

6

sigma(10,0)

4

Mod(5, 12)

5

a= mod(3,7)

a

3

a^100

4

3^100

515377520732011331036461129765621272702107522001

1/a

5

b=mod(-1,17)

sqrt(b)

4

c=mod(3,12)

1/c

Error in lines 1-1 Traceback (most recent call last): File "/projects/sage/sage-7.5/local/lib/python2.7/site-packages/smc_sagews/sage_server.py", line 982, in execute exec compile(block+'\n', '', 'single') in namespace, locals File "", line 1, in <module> File "sage/rings/integer.pyx", line 1858, in sage.rings.integer.Integer.__div__ (/projects/sage/sage-7.5/src/build/cythonized/sage/rings/integer.c:13046) return coercion_model.bin_op(left, right, operator.div) File "sage/structure/coerce.pyx", line 1055, in sage.structure.coerce.CoercionModel_cache_maps.bin_op (/projects/sage/sage-7.5/src/build/cythonized/sage/structure/coerce.c:9336) return PyObject_CallObject(op, xy) File "sage/structure/element.pyx", line 1638, in sage.structure.element.Element.__div__ (/projects/sage/sage-7.5/src/build/cythonized/sage/structure/element.c:12978) return (<Element>left)._div_(right) File "sage/rings/finite_rings/integer_mod.pyx", line 2430, in sage.rings.finite_rings.integer_mod.IntegerMod_int._div_ (/projects/sage/sage-7.5/src/build/cythonized/sage/rings/finite_rings/integer_mod.c:28670) raise ZeroDivisionError("Inverse does not exist.") ZeroDivisionError: Inverse does not exist.

a=mod(2, 645)

a^644

1

factor(644)

2^2 * 7 * 23

a^14

259

a^7

128

a^28

1

euler_phi(645)

336

a^336

1

factor(336)

2^4 * 3 * 7

xgcd(27720,13)

(1, -3, 6397)

13*6397

83161

mod(83161,27720)

1

-27720*3

-83160

mod(-83160,27720*13)

277200

mod(277200,13)

1

crt(1,2,3,4)

10

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

is_prime(219)

False

factor(219)

3 * 73

crt(32,-5,30,7)

2

210*11+2

2312

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

888 2^3 * 3 * 37 889 7 * 127 890 2 * 5 * 89 891 3^4 * 11 892 2^2 * 223 893 19 * 47 894 2 * 3 * 149 895 5 * 179 896 2^7 * 7 897 3 * 13 * 23 898 2 * 449 899 29 * 31 900 2^2 * 3^2 * 5^2 901 17 * 53 902 2 * 11 * 41 903 3 * 7 * 43 904 2^3 * 113 905 5 * 181 906 2 * 3 * 151 907 907

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

9552

crt(11, 19, 17,50)    

419

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

x |--> 262143/262144*x^18 + 131071/131072*x^17 + 65535/65536*x^16 + 32767/32768*x^15 + 16383/16384*x^14 + 8191/8192*x^13 + 4095/4096*x^12 + 2047/2048*x^11 + 1023/1024*x^10 + 511/512*x^9 + 255/256*x^8 + 127/128*x^7 + 63/64*x^6 + 31/32*x^5 + 15/16*x^4 + 7/8*x^3 + 3/4*x^2 + 1/2*x

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

89.00000000000004

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

[x == -1/2*(1/9*sqrt(11)*sqrt(3) + 17/27)^(1/3)*(I*sqrt(3) + 1) + 1/9*(-I*sqrt(3) + 1)/(1/9*sqrt(11)*sqrt(3) + 17/27)^(1/3) - 1/3, x == -1/2*(1/9*sqrt(11)*sqrt(3) + 17/27)^(1/3)*(-I*sqrt(3) + 1) - 1/9*(-I*sqrt(3) - 1)/(1/9*sqrt(11)*sqrt(3) + 17/27)^(1/3) - 1/3, x == (1/9*sqrt(11)*sqrt(3) + 17/27)^(1/3) - 2/9/(1/9*sqrt(11)*sqrt(3) + 17/27)^(1/3) - 1/3]

mod(24^37,77)

73

mod(73^13,77)

24