A tatuagem do programador

Se fosse obrigatório fazer tatuagem, ao invés de tatuar aquelas asneiras de jogador de futebol, todo programador deveria tatuar seu algoritmo preferido.

tattoo

Este seria o meu, tal é a beleza e simplicidade do algoritmo concebido por Euclides para achar o máximo divisor comum entre os números inteiros x e y.

O mestre Donald Knuth diz no volume 2 da bíblia:

“[The Euclidean algorithm] is the granddaddy of all algorithms, because it is the oldest nontrivial algorithm that has survived to the present day.”

E qual o algoritmo que você tatuaria?

Sugestão de Leitura: Prime Obsession

Existem livros que apesar de falarem de coisas bastante sofisticadas são muito agradáveis de ler e Prime Obsession é um deles.

Neste livro, John Derbyshare, conta a saga inacabada da busca do que é talvez o maior problema aberto na matemática e para o qual há um prêmio de US$ 1.000.000,00 para quem resolver: a Hipótese de Riemann.

Eu confesso que antes de ler o livro sabia muito pouca coisa sobre a Hipótese de Riemann (HR), mas o autor consegue com muita didática ir conduzindo o leitor passo a passo pelos labirintos matemáticos para que qualquer pessoa, ainda que sem formação matemática possa compreender a natureza do problema.

A HR diz que a função zeta de Riemann tem todos os zeros não triviais com parte inteira igual a 1/2.

E eu com isso, poderia dizer o leitor. A hipótese de Riemann além de ser uma excelente ferramenta para entender o comportamento dos números primos é o fundamento de muitas descobertas matemáticas, prová-las significa provar muitas outras coisas em consequência e refutá-la significa derrubar por terra vários outros resultados.

Além de agradável leitura, Prime Obsession leva o leitor a descobrir a era romântica da matemática, saindo de Leonard Euler, passando por Gauss, Riemann e vários gigantes que revolucionaram o campo nos séculos XVIII, XIX e XX.

Se você quiser dar um voto de confiança para este escriba aqui, leia este livro. Garanto que não se arrependerá.

Na Amazon: Prime Obsession: Bernhard Riemann and the Greatest Unsolved Problem in Mathematics

Na Livraria Cultura: Prime Obsession

Prime Obsession: Resenha

 

NOTA: A função zeta de Riemann é:

ζ(s) = 1 + 1/2s + 1/3s + 1/4s + ..

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)

Números Lychrel e reminiscências do ensino fundamental

Era um sábado de calor no Rio de Janeiro no ano de 1987. Estávamos na sexta série do Colégio São Bento e a prof. Sandra Carelli percebendo que naquele dia não iríamos conseguir desvendar os mistérios da matemática de Papy, era o que se ensinava lá, mostrou-nos uma curiosidade matemática que nunca esqueci.

Suponha que você tem um número, por exemplo 124. Se você inverter os dígitos do número e somar com o número original é grande a possibilidade de obter um número palíndromo.

124 + 421 = 565

No entanto para alguns números a operação precisa ser repetida algumas vezes para obter um palíndromo: 279 por exemplo:

279 + 972 = 1251; 1251 + 1521 = 2772

Tendo captado nossa atenção ela prometeu um chocolate para quem achasse um número para o qual isto não funcionasse. Depois de muito tentar ela nos falou que um número para o qual isto não teria, a princípio, como se formar um palíndromo seria 196.

Quase 25 anos depois eu fui procurar sobre o assunto na Internet e descobri os chamados números de Lychrel. Os números de Lychrel seriam os que tem a mesma propriedade que 196. Embora nem mesmo 196 seja comprovadamente um número de Lychrel há outros candidatos como por exemplo: 295, 394, 493, 592, 689, 691, 788, 790, 879, 887, 978, 986

Abaixo, como homenagem a minha professora Sandra, vai um programa em Python para encontrar números de Lychrel.

#!/usr/bin/env python
def reverseNum(n):
        st = str(n)
        return int("".join([st[i] for i in xrange(len(st)-1,-1,-1)]))
 
def isPalindrome (n):
        st = str(n)
        rev = str(reverseNum(st))
        return st==rev
 
def isLychrel (n, num_interations):
        p = n
        for i in xrange(num_interations):
                if isPalindrome(p):
                        return i
                p = p + reverseNum(p)
        return -1
 
for i in xrange(1000):
        p = isLychrel(i,100)
        if (p < 0):
                print i,p

Para saber mais sobre os números de Lychrel e a busca por um palíndromo para 196 veja: http://www.p196.org

Números primos e algoritmos – Parte III (final)

Hoje mostraremos algoritmos que são rápidos e determinísticos (ou seja, sempre indicam correto se o número é primo ou não), alguns destes foram desenvolvidos no ano de 1983, no entanto em 2002 foi mostrado pela primeira vez um algoritmo, chamado AKS, para encontrar números primos que reunia as seguintes características:

  • Geral: para qualquer tipo de número primo e não apenas para alguns subconjuntos de primos.
  • Tempo Polinomial: com relação ao número de dígitos.
  • Determinístico
  • Incondícional: Ou seja, foi provado que o algoritmo não depende de nenhuma conjectura.

O teorema sobre o qual o algoritmo AKS se sustenta é uma generalização do pequeno teorema de Fermat que mostramos na semana passada.

pow ((x-a),n) % pow(x,n)-a = n

Onde n é o número primo e a é um primo em comum com n.

Abaixo vai uma das possíveis implementações dele:

#include <iostream>
using namespace std;
#include <math.h>
 
#ifdef _M_IX86
 
#define umulrem(z, x, y, m)	\
__asm	mov		eax, x	\
__asm	mul		y	\
__asm	div		m	\
__asm	mov		z, edx
 
#define umuladdrem(z, x, y, a, m)	\
__asm	mov		eax, x	\
__asm	mul		y	\
__asm	add		eax, a	\
__asm	adc		edx, 0	\
__asm	div		m	\
__asm	mov		z, edx
 
#else
 
#ifdef _MSC_VER
typedef unsigned __int64	Tu64;
#else
typedef unsigned long long	Tu64;
#endif
 
#define umulrem(z, x, y, m)	\
	{	\
	z = (unsigned int)(x * (Tu64)y % m);	\
	}
 
#define umuladdrem(z, x, y, a, m)	\
	{	\
	z = (unsigned int)((x * (Tu64)y + a) % m);	\
	}
 
#endif
 
static bool IsPrime(unsigned int n)
{
	if (n < 2) return false;
	if (n < 4) return true;
	if (n % 2 == 0) return false;
 
	const unsigned int iMax = (int)sqrt(n) + 1;
	unsigned int i;
	for (i = 3; i <= iMax; i += 2)
		if (n % i == 0)
			return false;
 
	return true;
}
 
static unsigned int LargestPrimeFactor(unsigned int n)
{
	if (n < 2) return 1;
 
	unsigned int r = n, p;
	if (r % 2 == 0)
	{
		p = 2;
		do { r /= 2; } while (r % 2 == 0);
	}
	unsigned int i;
	for (i = 3; i <= r; i += 2)
	{
		if (r % i == 0)
		{
			p = i;
			do { r /= i; } while (r % i == 0);
		}
	}
	return p;
}
 
static unsigned int Powm(unsigned int n, unsigned int e, unsigned int m)
{
	unsigned int r = 1;
	unsigned int t = n % m;
	unsigned int i;
	for (i = e; i != 0; i /= 2)
	{
		if (i % 2 != 0)
		{
			umulrem(r, r, t, m);
		}
		umulrem(t, t, t, m);
	}
	return r;
}
 
static unsigned int Inv(unsigned int n, unsigned int m)
{
	unsigned int a = n, b = m;
	int u = 1, v = 0;
	do
	{
		const unsigned int q = a / b;
 
		const unsigned int t1 = a - q*b;
		a = b;
		b = t1;
 
		const int t2 = u - (int)q*v;
		u = v;
		v = t2;
	} while (b != 0);
	if (a != 1) u = 0;
	if (u < 0) u += m;
	return u;
}
 
class CPolyMod
{
protected:
	// (mod x^r - 1, n)
	const unsigned int m_r;
	const unsigned int m_n;
	unsigned int m_deg;
	unsigned int * mp_a;
 
private:
	CPolyMod():m_r(0), m_n(0) { mp_a = NULL; };
 
public:
	// default value is x
	CPolyMod(unsigned int r, unsigned int n)
		: m_r(r), m_n(n)
	{
		m_deg = 1;
		mp_a = new unsigned int [2];
		mp_a[0] = 0; mp_a[1] = 1;
	}
 
	CPolyMod(const CPolyMod & p)
		: m_r(p.m_r), m_n(p.m_n)
	{
		m_deg = p.m_deg;
		mp_a = new unsigned int [p.m_deg + 1];
		unsigned int i;
		for (i = 0; i <= p.m_deg; ++i)
			mp_a[i] = p.mp_a[i];
	}
 
	virtual ~CPolyMod()
	{
		if (mp_a != NULL)
			delete [] mp_a;
	}
 
private:
	void _polySquare()
	{
		const unsigned int deg = m_deg;
		const unsigned int n = m_n;
		const unsigned int * const p_a = mp_a;
 
		const unsigned int degr = deg + deg;
		unsigned int * const p_ar = new unsigned int [degr + 1];
		unsigned int k;
		for (k = 0; k <= degr; ++k)
			p_ar[k] = 0;
 
		unsigned int j;
		for (j = 1; j <= deg; ++j)
		{
			const unsigned int x = p_a[j];
			if (x != 0)
			{
				unsigned int i;
				for (i = 0; i < j; ++i)
				{
					const unsigned int y = 2 * p_a[i];
					unsigned int t = p_ar[j + i];
					umuladdrem(t, x, y, t, n);
					p_ar[j + i] = t;
				}
			}
		}
		unsigned int i;
		for (i = 0; i <= deg; ++i)
		{
			const unsigned int x = p_a[i];
			unsigned int t = p_ar[2 * i];
			umuladdrem(t, x, x, t, n);
			p_ar[2 * i] = t;
		}
 
		m_deg = degr;
		delete [] mp_a;
		mp_a = p_ar;
	}
 
	void _polyMul(const CPolyMod & p)
	{
		const unsigned int deg = m_deg;
		const unsigned int n = m_n;
		const unsigned int * const p_a = mp_a;
 
		const unsigned int degr = deg + p.m_deg;
		unsigned int * const p_ar = new unsigned int [degr + 1];
		unsigned int k;
		for (k = 0; k <= degr; ++k)
			p_ar[k] = 0;
 
		unsigned int j;
		for (j = 0; j <= p.m_deg; ++j)
		{
			const unsigned int x = p.mp_a[j];
			if (x != 0)
			{
				unsigned int i;
				for (i = 0; i <= deg; ++i)
				{
					const unsigned int y = p_a[i];
					unsigned int t = p_ar[j + i];
					umuladdrem(t, x, y, t, n);
					p_ar[j + i] = t;
				}
			}
		}
 
		m_deg = degr;
		delete [] mp_a;
		mp_a = p_ar;
	}
 
	void _Mod()
	{
		unsigned int deg = m_deg;
		unsigned int * const p_a = mp_a;
		while (deg >= m_r)
		{
			p_a[deg - m_r] += p_a[deg];
			if (p_a[deg - m_r] >= m_n) p_a[deg - m_r] -= m_n;
			--deg;
 
			while (p_a[deg] == 0) --deg;
		}
		m_deg = deg;
	}
 
	void _Norm()
	{
		const unsigned int deg = m_deg;
		const unsigned int n = m_n;
		unsigned int * const p_a = mp_a;
		if (p_a[deg] != 1)
		{
			const unsigned int y = Inv(p_a[deg], m_n);
			unsigned int i;
			for (i = 0; i <= deg; ++i)
			{
				unsigned int t = p_a[i];
				umulrem(t, t, y, n);
				p_a[i] = t;
			}
		}
	}
 
public:
	CPolyMod & operator = (const CPolyMod & p)
	{
		if (&p == this) return *this;
		m_deg = p.m_deg;
		delete [] mp_a;
		mp_a = new unsigned int [p.m_deg + 1];
		unsigned int i;
		for (i = 0; i <= p.m_deg; ++i)
			mp_a[i] = p.mp_a[i];
		return *this;
	}
 
	int operator != (const CPolyMod & p) const
	{
		if (m_deg != p.m_deg)
			return true;
		unsigned int i;
		for (i = 0; i <= m_deg; ++i)
			if (mp_a[i] != p.mp_a[i])
				return true;
		return false;
	}
 
	CPolyMod & operator += (unsigned int i)
	{
		const unsigned int t = i % m_n;
		mp_a[0] += t;
		if (mp_a[0] >= m_n) mp_a[0] -= m_n;
		return *this;
	}
 
	CPolyMod & operator -= (unsigned int i)
	{
		const unsigned int t = m_n - i % m_n;
		mp_a[0] += t;
		if (mp_a[0] >= m_n) mp_a[0] -= m_n;
		return *this;
	}
 
	CPolyMod Pow(unsigned int e) const
	{
		unsigned int er = 1;
		unsigned int j;
		for (j = e; j != 1; j /= 2)
		{
			er = 2 * er + (j % 2);
		}
 
		CPolyMod t(*this);
		unsigned int i;
		for (i = er; i != 1; i /= 2)
		{
			t._polySquare();
			t._Mod();
			if (i % 2 != 0)
			{
				t._polyMul(*this);
				t._Mod();
			}
		}
		t._Norm();
		return t;
	}
};

Com isto encerramos nossa série de introdução aos números primos. Se você quiser ler mais sobre o assunto eu recomendo este artigo avançado:
http://www.cse.iitk.ac.in/users/manindra/algebra/primality_v6.pdf
e este artigo em termos mais leigos:
http://www.flonnet.com/fl1917/19171290.htm

Números primos e algoritmos – Parte II

Nesta segunda parte vamos mostrar testes que rapidamente dizem se um numero é provavelmente primo, provavelmente porque como veremos abaixo eles checam rapidamente se não é primo, mas quando dizem que é primo … Aí tem alta probabilidade de ser primo e se você quiser ter certeza precisa usar um algoritmo mais lento.

Este que vou mostrar se baseia no Pequeno Teorema de Fermat. Este teorema é bem simples de enunciar:

Se um número p é primo então para todos os números a entre 1 e p-1 a equação abaixo tem que ser respeitada (em notação C++):

pow (a,p-1) % p = 1

Assim fazemos o teste para alguns números a entre 1 e p-1 e caso a equação seja verdadeira para eles podemos dizer que é provavelmente primo. Se falhar para algum número não é primo com certeza.

#define NUMTESTES 10 // vamos testar 10 vezes, pode ser um número maior ...
 
int isPrimo (int n) {
   if (n==1)
       return 0; // dica de um leitor da semana passada
   if (n<4)
       return 1;
   for (int i=0; i<NUMTESTES; i++) {
       if (i+1 > n-1)
            break;
       // aqui estou supondo que há uma biblioteca 
       // de precisão arbitrária como a mathX
       // ou usando uma técnica de obter modulo 
       // em exponenciação como a 
       // Fast Modular Exponentiation
       if (pow(i+1,n-1)%p != 1)  
            return 0;
   }
   return 1;
}

Caso o resultado seja 0 você pode estar certo de que não é primo, se for 1, você pode ter uma enorme probabilidade de ser primo.

Para algumas aplicações é preciso ter certeza que o número é primo e como existem os números de Carmichael que se comportam como primos nestes testes probabilísticos vamos ver na próxima parte como testar rapidamente e com certeza (já adianto que é bem complicado do ponto de vista matemático…)

Pierre Fermat ajuda você a encontrar números primos

Como sempre, se houver erros favor apontar nos comentários.