{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Real Double Field\n", "Real Field with 53 bits of precision\n", "3.141592653589793 \n", "3.14159265358979 \n", "3.14159265359\n" ] } ], "source": [ "# 1. Explain the difference between RDF and RR.\n", "# Illustrate with an example.\n", "# Answer: RDF corresponds to the hardware doubles,\n", "# whereas RR has a default precision of 53 bits.\n", "print RDF\n", "print RR\n", "a = RDF(pi)\n", "print a, type(a)\n", "b = RR(pi)\n", "print b, type(b)\n", "print float(pi)" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[26/15, 265/153, 362/209, 989/571, 8733/5042, 18817/10864, 70226/40545, 191861/110771, 716035/413403, 1694157/978122]\n" ] } ], "source": [ "# 2. Give a list of ten rational approximations for sqrt(3),\n", "# accurate with 2, 3, up to 11 decimal places.\n", "reset()\n", "x = sqrt(3)\n", "print [QQ(x.n(digits=k)) for k in range(2,12)]" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30]\n", "9 1\n", "15*x^5 + 4*x^4 + 23*x^3 + 26*x^2 + 6*x + 1\n", "\n", "(15) * (x^5 + 23*x^4 + 16*x^3 + 10*x^2 + 19*x + 29)\n", "True\n" ] } ], "source": [ "# 3. Consider ZZ_31.\n", "# What is the multiplicative inverse of 7?\n", "# Show that 15*x^5 + 4*x^4 + 23*x^3 + 26*x^2 + 6*x + 1\n", "# is irreducible over ZZ_31.\n", "reset()\n", "Z31 = GF(31)\n", "print list(Z31)\n", "i7 = Z31(1/7)\n", "print i7, Z31(7*i7)\n", "x = polygen(Z31)\n", "p = 15*x^5 + 4*x^4 + 23*x^3 + 26*x^2 + 6*x + 1\n", "print p\n", "print type(p)\n", "print p.factor()\n", "print p.is_irreducible()" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "15*x^5 + 4*x^4 + 23*x^3 + 26*x^2 + 6*x + 1\n", "\n", "(15) * (x^5 + 23*x^4 + 16*x^3 + 10*x^2 + 19*x + 29)\n" ] } ], "source": [ "# Instead of polygen(GF(31)), we can also use PolynomialRing,\n", "# which gives the same type of polynomial.\n", "reset()\n", "R. = PolynomialRing(GF(31))\n", "p = 15*x^5 + 4*x^4 + 23*x^3 + 26*x^2 + 6*x + 1\n", "print p\n", "print type(p)\n", "print factor(p)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "cos(2*pi/k) + I*sin(2*pi/k)\n" ] } ], "source": [ "# 4. Expand exp(I*2*Pi/k) as cos(2*Pi/k) + I*sin(2*Pi/k).\n", "reset()\n", "k = SR.symbol('k', domain='real')\n", "t = exp(I*2*pi/k)\n", "print t.real() + I*t.imag()" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "-z^7 + cos(x^3 - 1) + 3*sin(y)\n", "[-z^7, cos(x^3 - 1), 3*sin(y)]\n", "\n", "q0 : [z^7, -1]\n", "q1 : cos [x^3 - 1]\n", "q2 : [sin(y), 3]\n", "q00 : [z, 7]\n" ] } ], "source": [ "# 5. Consider q = cos(x^3 - 1) + 3*sin(y) - z^7. Draw the expression tree of q.\n", "reset()\n", "x, y, z = var('x,y,z')\n", "q = cos(x^3 - 1) + 3*sin(y) - z^7\n", "print q\n", "print q.operands()\n", "print q.operator()\n", "(q0, q1, q2) = q.operands()\n", "print 'q0 :', q0.operator(), q0.operands()\n", "print 'q1 :', q1.operator(), q1.operands()\n", "print 'q2 :', q2.operator(), q2.operands()\n", "q00 = q0.operands()[0]\n", "print 'q00 :', q00.operator(), q00.operands()" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "# Now we have enough information to draw the expression tree.\n", "#\n", "# add\n", "# +-- mul\n", "# | +-- pow\n", "# | | +-- z\n", "# | | +-- 7\n", "# | +-- -1\n", "# +-- cos\n", "# | +-- add\n", "# | +-- pow\n", "# | | +-- x\n", "# | | +-- 3\n", "# | +-- -1\n", "# +-- mul\n", "# +-- sin\n", "# | +-- y\n", "# +-- 3" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "# Let us start representing all leaves.\n", "Lx = LabelledOrderedTree([], label='x')\n", "Ly = LabelledOrderedTree([], label='y')\n", "Lz = LabelledOrderedTree([], label='z')\n", "L7 = LabelledOrderedTree([], label='7')\n", "Lm1 = LabelledOrderedTree([], label='-1')\n", "L3 = LabelledOrderedTree([], label='3')" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ " _mul\n", " / / \n", " pow -1\n", " / / \n", "z 7 " ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# We start with the - z^7.\n", "Lz7 = LabelledOrderedTree([Lz,L7], label='pow')\n", "Lop1 = LabelledOrderedTree([Lz7,Lm1], label='mul')\n", "ascii_art(Lop1)" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ " cos\n", " |\n", " _add\n", " / / \n", " pow -1\n", " / / \n", "x 3 " ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# The second operand is cos(x^3-1)\n", "Lx3 = LabelledOrderedTree([Lx,L3], label='pow')\n", "Lx3m1 = LabelledOrderedTree([Lx3,Lm1], label='add')\n", "Lop2 = LabelledOrderedTree([Lx3m1], label='cos')\n", "ascii_art(Lop2)" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ " mul\n", " / / \n", "3 sin\n", " |\n", " y" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Then the third operand is 3*sin(y)\n", "Lsin = LabelledOrderedTree([Ly], label='sin')\n", "Lop3 = LabelledOrderedTree([L3,Lsin], label='mul')\n", "ascii_art(Lop3)" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ " _______add_______\n", " / / / \n", " _mul cos mul\n", " / / | / / \n", " pow -1 _add 3 sin\n", " / / / / |\n", "z 7 pow -1 y\n", " / / \n", " x 3 " ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Now we can assemble the tree.\n", "tree = LabelledOrderedTree([Lop1,Lop2,Lop3], label='add')\n", "ascii_art(tree)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[ 0.833904329729911 -0.656545663082579 0.766892903499781 -0.485996489420197]\n", "[-0.518898657510414 0.756619263486535 -0.537707140548235 0.846783910822729]\n", "[-0.387030583921130 0.403675706277454 0.572595077475313 0.806936452267780]\n", "M\n", "[ 0.833904329729911 -0.656545663082579 0.766892903499781 -0.485996489420197]\n", "[-0.518898657510414 0.756619263486535 -0.537707140548235 0.846783910822729]\n", "[-0.387030583921130 0.403675706277454 0.572595077475313 0.806936452267780]\n" ] } ], "source": [ "# 6. Save a random matrix to file, destroy its reference,\n", "# and load then the saved matrix back into SageMath.\n", "reset()\n", "M = Matrix(3,4,lambda i,j: RR.random_element())\n", "print M\n", "M.save('/tmp/ourmatrix')\n", "var('M')\n", "print M\n", "M = load('/tmp/ourmatrix.sobj')\n", "print M" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "79*x^298 + 56*x^205 + 49*x^164 + 63*x^121 + 57*x^119 - 59*x^42\n", "[('load_arg', 0), ('ipow', 298), ('load_const', 79), 'mul', ('load_arg', 0), ('ipow', 205), ('load_const', 56), 'mul', 'add', ('load_arg', 0), ('ipow', 164), ('load_const', 49), 'mul', 'add', ('load_arg', 0), ('ipow', 121), ('load_const', 63), 'mul', 'add', ('load_arg', 0), ('ipow', 119), ('load_const', 57), 'mul', 'add', ('load_arg', 0), ('ipow', 42), ('load_const', -59), 'mul', 'add', 'return']\n", "number of operations : 30\n" ] } ], "source": [ "# 7. Consider p = 79*x^298+56*x^205+49*x^164+63*x^121+57*x^119−59*x^42.\n", "# Generate optimized code to evaluate p.\n", "# How many operations are needed to evaluate p?\n", "# Compare with the cost of the direct evaluation of p.\n", "reset()\n", "x = var('x')\n", "p = 79*x^298 + 56*x^205 + 49*x^164 + 63*x^121 + 57*x^119 - 59*x^42\n", "print p\n", "q = fast_callable(p,vars=['x'])\n", "print q.op_list()\n", "print 'number of operations :', len(q.op_list())" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[298, 205, 164, 121, 119, 42]\n", "number of multiplications: 949\n" ] } ], "source": [ "# In the fast_callable evaluation, we have six exponentiations (ipow),\n", "# six multiplications with a coefficient, and five additions.\n", "# In the direct evaluation, without ipow, the number of multiplications is\n", "# the sum of the degrees of the operands of p.\n", "pdegs = [t.degree(x) for t in p.operands()]\n", "print pdegs\n", "print 'number of multiplications:', sum(pdegs)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "[x - y, x, y]\n", "\n", "[x, -y]\n", "\n", "[y, -1]\n" ] } ], "source": [ "# 8. Draw the expression tree of x*y*(x-y).\n", "reset()\n", "var('x,y')\n", "p = x*y*(x-y)\n", "print p.operator()\n", "print p.operands()\n", "p0 = p.operands()[0]\n", "print p0.operator()\n", "print p0.operands()\n", "p1 = p0.operands()[1]\n", "print p1.operator()\n", "print p1.operands()" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "# So the expression tree becomes.\n", "#\n", "# mul\n", "# +-- add\n", "# | +-- x\n", "# | +-- mul\n", "# | +-- y\n", "# | +-- -1\n", "# +-- x\n", "# +-- y" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/plain": [ " _add\n", " / / \n", "x mul\n", " / / \n", " y -1" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# To represent the tree via LabelledOrderedTree,\n", "# we can start defining the leaves: x, y, and -1.\n", "Lx = LabelledOrderedTree([], label='x')\n", "Ly = LabelledOrderedTree([], label='y')\n", "Lm1 = LabelledOrderedTree([], label='-1')\n", "# The x - y is represented as x + (-1)*y\n", "Lmy = LabelledOrderedTree([Ly,Lm1],label='mul')\n", "Lop1 = LabelledOrderedTree([Lx,Lmy],label='add')\n", "ascii_art(Lop1)" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/plain": [ " ___mul__\n", " / / /\n", " _add x y\n", " / / \n", "x mul \n", " / / \n", " y -1 " ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Then the expression tree has Lop1 as the first child,\n", "# the other two children are Lx and Ly.\n", "tree = LabelledOrderedTree([Lop1, Lx, Ly], label='mul')\n", "ascii_art(tree)" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "the exact factorization :\n", "1/108*(3*(1/9*sqrt(26)*sqrt(3) + 1)^(1/3)*(I*sqrt(3) + 1) + 6*x + (-I*sqrt(3) + 1)/(1/9*sqrt(26)*sqrt(3) + 1)^(1/3))*(3*(1/9*sqrt(26)*sqrt(3) + 1)^(1/3)*(-I*sqrt(3) + 1) + 6*x + (I*sqrt(3) + 1)/(1/9*sqrt(26)*sqrt(3) + 1)^(1/3))*(3*x - 3*(1/9*sqrt(26)*sqrt(3) + 1)^(1/3) - 1/(1/9*sqrt(26)*sqrt(3) + 1)^(1/3))\n", "the numerical factorization\n", "(x - 1.52137970680457) * (x + 0.760689853402284 - 0.857873626595179*I) * (x + 0.760689853402284 + 0.857873626595179*I)\n", "\n", "Number Field in a with defining polynomial x^3 - x - 2\n", "(x - a) * (x^2 + a*x + a^2 - 1)\n", "Number Field in b with defining polynomial x^2 + a*x + a^2 - 1 over its base field\n", "(x + b + a) * (x - b) * (x - a)\n" ] } ], "source": [ "# 9. Consider p = x^3 - x - 2 and give all Sage commands for an exact,\n", "# numeric, and symbolic factorization.\n", "# In particular, answer the following questions.\n", "# 1. How can you write p as an *exact* product of\n", "# linear factors, with *exact* complex numbers?\n", "# 2. Compute a *numerical* factorization\n", "# of p over the complex numbers.\n", "# 3. Define a *symbolic* (i.e.: formal) factorization of p,\n", "# declaring sufficiently many roots.\n", "reset()\n", "x = var('x')\n", "p = x^3 - x - 2\n", "r = p.roots()\n", "print 'the exact factorization :'\n", "print prod([x - root[0] for root in r])\n", "print 'the numerical factorization'\n", "q = CC['x'](p)\n", "print q.factor()\n", "print ''\n", "K. = NumberField(p); print K\n", "kp = K['x'](p)\n", "print kp.factor()\n", "kpb = list(kp.factor())[1]\n", "L. = NumberField(kpb[0]); print L\n", "lkp = L['x'](kp)\n", "print lkp.factor()" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(x + y)*(x - y)\n", "(x - y)*z\n", "(x + y)*x - (x + y)*y\n" ] } ], "source": [ "# 10. Give all the sage commands to convert (x-y)*(x+y) into (x+y)*x = (x+y)*y.\n", "reset()\n", "x,y,z = var('x,y,z')\n", "p = (x-y)*(x+y)\n", "print p\n", "q = p.subs({x+y:z})\n", "print q\n", "r = q.expand()\n", "print r.subs(z=x+y)" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "(79.0000000000000*x^5 + 56.0000000000000*x^4 + 49.0000000000000*x^3 + 63.0000000000000*x^2 + 57.0000000000000*x - 59.0000000000000)/(45.0000000000000*x^5 - 8.00000000000000*x^4 - 93.0000000000000*x^2 + 43.0000000000000*x - 62.0000000000000)" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "(1.75555555555556, [1.73033923026938/(x - 1.34613302859424), (-0.209546357895326 - 0.317866843685468*I)/(x - 0.243395948193229 - 0.710544938063163*I), (-0.209546357895326 + 0.317866843685468*I)/(x - 0.243395948193229 + 0.710544938063163*I), (0.122648347698910 + 0.522394270600454*I)/(x + 0.827573573601461 - 1.06277001016849*I), (0.122648347698910 - 0.522394270600454*I)/(x + 0.827573573601461 + 1.06277001016849*I)])\n" ] } ], "source": [ "# 11. Consider r = (79*x^5 + 56*x^4 + 49*x^3 + 63*x^2 + 57*x - 59)/(45*x^5 \n", "# - 8*x^4 - 93*x^2 + 43*x - 62). \n", "# How many operations does it take to evaluate r\n", "# with a partial fraction decomposition? \n", "# Compare this number of operations with the Horner forms\n", "# of numerator and denominator for r.\n", "reset()\n", "R. = PolynomialRing(CC)\n", "r = (79*x^5 + 56*x^4 + 49*x^3 + 63*x^2 + 57*x - 59)/(45*x^5 - 8*x^4 - 93*x^2 + 43*x - 62)\n", "show(r)\n", "q = r.partial_fraction_decomposition()\n", "print q" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "# If we would use the Horner form to evaluate numerator and denominator,\n", "# then we need 5 multiplications and 5 additions, with one division.\n", "# The partial fraction decomposition gives 5 divisions and 5 (complex) additions." ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(x^100 - 1)/(x - 1)\n", "x^99 + x^98 + x^97 + x^96 + x^95 + x^94 + x^93 + x^92 + x^91 + x^90 + x^89 + x^88 + x^87 + x^86 + x^85 + x^84 + x^83 + x^82 + x^81 + x^80 + x^79 + x^78 + x^77 + x^76 + x^75 + x^74 + x^73 + x^72 + x^71 + x^70 + x^69 + x^68 + x^67 + x^66 + x^65 + x^64 + x^63 + x^62 + x^61 + x^60 + x^59 + x^58 + x^57 + x^56 + x^55 + x^54 + x^53 + x^52 + x^51 + x^50 + x^49 + x^48 + x^47 + x^46 + x^45 + x^44 + x^43 + x^42 + x^41 + x^40 + x^39 + x^38 + x^37 + x^36 + x^35 + x^34 + x^33 + x^32 + x^31 + x^30 + x^29 + x^28 + x^27 + x^26 + x^25 + x^24 + x^23 + x^22 + x^21 + x^20 + x^19 + x^18 + x^17 + x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1\n", "2.75666137216650\n" ] } ], "source": [ "# 12. Explain why normal forms are so important to symbolic computation. \n", "# What can we do if a normal form is too expensive to compute? Illustrate with a good example.\n", "# Answer: with a normal form we can decide whether two expressions are equal to each other.\n", "# The normalization for rational expression is the removal of common factors in the numerator\n", "# and the denominator. This may lead to expression swell, as the example below shows.\n", "# We can apply a probability one numerical evaluation test.\n", "# We generate a random number and evaluate the expression to see if the expression is zero or not.\n", "reset()\n", "x = var('x')\n", "p = (x^100 - 1)/(x - 1)\n", "print p\n", "print p.numerator()\n", "print p(x=RR.random_element())" ] } ], "metadata": { "kernelspec": { "display_name": "SageMath 8.3", "language": "", "name": "sagemath" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.15" } }, "nbformat": 4, "nbformat_minor": 2 }