{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\def\\CC{\\bf C}\n",
    "\\def\\QQ{\\bf Q}\n",
    "\\def\\RR{\\bf R}\n",
    "\\def\\ZZ{\\bf Z}\n",
    "\\def\\NN{\\bf N}\n",
    "$$\n",
    "# Exact and float point computations - Stability and unstability in dynamics\n",
    "\n",
    "There are a lot of different numbers in Sage. You will need to choose\n",
    "which kind of numbers you want to use depending on your computations. In\n",
    "this worksheet we will consider the following types of numbers:\n",
    "\n",
    "-   integers and rationals\n",
    "-   floating point numbers\n",
    "-   symbolic expressions\n",
    "-   algebraic numbers\n",
    "\n",
    "## Integers, rationals and floating point\n",
    "\n",
    "To create an integer or a rational number, you just write it as you\n",
    "would do on a sheet of paper"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "245"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "2 + 3^5    # an integer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "23/45"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "23 / 45    # a rational number"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To create a floating point number, you need to add a dot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.00000000000000"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "1.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "25.4300000000000"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "3.25 + 22.18"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Contrarily to integers and rationals, a floating point number has\n",
    "limited precision. Compare the results of the two following cells"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1024"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "2^100 + 2^10 - 2^100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.000000000000000"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "2.0^100.0 + 2.0^10.0 - 2.0^100.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Can you explain the above behavior?\n",
    "\n",
    "**Exercise:** Find the smallest power `n` so that `2^n + 1 - 2^n` and\n",
    "`2.0^n + 1 - 2.0^n` provide different answers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "# enter your commands here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you have an object `a` and want to know what kind of number it is,\n",
    "you can use one of the functions `parent` or `type`. The former returns\n",
    "the set in which your object belongs while the second one returns the\n",
    "type of the object (in a computer programming sense)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = 2\n",
    "b = 3/2\n",
    "c = 5.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Integer Ring\n",
      "Rational Field\n",
      "Real Field with 53 bits of precision\n"
     ]
    }
   ],
   "source": [
    "print parent(a)\n",
    "print parent(b)\n",
    "print parent(c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<type 'sage.rings.integer.Integer'>\n",
      "<type 'sage.rings.rational.Rational'>\n",
      "<type 'sage.rings.real_mpfr.RealLiteral'>\n"
     ]
    }
   ],
   "source": [
    "print type(a)\n",
    "print type(b)\n",
    "print type(c)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Maps of the interval: fixed points and iteration\n",
    "\n",
    "Let us consider the map $f_4(x) = 4 x (1 - x)$ from the interval $[0,1]$\n",
    "to itself.\n",
    "\n",
    "**Exercise:** Prove that $f_4$ is surjective and plot it"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# enter your commands here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Show that $f_4(3/4) = 3/4$ (in other words, the point $3/4$ is a fixed\n",
    "point of $f_4$). What do you expect from the following two commands (you\n",
    "have to guess whether the answer will be `True` or `False` and you can\n",
    "then check your answer by executing the cell):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "s = 3.0 / 4.0\n",
    "4 * s * (1 - s) == s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let $f_{7/2}(x) = \\frac{7}{2} x (1 - x)$.\n",
    "\n",
    "Prove that $f_{7/2}$ is not surjective on $[0,1]$ and that $5/7$ is a\n",
    "fixed point of $f_{7/2}$.\n",
    "\n",
    "**Exercise:** Is the following `True` or `False`? (you have to guess\n",
    "before executing the cell)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "s = 5.0 / 7.0\n",
    "7.0 / 2.0 * s * (1.0 - s) == s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Perform the same two computations as above with rational numbers instead\n",
    "of floating point."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# enter your commands here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On a computer a floating point number is a number of the form $m\\, 2^n$\n",
    "where $m$ (the mantissa) and $n$ (the exponent) have some fixed bounds.\n",
    "In particular, floating point numbers have *finite* precision.\n",
    "Computations with floating points numbers are inaccurate but very\n",
    "efficient.\n",
    "\n",
    "Compare the following computation with rationals:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "45842869094724092282256664559525216692173091162526497018856817207453767127095354052418072639857700888691134361639497243771528176716782876457748796945284249594870200383920805322146999777923783319507727283850422831309915607777814766948159289963355518294865659883896712136501856004613518112550490551045680741845733583765205748593300762225277706242970719821414680060203467479543923750220952764417/32415803605926610717906209990215925889160576452720547724429585547651493040781333832362975891811042370960101005436793459936444123023843707180863927766463725709098731479142615931281199106042945274354413006295825020803450451201071934075872286636802447240751072096022370680513420925571815159954608073162324835009657121509219599603167122882578517323407726312472428238351392267283960361495336307712\n",
      "1.41421356237310\n"
     ]
    }
   ],
   "source": [
    "s = 1\n",
    "for i in range(10):\n",
    "    s = (s + 2/s) / 2\n",
    "print s\n",
    "print s.numerical_approx()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and the same computation with floating point numbers:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.41421356237309\n"
     ]
    }
   ],
   "source": [
    "s = 1.0\n",
    "for i in range(10):\n",
    "    s = (s + 2.0 / s) / 2.0\n",
    "print s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What can you say? What are these loops computing?\n",
    "\n",
    "**Exercise:** Let us consider go back to the function\n",
    "$f_{7/2}(x) = {7/2} x (1 - x)$ from the interval \\[0,1\\] to itself.\n",
    "Starting from $x_0 = 5/7$ as a rational number, compute 100 iterations\n",
    "of the map $f$. Print the result as a rational number and get an\n",
    "approximation using the class function `numerical_approx`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# enter your commands here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Redo the same iteration starting with the floating point number\n",
    "`5.0 / 7.0` instead. Print the result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# enter your commands here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What do you conclude?\n",
    "\n",
    "## Symbolic expressions\n",
    "\n",
    "A symbolic expression is created anytime you invoke a symbolic function\n",
    "on exact input. For example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sqrt(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** What are the parent and type of the two above examples?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# enter your commands here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As for integers and rationals, you can use `numerical_approx` to obtain\n",
    "an approximation of your number"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print pi.numerical_approx()\n",
    "print sqrt(2).numerical_approx()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Iterating a map with symbolic expressions will give you more complicated\n",
    "expressions.\n",
    "\n",
    "**Exercise:** Startinf from `x_0 = sqrt(2) - 1` as a symbolic expression\n",
    "apply 10 times the map $x \\maspto 4 x (1 - x)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# enter your commands here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How many characters are there in this expression?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# enter your commands here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Symbolic expressions are useful to manipulate expression trees and apply\n",
    "simplification rules. However, most of the time this is *not* what you\n",
    "want to use.\n",
    "\n",
    "## Algebraic numbers\n",
    "\n",
    "You might want to perform exact computations on real numbers but\n",
    "integers and rationals are not enough. A field in between the rationals\n",
    "and the set of real numbers is the set of algebraic numbers. In Sage it\n",
    "is called `AA`. Elements of `AA` might be recognized because when not an\n",
    "exact rational they appear with a question mark at their right end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "AA(2)\n",
    "AA(2).sqrt()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** Check that the two numbers above are indeed elements of\n",
    "`AA`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# enter your commands here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can compare elements"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = AA(2).sqrt()\n",
    "b = AA(3).sqrt()\n",
    "213 * a < 174 * b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Taking square roots is not the only way to build elements from `AA`. The\n",
    "most universal way is to construct roots of polynomials (with\n",
    "coefficients in `QQ` or `AA`)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = polygen(QQ)\n",
    "p1 = x^3 - 3*x^2 + x - 1\n",
    "p2 = x^3 - x - 1\n",
    "r1 = p1.roots(AA)\n",
    "r2 = p2.roots(AA)\n",
    "print r1\n",
    "print r2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = r1[0][0]\n",
    "b = r2[0][0]\n",
    "y = polygen(AA)\n",
    "p = a * x^3 - b * x^2 + x - a*b\n",
    "p.roots(AA)[0][0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** Construct the number `(2/5)^(1/3)` as an element of `AA`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# enter your commands here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The advantage of `AA` is that it is very flexible. On the other hand it\n",
    "might be slow (even very slow). You can have faster numbers using number\n",
    "fields that are intermediate between rationals and real numbers. To\n",
    "construct a number field the basic syntax is as follows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = polygen(QQ)\n",
    "K = NumberField(x^3 - 2, 'cbrt2', embedding=AA(2)^(1/3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that you have to explicitely embbed the number field in `AA` in\n",
    "order for comparison to work properly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = K.gen()    # generator of the number field\n",
    "print cbrt2 > 1\n",
    "print cbrt2 < 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "Authors  \n",
    "-   Vincent Delecroix\n",
    "\n",
    "License  \n",
    "CC BY-SA 3.0"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 8.2",
   "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
}