13-computacao-numerica.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Computação numérica\n",
    "=====================\n",
    "\n",
    "Números e números\n",
    "-------------------\n",
    "\n",
    "Já vimos que o Python reconhece diferentes *tipos* de números:\n",
    "\n",
    "- números em ponto flutuante (`float`), como 3.14\n",
    "\n",
    "- inteiros (`int`) como 42\n",
    "\n",
    "- números complexos (`complex`), como 3.14 + 1*j\n",
    "\n",
    "### Limitações dos tipos de números\n",
    "\n",
    "#### Limitações de `int`s\n",
    "\n",
    "A matemática fornece o conjunto infinito dos números naturais $\\mathbb{ℕ} = \\{1, 2, 3,\\dots \\}$. Como o computador opera de modo *finito*, é impossível que ele represente todos esses números. Em vez disso, apenas um pequeno subconjunto de números é representado.\n",
    "\n",
    "O tipo `int` (geralmente [3]) representa números entre -2147483648 e +2147483647, e corresponde a 4 *bytes* (isto é 4\\*8 *bits* e 2<sup>32</sup> = 4294967296, que é o intervalo de -2147483648 a +2147483647).\n",
    "\n",
    "Você pode imaginar que o hardware usa uma tabela como esta para codificar inteiros usando bits (suponha, por simplicidade, que usemos apenas 8 *bits*):\n",
    "\n",
    "| número natural |representação binária|\n",
    "|---------------:|-------------------:|\n",
    "|               0|            00000000|\n",
    "|               1|            00000001|\n",
    "|               2|            00000010|\n",
    "|               3|            00000011|\n",
    "|               4|            00000100|\n",
    "|               5|            00000101|\n",
    "|               ⋮|                   ⋮|\n",
    "|             254|            11111110|\n",
    "|             255|            11111111|\n",
    "\n",
    "Usando 8 *bits*, podemos representar 256 números naturais (por exemplo, de 0 a 255), pois temos 2<sup>8</sup> = 256 formas diferentes de combinar oito 0s e 1s.\n",
    "\n",
    "Poderíamos também usar uma tabela ligeiramente diferente para descrever 256 números inteiros que variem, por exemplo, de -127 a +128.\n",
    "\n",
    "Assim é, *em princípio*, como os números inteiros são representados no computador. Dependendo do número de *bytes* utilizados, apenas números inteiros entre um valor mínimo e um máximo podem ser representados. Nos hardwares de hoje, é comum usar 4 ou 8 *bytes* para representar um número inteiro, o que leva exatamente aos valores mínimo e máximo de -2147483648 e +2147483647, como mostrado acima para 4 *bytes*, e +9223372036854775807 como o inteiro máximo para 8 *bytes* (isto é ≈ 9.2⋅10<sup>18</sup>).\n",
    "\n",
    "#### Limitações de `float`s\n",
    "\n",
    "Os números em ponto flutuante em um computador não são os mesmos que os números em ponto flutuante da matemática. (Isto é exatamente o mesmo que dize que os números inteiros (matemáticos) não são os mesmos em um computador: apenas um *subconjunto* do conjunto infinito dos números inteiros pode ser representado pelo tipo `int`, como mostrado anteriormente). Então, como os números em ponto flutuante são representados no computador?\n",
    "\n",
    "- Qualquer número real $x$ pode ser escrito como\n",
    "     $$x  = a \\cdot 10^b,$$\n",
    "\n",
    "onde $a$ é a mantissa e $b$ o expoente.\n",
    "      \n",
    "- Exemplos: \n",
    "\n",
    "| $x$                               | $a$     | $b$|\n",
    "|-----------------------------------|---------|----|\n",
    "| 123.45 = 1.23456 ⋅ 10<sup>2</sup> | 1.23456 |  2 |\n",
    "| 1000000 = 1.0 ⋅ 10<sup>6</sup>    | 1.00000 |  6 |\n",
    "| 0.0000024 = 2.4 ⋅ 10<sup>-6</sup> | 2.40000 | -6 |\n",
    "\n",
    "- Portanto, podemos usar 2 inteiros para codificar um número em ponto flutuante!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Seguindo (aproximadamente) o padrão IEEE-754, usa-se 8 *bytes* para um `float` $x$: estes 64 bits são divididos entre\n",
    "\n",
    "     - 10 bits para o expoente $b$ e\n",
    "\n",
    "     - 54 bits para a mantissa $a$.\n",
    "\n",
    "Isto resulta em\n",
    "\n",
    "- maior número em ponto flutuante possível ≈ 10 <sup>308</sup> (medida de qualidade para *b*)\n",
    "\n",
    "- menor número em ponto flutuante possível (positivo) ≈ 10 <sup>-308</sup> (medida de qualidade para $b$)\n",
    "\n",
    "- distância entre 1.0 e o número imediatamente superior ≈ 10<sup>-16</sup> (medida de qualidade para $a$)\n",
    "\n",
    "Observe que, *a priori*, é assim como os números em ponto flutuante são armazenados (na realidade, é um pouco mais complicado).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Limitações dos números complexos\n",
    "\n",
    "O tipo `complex` tem essencialmente as mesmas limitações que o tipo `float` porque um número complexo consiste de dois números `float`: um representa a parte real, o outro a parte imaginária.\n",
    "\n",
    "#### ... são esses tipos de números de valor prático?\n",
    "\n",
    "Na prática cotidiana, geralmente não encontramos números que excedam 10<sup>300</sup> (este é um número com 300 zeros!). Portanto, os números em ponto flutuante cobrem o intervalo de números de que geralmente precisamos.\n",
    "\n",
    "No entanto, você precisa estar ciente de que, em computação científica, números pequenos e grandes são usados, os quais podem (geralmente em resultados intermediários) exceder o intervalo de números em ponto flutuante.\n",
    "\n",
    "- Imagine, por exemplo, que tenhamos de calcular a quarta potência da constante ℏ = 1.0545716 ⋅ 10 <sup>-34</sup> *kg. m<sup>2</sup>/s*: \n",
    "\n",
    "- ℏ <sup>4</sup> = 1.2368136958909421 ⋅ 10 <sup>-136</sup> *kg<sup>4</sup>.m<sup>8</sup>/s<sup>4</sup>*, que está a \"meio caminho\" de nosso menor `float` positivo representável da ordem de 10<sup>-308</sup>.\n",
    "\n",
    "Se houver algum perigo de que possamos exceder os limites dos números em ponto flutuante, teremos que *escalonar* nossas equações de modo que (idealmente) todos os números sejam da ordem de 1. Escalonar nossas equações para que todos os números relevantes sejam próximos de 1 (normalização) também é útil para a depuração de código: se números muito maiores ou menores do que 1 surgirem, isso pode ser uma indicação de erro.\n",
    "\n",
    "### Usando números de ponto flutuante (sem o devido cuidado)\n",
    "\n",
    "Já sabemos que precisamos ter cuidado para que nossos valores em ponto flutuante não excedam os limites que podem ser representados no computador (*underflow* e *overflow*). Há outra complicação devido à forma como os números em ponto flutuante têm de ser representados internamente: nem todos eles podem ser representados exatamente no computador. O número 1.0 pode ser representado exatamente, mas os números 0.1, 0.2 e 0.3 não:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'1.00000000000000000000'"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "'%.20f' % 1.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'0.10000000000000000555'"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "'%.20f' % 0.1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'0.20000000000000001110'"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "'%.20f' % 0.2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'0.29999999999999998890'"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "'%.20f' % 0.3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Em vez disso, o número em ponto flutuante \"mais próximo\" do número real é escolhido.\n",
    "\n",
    "Isso pode causar problemas. Suponha que precisemos de um laço onde `x` tome os valores 0,1, 0,2, 0,3, ..., 0,9, 1,0. Podemos ser tentados a escrever algo como:\n",
    "\n",
    "```python\n",
    "x = 0.0\n",
    "while not x == 1.0:\n",
    "    x = x + 0.1\n",
    "    print (\" x = %19.17f\" % ( x ))\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Entretanto, este laço nunca terminará. Aqui estão as primeiras 11 linhas de saída do programa:\n",
    "\n",
    "    x=0.10000000000000001\n",
    "    x=0.20000000000000001\n",
    "    x=0.30000000000000004\n",
    "    x=0.40000000000000002\n",
    "    x=                0.5\n",
    "    x=0.59999999999999998\n",
    "    x=0.69999999999999996\n",
    "    x=0.79999999999999993\n",
    "    x=0.89999999999999991\n",
    "    x=0.99999999999999989\n",
    "    x=1.09999999999999987"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Como a variável `x` nunca assume exatamente o valor 1.0, o laço `while` continuará infinitamente.\n",
    "\n",
    "Portanto: **Nunca compare dois números em ponto flutuante por igualdade!**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Usando números de ponto flutuante sem o devido cuidado - 1\n",
    "\n",
    "Há várias alternativas para resolver este problema. Por exemplo, podemos comparar a distância entre dois números em ponto flutuante:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " x =0.10000000000000001\n",
      " x =0.20000000000000001\n",
      " x =0.30000000000000004\n",
      " x =0.40000000000000002\n",
      " x =0.50000000000000000\n",
      " x =0.59999999999999998\n",
      " x =0.69999999999999996\n",
      " x =0.79999999999999993\n",
      " x =0.89999999999999991\n",
      " x =0.99999999999999989\n"
     ]
    }
   ],
   "source": [
    "x = 0.0\n",
    "while abs(x - 1.0) > 1e-8:\n",
    "    x = x + 0.1\n",
    "    print ( \" x =%19.17f\" % ( x ))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Usando números de ponto flutuante sem o devido cuidado - 2\n",
    "\n",
    "Alternativamente, podemos iterar sobre uma sequência de inteiros e calcular o número em ponto flutuante a partir do inteiro:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " x =0.10000000000000001\n",
      " x =0.20000000000000001\n",
      " x =0.30000000000000004\n",
      " x =0.40000000000000002\n",
      " x =0.50000000000000000\n",
      " x =0.60000000000000009\n",
      " x =0.70000000000000007\n",
      " x =0.80000000000000004\n",
      " x =0.90000000000000002\n",
      " x =1.00000000000000000\n"
     ]
    }
   ],
   "source": [
    "for i in range (1 , 11):\n",
    "    x = i * 0.1\n",
    "    print(\" x =%19.17f\" % ( x ))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "x=0.10000000000000001\n",
    "x=0.20000000000000001\n",
    "x=0.30000000000000004\n",
    "x=0.40000000000000002\n",
    "x=                0.5\n",
    "x=0.60000000000000009\n",
    "x=0.70000000000000007\n",
    "x=0.80000000000000004\n",
    "x=0.90000000000000002\n",
    "x=                  1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Se compararmos esta saída com aquela do programa do caso 1, veremos que os números em ponto flutuante diferem. Isto significa que - em cálculo numérico - não é verdade que 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 = 1.0."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cálculo simbólico\n",
    "\n",
    "Usando o pacote SymPy, temos precisão arbitrária. Usando `sympy.Rational`, podemos definir a fração 1/10 exatamente com o cálculo simbólico. Adicionando a fração 10 vezes, seremos levados exatamente ao valor 1, conforme demonstrado por este script:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Atual x=1/10 = 0.1 \n",
      " Atingido x=1/10 \n",
      "Atual x= 1/5 = 0.2 \n",
      " Atingido x=1/5 \n",
      "Atual x=3/10 = 0.3 \n",
      " Atingido x=3/10 \n",
      "Atual x= 2/5 = 0.4 \n",
      " Atingido x=2/5 \n",
      "Atual x= 1/2 = 0.5 \n",
      " Atingido x=1/2 \n",
      "Atual x= 3/5 = 0.6 \n",
      " Atingido x=3/5 \n",
      "Atual x=7/10 = 0.7 \n",
      " Atingido x=7/10 \n",
      "Atual x= 4/5 = 0.8 \n",
      " Atingido x=4/5 \n",
      "Atual x=9/10 = 0.9 \n",
      " Atingido x=9/10 \n",
      "Atual x=   1 = 1.0 \n",
      " Atingido x=1 \n"
     ]
    }
   ],
   "source": [
    "from sympy import Rational\n",
    "dx = Rational (1 ,10)\n",
    "x = 0\n",
    "while x != 1.0:\n",
    "    x = x + dx\n",
    "    print(\"Atual x=%4s = %3.1f \" % (x , x . evalf ()))\n",
    "    print(\" Atingido x=%s \" % x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "No entanto, esse cálculo simbólico é muito mais lento, já que é feito através de software em vez de operações de ponto flutuante baseadas em CPU. O próximo programa aproxima os desempenhos relativos:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "laco usando float dx:\n",
      " desvio eh -1.88483681995422e-08\n",
      "laco float n=100000 leva 0.00944 segundos\n",
      "laco usando sympy simbolico dx:\n",
      " desvio eh                0\n",
      "laco sympy n = 100000 leva 1.80431 segundos\n",
      "laco simbolico eh um fator 191.1 mais lento.\n"
     ]
    }
   ],
   "source": [
    "from sympy import Rational\n",
    "dx_simbolico = Rational (1 ,10)\n",
    "dx = 0.1\n",
    "\n",
    "def loop_sympy (n):\n",
    "    x = 0\n",
    "    for i in range(n):\n",
    "        x = x + dx_simbolico\n",
    "    return x\n",
    "\n",
    "def loop_float(n):\n",
    "    x =0\n",
    "    for i in range(n):\n",
    "        x = x + dx\n",
    "    return x\n",
    "\n",
    "def medir_tempo(f, n):\n",
    "    import time\n",
    "    tempo_inicial = time.time()\n",
    "    resultado = f(n)\n",
    "    tempo_final = time.time()\n",
    "    print(\" desvio eh %16.15g\" % ( n * dx_simbolico - resultado ))\n",
    "    return tempo_final - tempo_inicial\n",
    "\n",
    "n = 100000\n",
    "print(\"laço usando float dx:\")\n",
    "tempo_float = medir_tempo(loop_float, n)\n",
    "print(\"laço float n=%d leva %6.5f segundos\" % (n, tempo_float))\n",
    "print(\"laço usando sympy simbólico dx:\")\n",
    "tempo_sympy = medir_tempo(loop_sympy, n)\n",
    "print(\"laço sympy n =% d leva %6.5f segundos\" % (n , tempo_sympy ))\n",
    "print(\"laço simbólico é um fator %.1f mais lento.\" % ( tempo_sympy / tempo_float ))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Este é, naturalmente, um exemplo artificial: adicionamos o código simbólico para demonstrar que esses erros de arredondamento originam-se da representação aproximada de números em ponto flutuante no hardware (e, portanto, linguagens de programação). Podemos, em princípio, evitar essas complicações através de cálculo com expressões simbólicas, mas, na prática, isso é muito lento. [4]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Resumo\n",
    "\n",
    "Em resumo, aprendemos que\n",
    "\n",
    "- números em ponto flutuante e inteiros utilizados em computação numérica geralmente são bastante diferentes dos \"números matemáticos\" (os cálculos simbólicos são exatos e usam os \"números matemáticos\"):\n",
    "\n",
    "    - há um número máximo e um número mínimo que podem  ser representados (tanto para números inteiros como para ponto flutuante)\n",
    "\n",
    "    - dentro desse intervalo, nem todo número em ponto flutuante pode ser representado no computador.\n",
    "\n",
    "- lidamos com essa limitação:\n",
    "\n",
    "    - nunca comparando dois números em ponto flutuante por igualdade (em vez disso, calculamos o valor absoluto da diferença)\n",
    "\n",
    "    - uso de algoritmos que são *estáveis* (isto significa que pequenos desvios de números corretos podem ser corrigidos pelo algoritmo. Ainda não mostramos nenhum exemplo desse tipo neste documento.)\n",
    "\n",
    "- Observe que há muito mais a dizer sobre métodos numéricos e técnicas algorítmicas para tornar a computação numérica tão precisa quanto possível, mas isto está fora do escopo desta seção."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercício: laço finito ou infinito\n",
    "\n",
    "1.  O que o seguinte fragmento de código computa? O laço terminará? Por quê?\n",
    "\n",
    "```python\n",
    "eps = 1.0\n",
    "while 1.0 + eps > 1.0:\n",
    "    eps = eps / 2.0\n",
    "print(eps)\n",
    "```"
   ]
  }
 ],
 "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}