

<h2>Fibonacci numbers</h2>
<p>Let us experiment a little bit with the <strong>Fibonacci</strong> <strong>numbers</strong>. Recall that they are defined by the inital condition f_1=f_2=1 and the recursion f_{n+2} = f_{n+1} + f_n</p>

{{{id=26|
F = [fibonacci(n) for n in range(1,15)]; F
///
[1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377]
}}}

{{{id=28|
[m for m in F if is_even(m)]
///
[2, 8, 34, 144]
}}}

{{{id=29|
sum(fibonacci(i) for i in range(1,5))
///
7
}}}

<p>Defining our first function:</p>

{{{id=30|
def sum_fib(n):
    return sum(fibonacci(i) for i in range(1,n+1))
///
}}}

{{{id=31|
[sum_fib(n) for n in range(1,13)]
///
[1, 2, 4, 7, 12, 20, 33, 54, 88, 143, 232, 376]
}}}

<p>Now we might conjecture that the sum of the first n Fibonacci numbers is equal to f_{n+2}-1. Let us test this conjecture!</p>

{{{id=32|
for n in range(1,100):
    assert( sum_fib(n) == fibonacci(n+2)-1)
print 'Done!'
///
Done!
}}}

<p>What would have happened if we had tried a wrong conjecture?</p>

{{{id=36|
for n in range(1,100):
    assert( sum_fib(n) == fibonacci(n+2))
print 'Done!'
///
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "_sage_input_59.py", line 10, in <module>
    exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("Zm9yIG4gaW4gcmFuZ2UoMSwxMDApOgogICAgYXNzZXJ0KCBzdW1fZmliKG4pID09IGZpYm9uYWNjaShuKzIpKQpwcmludCAnRG9uZSEn"),globals())+"\\n"); execfile(os.path.abspath("___code___.py"))' + '\n', '', 'single')
  File "", line 1, in <module>
    
  File "/private/tmp/tmpCmyTNL/___code___.py", line 4, in <module>
    assert( sum_fib(n) == fibonacci(n+_sage_const_2 ))
AssertionError
}}}

<h2>Partially ordered sets</h2>
<p>The set of all divisors of a number n form a <strong>poset</strong> by defining a&lt;=b if a divides b.</p>

{{{id=39|
n = 30
n.divisors()
///
[1, 2, 3, 5, 6, 10, 15, 30]
}}}

{{{id=40|
def order_relation(a,b):
    return b%a == 0
///
}}}

{{{id=60|
order_relation(15,30)
///
True
}}}

{{{id=61|
order_relation(6,15)
///
False
}}}

{{{id=41|
n = 30
P = Poset((n.divisors(), order_relation), cover_relations = False)
P.plot()
///
<html><font color='black'><img src='cell://sage0.png'></font></html>
}}}

{{{id=5|
B = Posets.BooleanLattice(3)
B.plot()
///
<html><font color='black'><img src='cell://sage0.png'></font></html>
}}}

{{{id=12|
n = 60
P = Poset((n.divisors(), order_relation), cover_relations = False)
P.plot()
///
<html><font color='black'><img src='cell://sage0.png'></font></html>
}}}

<p>We can ask about some properties of the poset:</p>

{{{id=13|
P.is_graded()
///
True
}}}

{{{id=18|
P.is_meet_semilattice()
///
True
}}}

{{{id=42|
P.is_join_semilattice()
///
True
}}}

{{{id=14|
P.minimal_elements()
///
[1]
}}}

{{{id=19|
P.maximal_elements()
///
[60]
}}}

{{{id=20|
P.mobius_function(P(1),P(30))
///
-1
}}}

{{{id=59|
P.mobius_function(P(1),P(4))
///
0
}}}

{{{id=43|
sum([P.mobius_function(P(1),v) for v in P])
///
0
}}}

<h2>Mersenne numbers</h2>
<p>Let us test some results on <strong>Mersenne numbers:</strong></p>

{{{id=2|
def mersenne(p):
    return 2^p-1
///
}}}

<p>Is p always prime when the Mersenne number is a prime?</p>

{{{id=46|
[is_prime(p) for p in range(1000) if is_prime(mersenne(p))]
///
[True, True, True, True, True, True, True, True, True, True, True, True, True, True]
}}}

<p>Is the converse also true?</p>

{{{id=48|
all( is_prime(mersenne(p)) for p in range(1000) if is_prime(p) )
///
False
}}}

<p>Finding the smallest counterexample:</p>

{{{id=50|
exists( (p for p in range(1000) if is_prime(p)), lambda p: not is_prime(mersenne(p)) )
///
(True, 11)
}}}

<h2>Game of poker</h2>
<p>&nbsp;We have 52=4x13 cards with suits Club, Spade, Heart and Diamond and values 2,3,...,Ace</p>

{{{id=64|
Suits = Set( ['Club', 'Spade', 'Heart', 'Diamond'] )
Values = Set( [2, 3, 4, 5, 6, 7, 8, 9, 10, 'Jack', 'Queen', 'King', 'Ace'] )
///
}}}

<p>The cards are the Cartesian product:</p>
<p>Cards = Suits x Values = {(s,v) | s in Suits, v in Values}</p>

{{{id=65|
Cards = CartesianProduct(Suits, Values)
///
}}}

{{{id=67|
Suits.cardinality()
///
4
}}}

{{{id=68|
Values.cardinality()
///
13
}}}

{{{id=69|
Cards.cardinality()
///
52
}}}

<p>We can deal a random element:</p>

{{{id=71|
Cards.random_element()
///
['Heart', 'Jack']
}}}

<p>Now we need a brief technical digression:</p>

{{{id=73|
type(Cards.random_element())
///
<type 'list'>
}}}

<p>Python lists are not immutable, hence they cannot be used to put elements in a set. This can be fixed by transforming the lists into tuples:</p>

{{{id=74|
Cards = CartesianProduct(Suits, Values).map(tuple)
Cards.random_element()
///
('Spade', 7)
}}}

{{{id=76|
Set( [ Cards.random_element(), Cards.random_element() ] )
///
{('Club', 6), ('Diamond', 5)}
}}}

<p>A poker hand consists of 5 cards. All cards are distinguishable and the order of the cards does not matter for a hand. A hand is hence a 5-subset&nbsp; of the set of all cards.</p>

{{{id=77|
Hands = Subsets(Cards,5)
Hands.random_element()
///
{('Club', 9), ('Heart', 'Jack'), ('Diamond', 'King'), ('Diamond', 4), ('Heart', 'Ace')}
}}}

<p>Total number of hands:</p>

{{{id=80|
binomial(52,5)
///
2598960
}}}

{{{id=81|
Hands.cardinality()
///
2598960
}}}

<p>Let us now calculate the probability for a <strong>flush,&nbsp;</strong> that is, all cards in a given hand have the same suit:</p>

{{{id=83|
Colors = CartesianProduct(Suits, Subsets(Values,5))
Colors.cardinality()
///
5148
}}}

{{{id=84|
100.0 * Colors.cardinality() / Hands.cardinality()
///
0.198079231692677
}}}

<p>We can test our reasoning by doing a little simulation: we produce 10000 random hands and count how many of them have the same color.</p>

{{{id=87|
def is_same_suit(hand):
        return len(set(suit for (suit,value) in hand)) == 1
///
}}}

{{{id=88|
n = 10000
nsuits = 0
for i in range(n):
    hand = Hands.random_element()
    if is_same_suit(hand):
        nsuits += 1
print n, nsuits
print 100.0 * nsuits / n
///
10000 20
0.200000000000000
}}}

{{{id=95|

///
}}}