CoCalc Public Filesdensity.sagews
Author: Paul Zeitz
Views : 33
#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.

EXAMPLES:

sage: CRT_list([2,3,2], [3,5,7])
23
sage: x = polygen(QQ)
sage: c = CRT_list([3], [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)]
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.

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