 SharedNumber Theory 17S Scratch Paper.sagewsOpen in CoCalc
Author: Paul Zeitz
Views : 1
# 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