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

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)

prime_pi(10^6)  

7*8

5667^79

euler_phi(9)

mod(6^16,17)

sum=0
for k in [1..1000]:
    sum += k^3

sum

n=10

    print k

n=100
for k in [0..n]:
    print k, 2^k

b=2^12
b

b.digits()[-1]

mylist = [1,2,3,4,5,6,77]

mylist

mylist[-1]

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]]

sum = 0
vals =[]
n = 2016
for k in range(n):
    sum += moebius(k)
    vals.append(sum)

list_plot(vals,plotjoined=true)

CRT_list?

CRT_list([1,2,3,4,5,6,7,8,9],[2,3,4,5,6,7,8,9,10])

michelle = Primes()

michelle

michelle.next(9)

michelle.unrank(9)

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

modList

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)
        

n = 100000
sum = 0
for k in [1..n]:
    if moebius(k) != 0:
        sum += 1
sum

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