 CoCalc Public FilesNumber Theory 17S Scratch Paper.sagews
Author: Paul Zeitz
Views : 31
# 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 = #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