Primos na Espiral de Ulam

Uma das razões dos números primos serem tão legais se deve ao fato deles se comportarem de forma estranha. Eles parecem aleatórios.

Algumas vezes você tem longos espaços entre dois números primos e, de repente, como os ônibus, vêm dois de uma vez só. No entanto, no fundo no fundo, eles não são completamente aleatórios.

Um matemático Polonês chamado Stanislaw Ulam, que foi para os EUA pouco antes da Segunda Guerra Mundial, estava, depois da guerra, em 1963, assistindo uma apresentação de um trabalho chatíssimo e longo com um papel e uma caneta na mão e resolveu fazer a seguinte brincadeira para se distrair:

No centro do papel colocou o número 1 e foi fazendo uma espiral quadrada com a sequência, conforme a figura abaixo:
ulam1

Depois disso, começou a circular os números primos nesta espiral, conforme a figura a seguir:
ulam2

Ele ficou surpreso pelo fato dos números primos caírem em diagonais. Como todos os números primos, exceto o 2, são números ímpares e como as diagonais nesta espiral alternam entre pares e ímpares, não é surpresa que os números caiam em diagonais alternadas, mas que algumas diagonais tenham mais primos do que outras.

Pouco tempo depois ele resolveu fazer um programa usando o computador MANIAC II para imprimir pixels nos pontos primos (exatamente como eu fiz em vermelho no excel) e conseguiu fazer uma imagem com os números até 65.025 (255 x 255).

Aqui abro um parêntese para falar do MANIAC II.

Ele era um computador criado em 1957 com 4096 words de 48bits (24kbytes) de memória RAM em Magnetic-core Memory e 12288 words de 48bits (576kbytes) em Williams tubes. Em média, uma multiplicação neste computador demorava 180 microsegundos e uma divisão 300 microsegundos. Uma verdadeira eternidade para os tempos de hoje.

Voltando à espiral do Ulam, o que ele conseguiu com a impressão da matriz de 255 x 255 de pixels foi a confirmação do que ele havia visto com 100 números. Realmente há um padrão no qual os números primos aparecem em diagonais, com intervalos, é claro, e algumas diagonais parecem ter mais primos do que outras.

ulam-255x255

Todas as linhas nesta espiral obedecem à seguinte equação quadrática: 4x^2+bx+c

Por exemplo, a diagonal que começa no número 3 tem a seguinte equação: 4x^2-2x+1
3, 13, 31, 57, 91, … (Confira na imagem acima)

Na prática, há uma hipótese de que estas diagonais podem servir para procurarmos números primos grandes, já que algumas diagonais têm mais primos do que outras diagonais. Melhor dizendo, algumas equações quadráticas têm mais chance de retornar números primos do que outras.

Um exemplo de diagonal com 40 números primos em sequência é a seguinte: x^2-x+41 que gera a seguinte sequência de 40 números primos:

41, 43, 47, 53, 61, 71, 83, 97, 113, 131, 151, 173, 197, 223, 251, 281, 313, 347, 383, 421, 461, 503, 547, 593, 641, 691, 743, 797, 853, 911, 971, 1033, 1097, 1163, 1231, 1301, 1373, 1447, 1523, 1601 (O 41o número é igual a 41^2 e, portanto, não é primo)

As diagonais com as maiores densidades de números primos conhecidas são a belezuras abaixo:

x^2 + x + 3399714628553118047

e (desculpe, não cabe na tela)

x^2 + x + 332518109806968781031500852571295088573128477514981900349983874538507313

O que interessa nisto tudo, é que, parece que, há fórmulas com mais densidade de primos do que outras (isso ainda não foi provado) e isto pode ajudar a resolver outros problemas, como a hipótese de Goldbach (que diz que todos os números pares maiores do que 2 podem ser expressados pela soma de dois números primos) ou a hipótese da existência de infinitos números primos gêmeos.

Como não poderia faltar, eu fiz uma implementação em javascript para mostrar esse grid, só que em vez de 255 x 255, o grid que eu fiz é de 1000×1000. O resultado está neste link. Obviamente você pode aumentar o tamanho do canvas para gerar coisas grandes, como essa de 25.000.000 de números (5000 x 5000) que eu fiz usando o mesmo código (Clica que aumenta).

ulam-5000x5000

 

Edição das 11:50 (Não podia faltar o código Python. O @jbvsmo depois dirá que está lento….)

#!/usr/bin/env python
import sys,math
import Image
 
size = int(sys.argv[1])
 
def sieveGen(siz):
    l = [2,3,5,7,11,13,17,19,23,29]
    if (siz < 31):
        return l
    for i in xrange(31,siz,2):
        isP = True
        rT = math.sqrt(i)
        for j in l:
            if i%j == 0:
                isP = False
                break
            if j>rT:
                break
        if isP:
            l.append(i)
    return l
 
mySieve = sieveGen(size)
 
 
def isPrime(n):
    if (n == 1):
        return False
    if (n < size):
        return (n in mySieve)
    isP = True
    sqrtN = math.sqrt(n)
    for j in mySieve:
        if (n%j)==0:
            return False
        if (j>sqrtN):
            return True
    return True
 
 
 
 
def spiral(N):
    im = Image.new("RGB", (N, N), "white")
    pix = im.load()
    red = (255,0,0)
    if(N%2):
        x = y = ((N-1)/2)
    else:
        x = y = (N/2)
    N2 = N*N
    dx = 1
    dy = 0
 
    val = 1
    amp = 1
    c = 0
 
    while (val <= N2):
        mvd = 0
        while ((mvd < amp) and (val <= N2)):
            if isPrime(val):
                pix[x-1,y-1] = red
            x += dx
            y += dy
            mvd += 1
            val += 1
 
        c += 1
        if (c == 2):
            c = 0
            amp += 1
 
        if (dx == 1):
            dx,dy = 0,1
        else:
            if (dy == 1):
                dx,dy = -1,0
            else: 
                if (dx == -1):
                    dx,dy = 0,-1
                else:
                    if (dy == -1):
                        dy,dx = 0,1
    im.transpose(Image.FLIP_TOP_BOTTOM).save("ulam.png")
 
spiral(size)

Portando para o Python 3

Em 2008 foi lançada a versão 3 da linguagem Python com diversas funcionalidades novas, mas parcialmente incompatível com os programas feitos para Python 2. Nestes 5 anos, a adoção ainda não é grande exatamente por causa deste fator e porque normalmente as pessoas não querem ter que aprender os detalhes que são necessários para essa mudança.

Hoje em dia existem diversas ferramentas que permitem desenvolver software que seja compatível com as duas versões, mas para isso é preciso entender o que mudou e adaptar. Este é será um guia rápido com dicas para migrar para Python 3 e aproveitar as melhorias que essa nova versão traz.

Parte I – Quando migrar

Por causa da incompatibilidade que existe, muitos desenvolvedores não podem começar a migração pois dependem de diversas bibliotecas que ainda não fizeram sua parte. Hoje em dia, grande parte das bibliotecas mais importantes da linguagem já são compatíveis com Python 3, o que torna bem mais plausível converter seu código.

Para computação númerica NumPy, SciPy, Matplotlib, IPython e Pandas são as ferramentas mais importantes e todas já foram portadas. Em desenvolvimento web, os frameworks mais importantes – Pyramid e Django – também funcionam perfeitamente em Python 3. Para criar GUI’s, PyQT4 e Tkinter podem ser utilizados sem problemas.

– O que falta então?

Normalmente são as pequenas bibliotecas que fazem trabalhos específicos e que já não são atualizadas há anos. Inclusive, este é um bom momento para se livrar delas. Procure no PyPI por algo semelhante que resolva seu problema e que seja compatível com Python 3.

– E se eu não quiser migrar?

Aproximadamente em 2017, o suporte para Python 2 vai terminar e este já não recebe novas funcionalidades desde 2010, com o lançamento do Python 2.7. Atualmente, só atualizações de segurança e outros bugs sérios são lançadas. Além disso, em breve o Python 3 será instalado por padrão na maioria das distribuições Linux, como já o é no Arch Linux e será no Ubuntu 13.04 ou 13.10.

– Mas mexer em código antigo que já funciona há anos é difícil…

Sim, isso é um grande impedimento. O ideal é começar a aprender com código novo e, quando surgir a oportunidade (ou necessidade) modificar os códigos antigos. A dica é, comece pelas suas bibliotecas para depois partir para as aplicações.

– Por onde eu começo?

O meu texto preferido sobre o tema é Dive into Python 3 que traz todos os detalhes sobre a linguagem e mostra o caso de migração de um software real. Aqui só faço uma introdução sobre o tema.

Independente de querer migrar ou não o seu código, é interessante conhecer o Python 3 e passar suportá-lo em código novo. Mesmo que você esteja escrevendo código com bibliotecas obsoletas, existem boas práticas que deixam permitem desenvolver já pensando no futuro.

Continuar lendo Portando para o Python 3

É preciso saber programar para conferir título de capitalização

Em primeiro lugar eu queria dizer que títulos de capitalização são a maior safadeza já inventada. É a típica parceria Caracu, onde o Banco entra com a cara.

Tendo dito isto, talvez você não saiba como é feito o sorteio dos tais prêmios. Um infeliz amigo que foi obrigado pelo seu gerente a fazer um título de capitalização pediu um help ontem e coloco abaixo a forma como resolvi.

Vamos às regras do título de capitalização do Banco Maldito (você sabe qual é):

  1. Só valem os sorteios da Loteria Federal feitos no último sábado do mês.
  2. Se o mês for Março, Junho, Setembro ou Dezembro vão ser escolhidos dois números. Um tradicional e outro especial.
  3. O número tradicional é formado pela dezena simples e unidade simples do primeiro prêmio da loteria federal e pela unidade simples do 2o, 3o, 4o e 5o prêmios da loteria federal.
  4. No caso do sorteio especial o segundo número é formado pela centena simples do primeiro prêmio da loteria federal e pela dezena simples do 1o, 2o, 3o, 4o e 5o prêmio.

Para resolver você precisa de uma tabela de resultados da loteria federal que pode ser conseguida aqui: http://www1.caixa.gov.br/loterias/loterias/federal/download.asp

Com isto em mãos e o programa em Python que vai abaixo você resolve seu problema. (obviamente os mestres de Python que leem o blog terão soluções melhores).

Mas o melhor conselho é *nunca* fazer título de capitalização.

Título de Capitalização é safadeza

#!/usr/bin/env python
import urllib
import time
import calendar
from bs4 import BeautifulSoup
 
fhtml = open("D_LOTFED.HTM").read()
soup = BeautifulSoup(fhtml)
meus_numeros = [
394465,558487,
572418,640294,
439592,329068,
368570,765895,
023206,847826]
 
l = []
for i in soup.table.tbody.findAll("tr"):
        a = []
        for k in i.findAll("td"):
                a.append(k.text)
        l.append(a)
 
def modalidadeMensal (a1):
        return int(a1[2][3:]+a1[3][4]+a1[4][4]+a1[5][4]+a1[6][4])
 
def modalidadeEspecial(a1):
        return int(a1[2][2:3]+a1[3][3]+a1[4][3]+a1[5][3]+a1[6][3])
 
def isLastSaturday(st):
        tupl = time.strptime(st,"%d/%m/%Y")
        if tupl[6]==5:
                if ((tupl[2])+7) > calendar.monthrange(tupl[0],tupl[1])[1]:
                        return True
        return False
 
def isEspecial(st):
        if (isLastSaturday(st)):
                tupl = time.strptime(st,"%d/%m/%Y")
                if tupl[1] in [3,6,9,12]:
                        return True
        return False
 
for i in l[1:]:
        if isLastSaturday(i[1]):
                print i[0],i[1],modalidadeMensal(i)
                if modalidadeMensal(i) in meus_numeros: print i
        if isEspecial(i[1]):
            print "*",i[0],i[1],%modalidadeEspecial(i)
            if modalidadeEspecial(i) in meus_numeros: print i

Quebrando o pau, matematicamente

Quando estudei na PUC-Rio havia uma matéria chamada Modelos Probabilisticos ministrada pelo pessoal da Telecom, que era conhecida por suas provas criativas.

Numa destas provas havia uma questão cujo enunciado é o que vai abaixo e serve como nossa leitura de final de semanamatem:

“Suponha que você quebre um palito num ponto aleatório do comprimento dele, escolha o maior pedaço e quebre novamente de forma aleatória. Qual é a probabilidade de você formar um triângulo com os três pedaços?”

Antes de sair procurando alopradamente no Google, não havia o Google ainda em 1996 quando fazia esta matéria, note que este problema é ligeiramente distinto do “three stick triangle problem”, neste nosso você escolhe para a segunda quebra o maior pedaço o que é diferente de quebrar um pedaço de pau em 3 pedaços diretamente.

Há duas formas de resolver este problema, uma forma computacional, que não era o caso na época, e uma forma analítica. A forma analítica é interessante porém exige um conhecimento matemático mais aprofundado. Já a forma computacional está disponível para qualquer um que saiba programar como mostramos no script Python abaixo.

from random import random
 
def isTriangle (a,b,c):
        lista = [a,b,c]
        lista.sort()
        if (lista[2] < (lista[1]+lista[0])):
                return True
        else:
                return False
 
def breakRule():
        p = random()
        if (p < 0.5):
                a = p
                p1 = random()
                b = (1-p) * p1
                c = (1-p) - b
        else:
                a = (1.0 - p)
                p1 = random()
                b = p * p1
                c = p - b
        return isTriangle (a,b,c)
 
sim = 0
nao = 0
 
for i in xrange(100000):
        if breakRule():
                sim += 1
        else:
                nao += 1
 
print "Número de vezes que deu certo: %d"%sim
print "Número de vezes que deu errado: %d"%nao
print "Probabilidade: %f"%(float(sim))/(sim+nao)

Agora a forma analítica de resolver.

Considerando a régua como tendo comprimento 1, vamos dividir o problema em duas partes que são simétricas. Na primeira parte a quebra ocorreu no ponto n (0 < n < 0.5) até 50% do comprimento da régua. Desta forma a segunda quebra não pode ocorrer, para respeitar a desiguldade triangular:

 

A < B + C

onde A é o comprimento do maior pedaço, nos trechos:

[n;\frac{1}{2}] ; [1 - (\frac{1}{2} - n);1]

Como pode-se ver na figura abaixo:

Desta forma, ao quebrar uma vez a régua no ponto n < 0.5 a probabilidade de formar um triângulo é dada por:

P_t(n) = \frac{n}{1-n}

Portanto a probabilidade de formar um triângulo é:

\int\limits_0^\frac{1}{2} P_t(n)dn = \int\limits_0^\frac{1}{2} \frac{n}{1-n}dn

Resolvendo a integral temos:

\log{2} - \frac{1}{2} = 0.193147

Como falamos que eram dois problemas simétricos a probabilidade de formar um triângulo é:

2 (\log{2} - \frac{1}{2}) = 0.386294

Que é o valor obtido pela nossa simulação.

Twitter – contando seguidores e comparando

Um amigo pediu para eu fazer um script que mostrasse a razão entre os seguidores de Jose Serra e de Dilma no twitter.

Serve para outros acompanhamentos e para obter outros dados de usuário.

import urllib
from BeautifulSoup import BeautifulStoneSoup
 
def usercnt(user):
        xml = urllib.urlopen ("http://twitter.com/users/show/%s"%user).read()
        soup = BeautifulStoneSoup(xml)
        return int(soup.user.followers_count.text)
 
def joseserra():
        return usercnt("joseserra_")
 
def dilma():
        return usercnt("dilmabr")
 
print joseserra()/float(dilma())

Neste momento a resposta do script é 1,9326

A resposta a chamada joseserra() é 505523
A resposta a chamada dilma() é 261595
Dá para fazer uns belos gráficos comparativos com esta ferramenta.

Ruido sonoro tem cor. E uma nova seção do Blog

http://en.wikipedia.org/wiki/Colors_of_noise

Como vocês podem ver no link acima ruído sonoro tem cor, pelo menos uma sensação de cor.

O ruído é útil para pessoas que sofrem de zumbido no ouvido e também ajuda a se concentrar mascarando sons externos.

Como precisava de uma fonte de ruído para fazer uma app para o IPhone escrevi um script em python e veio a idéia de criar uma seção com dicas de Python que agora fica no menu horizontal do Zeletron ou na página: http://www.zeletron.com.br/dicas-python

Quem tiver alguma receita de python interessante pode mandar para inserirmos lá.

MedCalc versao 0.6.0

Via BugSa

Devido a uma série de fatores eu fiquei sem celular da Nokia e por isso o MedCalc parou temporariamente.

Como muitos leitores, não sabia que tinha tantos … Mandaram e-mail pedindo um update estou postando outra versão (0.6.0) com suporte a Python 2.0 para S60 – esta versão foi empacotada pelo Marcelo Barros a quem agradeço publicamente.

Link para download: http://sites.google.com/site/medcalculator

Instruções:
1) Caso você não tenha o Python para S60 instale do link: https://garage.maemo.org/frs/download.php/7611/PyS60_binaries_certificate_error_fixed.zip
2) Instale o MedCalc

Se alguém tiver um celular S60 sobrando para emprestar eu posso fazer outras versões mais rapidamente. Quem quiser emprestar ou doar entre em contato comigo 🙂

Programa para descobrir sua idade auditiva

A audição é um sentido fascinante. Tipicamente o ser humano escuta frequências que variam de 20Hz a 20000Hz, no entanto, com o passar dos anos as pessoas vão perdendo audição na faixa mais próxima a 20000Hz, este fenômeno é conhecido como Presbiacusia.

Isto já foi explorado antes naquele toque de celular que só pessoas jovens escutam, no entanto ficou como uma curiosidade e há usos interessantes…

Vamos ensiná-lo a fazer um programa que teste a sua audição usando um computador e poucas linhas de código.

Você só precisa instalar a linguagem de programação Python.

Testando se escuta 10000Hz:

Abra um editor, digite as linhas abaixo e salve como teste.py; a seguir execute o arquivo teste.py – você ouvirá um som agudo por um segundo.

import winsound
winsound.Beep(10000,1000)

Como podemos inferir do código acima, a função winsound executa uma frequência determinada em hertz (de 37Hz a 32767Hz) durante um tempo especificado em milisegundos (no nosso caso 1000 milisegundos)

E 15000Hz? Fácil …

import winsound
winsound.Beep(15000,2000)

Note que neste caso aumentamos a duração para 2000ms ou 2 segundos.

Vamos fazer agora um programa para testar aproximadamente a frequência de corte da sua audição (vamos fazer incrementos de 500Hz e iniciar em 10000 Hz)

import winsound
for i in xrange(21):
	freq = i*500+10000 #10000, 10500, 11000, etc
	print "Tocando %dHz"%freq
	winsound.Beep(freq,3000) 
	#note que voce escuta um clique no comeco e fim do som
	#isto é a inércia do alto-falante

Agora que você descobriu a máxima frequência que você escuta temos um excelente meio de espantar as crianças de perto do computador, sem lhe incomodar…

import winsound
while 1:
       winsound.Beep(17500,200000)

Caso você queira espantar o gato e não as crianças use o seguinte código:

import winsound
while 1:
       winsound.Beep(23000,200000)

Utilize uma caixa de som bem potente 🙂

Artigo da Semana no Forum Nokia

Ontem fui acordado por um SMS do Marcelo Barros me dando uma ótima notícia!

O artigo “How to display transparent PNG on canvas with masks“, que eu escrevi em abril, foi escolhido pelo pessoal do Forum Nokia para ser o “Featured Article of Week“.


Para os desenvolvedores Python de plantão, deixo os links para o artigo em Inglês e em Português:

Agradeço mais uma vez ao pessoal do Forum Nokia pela escolha!

Livro de Python Grátis em Português

O Marcelo Barros colocou no blog dele uma dica ótima para os desenvolvedores Python; um link para baixar um livro de Python gratuito e em português. O livro tem licença Creative Commons.

Visite o blog dele para baixar o livro.

http://jedizone.wordpress.com/2009/05/27/python-para-desenvolvedores-novo-livro-em-creative-commons/

Boa leitura.

Ovi Store começa a aparecer hoje

Hoje em meu aplicativo de Download vi um icone Ovi Store. Pensei que era algum erro mas agora acabei de ler:

“The buzz has started with the Ovi Store showing up in folks’ Download! folder on certain operators in Asia (see report on AAS). We knew this would slip out, since so many of you are clever and eager to check out the store. What’s happening is that the Ovi Store is moving to production servers, meaning that the Ovi […]” – http://conversations.nokia.com/2009/05/25/quick-note-ovi-store-almost-here/

e aqui também:

http://info.abril.com.br/noticias/tecnologia-pessoal/ovi-store-da-nokia-comeca-a-funcionar-25052009-17.shl

Já ia esquecendo: com um dia de atraso saiu o novo PythonS60: https://garage.maemo.org/frs/?group_id=854

http://www.zeletron.com.br/2009/05/ei-facam-o-favor-de-aceitar-aplicativos.html