Aldeia Numaboa UM PORTAL DIFERENTE EM PORTUGUÊS DO BRASIL

Tutorial PARI/GP (I)

Detalhes

Categoria: Matemática Numaboa

Atualização: Domingo, 12 Abril 2009 16:16

Autor: vovó Vicki

Acessos: 13370

A Teoria dos Números, em última análise, é a ciência das propriedades algébricas elementares dos números inteiros. É um dos ramos da Matemática mais antigos e um dos preferidos pelos matemáticos profissionais e amadores. Muitos dos padrões mais importantes na teoria dos números foram descobertos através da análise de cálculos minuciosos e isto, até há pouco tempo, exigia um imenso talento e uma paciência apreciável. Atualmente, com o advento da sofisticada tecnologia da computação, a pesquisa de cálculos experimentais pode ser desenvolvida por praticamente todos os interessados. Só para se ter uma idéia, vários dos elementos chave que possibilitaram provar o Último Teorema de Fermat foram padrões descobertos com a ajuda do computador.

Em 1980, uma equipe de matemáticos franceses (Batut, Bernardi, Cohen, Olivier e Belabas) começaram a desenvolver um software para facilitar os cálculos da teoria dos números. Criaram o sistema PARI/GP, hoje um dos preferidos dos pesquisadores e, para coroar o sucesso, o software é distribuído gratuitamente.

PARI/GP é um sistema que permite realizar, entre outras coisas, cálculos de fatorização, teoria de números algébricos e curvas elípticas com enorme rapidez. Além disso, possui uma gama enorme de funções para cálculos de entidades matemáticas como matrizes, polinômios, séries de potências, números algébricos, etc, e muitas funções transcendentais. PARI é a biblioteca de funções e GP é a linguagem que permite acessá-las.

Se você tiver interesse, sugiro que faça o download do PARI/GP diretamente do site dos desenvolvedores para obter a versão mais atualizada. Instale o programa e, para que você possa começar a se divertir, aqui vai este pequeno tutorial baseado no tutorial dos autores. Neste texto, tudo o que é sugerido para você digitar está destacado em verde.

• Saudações!

Imagino que você já tenha instalado seu PARI/GP, que ele já esteja rodando e que a tela mostre gp

Experimente digitar 2 + 2, seguido de [ENTER]. O resultado não é nada além do esperado: deve aparecer %1 = 4, indicando que sua operação número 1 resultou em 4. Caso você tenha o hábito de digitar ponto e vírgula após os comandos, ou seja, "2 + 2;", o resultado é decepcionante: não aparece nada. Isto é porque a GP interpreta o ponto e vírgula como um comando para NÃO imprimir e o resultado é armazenado apenas na história da sessão. Experimente outras operações básicas, como subtração, multiplicação e divisão. A GP se comporta como uma calculadora, só que muito mais rápida!

Peraí! Na divisão, a GP parece estar brincando. Se digitamos 1/7, ela responde 1/7. Se digitamos 1/3, ela responde 1/3. O que é isto? É que a GP, calibrada para trabalhar com números inteiros, sempre tenta responder com o resultado mais preciso possível e o resultado mais preciso de 1/7 é 1/7 Caso você queira a expansão desta divisão, experimente digitar 1./7, 1/7., 1./7. ou 0. + 1/7. Neste caso, a maioria das máquinas devolve um valor fracionário de 28 casas (algumas dão 38 casas), ou seja, 0.142857... Neste caso, o ponto após qualquer um dos valores indica que são números reais imprecisos e a GP não se sente obrigada a responder com o valor mais exato.

:( Este negócio está ficando chato. Você deve estar pensando "Quer dizer que eu entupi meu HD com alguns mega de software só pra fazer umas continhas?". Se quiser desistir, então digite \q para sair do programa. Se achar que vale a pena continuar, experimente digitar exp(1). Melhor ainda, digite log(exp(1)).

Bem, se você chegou até aqui, é sinal de que está gostando do tutorial. Então está na hora de experimentar algumas novidades. Digite pi. O que? A resposta foi "pi" ao invés de 3.1415... Então pergunte o que está acontecendo com ?pi. A GP responde com "pi: user defined variable", ou seja, está mais perdida do que nós e informa que pi é uma variável definida pelo usuário. Tente digitar Pi (desta vez, com letra maiúscula). Ooooba, agora sim. Recebemos como resposta o valor de Pi com a precisão que nosso equipamento permite. Se perguntarmos ?Pi, a resposta vem rapidinho "Pi=Pi(): the constant Pi with current precision".

A coisa está ficando quente... digite exp(Pi * sqrt(163)). Hmmmm, este monte de zeros depois da vírgula e um 1 no fim não está com cara de ser um inteiro. É chegada a hora de mexer com a precisão. Digite \p 50 e depois repita a pergunta exp(Pi * sqrt(163)). Ahá! Aquele monte de zeros se transformou numa porção de noves e mais alguma coisa mostrando que o resultado não era mesmo um número inteiro.

A GP lida com casas decimais tão bem quanto com dígitos de números inteiros. Experimente calcular 1000!. Apesar da resposta do fatorial de 1000 ser um baita de um inteiro, ela aparece rapidinho na tela.

Agora seria bom voltar para a precisão de 28 dígitos, ou seja, dar o comando \p 28. Digite exp(75) e note que a GP mostra o resultado na notação científica 3.7332419667990011640254908317 E32. Pelo jeito, a precisão não está permitindo que o resultado seja mostrado com todos os números. Então é fácil. Digitamos, por exemplo, \p 40 e depois chamamos a última operação efetuada com % ou %número da operação. Ué, não mudou nada? É que a operação foi efetuada com precisão 28 e a GP não se esquece disso. Como agora estamos trabalhando com uma precisão de 40, é preciso repetir a operação para ver que exp(75) fornece 3733241966799001164025490831726470.0143428

Será que precisamos guardar de cabeça a precisão com a qual estamos trabalhando? Nada disso! Digite apenas \p e a GP devolve a informação. Além desta variável interna, a GP possui uma porção de outras. Digite default para obter a lista do que pode ser alterado de acordo com a conveniência do momento.

• Esquentando os motores

Antes de continuar, digite \p 28 para voltar para a precisão padrão. Agora, peça para a GP calcular 1/0. A GP responde com uma mensagem de erro óbvia: *** division by zero. Tente obter floor(exp(75)). Agora a resposta é *** floor: precision too low in truncr (precision loss in truncation), ou seja, precisão muito baixa no truncr (perda de precisão ao truncar). Como a função floor pede a parte inteira e esta é maior do que a precisão permite, a GP avisa com a mensagem de erro correspondente.

Ainda falando de erros, experimente calcular log(x). Creeeedo, a mensagem de erro diz *** log: log is not meromorphic at 0, ou seja, log não é meromórfico em 0 8) Para tirar a dúvida se a GP não sabe o que é log, digite ?log. A resposta é log(x): natural logarithm of x. Se a função é conhecida, então o erro é do usuário que não informou o valor de x. Elementar meu caro Watson ;)

Calcule a raiz quadrada de -1 com sqrt(-1). Você estava esperando uma nova mensagem de erro? Hahaha, a GP foi mais esperta e devolveu um número complexo! Experimente log(-2), exp(I*Pi) e I^I. Não custa repetir que a GP faz diferença entre letras maiúsculas e minúsculas (compare os resultados de I^2 e i^2). Só de brincadeira, experimente calcular 6*zeta(2)/Pi^2. O resultado está muito próximo, você não acha também?

A GP, apesar de chamar a função log, parece não ter entendido log(x). Que tal ver o que acontece com exp(x)? Eu fiquei embasbacada com o resultado 8) Esperava exp(x) e apareceu 1 + x + 1/2*x^2 + 1/6*x^3... etc e tal. É que a GP devolveu a Expansão de Taylor de 16 termos da função ao redor de x=0. 16 é o valor default para seriesprecision e que pode ser alterado digitando-se \ps n ou default(seriesprecision, n), onde n é o número de termos desejados na série de potências. Estas séries terminam sempre com 0(x^n), notação conhecida como "big-oh" e marca registrada das séries de potências.

A expansão de Taylor pode ser obtida para qualquer função que possa ser expandida ao redor de zero (o que, aliás, explica porque log(x) não funcionou). Experimente cos(x), cos(x)^2 + sin(x)^2, exp(cos(x)), gamma(1+x), exp(exp(x) - 1) e 1/tan(x). Aproveite para testar diversos comprimentos de série. Mas, e se quisermos uma expansão ao redor de um ponto diferente de zero? Basta informar a GP. Por exemplo, exp(x+5), cos(x+2) e até log(x+3).

Vamos testar uma coisa diferente. Digite (1 + x)^3 para obter o resultado x^3 + 3*x^2 + 3*x + 1. O resultado é um polinômio, novamente o resultado mais exato possível! Ué, e onde foi parar a expansão de Taylor com o "big-oh"? Se quisermos forçar a obtenção de uma série de potências é preciso usar a função Ser(), ou seja, Ser((1 + x)^3).

Resumindo: até este ponto vimos que, além de números inteiros, reais e racionais, o PARI/GP é capaz de lidar com números complexos, polinômios, funções racionais e séries de potências. O número de funções que ele oferece para trabalhar com estes tipos é muito grande. Neste tutorial serão vistos apenas alguns deles.

• Outros tipos na biblioteca PARI

Digite p = x * exp(-x). Como seria de se esperar, à variável p é atribuída a expansão para uma série de potências de 16 termos (se você não tiver alterado o valor default). Agora digite pr = serreverse(p) para pedir a reversão da série de potências p, ou seja, a função inversa. Isto só é possível para séries de potências cujo primeiro coeficiente não zero seja o de x¹. Para conferir se o resultado está correto, use subst(p, x, pr) ou subst(pr, x, p) para voltar ao x + 0(x^18).

Os coeficientes de pr seguem uma fórmula muito simples: multiplicando-se os coeficientes x^n por n! eles se transformam em inteiros. É exatamente isto o que a função serlaplace() faz. Digite ps = serlaplace(pr) para observar o resultado. Os coeficientes xⁿ agora são iguais a n^{n - 1}.

Vetores e Matrizes

A PARI não poderia deixar os vetores e as matrizes de fora. Por exemplo, digite [1, 2, 3, 4]. O resultado é o vetor da linha cujas coordenadas são 1, 2, 3 e 4. Se quisermos o vetor da coluna, basta pedir [1, 2, 3, 4]~, onde o til significa uma transposição. A não ser pelo til, o resultado não é muito diferente: obtemos [1, 2, 3, 4] e [1, 2, 3, 4]~. Então, qual é a vantagem que Maria leva? Experimente digitar [1, 2, 3, 4]~ e depois \b. A coluna agora é apresentada como

    [1]
    [2]
    [3]
    [4]

Agora peça m = [a, b, c; d, e, f]. Você acaba de definir uma matriz com duas linhas e três colunas. As matrizes são definidas por linhas separadas por ponto e vírgula ";". A matriz é mostrada em forma retangular, ou seja

    [a b c]
    [d e f]

Se você quiser vê-la na forma original, como foi definida, digite \a. Para manter este padrão de apresentação, digite default(output, 0). Se preferir a forma retangular, digite default(output, 1).

Uma vez que a matriz m está definida, experimente chamar alguns elementos. Abaixo estão alguns exemplos e os respectivos resultados:

    m[1,2]   .................   b
    m[1,]    .................   [a, b, c]
    m[,2]    .................   [b, e]~

A primeira referência é sempre a da linha e, a segunda, a da coluna. O primeiro índice é sempre 1 (e não o zero) e este valor default não pode ser alterado. Como temos acesso a todos os elementos da matriz, é lícito supor que estes elementos também possam ser modificados. Digite m[1, 2] = 5;. O ponto e vírgula, como já vimos, impede que o resultado seja mostrado na tela. Mas ele tem outras utilidades como, por exemplo, poder digitar vários comandos em uma só linha. Agora digite novamente m[1,] para obter [a, 5, c]. Nada impede que alteremos um vetor inteiro, por exemplo, m[1,] = [15, -17, 8]. Finalmente, tente m[,2] = [j, k]. Neste caso, obtemos uma mensagem de erro porque digitamos um vetor de linha para o vetor de coluna m[,2]. Corrija para m[,2] = [j, k]~ e deixe o til se encarregar da transposição.

A PARI possui várias matrizes "pre-calculadas" (veremos a construção mais genérica logo adiante), dentre elas a "Matriz de Hilbert" cujo coeficiente da linha i e coluna j é igual a (i + j - 1)^-1. Chame, por exemplo, h = mathilbert(20) para ver um exemplo.

O que é interessante a respeito das matrizes de Hilbert é que seus inversos e determinantes podem ser calculados explicitamente - e o inverso possui coeficientes inteiros. Além disso, sabe-se que são matrizes numericamente muito instáveis... menos na PARI. Como os coeficientes são dados como números racionais, o cálculo será efetuado com exatidão de modo que não ocorra nenhum erro. Confira. Digite d = matdet(h) para obter como resultado um número racional cujo numerador é 1 e com um denominador de 226 dígitos. Não perca tempo contando cada um dos dígitos do denominador, use a função sizedigit(1/d) ou transforme o resultado numa string e peça seu comprimento com #Str(1/d).

Agora digite hr = 1. * h; (o ponto e vírgula é pra não entupir novamente a sua tela) e, em seguida, dr = matdet(hr). Aqui é preciso esclarecer duas coisas. A primeira é que o cálculo é muito mais rápido do que no caso dos racionais. A explicação para isto é que a PARI está lidando com números reais com precisão 28 enquanto que, no caso dos racionais, ele está lidando com inteiros com 226 dígitos.

A segunda coisa, mais importante que a primeira, é que o resultado está completamente errado! Compare os resultados obtidos com 1. * d e 1. * h para ver que a coisa está realmente feia. Esta é, justamente, a catastrófica instabilidade das matrizes de Hilbert mencionada acima. Na verdade, a situação é ainda pior. Digite norm12(1/h - 1/hr). O resultado é maior do que 10⁵⁰, mostrando que alguns coeficientes de 1/hr mostram erros de até 10²⁴ (na verdade, o maior erro é igual a 4224*10²⁴ para o coeficiente da linha 15, coluna 15, que é um inteiro de 28 dígitos). Para obter o resultado correto, após arredondar o inverso, é preciso usar uma precisão de 56 dígitos (experimente). Ah, e a função norm12 dá o quadrado da norma L², isto é, a soma dos quadrados dos coeficientes.

Apesar dos vetores poderem ser definidos manualmente, digitando explicitamente seus componentes, muitas vezes estes elementos satisfazem uma lei simples, o que permite usar uma sintaxe diferente. Por exemplo, assuma que você queira um vetor cuja i -ésima coordenada seja igual a i ². Sem problemas. Neste caso entre com vector(10, i, i^2) se quiser um vetor de 10 elementos. De forma semelhante, se você digitar matrix(5, 5, i, j, 1/(i+j-1)) você obterá uma matriz de Hilbert de ordem 5 (ou seja, a função mathilbert é redundante). As variáveis i e j são variáveis que podem ser desprezadas. São usadas apenas para numerar, respectivamente, as linhas e as colunas (naturalmente, no caso de vetores, apenas uma delas estará presente). Não se esqueça de, juntamente com as dimensões do vetor ou matriz, indicar explicitamente o nome destas variáveis. Também é possível omitir as variáveis e a expressão final para se obter zero entradas, como em matrix(10, 20).

Observação: a letra I é reservada para o número complexo igual à raiz quadrada de -1. Portanto, é proibido usá-lo como variável. Existem mais dois nomes de variáveis reservados: Pi e Euler. Todos os nomes das funções também são proibidos. Por outro lado, não existem restrições em relação a i, pi ou euler.

Ao criar vetores ou matrizes, muitas vezes é útil usar operadores booleanos e a declaração if(). Na verdade, uma expressão if possui um valor que é igual à parte avaliada do if. Por exemplo, você pode digitar matrix(8, 8, i, j, if ((i-j)%2, 1, 0)) para obter uma matriz do tipo tabuleiro de 0 e 1. Mas, antes de usar uma matriz, esta precisa ser criada. Se, por exemplo, você digitar for (i = 1, 5, v[i] = 1/i) antes de criar o vetor v, vai obter uma mensagem de erro. O vetor v pode ser algo como v = vector(5).

Outro modo muito prático de criar vetores e matrizes é extraindo-os de outros maiores usando vecextract(). Por exemplo, se h é a matriz de Hilbert citada acima, vecextract(h, "11..20", "11.20") é seu quadrante inferior direito.

O tutorial de PARI/GP, apesar de ser apenas um resumo do que é possível realizar com este sistema, é bastante longo. Por este motivo transformei-o dois módulos. Se quiser continuar a leitura (e a aventura com o PARI/GP), continue com o Tutorial PARI/GP (II).

Grande abraço

vovó Vicki