Operações elementares no computador

......

Transcript

1. LNCC UFRJ Operações elementares no computador Prof. Paulo R. G. Bordoni
2. Flores lindas! Veremos a seguir que nem tudo são flores nos jardins da Análise Numérica.
3. A conversão base 10  base 2 produz ”bízimas” periódicas. Bízimas?
4. O algoritmo para obtenção da representação na base 2 de um número entre 0 e 1, escrito na base 10, pode ser descrito por: Multiplique o número dado sucessivamente por 2 e vá coletando os “vai 0” ou “vai 1”.
5. Nesse exemplo eu mostro como converter 0.3 para a base 2. 0.3 = 0. 0.3 ∗ 2 0.6 0 0.6 ∗ 2 1.2 1 0.2 ∗ 2 0.4 0 0.4 ∗ 2 0.8 0 0.8 ∗ 2 1.6 1 0.6 ∗ 2 1.2 1 1 0 0 . . .
6. Para números racionais, o algoritmo de conversão só termina em um número finito de passos quando o número é uma potência negativa de 2 ou uma soma de potências negativas de 2. Em todos os outros casos, a representação do número racional, na base 2, é uma seqüência periódica – uma “bí”zima periódica
7. Como as representações IEEE 754 possuem apenas uma quantidade finita de bits na fração, precisaremos truncá-los. O que introduzirá um erro na representação IEEE do número!.
8. Se as bízimas forem poucas não haverá muito problema, Mestres! Vou mostrar a você que, infelizmente, serão muitas bízimas. Na verdade a grande maioria!
9. 80% dos números decimais com apenas 1 dígito após a vírgula geram bízimas periódicas na base 2: Surfista! Construa um programa Python que recebe um número entre zero e um, escrito na base 10, com apenas uma casa decimal após a vírgula, e o converta para a base dois. 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
10. 0.00 0.01 ... 0.24 0.25 0.26 ... 0.49 0.50 0.51 ... 0.74 0.75 0.76 ... 0.99 Confira que 96% dos números decimais com apenas 2 dígito após a vírgula geram dízimas periódicas na base 2:
11. Repita isso para números com 3, 4, 5, 6, 7 casas decimais após a vírgula ... Depois dê uma olhada na próxima transparência!
12. 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 80% 0.00 0.01 ... 0.24 0.25 0.26 ... 0.49 0.50 0.51 ... 0.74 0.75 0.76 ... 0.99 96% 0.000 0.001 ... 0.124 0.125 0.126 ... 0.499 0.500 0.501 ... 0.874 0.875 0.876 ... 0.999 99.2% 0.0000000 0.0000001 ... 0.xxxxxxx 0.xxxxxxx 0.xxxxxxx ... 0.4999999 0.5000000 0.5000001 ... 0.xxxxxxx 0.xxxxxxx 0.xxxxxxx ... 0.9999999 99.9936% A grande maioria é bízima periódica. No padrão Single: 99,99 %.
13. A consequência poderá ser funesta ao calcularmos valores de funções ... São apenas 274.877, de números corretos do total de 4.294.967.296, no caso dos Single. Não chega a 6 milésimos %.
14. Não minha Loirinha querida! É que no padrão Single, temos 24 = 23 + 1 bits de precisão. O Mestre parou em sete dígitos após a vírgula. Por quê? Sete é algum número mágico Mestra?!
15. É só responder à pergunta: “24 casas binárias correspondem a quantas casas decimais?” Em matemática: qual o valor de k para o qual 10 = 224 Resolvendo a equação acima encontramos = 24 ∗ 10 (2) = 24 ∗ 0.301030 = 7.224719 ≅ 7 Continuo não entendendo!
16. Mestra, como eu faço para buscar informações desse tipo na Numpy? Pela “Search page”. Marquei o endereço ali em cima.
17. Qual é a informação que você quer buscar Loirinha? Precisa especificar com clareza! Ora Mestra, sobre ponto flutuante!
18. Fui buscar por “floating point information”. Veio muita coisa!
19. Pois é... Agora vocês precisam separar o joio do trigo. Eliminei algumas que não tem nada a ver!
20. Não risquei as “Data type ...” porque float é um tipo de dado. Explique porque você deixou as outras, Surfista. Acho importante esclarecer seu critério de busca!
21. Vejam, achamos! Inclusive o “helper” marca em amarelo para ajudar. Vou clicar em finfo.
22. Informações sobre os tipos de ponto flutuante e dois “links” importantes.
23. É mesmo Surfista. Inclusive tem tudo que vimos em “Números reais no computador ou IEEE 754”
24. Maravilha, mas não entendi direito como usar essa finfo.
25. O Surfista tem razão, Mestre. Dê uma mãozinha, pois o “help” não está muito claro! É só instanciar o parâmetro dtype com o tipo que desejamos obter as informações. Por exemplo: dtype = numpy.float32.
26. Vejam meu programinha:
27. O resultado: Inclusive notem que o tipo float16 é quase igual ao do Toy System.
28. Façam um programa que calcula a ULP(x) para x dado. Como exercício, procurem as informações nos “links” indicados como importantes pelo Surfista: MachAr e iinfo.
29. O Eps já sabemos onde obter. Só falta descobrir como calcular exp(). Essa é fácil, = ∗ 2exp().
30. Dei uma busca em “Floating point routines”. Como o resultado é quase igual à nossa busca anterior, já eliminei as 3 primeiras.
31. A “Data type routines” também risco – foi a anterior. As duas seguintes não tem nada a ver. Vamos tentar “Mathematical functions”.
32. Acho que já achamos, Surfista! Veja onde marquei.
33. Isso mesmo Loirinha, ela recebe e retorna a () e () na base 2, que usaremos para calcular (). Vou lá ver!
34. Os detalhes dela!
35. Nosso programa, Mestre. Caprichamos! Inclusive recalculamos o valor de x para conferir!
36. Executamos para diversos valores de x e do tipo de x.
37. Mestres, vocês encerraram a aula passada falando sobre número de condicionamento e fatores de propagação de perturbações ∆. Mas não fizeram nenhum programa.
38. Sim! Além disso, qual a utilidade prática do conceito?
39. Bem, minha querida, escolhendo ∆ = (0 ), você terá uma estimativa do erro propagado pela f ao aproximar 0 por (0 ): (∆ )(0 ) ≅ , 0 (0 ) (0 ) (∆ )(0 ) 0 (0 )
40. Calculei os para algumas funções: () (, ) (, ) ||−1 n () () || () 1 1 () () | | | | () | | | ()| () 2() | ( + ())| () 1 1 − 2 | | 1 − 2
41. E devolverá os valores de , e , . Ele estimará também os erros absoluto e relativo propagados por f em a partir de () e Eps. Finalmente exibirá os gráficos da f e da reta tangente ao gráfico de f por numa pequena vizinhança de e a propagação de uma perturbação ∆. Nosso programa receberá o valor de , a expressão () da função f e a expressão de sua derivada ().
42. Executamos o programa para = 3 e = 247.4, em precisão simples (float32). Atenção!
43. Escolhemos próximo de /2, ponto de descontinuidade de tan(). Vejam os valores dos números de condicionamento. Notem que perdemos quase 2 dígitos na precisão. Atenção!
44. Atenção! Executamos o programa para = arcsin() e = 0.996. Perdemos praticamente 1 dígito.
45. Antes de mostrar o código do programa, chamo a atenção para o fato que, na grande maioria dos exemplos, as escalas vertical e horizontal diferirem em várias ordens de grandeza. Confiram! Não há como fazer diferente. Atenção! Atenção! Atenção!
46. Este é o programa prometido.
47. A continuação.
48. Os sistemas IEEE 754 de ponto flutuante NÃO foram criados para aproximar números reais no computador.
49. Nossa análise até agora mostrou que: Mais de 99,99% dos números ∗ = () do padrão IEEE, Single ou Double apresentam: • Erro absoluto limitado pela (). • Erro relativo limitado pelo “épsilon da máquina”, Eps.
50. = − ∗ < = − ∗ < Insistindo no aspecto fundamental: ∈ ℝ mas ∗ = ∈ ℚ, com: Mas você está focando apenas no aspecto estático das coisas. E a dinâmica dos processos?
51. De fato Professora, a principal razão de existência dos números reside na possibilidade de computar com eles. “C’est la raison d’etre” dos computadores. Como se comportarão as operações algébricas perante tais erros?
52. Mestres, Existem conhecimentos subjacentes ao entrelaçamento operações algébricas ↔ IEEE 754 que precisam ser aprofundados. Sim Filósofo, precisaremos entender o comportamento das operações algébricas no seio do IEEE 754.
53. As operações algébricas são funções definidas para números reais. A adição recebe um par de números reais e devolve um número real. Isto também vale para os racionais e os complexos. + Adição z = x + y ℝ (x,y) y x ℝ ℝ
54. ∗ Multiplicação (x,y) y x ℝ ℝ z = x ∗ y ℝ O mesmo vale para a multiplicação em ℚ, ℝ e ℂ.
55. Vocês já devem ter ouvido falar que ℚ, ℝ e ℂ com as operações de adição + e multiplicação ∙ constituem estruturas algébricas com propriedades os tornam um corpo. Em outras palavras, são grupos abelianos para a adição e a multiplicação (tirando o 0) e que satisfazem também a distributividade da multiplicação sobre a adição.
56. Só para recordar, um grupo abeliano é um conjunto com uma operação ∘∶ × → satisfazendo as propriedades básicas: 1. Comutativa: ∘ = ∘ , , ∈ . 2. Associativa: ( ∘ ) ∘ = ∘ ( ∘ ), , , ∈ . 3. Existe um elemento, ∈ , nominado identidade, para o qual ∘ = 4. Cada ∈ possui seu inverso −1, solução de ∘ −1 = .
57. Quando a operação é a adição, a identidade é o zero, 0. Quando é a multiplicação a identidade é o 1. Na multiplicação de matrizes é a matriz identidade I. Na composição de funções, é a função identidade, : → .
58. Além disso ℚ e ℝ, é são corpos ordenados. Isto é, para quaisquer , , ∈ ℚ, ℝ valem as afirmações: • > ⟹ + > + • > 0 > 0 ⟹ > 0 Aliás, não existe uma relação de ordem natural em ℂ. Não há como dizer que um número complexo é maior que outro.
59. Finalmente vamos mergulhar fundo na continuidade dos números reais. Algo tão importante que nossos Mestres vestirão roupas de gala!
60. O axioma da existência do supremo: “Todo subconjunto (≠ ∅) limitado superiormente de ℝ possui um supremo”. É ele quer possibilita falarmos em continuidade, derivadas, integrais, ... Todo o Cálculo Diferencial e Integral.
61. Toda a “matemática do contínuo”, caros aluninhos. Os processos limite e as infinitas casa decimais em números como , 2, etc.
62. Um subconjunto ⊂ ℝ é limitado superiormente quando: ∃ ∈ ℝ tal que ≥ , ∀ ∈ ℝ. Mestre, você está um gatão! Mas, o que significa “conjunto limitado superiormente”?
63. O supremo s de um conjunto limitado superiormente X ⊂ ℝ, anotado (), é o menor de seus limitantes superiores. Em símbolos: = () ⟺ 1) ∀ ∈ , ≤ 2) ∀ > 0, ∃ ∈ tal que − < ≤
64. • 2 = ∈ ℚ . 2 < 2 • 1 = ∈ ℝ . < 1 • 1 = ∈ ℝ . ≤ 1 • 2 = 1 + 1 2 , 1 + 1 2 + 1 4 , 1 + 1 2 + 1 4 + 1/8, ⋯ • 9 = 0, 1, 2, 3, 4, 5, 6, 7, 8, 9 • 2 = , é í í ê 1 Cruz-credo, Mestra! Dê um exemplo, pelo amor de Deus Veja alguns:
65. Limitantes superiores s A 1ª parte da definição, ∀ ∈ , ≤ , garante que s é um limitante superior de X.
66. − ∈ s A 2ª parte da definição, ∀ > 0, ∃ ∈ tal que − < ≤ , garante que s é o menor de todos os limitantes superiores.
67. Nessa 2ª figura vocês não mostraram os elementos do conjunto X. Por quê? Só para descobrir se você estava atento! Imagine, p/ex., X como as aproximações decimais por falta para 2: = { 1.4, 1.41, 1.414, 1.4142, … }
68. Um ponto a ser destacado a respeito do supremo s de X é que ele pode ser um elemento de X, mas isto não é uma exigência. Quando ∈ , ele é o máximo de X, anotado max()
69. Surfistas e Loirinhas, neste exato momento, interrompam a leitura deste conjunto de transparências e dirijam-se à aba “Leituras adicionais” do site do Mestre para ler os textos indicados a seguir.
70. O primeiro dos textos é:
71. O segundo:
72. E o último:
73. Mestres, os sistemas IEEE de ponto flutuante são constituídos apenas por números racionais. Para eles não vale a propriedade de existência do supremo! Exatamente, Surfista. Esta é a maior perda. Ela é terrível!
74. Mas, então como faremos as coisas de Cálculo no computador? A palavra chave é engenhosidade Loirinha, criatividade e conhecimento.
75. Observem que para obter o oposto, −, e o inverso, 1/, de ∈ ℝ temos que resolver os problemas inversos: • Qual ∈ ℝ satisfaz a equação + = 0. • Qual ∈ ℝ satisfaz a equação ∗ = 1 São as operações de subtração e divisão, − e / que fornecem a solução para esses problemas.
76. Assim, operacionalmente, não basta falar em adição e multiplicação. Precisaremos das quatro operações elementares +, −, ∗ e /. E, também da igualdade =, ≠ , das operações lógicas ( ∧,∨, →, ↔ ) e da relação de ordem >, ≥ .
77. É por esse motivo operacional que os computadores possuem uma unidade de aritmética e lógica na qual essas operações foram implementadas de acordo com os padrões do IEEE 754. Eu nunca tinha me detido em pensar nessas questões, Filósofo!
78. As linguagens C e Fortran padrão usam a unidade de aritmética e lógica. Linguagens como Python, Java e as .net trabalham com máquinas virtuais que obedecem o padrão IEEE 754 em suas implementações. Assim, todas as operações lógicas, comparações e as operações elementares +, −, ∗ , / e já satisfazem o padrão IEEE 754.
79. Mas será que, operacionalmente, a aritmética do IEEE 754 manterá as propriedades de corpo? Não será possível! Mas essa tua pergunta, desencadeará uma série de outras perguntas.
80. Além da perda do contínuo, Loirinha, como os números do IEEE 754 estão limitados à faixa (−, +), perdemos também a propriedade do fechamento. Isto não será tão problemático pois os sistemas IEEE 754 nos informarão com mensagens do tipo NaN (Not A Number).
81. Este exemplo simples, com o Toy System, que obedece o IEEE 754, confirma que a adição não é fechada. Vejam, a soma de = 0.125 e = 1.0, ambos números do Toy, é + = 1.125, que não está no Toy System. +
82. Bem Loirinha, o passo seguinte é obrigar “a ovelha desgarrada” retornar ao lar. Mas e aí, o que eu faço com o resultado? +
83. Em outras palavras, precisamos transformar o resultado num número do tipo IEEE 754. Sim, e esse “modus operandi” se aplica a todos os tipos float do IEEE 754.
84. Não Loirinha, quando somamos dois números, o computador (o IEEE 754) já dá o resultado diretamente! Mas, depois de calcular + , eu sempre vou ter que aproximar o resultado, isto é calcular ( + ) ?
85. IEEE 754 Operações com IEEE 754 ℚ IEEE 754 IEEE 754 A imagem a gravar é esta, para todas as operações elementares e a . ∘ = +,∗, −,/, ∗ = ( ∘ ) = ∘ ∘ (, ) ° *
86. Uma nova operação, ? ° * Na realidade é uma das operações ∘ seguida de fl( ). ° *
87. Inclusive, conforme já adiantamos, o IEEE 754 já devolve informações do tipo +Inf, -Inf ou NaN nos casos excepcionais. Portanto, para nossos programas não travarem, temos que incluir rotinas para tratamento desses tipos de erro.
88. As operações aritméticas elementares e a foram construídas na IEEE 754 de forma a satisfazer o teorema fundamental da representação de ponto flutuante:
89. Se , ∈ ℝ e , e ∘ estão na faixa normal de um dos sistemas IEEE 754, então ∘ = ( ∘ )(1 + ) para algum satisfazendo < . Substitua ∘ por +, −, ∗, / ou .
90. • Erro relativo: ∘ − ∘ | ∘ | < , ∘ ≠ 0 • Erro absoluto: ∘ − ∘ < ( ∘ ) Em outras palavras, são observadas as limitações:
91. Iniciamos com flores lindas! Agora os espinhos. Começaremos calculando o número de condicionamento das operações elementares (Lembrem-se, elas são funções).
92. Mestre qual é o número de condicionamento da adição/subtração? Para simplificar trabalharemos com ↦ = ± , com a > 0, fixo. Então = ±1 de modo que = = 1, e o condicionamento absoluto é ótimo.
93. Já o condicionamento relativo é dado por = | | = = ± ∙ 1 Em particular, na subtração, para = + , com > 0, teremos ↦ = − = − + = −. Logo = + = + 1
94. Portanto, quanto menor o valor de tanto maior será . Por exemplo, para = 1 e = 0.001, teremos = 1 0.001 + 1 = 1001.
95. E, quanto menor o valor de tanto maior será e a ampliação do erro. Portanto o erro relativo será ampliado em 1.000 vezes.
96. Esta situação conduz a um problema clássico, conhecido como cancelamento catastrófico. Não entendi o porque do nome cancelamento catastrófico!
97. f Quando os números são aproximadamente iguais, todos os bits mais significativos da fração são iguais e, na subtração, eles se cancelam, restanto os últimos. Como, com 99 ,99% de certeza, os floats32 apresentam erro na 24ª casa binária da fração, após a normalização, esse erro é então promovido para as primeiras casas. 0 0 0 0 X X
98. Usando o programa feito pelos Mestres poderemos ver o que acontece na prática! Sim, como exemplo usaremos ↦ = − 294. e = 294,001.
99. Agora entendi perfeitamente o nome cancelamento catastrófico! Perdemos 5.4 dígitos dos 6.9 de float32. E o erro relativo é quase 300.000 vezes o Eps.
100. Nossa! Escolhendo float16, ficamos devendo dígitos na precisão relativa. Mas a parte gráfica é idêntica. Entretanto o erro absoluto também é grande.
101. Trabalhando com precisão dupla (float64) conseguimos uma resposta correta, posto que, para = 294.001, = − 294. = 0.00099 ⋯ 76 ≅ 0.001. Dos 15.7 dígitos, perdemos 1/3, ficando com 10.2 – uma boa margem de confiança. Claro que, quanto maior a quantidade de bits na fração, melhor a aproximação. É o motivo fundamental de Python, dos Físicos e Engenheiros usarem precisão dupla.
102. Entendeu agora, Loirinha, o motivo do nome cancelamento catastrófico. Realmente uma catástrofe. Nunca imaginei que pudessem ocorrer erros absurdos assim no computador!
103. Então Loirinha, a multiplicação por a, é dada por = ∗ , Portanto = . Assim temos = = . Mestre qual é o condicionamento da multiplicação?
104. Portanto multiplicar por valores grandes amplifica o erro absoluto. Maior o multiplicador, maior é o erro absoluto no produto Observem que lim →∞ () = + ∞
105. Já o condicionamento relativo é dado por = | | = = ∗ = 1 ∙ Assim, o condicionamento relativo da multiplicação é ótimo.
106. Atenção! Rodei o programa de vocês para a multiplicação. Escolhi um multiplicador grande.
107. Mestre e no caso da divisão? Então Loirinha, a divisão de a por é dada por = /, Portanto = −/2. Assim temos = = /2.
108. Observem que lim →0 () = + ∞ Assim, dividir por valores pequenos propaga muito o erro absoluto.
109. Já o condicionamento relativo da divisão é dado por = | | = = / /2 = 1 ∙ Assim, o condicionamento relativo da divisão também é ótimo.
110. Atenção! Escolhi um divisor pequeno, = 0.001. Prestem atenção no condicionamento e no erro absolutos.
111. Mestre, a raiz quadrada, também é uma operação implementada como elementar no IEEE 754. Como é o condicionamento de ? Então Loirinha, seja r a função radiciação, = , para ≥ 0. Sabemos que = 1 2 , ≠ 0 . Assim temos, para > 0: = = 1 2 .
112. Observem que lim →0 , = + ∞ Assim, extrair a raiz quadrada de valores pequenos de propaga muito o erro absoluto.
113. Já o condicionamento relativo da radiciação é dado, para > 0, por = || | | = = 1/2 = 1/2 ∙ Assim, o condicionamento relativo da radiciação também é ótimo.
114. Escolhi = 2.∗ 10−5, um valor bem pequeno. Pelos resultados parece que não teremos muitos problemas na radiciação. Atenção!
115. Antes de encerrar, recomendamos fortemente a leitura deste artigo, colocado em “Leituras adicionais”, no site do Mestre.
116. Tchau, até a próxima aula!