The number of $3$-sets chosen from a $10$-set, ${10 \choose 3}$.
{{{id=1| binomial(10, 3) /// }}}The coefficients of an expansion of $(a+b)^n$.
{{{id=3| var('a, b') expr = (a+b)^10 expr.expand() /// }}}All of the coefficients.
{{{id=5| bc = binomial_coefficients(10) bc /// }}}A Python dictionary, indexed by powers of the two variables in the expansion.
{{{id=7| bc[(3, 7)] /// }}}Tote up all of these binomial coefficients, to get $2^{10}$. (Size of the power set, or the result of setting $a=1$ and $b=1$).
{{{id=9| sum(bc.values()) /// }}}Actual subsets of size 3 from a $10$-set; one way to understand a binomial coefficient.
{{{id=11| S = list("U-VIctoria") S /// }}} {{{id=12| sub = Subsets(S, 3) sub /// }}}sub is a “generator”. We can list the possibilities.
{{{id=14| sub.list() /// }}}We can iterate over sub.
{{{id=16| for three in sub: print three /// }}}$C_n = \frac{1}{n+1}{2n \choose n} = \frac{(2n)!}{(n+1)!\,n!}$
{{{id=20| (1/9)*binomial(16, 8) /// }}}In honor of Eric Temple Bell.
{{{id=22| bell_number(6) /// }}}Number of partitions of a set into disjoint non-empty sets.
{{{id=24| S = list('campus') part = SetPartitions(S) part /// }}} {{{id=25| part.list() /// }}} {{{id=26| part[34] /// }}} {{{id=27| len(part) /// }}} {{{id=28| part.cardinality() /// }}}Stirling numbers come in two flavors, “first” and “second”, or “cycle” and “subset”. We'll demonstrate the first.
{{{id=30| stirling_number1(6, 3) /// }}}The number of permutations on $n$ symbols (in cycle notation) having exactly $k$ cycles, $\left\{n\atop k\right\}$.
{{{id=32| perm = Permutations(6) a = perm[134] a /// }}} {{{id=33| a.cycle_string() /// }}}Now we get the trivial cycles. List length is what we want.
{{{id=35| a.cycle_tuples() /// }}}Collect all permutations with 3 cycles.
{{{id=37| three = [p for p in perm if len(p.cycle_tuples())==3] three /// }}}How many?
{{{id=39| len(three) /// }}}Sage was born of necessity to do number theory.
{{{id=41| p = next_prime(10^25) q = next_prime(10^25+5*10^24+10^12) m = p*q print p, q print m /// }}}Factor 50-digit number ($\approx$166 bits).
{{{id=43| m.factor() /// }}}Euler $\phi$ function. (“totient” function.)
{{{id=45| euler_phi(100) /// }}}Integers less than $100$ and relatively prime to $100$. (Note the srange function to generate Sage integers.)
{{{id=47| relp = [x for x in srange(100) if gcd(x, 100) == 1] relp /// }}} {{{id=48| len(relp) /// }}}Fact: $\displaystyle\sum_{d\mid n}\,\phi(d) = n$
Proof: Group the fractions, $\displaystyle\frac{i}{n}$, $0\leq i\leq n-1$, by denominators once written in reduced terms.
{{{id=50| n = 100 sum([euler_phi(d) for d in divisors(n)]) == n /// }}}You need to declare symbolic variables (except x comes pre-defined). That done, summations simplify as expected.
{{{id=52| var('i, n') expr = sum(i^2, i, 0, n) expr /// }}}We'll recognize this result if we factor.
{{{id=54| expr.factor() /// }}}Arbitrarily complicated polynomials as summands can be simplified.
{{{id=56| var('i, n') expr = sum(2*i^5 - 6*i^4 +7*i^2 - 8, i, 0, n) expr /// }}}We can convert this symbolic expression to a callable function.
{{{id=58| var('t') g(t) = expr.subs(n=t) g /// }}}And call it --- thus making $n$ concrete.
{{{id=60| g(10) /// }}}Straightforward to plot a discrete function, we will plot using an expression.
{{{id=62| var('i, n') expr = sum(2*i^3 - 12*i^2, i, 0, n) points = [(k, expr.subs(n=k)) for k in range(10)] list_plot(points) /// }}} {{{id=63| list_plot(points, size=200, color='purple') /// }}} {{{id=64| list_plot(points, color='red', plotjoined=True) /// }}}Numbers of certain objects can sometimes be counted by recurrence relations. We would like closed-form expressions for terms of sequences defined this way.
Perrin Sequence:
$p(0) = 3$; $p(1) = 0$; $p(2)=2$
$p(n) = p(n-2) + p(n-3)$
Looks like the Fibonacci sequence, but “skips back” two terms, not one.
Compute by hand: $3, 0, 2, 3, 2, 5, 5, 7, 10, 12, 17, \dots$
This is in Sloane's Online Encyclopedia of Integer Sequences as sequence number A001608.
A brute-force approach with a Python function. Impractical above about $n=50$.
{{{id=66| def perrin(n): if n == 0: return 3 elif n == 1: return 0 elif n == 2: return 2 else: return perrin(n-2) + perrin(n-3) /// }}} {{{id=67| perrin(10) /// }}} {{{id=68| perrin(20) /// }}} {{{id=69| perrin(23).factor() /// }}}Fact: If $q$ is prime, then $q$ divides $p(q)$.
(First composite number that behaves this way is $521^2$.)
Generating function: $f(x)=\displaystyle\sum_{i=0}^\infty\,p_i\,x^i$
Theory gives easy computation for Perrin sequence, denominator comes from recurrence relation, numerator is simple polynomial multiplication. $$f(x) = \frac{3-x^2}{1-x^2-x^3}$$
We can expand $f$ as a Taylor series.
{{{id=71| var('x') f=(3-x^2)/(1-x^2-x^3) f /// }}} {{{id=72| series = f.taylor(x, 0, 20) series /// }}} {{{id=73| [series.coeff(x, i) for i in range(20)] /// }}}Partial Fractions can simplify a rational generating function.
$a(0)=7$; $a(1)=41$; $a(2)=204$
$a(n) = 7a(n-1) - 12a(n-2) + 10a(n-3)$
Generating function --- the rational function:
{{{id=75| h = (x^2 - 8*x + 7)/(1 - 7*x + 12*x^2 -10*x^3) h /// }}}Check $h(3)$:
{{{id=77| 7*204 -12*41 + 10*7 /// }}} {{{id=78| h.taylor(x, 0, 3) /// }}}Create partial fraction decomposition and examine the pieces:
{{{id=80| h.partial_fraction() /// }}} {{{id=81| denom1 = 1/(2*x^2 - 2*x + 1) denom1.taylor(x, 0, 30) /// }}} {{{id=82| denom2 = 1/(5*x - 1) denom2.taylor(x, 0, 8) /// }}}SymPy is a pure Python package included in Sage, but not integrated with Sage.
http://docs.sympy.org/dev/modules/solvers/solvers.html#recurrence-equtions (sic)
Import pieces of the SymPy library.
{{{id=84| from sympy import Function, rsolve from sympy.abc import n y = Function('y') /// }}}Define the recurrence as an expression in $a(\cdot)$ that equals zero.
{{{id=86| k = y(n+3)- 7*y(n+2) + 12*y(n+1) - 10*y(n) /// }}}And solve:
{{{id=88| rsolve(k, y(n)) /// }}} {{{id=89| rsolve(k, a(n), {a(0):7, a(1):41, a(2):204}) /// }}}DEMONSTRATION: Modifying Sage source code.
This worksheet available at: http://buzzard.ups.edu/talks.html
{{{id=91| /// }}}