Author: Paul Zeitz
Views : 5
#density estimator n=10000 #check from 1...n sum = 0 for k in [1..n]: if moebius(k)!=0: #test for k's membership sum += 1 float(sum/n)
0.6083
r=[0,1] n=10 for k in [1..n]: r.append(.5*(r[k]+2./(k+1))) print 2^k, r[k], r[k]*2^k,prime_pi(2^k)
2 1 2 1 4 1.00000000000000 4.00000000000000 2 8 0.833333333333333 6.66666666666667 4 16 0.666666666666667 10.6666666666667 6 32 0.533333333333333 17.0666666666667 11 64 0.433333333333333 27.7333333333333 18 128 0.359523809523810 46.0190476190476 31 256 0.304761904761905 78.0190476190476 54 512 0.263492063492063 134.907936507937 97 1024 0.231746031746032 237.307936507937 172
prime_pi(10^6)
78498
7*8
56
5667^79
327268193105661183797371787183788407076350243721782765059766772115276680877327201409266959755406511909524212919139023898815457959703176202017356583363286595211739210796061985288573900187007140725411329731580142846117411835707993405458437717490490789889678522564843345607290152996427337407660016203
euler_phi(9)
6
mod(6^16,17)
1
sum=0 for k in [1..1000]: sum += k^3
sum
250500250000
n=10 print k 7fb34dc5-9c66-4672-8395-c24770f5fdd3 n=100 for k in [0..n]: print k, 2^k
0 1 1 2 2 4 3 8 4 16 5 32 6 64 7 128 8 256 9 512 10 1024 11 2048 12 4096 13 8192 14 16384 15 32768 16 65536 17 131072 18 262144 19 524288 20 1048576 21 2097152 22 4194304 23 8388608 24 16777216 25 33554432 26 67108864 27 134217728 28 268435456 29 536870912 30 1073741824 31 2147483648 32 4294967296 33 8589934592 34 17179869184 35 34359738368 36 68719476736 37 137438953472 38 274877906944 39 549755813888 40 1099511627776 41 2199023255552 42 4398046511104 43 8796093022208 44 17592186044416 45 35184372088832 46 70368744177664 47 140737488355328 48 281474976710656 49 562949953421312 50 1125899906842624 51 2251799813685248 52 4503599627370496 53 9007199254740992 54 18014398509481984 55 36028797018963968 56 72057594037927936 57 144115188075855872 58 288230376151711744 59 576460752303423488 60 1152921504606846976 61 2305843009213693952 62 4611686018427387904 63 9223372036854775808 64 18446744073709551616 65 36893488147419103232 66 73786976294838206464 67 147573952589676412928 68 295147905179352825856 69 590295810358705651712 70 1180591620717411303424 71 2361183241434822606848 72 4722366482869645213696 73 9444732965739290427392 74 18889465931478580854784 75 37778931862957161709568 76 75557863725914323419136 77 151115727451828646838272 78 302231454903657293676544 79 604462909807314587353088 80 1208925819614629174706176 81 2417851639229258349412352 82 4835703278458516698824704 83 9671406556917033397649408 84 19342813113834066795298816 85 38685626227668133590597632 86 77371252455336267181195264 87 154742504910672534362390528 88 309485009821345068724781056 89 618970019642690137449562112 90 1237940039285380274899124224 91 2475880078570760549798248448 92 4951760157141521099596496896 93 9903520314283042199192993792 94 19807040628566084398385987584 95 39614081257132168796771975168 96 79228162514264337593543950336 97 158456325028528675187087900672 98 316912650057057350374175801344 99 633825300114114700748351602688 100 1267650600228229401496703205376
b=2^12 b
4096
b.digits()[-1]
4
mylist = [1,2,3,4,5,6,77]
mylist
[1, 2, 3, 4, 5, 6, 77]
mylist[-1]
77
a=2 n=10000 powers = [a^k for k in range(n)] lastDigits = [powers[i].digits()[-1] for i in range(n)] [float(lastDigits.count(j)/n) for j in [1..9]]
[0.301, 0.1761, 0.1249, 0.097, 0.0791, 0.067, 0.0579, 0.0512, 0.0458]
sum = 0 vals =[] n = 2016 for k in range(n): sum += moebius(k) vals.append(sum)
list_plot(vals,plotjoined=true) CRT_list? File: /projects/sage/sage-6.7/local/lib/python2.7/site-packages/sage/rings/arith.py Signature : CRT_list() Docstring : Given a list "v" of elements and a list of corresponding "moduli", find a single element that reduces to each element of "v" modulo the corresponding moduli. See also: * "crt()" EXAMPLES: sage: CRT_list([2,3,2], [3,5,7]) 23 sage: x = polygen(QQ) sage: c = CRT_list(, [x]); c 3 sage: c.parent() Univariate Polynomial Ring in x over Rational Field It also works if the moduli are not coprime: sage: CRT_list([32,2,2],[60,90,150]) 452 But with non coprime moduli there is not always a solution: sage: CRT_list([32,2,1],[60,90,150]) Traceback (most recent call last): ... ValueError: No solution to crt problem since gcd(180,150) does not divide 92-1 The arguments must be lists: sage: CRT_list([1,2,3],"not a list") Traceback (most recent call last): ... ValueError: Arguments to CRT_list should be lists sage: CRT_list("not a list",[2,3]) Traceback (most recent call last): ... ValueError: Arguments to CRT_list should be lists The list of moduli must have the same length as the list of elements: sage: CRT_list([1,2,3],[2,3,5]) 23 sage: CRT_list([1,2,3],[2,3]) Traceback (most recent call last): ... ValueError: Arguments to CRT_list should be lists of the same length sage: CRT_list([1,2,3],[2,3,5,7]) Traceback (most recent call last): ... ValueError: Arguments to CRT_list should be lists of the same length TESTS: sage: CRT([32r,2r,2r],[60r,90r,150r]) 452
CRT_list([1,2,3,4,5,6,7,8,9],[2,3,4,5,6,7,8,9,10])
2519
michelle = Primes()
michelle
Set of all prime numbers: 2, 3, 5, 7, ...
michelle.next(9)
11
michelle.unrank(9)
29
4b57a35c-456d-4a7c-a0ae-b1ebfddd20e1 l = 2015 valList = [-k for k in range(l)] modList = [michelle.unrank(k)^2 for k in range(l)] n =CRT_list(valList, modList) n [moebius(n+k) for k in range(l)] 0c6457ad-e1a5-439a-bd8d-db02ed8cc8cd int(log(n,10))
15091
modList
[4, 9, 25, 49, 121, 169, 289, 361, 529, 841]
n=400 relPrimeArray=[[gcd(x,y)==1 for x in [1..n]] for y in [1..n]] matrix_plot(relPrimeArray,origin='lower',cmap='winter') sum =0 #let's compute the number of green pixels for x in [1..n]: for y in [1..n]: if gcd(x,y)==1: sum += 1 print float(sum/n^2) 0.60846875
n = 100000 sum = 0 for k in [1..n]: if moebius(k) != 0: sum += 1 sum
60794
######################################################################################### # This is an INTERACTIVE simulation which allows you to change the price per game # and the total number of games. It also prints out the maximum payout, and # when this happened, the final amount of money on hand, the minimum amount, and the maximum, # and the average payout (which, theoretically, is INFINITY). # # Note the use of the @interactive command, followed by a function definition (the # name of the function can be anything, but we choose the generic "_" since we will never use the # function's name. You can look up the use of Sage interactive online, or just study this example. # It's pretty flexible, and easier to use than Mathematica's Manipulate command, although less slick. ############################################################################################## # We include the already-defined first_head function, so that this is a stand-alone cell. def first_head(): coin=0 # 0 for tails, 1 for heads. Start by setting the coin to tails count=0 # Keep track of the number of flips before first head while coin == 0: # the two lines in this "while loop" will be executed as long as the coin value is 0 # but once it equals 1, the loop ends and we move to the final line of the function def. coin = randint(0,1) # set coin equal to 0 or 1 randomly count += 1 # increment the count return count # this is the count value when the coin first equals 1 (and thus, the while loop ends). ############################################################################################################ # interactive part starts here ############################################################################################################ @interact def _(price= [5,10,15,20,100], nGames = [100,1000,10000,100000]): # note that the selection boxes for price and nGames are defined in the function defn above. accumulation = 0 #player's profit sequence = [] #keep track of the first_head() outputs daily = [] #accumlation on each day; to be built for _ in range(nGames): fh = first_head() sequence.append(fh) accumulation += 2^fh-price daily.append(accumulation) # Compute average payout totalPayout = 0 for i in range(len(sequence)): totalPayout += 2^sequence[i] avg = float(totalPayout/nGames) # Compute longest run of tails, max min, etc. longestRun=max(sequence) longestRunIndex=sequence.index(longestRun) print('Luckiest: game # %s, first head at position %s, winning \$%s.'%(longestRunIndex, longestRun, 2^longestRun)) print('Lowest: %s'%min(daily)) print('Highest: %s'%max(daily)) print('Final: %s'%daily[nGames-1]) print('Average payout: %s.'%avg) # Plot the money on hand. l=list_plot(daily,plotjoined = True,gridlines=True) show(l) #in interactive mode, the plot doesn't display with just the list_plot() command. #instead, we use the show() function.
Interact: please open in CoCalc