#!/usr/bin/env python
# coding: utf-8

# # Инфраструктура Python. CAS

# Символьная (и точная) математика в Python делается через библиотеку `sympy` (входит в Anaconda). Она несколько уступает CAS с многолетней историей вроде Maple или Mathematica, зато тесно интегрирована с Python http://docs.sympy.org/latest/index.html
# 
# Попробовать `sympy` можно онлайн http://gamma.sympy.org/

# In[1]:


import networkx, PIL
from io import BytesIO
import itertools
import os

get_ipython().run_line_magic('pylab', 'inline')


# In[2]:


import sympy
from sympy import *
from sympy.solvers.diophantine import diophantine
from sympy.ntheory.factor_ import digits

from sympy.abc import h, k, n, t, x, y, z


# In[3]:


sympy.__version__


# ## Символьные вычисления

# Библиотека `sympy` предназначена для символьных вычислений. Она представляет математические выражения (переопределяя сотни имен `numpy`) в виде деревьев, над которыми может проводить различные операции, например, интегрирование.

# In[4]:


expr = integrate(tanh(x), x)
expr


# In[5]:


print expr
print srepr(expr)


# In[6]:


sympy.printing.tree.print_tree(expr)


# In[7]:


import sympy.printing.dot
dotstr = sympy.printing.dot.dotprint(expr)
graph = networkx.nx_pydot.read_dot(BytesIO(dotstr))

print expr
print
print '\r\n'.join(list(graph.nodes))

A = networkx.nx_agraph.to_agraph(graph)
for node in A.nodes():
    if not node.attr['fontname']:
        node.attr['fontsize'] = 8
        node.attr['fontname'] = 'Verdana'
        node.attr['shape'] = 'none'
        node.attr['style'] = 'filled'
        node.attr['color'] = '#F0F0F0'

io = BytesIO()
A.draw(io, format='png', prog='dot')
PIL.Image.open(io)


# Ввести собственные символы (или распарсить выражение) можно, используя `Symbol`, `var` (`symbols`), `S` (`sympify`).

# In[8]:


Symbol('x') ** 2


# In[9]:


sum(var('a b c'))


# In[10]:


x1, x2 = symbols('x_1 x_2')
solve(x1 ** 2 + x2 ** 2 - 1)


# In[11]:


[S('1/2') ** i for i in range(5)]


# In[12]:


solve(S('x1 ** 2 + x2 ** 2 - 1'))


# ## Отображение выражений

# Выражения `sympy` можно просматривать в различных формах, конвертировать в код на разных языках программирования, впрочем, не всегда успешно.

# In[13]:


alpha, beta, x = var('alpha beta x')
expr = (sin(alpha)**2 + 2*cos(beta)) / x
print latex(expr)
print ccode(expr)
print jscode(expr)
print mathematica_code(expr)


# In[14]:


print latex(divisor_sigma(n))
print
print ccode(divisor_sigma(n))
print
print jscode(divisor_sigma(n), n)
print
print mathematica_code(divisor_sigma(n)) # должно быть DivisorSigma[1,n]


# In[15]:


print expr


# In[16]:


pprint(expr)


# In[17]:


import IPython
display(IPython.display.Latex('$%s$' % latex(expr)))


# Оставить один из этих форматов вывода в качестве дефолтного в IPython можно одним из трех вызовов

# In[18]:


init_printing(use_latex=False, pretty_print=False)
expr


# In[19]:


init_printing(use_latex=False, pretty_print=True)
expr


# In[20]:


init_printing(use_latex='mathjax')
display(expr)
display([1, 2, 3])


# Просматривать встроенные типы в виде $\LaTeX$ может быть неудобно, так что стоит использовать

# In[21]:


init_printing(use_latex='mathjax', print_builtin=False)
display(expr)
display([1, 2, 3])


# ## Стандартные операции над выражениями

# In[22]:


print N(pi, 100)
print pi.evalf(100)


# In[23]:


print type(N(pi))
print type(float(pi))


# In[24]:


print N(I ** I)


# In[25]:


simplify((x**2 - 4)/((x+3)*(x-2)))


# In[26]:


div(x**2 - 4 + x, x-2)


# In[27]:


factor(x**4/2 + 5*x**3/12 - x**2/3)


# In[28]:


apart(x/((x+2)*(x+1)), x)


# In[29]:


together(1/(x+2) + 1/(x+1))


# In[30]:


S('(a+sqrt(b))**5').expand()


# In[31]:


expand((x+y)**5)


# In[32]:


expand_trig(sin(3*x))


# In[33]:


sin(3*x).rewrite(tan)


# In[34]:


((x+y)**5).subs(y, 2)


# In[35]:


((x+y)**5).subs([(x, 1), (y, sqrt(y))]).expand()


# In[36]:


p = expand(Poly('x + y') ** 10)
print p.coeffs()
print p.monoms()


# Выражение можно превратить в функцию

# In[37]:


f = Lambda(x, x ** 2)
print f(5)
print integrate(f(t), t)


# ## solve

# Функция `solve` это причина, по которой CAS вообще нужны

# In[38]:


print solve(x ** 2 - 5 * x + 6)


# In[39]:


solve(x ** 2 - 5 * x + 6 > 0)


# In[40]:


map(display, solve(x**3 + 4*x + 181, x));


# In[41]:


solve(sin(x) - S('1/sqrt(2)'))


# In[42]:


solveset(sin(x) - S('1/sqrt(2)'))


# In[43]:


solve(exp(x) + x - 2)[0]


# In[44]:


solve([x + y - 5, y - x - 6], (x,y))


# Увы, пока что `sympy` не умеет решать системы неравенств, в отличие от Sage.

# In[45]:


solve([x + y > 5, y - x < 6], (x,y))


# ## Матан

# Предел

# In[46]:


limit(sin(x) / x, x, 0)


# In[47]:


limit((cos(x) - 1) / x ** 2, x, 0)


# Дифференцирование

# In[48]:


expr = sin(x)**2 + 2*cos(x)


# In[49]:


expr.diff(x)


# In[50]:


diff(sin(x) / x)


# In[51]:


diff(x / (x**2 + 2*x + 1), x)


# In[52]:


[diff(sin(2 * x), x, u) for u in range(5)]


# Интегрирование

# In[53]:


expr.integrate(x)


# In[54]:


integrate(x / (x**2 + 2*x + 1), x)


# In[55]:


integrate(sin(x)/x, x)


# In[56]:


sigma = var('sigma')
integrate(exp(-x**2 / 2 / sigma ** 2) / sqrt(2 * pi) / sigma, x)


# In[57]:


integrate(sin(x)/x, (x, -oo, oo))


# Ряд Тейлора

# In[58]:


series(tan(x), x, 0, 15)


# Ряды Фурье

# In[59]:


ser = fourier_series(x**2, (x, -pi, pi))
ser


# In[60]:


fourier_transform(exp(-x ** 2), x, k)


# In[61]:


inverse_fourier_transform(sqrt(pi)*exp(-(pi*k)**2), k, x)


# In[62]:


serfn = lambdify(x, ser.truncate(3))

xx = linspace(-numpy.pi, numpy.pi)
matplotlib.pyplot.plot(xx, xx ** 2);
matplotlib.pyplot.plot(xx, serfn(xx));


# Векторные поля

# In[63]:


from sympy.vector import CoordSys3D, express, divergence, curl, gradient, scalar_potential


# In[64]:


R = CoordSys3D('R')
vector_field = R.x * R.y * R.i + R.y ** 2 * R.j + R.z ** 2 * R.k
display(divergence(vector_field))
display(curl(vector_field))


# In[65]:


scalar_field = R.x * R.y + R.y ** 2 + R.z ** 2
gr = gradient(scalar_field)
gr


# In[66]:


print sympy.vector.is_conservative(gr)
print sympy.vector.is_solenoidal(gr)


# In[67]:


scalar_potential(R.y * R.i + (R.x + 2 * R.y) * R.j + R.z * R.k, R)


# Дифференциальные уравнения

# In[68]:


dsolve(y(t).diff(t, t) - y(t) - exp(t), y(t))


# In[69]:


dsolve((Eq(x(t).diff(t), 2 * x(t) + 3 * y(t)), Eq(y(t).diff(t), 1 * x(t) + 4 * y(t))))


# Урчпы

# In[70]:


f = Function('f')
u = f(x, y)
u_x = u.diff(x)
u_y = u.diff(y)
pdsolve(Eq(1 + (2*(u_x/u)) + (3*(u_y/u))))


# ## Линал

# Для вещественных чисел в `sympy` определены функции `sign Min Max`, для комплексных `re im arg conjugate Abs`.

# Матрицы, транспонирование, арифметика, определитель, след, обратная, характеристический полином, собственные значения, жорданова форма, канонический базис

# In[71]:


Matrix(4,4, lambda i,j: (i+1)*(j+1))


# In[72]:


M = Matrix([[1,2],[3,4]])
display(M.T)
display(M * M)
display(M.multiply_elementwise(M))


# In[73]:


m = Matrix([[S('1/2'), S('1/2'), S('0')], [S('1/2'), S('0'), S('1/2')], [S('0'), S('1/2'), S('1/2')]])
m


# In[74]:


det(m), trace(m)


# In[75]:


m.inv()


# In[76]:


FiboM = Matrix([[1,1],[1,0]])
FiboM


# In[77]:


FiboM ** (-1)


# In[78]:


FiboM ** 8


# In[79]:


print FiboM.det(), FiboM.trace()


# In[80]:


FiboM.charpoly(var('lamda'))


# In[81]:


FiboM.eigenvals()


# In[82]:


FiboM.eigenvects()


# In[83]:


kb, jf = simplify(FiboM.jordan_form())
display(jf)
display(kb)


# In[84]:


simplify(kb * jf * kb.inv())


# In[85]:


simplify(kb * (jf ** k) * kb.inv())[1,0]


# In[86]:


x1, y1, z1, x2, y2, z2 = var('x1 y1 z1 x2 y2 z2')
cross([x1,y1,z1], [x2,y2,z2])


# In[87]:


Matrix(symarray('a', (3,))).norm()


# In[88]:


Matrix([x1,y1,z1]).norm()


# In[89]:


X = MatrixSymbol('X', 3, 3)
Y = MatrixSymbol('Y', 3, 3)
(X * Y)[1, 2]


# ## Конкретная математика

# Суммирование

# In[90]:


summation(k ** 2, (k, 1, n))


# In[91]:


summation(x ** k, (k, 0, oo))


# Не все суммы поддаются движку sympy

# In[92]:


A, B = var('A B', positive=True, integer=True)
summation(Abs(x - y), (x, 1, A), (y, 1, B))


# Решение рекуррентности (вывод формулы Бине)

# In[93]:


rsolve(y(n+2)-y(n+1)-y(n), y(n), {y(0): 0, y(1): 1})


# Цепная дробь

# In[94]:


list(continued_fraction_iterator(Rational(3, 8)))


# $$\frac{3}{8} = 0 + \frac{1}{2 + \frac{1}{1 + \frac{1}{2}}}$$

# In[95]:


list(itertools.islice(continued_fraction_convergents(continued_fraction_iterator(pi)), 0, 10))


# In[96]:


p, q, r = var('p q r')
k2 = to_dnf((p >> q) & r)
print k2
print list(satisfiable(k2, all_models=True))


# ## Теория чисел

# In[97]:


digits(13, 2)


# In[98]:


print primepi(1000)
print list(primerange(1, 1000))


# In[99]:


print prime(100)
print isprime(2 ** 127 - 1)
print nextprime(1000000000)


# In[100]:


print factorint(84759364534527)


# In[101]:


get_ipython().run_cell_magic('time', '', 'factorint(2 ** 67 - 1)\n')


# Профессор Коул к 1903 году потратил на это «все воскресенья в течение трех лет».

# In[102]:


print gcd(244, 68), gcdex(244, 68), -5 * 244 + 18 * 68


# In[103]:


print sum(divisors(284)[:-1]), sum(divisors(220)[:-1])


# In[104]:


print divisor_count(12), divisor_sigma(12)
print divisors(12)
print totient(12)


# Китайская теорема

# In[105]:


from sympy.ntheory.modular import crt
print crt([99, 97, 95], [49, 76, 65])
print 639985 % 99, 639985 % 97, 639985 % 95, 99 * 97 * 95


# Проверка обращения Мёбиуса $G(n) = \sum_{d|n} F(d) \Leftrightarrow F(n) = \sum_{d|n} \mu(d) G(\frac{n}{d})$

# In[106]:


j_divides_i = Matrix(fromfunction(lambda i,j: where((i+1) % (j+1) == 0, 1, 0), (20, 20)))
inv_j_divides_i = j_divides_i.inv()
inv_j_divides_i


# In[107]:


mo = [[int(mobius((i + 1) / (j + 1))) if (i + 1) % (j + 1) == 0 else 0 for j in range(20)] for i in range(20)]
Matrix(mo) == inv_j_divides_i


# In[108]:


diophantine(2*x + 3*y - 5)


# In[109]:


diophantine(x**2 - 4*x*y + 8*y**2 - 3*x + 7*y - 5)


# In[110]:


get_ipython().run_cell_magic('time', '', 'a = 30741415699351712223222592216948664852787225938672761888182663459401433940293299920472821794247792558438137266179003170954852009167996\nb = 5493768938913164268735898548345740984062864198582915532893452862692397723477349340213609803716310002722689784104429402408034399680657214203589114662020419403506380275921075865806868813167390548527002018378296772011002830464813992537445273710373363121247809539380530302504021135339366740636410963119289644130568327436166142653219502858785723203265988729334812946807271940413511096470449043813939925759131485224930805050595388961884769573807597931344321170355620526497042548348665782820474459300583235597864519400688289208986957832681459966933968204573143630154199268824952951642924625108969237666308633980889261559554117298872683137210357137988177192697560524562835477922442010038014370796719316949327823897102946508453004777622812045656957796013010907193455701053005247975317200627083352360327185314875288693488759973999715297338413318230762734322668133288476183336301663326505090120348966948964517423252467946014780068274519231694498484297\nprint max(list(diophantine(x**2 + a*(y**2) - b)))\n')


# ## Перестановки

# In[111]:


from sympy.combinatorics import Permutation


# In[112]:


p = Permutation([0, 3, 2, 1])
q = Permutation([1, 2, 3, 0])
print p, q


# In[113]:


print Permutation([0,1,2,4,3]).rank()
print Permutation.unrank_lex(5, 3).array_form


# In[114]:


print p.cycles, p.cycle_structure, q.cycle_structure


# In[115]:


print p.inversions(), p.is_even, p.parity()


# In[116]:


print p.support(), q.support()


# In[117]:


print p.order(), q.order()


# In[118]:


list(Permutation.from_sequence('eafbcd'))


# In[119]:


print p.rank(), p.next_lex().array_form


# In[120]:


print list(p * q), list(q * p)


# In[121]:


q.transpositions()


# In[122]:


pp = Permutation(0, 2)(1, 3)
print pp(range(4)), pp.list(), pp.cyclic_form, pp.array_form


# In[123]:


Permutation([[1,2]], size=5).list()


# In[124]:


print Permutation([[1,2]], size=5).cyclic_form, Permutation([[1,2]], size=5).full_cyclic_form


# In[125]:


Permutation(0, 4) == Permutation([4, 1, 2, 3, 0])


# In[126]:


p.length()


# In[127]:


[i ^ p ** 3 for i in range(4)]


# In[128]:


q(p(range(4))), (p * q)(range(4))


# ## Группы

# In[129]:


print list(sympy.combinatorics.dihedral(5))


# In[130]:


for gr in 'symmetric cyclic alternating dihedral'.split():
    print gr, len(list(getattr(sympy.combinatorics, gr)(5)))


# Мультипликативная группа, порядок элемента в $Z_n^*$, степень по модулю, корень по вычету, дискретный логарифм

# In[131]:


pow(2, 100, 1000)


# $2^{504} \equiv 1 \pmod{1009} \\
# 2^{503} \equiv 2^{-1} \equiv 505 \pmod{1009}$

# In[132]:


print pow(2, 504, 1009), n_order(2, 1009), mod_inverse(2, 1009), pow(2, 504-1, 1009)


# In[133]:


print nthroot_mod(505, 503, 1009), discrete_log(1009, 2, 505)


# ## Комбинаторика

# In[134]:


print [fibonacci(k) for k in range(10)]


# In[135]:


print factorial(100)
print subfactorial(100)


# In[136]:


print [binomial(10, k) for k in range(11)]


# Число Каталана — число правильных скобочных последовательностей длины $2n$

# In[137]:


print [catalan(k) for k in range(10)]


# Числа Стирлинга 1-го рода $\left[ n \atop k \right]$ — коэффициенты при $x^k$ разложения $x^{\underline{n}} = x (x-1) \ldots (x-n+1)$.
# 
# Числа Стирлинга 2-го рода  $\left\{ {n \atop k} \right\}$ — коэффициенты при $x^{\underline{k}}$ разложения $x^n$. Заодно число разбиений множества размера $n$ на $k$ непустых подмножеств.

# In[138]:


k = var('k')
expand(product(x-k, (k, 0, 9)))


# In[139]:


print [sympy.combinatorial.numbers.stirling(10, k, kind=1) for k in range(11)]


# In[140]:


print [sympy.combinatorial.numbers.stirling(10, k, kind=2) for k in range(11)]


# Число Белла — число разбиений множества размера $n$

# In[141]:


print [sympy.combinatorial.numbers.bell(k) for k in range(10)]


# Числа Бернулли известны прежде всего из-за их роли в формуле Фаульхабера https://en.wikipedia.org/wiki/Faulhaber%27s_formula о закрытой форме суммы $\sum_{k=1}^n k^p$.

# In[142]:


print [bernoulli(k) for k in range(10)]


# Число разбиений

# In[143]:


npartitions(50)


# In[144]:


import sympy.combinatorics
pa = sympy.combinatorics.partitions.IntegerPartition([1] * 6)
while True:
    print pa.partition
    if len(pa.partition) == 1: break
    pa = pa.next_lex()


# In[145]:


print list(sympy.combinatorics.subsets.GrayCode(4).generate_gray())


# In[146]:


ss = sympy.combinatorics.Subset([], 'abcd')
while True:
    print ss.rank_gray, ss.subset
    ss = ss.next_gray()
    if ss.subset == []: break


# In[147]:


ss = sympy.combinatorics.Subset([], 'abcd')
while True:
    print ss.rank_binary, ss.subset
    ss = ss.next_binary()
    if ss.subset == []:
        break


# In[148]:


list(sympy.combinatorics.subsets.ksubsets('abcd', 2))


# ## Геометрия

# In[149]:


print Line(Point(1,1), slope=S('2/5'))


# In[150]:


print Segment(Point(0,0), Point(1,1)).midpoint


# In[151]:


A = Point(0,0)
B = Point(3,0)
C = Point(0,4)
print Point.is_collinear(A, B, C)
print Line(A, B).is_perpendicular(Line(A, C))
print

t = Triangle(A, B, C)
print t
print t.area, t.perimeter
print t.angles
print

print 'Медианы'
print t.medians
print Line.are_concurrent(t.medians.values())
print t.centroid # intersection(t.medians.values()[0], t.medians.values()[1])
print

print 'Высоты'
print t.altitudes
print Line.are_concurrent(t.altitudes.values())
print intersection(t.altitudes.values()[0], t.altitudes.values()[1])
print

print 'Биссектрисы'
print t.bisectors()
print Line.are_concurrent(t.bisectors().values())
print intersection(t.bisectors().values()[0], t.bisectors().values()[1])
print 'Вписанная окружность', t.incenter, t.inradius, t.incircle
print

AB = Segment(A, B)
AC = Segment(A, C)
print AB.perpendicular_line(AB.midpoint).intersect(AC.perpendicular_line(AC.midpoint))
print 'Описанная окружность', t.circumcenter, t.circumradius, t.circumcircle
print

print t.bounds


# In[152]:


Line(A, B).angle_between(Line(B, C))


# In[153]:


print Line(B, C).arbitrary_point()


# In[154]:


Line(A, B).distance(Point(C))


# In[155]:


Segment(B, C).length


# In[156]:


Ray(A, B).distance(Point(-1,1))


# In[157]:


Segment(B, C).distance(A)


# In[158]:


print Segment(B, C).parallel_line(Point(5,5))
print Segment(B, C).perpendicular_line(Point(5,5))
print Segment(B, C).perpendicular_segment(Point(5,5))


# In[159]:


print Line(B, C).projection(Point(5,5))


# In[160]:


print Line(B, C).projection(Segment(Point(5,5), Point(6,6)))


# In[161]:


print AB.perpendicular_bisector()


# In[162]:


print t.rotate(pi/3)


# In[163]:


L = Line(C, B)
L2 = Line(A, B)
print L.equation(), L.coefficients
print L.intersect(L2)


# Задача с тестового тура ПВГ: В равнобедренном треугольнике известно основание $b=126$ и медиана $m_a = \frac{\sqrt{35997}}{2}$, проведенная к боковой стороне. Найдите расстояние между центрами вписанной и описанной окружностей.

# In[164]:


h = var('h')
A = Point(0, 0)
B = Point(63, h)
C = Point(126, 0)
T = Triangle(A, B, C)
hh = max(solve(T.medians[A].length - sqrt(35977) / 2))
print 'h =', hh
B = Point(63, hh)
T = Triangle(A, B, C)
print Segment(T.incircle.center, T.circumcircle.center).length


# In[165]:


el = Ellipse(Point(5,5), 1, 2)
el.equation()


# In[166]:


print el.rotate(pi/2, Point(0,0))


# In[167]:


print intersection(el, Line(Point(2,2), Point(6,7)))


# In[168]:


el.intersect(L)


# In[169]:


list(el.intersect(Line(Point(0,0), Point(5,5))))[0].x


# In[170]:


el.area


# In[171]:


el.circumference


# In[172]:


el.eccentricity


# In[173]:


print el.foci


# In[174]:


print el.tangent_lines(Point(5,3))


# In[175]:


print t


# In[176]:


Plane(Point(1,2,3), normal_vector=(1,1,1)).distance(Point(0,0,0))


# In[177]:


pi1 = Plane(Point(1,2,2), Point(2,5,4), Point(3,2,1))
pi1.equation()


# In[178]:


pi1.angle_between(Line(Point(1,2,3), Point(5,1,2)))


# In[179]:


pi2 = Plane(Point(1,2,2), Point(2,5,4), Point(3,2,1))
pi3 = Plane(Point(1,2,3), Point(5,1,2), Point(-1,2,0))
pi2.angle_between(pi3)


# In[180]:


print dot(pi2.normal_vector, pi3.normal_vector), pi2.normal_vector, pi3.normal_vector


# In[181]:


print pi2.projection(Point(0,0,0))


# In[182]:


numpy.random.seed(0)
xy = numpy.random.rand(12, 2)
points = [Point(x,y) for x,y in xy]

poly = convex_hull(*points)
print len(poly.vertices)

xx = list([float(pt.x) for pt in poly.vertices])
yy = list([float(pt.y) for pt in poly.vertices])

matplotlib.pyplot.gca().set_aspect('equal')
matplotlib.pyplot.scatter([pt.x for pt in points], [pt.y for pt in points]);
matplotlib.pyplot.plot(xx + xx[0:1], yy + yy[0:1]);


# ## Графики

# Библиотека `sympy` умеет рисовать графики, используя библиотеку `matplotlib`, немного оборачивая ее для автоматического выбора диапазона, адаптивного числа точек, перенесения осей в 0, а заодно и затрудняя использование возможностей `matplotlib`, для которых обертки нет. Манки-патч в функции `plot_plus` немного решает эту проблему для самых вопиющих вещей.

# In[183]:


def plot_plus(plot_fn, *args, **kwargs):
    parse_args = {'set_aspect', 'figsize'}
    if set(kwargs.keys()) & parse_args:
        kw = {k:v for k,v in kwargs.items() if k not in parse_args}
        kw['show'] = False
        p = plot_fn(*args, **kw)
        p._backend = p.backend(p)
        if 'figsize' in kwargs:
            p._backend.fig.set_size_inches(kwargs['figsize'])
        if 'set_aspect' in kwargs:
            p._backend.ax.set_aspect(kwargs['set_aspect'])
        p._backend.show()
        return p
    else:
        return plot_fn(*args, **kwargs)

import functools
for fn in 'plot plot_parametric plot3d plot3d_parametric_line plot3d_parametric_surface plot_implicit'.split():
    globals()[fn] = functools.partial(plot_plus, getattr(sympy.plotting, fn))


# In[184]:


x, y, t, u = var('x y t u')


# In[185]:


plot(x**2);


# In[186]:


plot(sin(x));


# In[187]:


plot(x, x**2, x**3, (x, -2, 2));


# In[188]:


plot(x ** 2, xlim=(-2,2), ylim=(0,2), set_aspect='equal', figsize=(12,6));


# In[189]:


plot(1/x, x, (x, -5, 5), ylim=(-5, 5), set_aspect='equal', figsize=(8,8));


# In[190]:


plot_parametric(cos(t * 5) * cos(t), cos(t * 5) * sin(t), set_aspect='equal');


# In[191]:


plot3d(sin(x**2 + y**2), (x, -2, 2), (y, -2, 2), nb_of_points_x=200, nb_of_points_y=200);


# In[192]:


plot3d_parametric_surface(cos(t), sin(t) * sin(u), sin(t) * cos(u), set_aspect='equal');


# In[193]:


plot_implicit(x + y < 2, line_color='yellow', set_aspect='equal');


# In[194]:


plot_implicit(x ** 2 + y ** 2 + x * y - 2, line_color='red', set_aspect='equal');


# ## Статистика

# In[195]:


import sympy.stats
from sympy.stats import Normal, Die, P, E, variance, std, cdf, density, sample


# In[196]:


p, t = symbols('p t', positive=True, real=True)
x = sympy.stats.Binomial('x', 5, p)
simplify(E(exp(t*x)))


# In[197]:


Twodices = Die('X', 6) + Die('Y', 6)


# In[198]:


throws = [int(sample(Twodices)) for i in range(1000)]


# In[199]:


P(Eq(Twodices, 6))


# In[200]:


twodices = density(Twodices)
twodices


# In[201]:


hist(throws, range=(0.5, 12.5), bins=12);
bar(range(13), [1000 * float(twodices.get(i, 0)) for i in range(13)], color='r', alpha=0.3);


# In[202]:


iq = sympy.stats.Normal('iq', 100, 15)
print E(iq), std(iq), variance(iq)
print N(P(iq > 150))
display(density(iq)(x))
display(simplify(cdf(iq)(x)))
print [sample(iq) for i in range(10)]


# In[203]:


float(E(iq, iq > 100))


# In[204]:


P(iq > 150).rewrite(Integral)


# ## Физика

# In[205]:


from sympy.physics import units


# In[206]:


print units.acceleration_due_to_gravity.convert_to(units.m / units.s ** 2)


# In[207]:


print units.year.convert_to(units.second)


# In[208]:


print units.convert_to(1 * units.pound * units.acceleration_due_to_gravity * 1 * units.foot, units.joule).n()


# Третий закон Кеплера

# In[209]:


T = symbols("T")
a = units.Quantity('venus_a', units.length, 108208000e3 * units.meter)
M = units.Quantity('solar_mass', units.mass, 1.9891e30 * units.kilogram)

G = units.gravitational_constant
eq = Eq(T**2 / a**3, 4*sympy.pi**2 / G / M)
q = solve(eq, T)[1]
print q
print units.convert_to(q, units.day).n()


# ## Интероп с Python

# Выражение, результат символьного вычисления, можно сунуть в питон-функцию

# In[210]:


def naive_sum_of_squares(n):
    result = 0
    for i in range(1, n + 1):
        result += i ** 2
    return result


# In[211]:


get_ipython().run_line_magic('timeit', 'naive_sum_of_squares(1000000)')


# In[212]:


def fast_naive_sum_of_squares(n):
    return sum(arange(1, n + 1) ** 2)


# In[213]:


get_ipython().run_line_magic('timeit', 'fast_naive_sum_of_squares(1000000)')


# In[214]:


k, n = var('k n')
sum_of_squares = lambdify(n, summation(k ** 2, (k, 1, n)))
print sum_of_squares(5)


# In[215]:


get_ipython().run_line_magic('timeit', 'sum_of_squares(1000000)')


# In[216]:


print lambdify(x, Si(x), modules='sympy')(5.0)


# ## Интероп с Maple

# Иногда вызов внешней CAS удобнее или быстрее дает ответ

# In[217]:


def eval_maple(expr):
    import os, delegator
    cmaple_dir = 'c:\\Program Files (x86)\\Maple 15\\bin.win'
    cmaple_exe = 'c:\\Program Files (x86)\\Maple 15\\bin.win\\cmaple.exe'
    cd = os.getcwdu()
    os.chdir(cmaple_dir)
    maple = delegator.run([cmaple_exe], block=False)
    maple.expect('>')
    maple.send('interface(prettyprint=0);')
    maple.expect('>')
    maple.send(expr)
    maple.expect('>')
    maple.send('quit();')
    os.chdir(cd)
    return '\r\n'.join(maple.out.splitlines()[3:-2])


# In[218]:


integrate(sin(exp(x)), x)


# In[219]:


S(eval_maple('int(sin(exp(x)), x);'))


# In[220]:


get_ipython().run_cell_magic('time', '', 'print sympy.factorint(987345823683476598732645172)\n')


# In[221]:


get_ipython().run_cell_magic('time', '', "print eval_maple('ifactor(987345823683476598732645172);')\n")


# In[222]:


def parse_ifactor_output(s):
    result = {}
    import re
    s = re.sub('\s+', '', s)
    for item in s.split('*'):
        m = re.match('^``\((\d+)\)(?:\^(\d+))?$', item)
        if not m:
            raise Exception('ifactor says ' + s + ', item = ' + item)
        base = int(m.groups()[0])
        power = int(m.groups()[1]) if m.groups()[1] is not None else 1
        result[base] = power
    return result

def maple_factorint(n):
    output = eval_maple('ifactor(%d);' % n)
    return parse_ifactor_output(output)


# In[223]:


get_ipython().run_cell_magic('time', '', 'print maple_factorint(987345823683476598732645172)\n')