{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# COMPIII - 2° Semestre de 2023\n",
    "# Nelson Kuhl - IME - USP (nmkuhl@usp.br)\n",
    "## Exemplos com Aritmética de Ponto Flutuante"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np       # sempre\n",
    "import sys               # útil para informações do sistema\n",
    "from sigfig import round # para testes com precisão escolhida"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Para especificações relacionadas a ponto flutuante, veja o tópico **sys.float_info** na página da biblioteca [sys](https://docs.python.org/3/library/sys.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alguns exemplos de arredondamento:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pi = np.pi            #Pi com a precisão do python\n",
    "pi_13 = round(pi, 13) # Pi com 13 algarismos\n",
    "pi_4 = round(pi, 4)   # Pi com 4 algarismos\n",
    "pi_2 = round(pi, 2)   # Pi com 2 algarismos\n",
    "\n",
    "print(\"Arredondamentos do número pi:\")\n",
    "print()\n",
    "print(\"precisão do python: \", pi)\n",
    "print(\"13 algarismos: \", pi_13, \"(o mesmo que 12)\")\n",
    "print(\"4 algarismos: \", pi_4)\n",
    "print(\"2 algarismos: \", pi_2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Frações decimais podem não ser exatas na base 2. Por exemplo,\n",
    "$$\n",
    "    (0.1)_{10} = (0.00011001100110011\\dots)_2\n",
    "$$\n",
    "Portanto, em uma máquina que usa a base 2, a fração 1/10 é arredondada. Vejamos alguns exemplos:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Números com representação exta na base 2:\")\n",
    "print(\"Número:\", 1.0, \"\\t Fração:\", 1.0.as_integer_ratio(), \"\\t Alocado na memória:\", format(1.0, '.18g') )\n",
    "print(\"Número:\", 0.5, \"\\t Fração:\", 0.5.as_integer_ratio(), \"\\t Alocado na memória:\", format(0.5, '.18g') )\n",
    "print(\"Número:\", 4.0, \"\\t Fração:\", 4.0.as_integer_ratio(), \"\\t Alocado na memória:\", format(4.0, '.18g') )\n",
    "print(\"Número:\", 0.125, \"\\t Fração:\", 0.125.as_integer_ratio(), \"\\t Alocado na memória:\", format(0.125, '.18g') )\n",
    "print(\"Número:\", 0.25, \"\\t Fração:\", 0.25.as_integer_ratio(), \"\\t Alocado na memória:\", format(0.25, '.18g') )\n",
    "print()\n",
    "\n",
    "print(\"Números com representação não finita na base 2:\")\n",
    "print(\"Número:\", 0.1, \" \\t Fração:\", 0.1.as_integer_ratio(), \"  Alocado na memória:\", format(0.1, '.18g') )\n",
    "print(\"Número:\", 0.05, \"\\t Fração:\", 0.05.as_integer_ratio(), \"  Alocado na memória:\", format(0.05, '.18g') )\n",
    "print(\"Número:\", 0.4, \" \\t Fração:\", 0.4.as_integer_ratio(), \"   Alocado na memória:\", format(0.4, '.18g') )\n",
    "print(\"Número:\", 0.0125, \"\\t Fração:\", 0.0125.as_integer_ratio(), \" Alocado na memória:\", format(0.0125, '.18g') )\n",
    "print(\"Número:\", 0.025, \"\\t Fração:\", 0.025.as_integer_ratio(), \" Alocado na memória:\", format(0.025, '.18g') )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Representação de números reais no Python:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"base: \", sys.float_info.radix)\n",
    "print(\"número de dígitos para a mantissa:\", sys.float_info.mant_dig, \"; na base 2 significa ...\")\n",
    "print(\"expoente máximo: \", sys.float_info.max_exp - 1)\n",
    "print(\"expoente mínimo: \", sys.float_info.min_exp - 1)\n",
    "print(\"máximo inteiro e tal que 10**e está em F: \", sys.float_info.max_10_exp)\n",
    "print(\"mínimo inteiro e tal que 10**e está em F: \", sys.float_info.min_10_exp)\n",
    "print(\"máximo float positivo representável: \", sys.float_info.max)\n",
    "print(\"mínimo float positivo representável: \", sys.float_info.min)\n",
    "print(\"número de dígitos decimais representados corretamente: \", sys.float_info.dig)\n",
    "print(\"epsilon de máquina: \", sys.float_info.epsilon)\n",
    "print(\"modo de arredondamento: \", sys.float_info.rounds, \"; (obs: arredondamento para o mais próximo é o modo 1)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Algumas consequências de arredondamento:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Somar 0.1 dez vezes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = 0.0\n",
    "for i in range(10):\n",
    "    x += 0.1\n",
    "    \n",
    "print(\"x = 1.0? \", x == 1.0)\n",
    "print(\"valor calculado: \", x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A adição e subtração não são associativas (mas a adição é comutativa). \n",
    "Exemplo com 3 algarismos significativos\n",
    "\n",
    "(ver Humes et al, *Noções de Cálculo Numérico*, McGraw-Hill, 1984)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = round(4.26 + 9.24, 3)\n",
    "b = round(a + 5.04, 3)\n",
    "print('(4.26 + 9.24) + 5.04 \"=\"', a, '+ 5.04 \"=\"', b)\n",
    "\n",
    "a = round(9.24 + 5.04, 3)\n",
    "b = round(4.26 + a, 3)\n",
    "print('4.26 +(9.24 + 5.04) \"=\"', '4.26 +', a,'\"=\"', b)\n",
    "print()\n",
    "\n",
    "a = round(4210 - 4.99, 3)\n",
    "b = round(a - 0.02, 3)\n",
    "print('(4210 - 4.99) - 0.02 \"=\"', a, '- 0.02 \"=\"', b)\n",
    "\n",
    "a = round(4.99 + 0.02, 3)\n",
    "b = round(4210 - a, 3)\n",
    "print('4210 - (4.99 + 0.02) \"=\"', '4210 -', a, '\"=\"', b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Idem para multiplicação/divisão.\n",
    "Exemplo com 3 algarismos significativos (Humes et al ...)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = round(0.123/7.97, 3)\n",
    "b = round(a * 84.9, 3)\n",
    "print('(0.123 / 7.97) * 84.9 \"=\"', a, '* 84.9 \"=\"', b)\n",
    "\n",
    "a = round(0.123 * 84.9, 3)\n",
    "b = round(a / 7.97, 3)\n",
    "print('(0.123 * 84.9) / 7.97 \"=\"', a, '/ 7.97 \"=\"', b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A propriedade distributiva também não é necessariamente válida.\n",
    "Exemplo com 3 algariemos significativos (Humes et al ...)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = round(4.99 + 0.02, 3)\n",
    "b = round(15.9 * a, 3)\n",
    "print('15.9 * (4.99 + 0.02) \"=\" 15.9 *', a, '\"=\"', b)\n",
    "\n",
    "a = round(15.9 * 4.99, 3)\n",
    "b = round(15.9 * 0.02, 3)\n",
    "c = round(a + b, 3)\n",
    "print('(15.9 * 4.99) + (15.9 * 0.02) \"=\"', a, '+', b, \"=\", c)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Equação quadrática:\n",
    "$$\n",
    "    ax^2 + bx + c = 0, \\quad, a \\ne 0.\n",
    "$$\n",
    "Raízes:\n",
    "$$\n",
    "    x_{1,2} = \\frac{-b \\pm\\sqrt{b^2 - 4ac}}{2a} .\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Exemplo\n",
    "\n",
    "(George E. Forsythe, *What is a satisfactory quadratic equation solver?* In In B. Dejon and P. Henrici,\n",
    "editors, Constructive Aspects of the Fundamental Theorem of Algebra, pages 53–61. Wiley-Interscience, London, 1969)\n",
    "\n",
    "\n",
    "$$\n",
    "    \\frac{1}{2}x² - 28x + \\frac{1}{2} = 0.\n",
    "$$\n",
    "Usando 5 algarismos significativos:\n",
    "$$\n",
    "    x_1 = 28 + \\sqrt{783} \"=\" 28 + 27.982 = 55.982\\\\\n",
    "    \\ \\ \\ x_2 = 28 - \\sqrt{783} \"=\" 28 - 27.982 = 0.018000\n",
    "$$\n",
    "\n",
    "O único arredondamento que ocorreu foi no cálculo de $\\sqrt{783}$. Note que apesar de usarmos 5 algarismos significativos, o valor calculado de $x_2$ tem apenas dois algarismos (CANCELAMENTO!)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "O cancelamento acima pode ser evitado lembrando-se que as soluções $x_1$ e $x_2$ de uma equação quadrática satisfazem\n",
    "\n",
    "$$\n",
    "    x_1*x_2 = \\frac{c}{a}.\n",
    "$$\n",
    "\n",
    "Logo, podemos evitar o cancelamento em $x_2$ acima calculando (com 5 algarismos significativos)\n",
    "\n",
    "$$\n",
    "    x_2 = \\frac{1}{x_1} = \\frac{1}{55.982} \"=\" 0.017863\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Observe que as raízes com uma precisão maior são iguais a (assumiremos como valores extatos):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_1 = 28 + np.sqrt(783)\n",
    "x_2 = 1/x_1\n",
    "print()\n",
    "print('x_1 =', x_1)\n",
    "print('x_2 =', x_2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Arredondando-se esses resultados para 5 algarismos significativos, obtemos os valores calculados anteriormente com a fórmula que evita o cancelamento."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "O erro relativo entre x_2 e o valor obtido com cancelamento é"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "erro_ruim = abs((x_2 - 0.018) / x_2)\n",
    "print(\"erro_ruim =\", format(erro_ruim, '.1e'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Os erros relativos quando usamos a fórmula melhor são:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "erro_x1 = abs((55.982 - x_1) / x_1)\n",
    "erro_x2 = abs((0.017863 - x_2) / x_2)\n",
    "print(\"erro em x1: \", format(erro_x1, '.1e'))\n",
    "print(\"erro em x2: \", format(erro_x2, '.1e'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}