Turbinando a fatoração de fermat

Podemos turbinar a fatoração de Fermat, evitando  a computação das raízes quadradas de de todos $a^2 - N$.

Dado um valor de m, podemos reduzir a quantidade de números que devem ser computados.

Por exemplo, os quadrados perfeitos módulo 20 terminam em 0,1,4,5,9 ou 16. Dessa maneira, Eliminamos 14 possibilidades módulo 20.

$b^2 = a^2 - N$ deve terminar em 0,1,4,5,9 ou 16.

Vamos calcular as possibilidades $a^2$ módulo 20 ser um quadrado perfeito dado que $b^2$ é um quadrado perfeito.  

a^2 = b^2 + N 

Vamos supor que  N módulo 20 = 3

$b^2 + N$ mod 20 
0 + 3 = 3
1 + 3 = 4
4 + 3 = 7
5 + 3 = 8
9 + 3 = 12
16+ 3 = 19

Dado que $b^2$ é um quadrado perfeito, então a^2 mod 20 = 1

Os valores de a tal que $a^2$ mod 20 = 1 são:
1,9,11,19

Precisamos apenar computar os valores de a que termina 1,9,11 e 19. Dessa maneira, eliminamos 16 possibilidade de 20. 

#include <stdio.h>
#include <stdlib.h>
#include <set>

using namespace std;

int main(){
 int n,m,i;
 set <int> squared;
 set <int> squaredA;
 set <int> a;
 set <int>::iterator it;
  
 scanf("Entre com m:");
 scanf("%d",&m);

 printf("Os seguintes numeros sao quadrado perfeitos modulo %d\n",m);
 for(i=0;i<m;i++)
  squared.insert((i*i)%m);
 
 for(it= squared.begin(); it != squared.end(); it++)
  printf("%d\n", *it);
 
 printf("Entre com n:");
 scanf("%d",&n);
 
 printf("%d modulo %d = %d\n",n,m,n%m);
 
 //a^2 - b^2 = N
 //a^2  = N + b^2
 
 printf("a^2 tem que terminar em\n");
 for(it= squared.begin(); it != squared.end(); it++)
  if( squared.find((*it + n%m)%m)!= squared.end() ){
   printf("%d\n", (*it + n%m)%m);
   squaredA.insert((*it + n%m)%m);  
  }
 
 printf("a tem que terminar em\n"); 
 for(i=0;i<m;i++){
  if( squaredA.find((i*i)%m)!= squaredA.end() ){
   a.insert(i%m);
  }
 }
  
 for(it= a.begin(); it != a.end(); it++)
  printf("%d\n",*it);
  
 system("PAUSE"); 

}

Saída


20
Os seguintes numeros sao quadrado perfeitos modulo 20
0
1
4
5
9
16
Entre com n:17
17 modulo 20 = 17
a^2 tem que terminar em
1
a tem que terminar em
1
9
11
19
Pressione qualquer tecla para continuar. . .

Fonte: http://en.wikipedia.org/wiki/Fermat's_factorization_method

Fatoração de Fermat

Este algoritmo de fatoração foi desenvolvido por  Pierre de Fermat. O algoritmo é baseado na representação do números inteiros ímpares como a diferença de dois  quadrados.

$N = a^{2} - b^{2}$

Essa diferença pode ser fatorada da seguinte maneira:

$a^{2} - b^{2} = (a+b)(a-b)$

Se N = cd é uma fatoração de N então
$N = a^{2} - b^{2}$
$c = a+b$
$d = a-b$
$c+d = 2a$
$a = \frac{c+d}{2}$
$c-d = 2b$
$b = \frac{c-d}{2}$
$N = (\frac{c+d}{2})^2 - (\frac{c-d}{2})^2$

Como a = $\sqrt{N + b^2}$,logo a >= $\sqrt{N}$ e a < N.

A idéia do algorimor é percorrer todos os valores de  a = |$\sqrt{N}$|...N. Para cada valor de a, verifique se b = $\sqrt{a^2 - N}$ é inteiro e se (a+b) (a-b) são fatores não triviais de N.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <map>

using namespace std;


int eh_quadrado(int x){
 int y = (int)(sqrt(x+0.5));
 return y*y == x;
} 

int fermat_factor(int N){
 int a,b,b2,c,d;
 for(a=(int)ceil(sqrt(N));a<=N;a++){
  b2 = a*a - N;
  if( eh_quadrado(b2) ){
   b = (int)sqrt(b2);
   c = a+b;
   if(c!=1 and c!=N)
    return c;
  }
 }
}

int primo(int n){ 
 if(n==1) return 1;
 else if( n< 4) return 1;
 else if( n%2 ==0) return 0;
 else if( n< 9) return 1;
 else if( n%3 ==0) return 0;
 else {
  int r = (int)sqrt(n);
  int f = 5;
  
  while( f<=r){
   if(n%f==0) return 0;
   if(n%(f+2)==0) return 0;
   f=f+6;
  }
 }
 return 1;
}


void factorization(int N, map <int, int> & fatores ){
 int f;
 
 if(primo(N)){
  fatores[N]++;
 }else{
  f = fermat_factor(N);
  factorization(f, fatores);
  factorization(N/f, fatores);
 }
}



void fermat_factorization(int N, map <int, int>  & fatores){
 while(N%2==0) { fatores[2]++; N=N/2;}
 factorization(N, fatores);
}

void trial_factorization(int N, map <int, int>  & fatores){
 
 int i = 3;
 int limite;
 
 while(N%2==0) { fatores[2]++; N=N/2;}
 
 limite = (int)sqrt(N);
 
 while( N!=1 && i <= limite ){   
   while(N%i==0) { fatores[i]++; N=N/i;}
   limite = (int)sqrt(N);  
   i = i+2;
 }
 
 if(N!=1){ fatores[N]++;}
 
}


int main(){
 
 /*
 for(int i=1;i<=10000;i++)
  if(eh_quadrado(i))
   printf("%d\n",i);
 */
 
 map <int, int> fatores;
 map <int, int>::iterator iter;
 
 printf("fatores:\n");
 fermat_factorization(675734124, fatores);
 
 printf("%d\n", fatores.size());
 
 for(iter = fatores.begin(); iter!= fatores.end(); iter++){
  printf("%d %d\n", iter->first , iter->second);
 }
 
 fatores.clear();
 
 trial_factorization(675734124, fatores);
 
 printf("fatores:\n");
 
 printf("%d\n", fatores.size());
 
 for(iter = fatores.begin(); iter!= fatores.end(); iter++){
  printf("%d %d\n", iter->first , iter->second);
 }
 return 0;

}

Saída:


fatores:
5
2 2
3 1
13 1
113 1
38333 1
fatores:
5
2 2
3 1
13 1
113 1
38333 1

Algoritmo Pollard's rho

O Algoritmo $\rho$ de Pollard é um algoritmo de fatorização desenvolvido por Pollard em 1975. O algoritmo $\rho$ Pollard é baseado em dois aspectos importantes.

1) O algoritmo utiliza uma função módulo n como um gerador de sequência pseudo-aleatória. 

$x_{n+1} = {x_{n}}^2 + a$ (mod n)

A sequência gera números pseudo-aleatórios  distintos  até cair em um ciclo. O tempo esperado até a sequência tornar-se cíclica e o tamanho esperado do ciclo são proporcionais a $\sqrt{n}$. A mesma observação realizada no paradoxo do aniversário garante esse fato

Podemos descobrir se dois números x e y, são congruentes módulo p,  realizamos o seguinte cálculo

$MDC( abs(x-y), n)  \leq n$

que é igual a p.


2) A detecção do ciclo na sequência é baseada na idéia atribuída a
Floyd conhecida como algoritmo da tartaruga e do coelho comparando a sequência $x_{i}$ com $x_{2i}$ para todo i. A sequência $x_{i}$ representa a tartaruga e a sequência $x_{2i}$ representa o coelho que move duas vezes mais rápido.

Algoritmo 

Entrada: n, the integer to be factored; and f(x), a pseudo-random function modulo n

Saída: a non-trivial factor of n, or failure.

1. x ← 2, y ← 2; d ← 1
2. While d = 1:
   1. x ← f(x)
   2. y ← f(f(y))
   3. d ← GCD(|x − y|, n)
3. If d = n, return failure.
4. Else, return d.

Fonte: http://en.wikipedia.org/wiki/Cycle_detection


No caso do algoritmo não encontrar um fator, nós vamos utilizar um f(x) diferente. O algoritmo não funciona
quando n é primo, uma vez que, d sempre será 1.

int mulmod(int x, int y, int n){
 return ( x%n * y%n )%n;
}

int addmod(int x, int y, int n){
 return ( x%n + y%n )%n;
}

//f(x) mod n  = x*x + c mod n
int f(int x, int c, int n){
 return addmod( mulmod(x,x, n) , c , n );
}

int abs(int x){
 if(x<0) return -x;
 else return x;
}

int gcd(int a, int b){
 int r;
 while(b!=0){
  r = a%b;
  a = b;
  b = r;
 }
 return a;
}

int primo(int n){
 
 if(n==1) return 0;
 else if( n< 4) return 1;
 else if( n%2 ==0) return 0;
 else if( n< 9) return 1;
 else if( n%3 ==0) return 0;
 else {
  int r = sqrt(n);
  int f = 5;
  
  while( f<=r){
   if(n%f==0) return 0;
   if(n%(f+2)==0) return 0;
   f=f+6;
  }
 }
 return 1;
 
}

int rho( int (*f)(int, int , int) ,int c, int n){
 int x,y,d;
 x = 2;
 y = 2;
 d = 1;
 
 if (primo(n)) return -1;
 
 while(d==1){
  x = f(x,c,n);
  y = f(f(y,c,n),c, n);
  d = gcd(abs(x-y), n);
    
 }
 
 if(d==n) return rho(f,++c, n);
 else return d;
 
}




Seja n = 8051 e f(x) = (x2 + 1 ) mod 8051.

i
xi
yi
GCD(|xi − yi|, 8051)
1
5
26
1
2
26
7474
1
3
677
871
97


Seja n = 8051 e f(x) = (x2 + 2 ) mod 8051.

i
xi
yi
GCD(|xi − yi|, 8051)
1
6
38
1
2
38
5709
1
3
1446
3607
1
3
5709
1227
83

Fonte:
http://en.wikipedia.org/wiki/Pollard's_rho_algorithm
http://comeoncodeon.wordpress.com/2010/09/18/pollard-rho-brent-integer-factorization/
http://introcs.cs.princeton.edu/java/78crypto/PollardRho.java.html

Logaritmo discreto

O problema do logaritmo discreto é encontrar o valor de x  dado a, b e n tal que

 $a^{x} = b (mod n)$

O algoritmo básico seria calcular todas as potências de  a,a^2,a^3,...,  até encontrar o valor b.

O algoritmo Shanks basea-se na reescrita de x como $x = im + j$, com $m = ceil (sqrt(n) )$ , $0<=i
$a^{im + j} = b$
$a^{im} a^{j} = b$
$a^{j}  = b(a^{-m})^{i}$

O algoritmo pré-computa $a^{j}$ para alguns valores de j. Para cada valor de i, ele testa se existe algum j que

$a^{j}        = b(a^{-m})^{i}$

Se sim, devolve o valor im+j
Se não, passa para o próximo valor de i

Pseudo-código

1. m ← Ceiling(√n)
2. para todo j onde 0 ≤ j < m:
   1. Compute α^j mod n and armazene o par (j, α^j) em uma tabela. 
3. Compute α^−m mod n
4. γ ← β. (set γ = β)
5. For i = 0 to (m − 1):
   1. verifique se existe γ como a segunda componente (α^j) de algum par na tabela.
   2. Se sim, devolva im + j.
   3. If não, γ ← γ • α^−m mod n

Código

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
using namespace std;

int n,a,b,m;
int g;

typedef struct {
 int index;
 int value;
} table; 

//função que calcula o mdc estendido ax + by = mdc(a,b)
int mdc(int  a, int b, int *x, int *y);

//calcula o inverso multiplicativo modulo n
//a*x = 1 mod n
int inverso(int a, int n);

//algoritmo de exponeciacao rapida modular
int fastexp(int a,int b, int n);

//busca binaria do valor g no vetor giant de tamanho n
int busca_binaria(int g, int n, table giant[]);

//funcao usada no qsort
int compara(const void *a, const void *b);

int shanks(int a, int b, int n){
 
 int g,m,r;
 int i,j;
 

 m = (int)ceil( sqrt(n) );
 
 printf("valor de m: %d\n",m);
 
 table giant[m];
 
 a = a%n;
 b = b%n;
 
 giant[0].index = 0;
 giant[0].value = 1;
 
 //Passo de Gigante
 for(j=1;j<m;j++){
  giant[j].index = j;
  giant[j].value = (giant[j-1].value*a)%n;
 }
 
 qsort(giant, m , sizeof(table), compara );
 
 //Passo de Gigante
 printf("Passo gigante\n");
 for(j=0;j<m;j++){
   printf("%d %d\n",giant[j].index, giant[j].value);
 }
 
 r = fastexp(inverso(a,n),m,n);
 g = b;
 
 //Passo de Baby
 printf("Passo de bebe\n");
 for(i=0;i<m;i++){
  printf("g: %d\n",g);
  j = busca_binaria(g,m,giant);
  if( j!=-1 ){
    printf("j: %d\n", giant[j].index); 
    return i*m + giant[j].index;
  }else{
   g = (g*r)%n;
  }
 }    

}

int main(){
 
 int j;
 int r;
 
 scanf("%d %d %d",&a,&b,&n);

 printf("%d\n",shanks(a,b,n));
 
 return 0; 

}

int mdc(int  a, int b, int *x, int *y) {
  int xx, yy, d;
  if(b==0) {
    *x=1; *y=0;
    return a;
  }

  d = mdc(b, a%b, &xx, &yy);
  *x = yy;
  *y = xx - a/b*yy;
  return d;
}


int inverso(int a, int n){
  int x,y,d;
  d = mdc(a,n,&x,&y);
  if(x<0){
    x = x+n;
  }
  return x;
}

int fastexp(int a,int b, int n){

  long long int x;

  if(b==0) return 1;
  if(b==1) return a;

  if(b%2==0){
    x = fastexp(a,b/2,n)%n;
    return (x*x)%n;
  }else{
    return (a*fastexp(a,b-1,n))%n;
  }

}

int compara(const void *a, const void *b){
 return ((table*)a)-> value - ((table*)b)->value;
}

int busca_binaria(int g, int n, table giant[]){
 int i,f;
 int m;
 
 i = 0;
 f = n-1;

 while(i<=f){
  m = (i+f)/2;
  if(giant[m].value == g) return m;
  else if(giant[m].value > g){
   f = m-1;
  }else {
   i = m+1;
  }
 }
 
 return -1;
  
}

Saída

5 315 317

valor de m: 18

Passo gigante

0 1

1 5

2 25

8 81

9 88

6 92

10 123

3 125

7 143

17 154

13 159

14 161

15 171

16 221

12 222

5 272

11 298

4 308

Passo de bebe

i 0 g: 315

i 1 g: 303

i 2 g: 219

i 3 g: 265

i 4 g: 270

i 5 g: 305

i 6 g: 233

i 7 g: 46

i 8 g: 5

j: 1

145

Pressione qualquer tecla para continuar. . .

$log_{5} 315 = mi + j = 18(8) + 1 = 145 $

A complexidade deste algoritmo é O($\sqrt(n) lg \sqrt(n)$)

Calculadora de logaritmos discretos
http://www.numbertheory.org/php/discrete_log.html


Referências:
http://en.wikipedia.org/wiki/Baby-step_giant-step
http://pastebin.com/F8Nin82x
http://pt.scribd.com/doc/63892535/38/Algoritmos-para-PLD

Inverso Multiplicativo

O objetivo deste post é apresentar diversas maneiras de encontrar o inverso multiplicativo módulo P. O inverso multiplicativo de um número a é um número x tal que
ax = 1 mod P
A forma mais direta de encontrar o inverso multiplicativo é buscar o x que satisfaz a propriedade acima:

/*
Encontrar x tal que ax = 1 mod P
*/

int inv(int a){
  int x;

  for(x=1;x<=P;x++){
    if((a*x)%P==1)
      return x;
  }

}

Durante a busca do inverso multiplicativo de a, podemos encontrar também o inverso multiplicativo de x.

int inverso[P+1];
int inv2(int a){
  int x;
  if(inverso[a]!=0) return inverso[a];
  else{
    for(x=1;x<=P;x++){
      if((a*x)%P==1){
        inverso[a] = x;
        inverso[x] = a;
      }
    }
    return inverso[a];
  }

}

Pequeno Teorema de Fermat

a^p-1= 1 mod p

Assim,

a*a^p-2=1 mod p

inv(a) = a^p-2 mod p

/*Pequeno Teorema de Fermat (a^P-1) = 1 mod P*/
int inv3(int a){
  int i;
  long long int x;
  x=1;
  for(i=1;i<=P-2;i++) x = (x*a)%P;
  return x;
}

Exponenciação Rápida

int fastexp(int a,int b){

  long long int x;

  if(b==0) return 1;
  if(b==1) return a;

  if(b%2==0){
    x = fastexp(a,b/2)%P;
    return (x*x)%P;
  }else{
    return (a*fastexp(a,b-1))%P;
  }

}


int inv4(int a){
  return fastexp(a,P-2);
}

Algoritmo de Euclides Estendido

/*
ax + by = 1
ax = 1 mod b
*/

int mdc(int  a, int b, int *x, int *y) {
  int xx, yy, d;
  if(b==0) {
    *x=1; *y=0;
    return a;
  }

  d = mdc(b, a%b, &xx, &yy);
  *x = yy;
  *y = xx - a/b*yy;
  return d;
}


int inv5(int a){
  int x,y,d;
  d = mdc(a,P,&x,&y);

  if(x<0){
    x = x+P;
  }

  return x;


}

Crivo de Eratóstenes

Embora não exista uma regra para descobrir todos os primos, é possível encontrar todos os primos menores que um dado número. Uma maneira de encontrar é através de um processo chamado de crivo que elimina sistematicamente os números compostos deixando passar pelo crivo apenas os números primos. Vamos apresentar o crivo desenvolvido por Eratóstenes (200 A.C), o terceiro bibliotecário-chefe da Biblioteca de Alexandria. 

Vamos aplicar o crivo de Eratóstenes para encontrar todos os números primos menores que 10.

Inicialmente, vamos considerar que todos os números de 2 até 10 são primos:
2 3 4 5 6 7 8 9 10

A cor preta representa os números não marcados. A cor vermelha representa os números eliminados. A cor azul representa os números que passaram pelo crivo, ou seja, os números primos. O primeiro número não marcado recebe a cor azul e todos os seus múltiplos recebem a cor vermelha.
2 3 4 5 6 7 8 9 10
2 3 4 5 6 7 8 9 10
2 3 4 5 6 7 8 9 10
2 3 4 5 6 7 8 9 10

Algumas modificações simples podem ser realizadas no algoritmo para tornar o algoritmo executa menos operações:

1. Inicialmente, vamos considerar que todos os números no intervalo são primos.
2. Elimine os múltiplos do número p que passou pelo crivo a partir de p*p.
3. Faça o crivo até 

 

Simulação do Crivo de Eratóstenes

http://www.cut-the-knot.org/Curriculum/Arithmetic/Eratosthenes.shtml

Código em C

#include <stdio.h>
#include <math.h>
#include <stdlib.h>

#define MAX 1000001
#define NUM 78498 
 
int main(){
  int i,j;
  int limite;
  char ehprimo[MAX];
  int cont=0;
  int primos[NUM];
  FILE *fp;
  fp = fopen("primos.txt","wt");

  for(i=2;i<MAX;i++) ehprimo[i]=1;
  limite = (int)sqrt(MAX);
  for(i=2;i<limite;i++){
    if(ehprimo[i]){
      for(j=i*i;j<MAX;j=j+i)
        ehprimo[j] = 0;
    }
  }

  for(i=2;i<MAX;i++){
    if(ehprimo[i]){
      fprintf(fp,"%d %d\n",cont,i);
      primos[cont]=i;
      cont++;
    }
  }
  printf("%d\n",cont);   
  system("PAUSE");

}

Esse programa gera um arquivo contendo todos os números primos menores que 1000000.

Código C++

#include <iostream>
#include <fstream>
#include <bitset> 
#include <vector>
#include <cmath>
#include <cstdlib>
#define MAX 1000001

using namespace std;

int main(){

  bitset<MAX> bs; //vetor de bits
  vector <int> primos; //vetor de primos
  int limite;
  ofstream outfile ("primos2.txt");

  bs.reset(); //seta todos os numeros para 0 
  bs.flip();  //seta todos os numeros para 1
  
  cout << "Inicio do Crivo" << endl;
  limite = (int)sqrt(MAX);

  for (int i = 2; i <= limite; i++){ 
    if (bs.test((size_t)i)) {
      for (int j = i * i; j < MAX; j += i) 
        bs.set((size_t)j, false); 
    }
  }
  
  for(int i=2; i< MAX; i++)
    if (bs.test((size_t)i)){
      outfile << i << endl;
      primos.push_back(i);
    } 
  system("PAUSE");
}

Fatoração em primos

O resultado mais importante sobre os números primos é o chamado Teorema Fundamental da Aritmética que garante que todo número maior que 1 pode ser decomposto em fatores primos. Vamos apresentar o algoritmo que executa mais operações, mas que é mais fácil de ser entendido. O algoritmo consiste em testar os números menores ou iguais   e ir dividindo a medida do possível:

#include <stdio.h>
#include <math.h>
#include <vector>
#include <iostream>
#include <stdlib.h>
using namespace std;

vector <int> fatora(int n){
  int i,limite;
  vector <int> primos;
  i=2;
  limite = (int)sqrt(n);
  while(n > 1 && i<=limite){
    //se i divide n entao i eh primo retira todos os fatores i de n
    while(n%i==0){
      primos.push_back(i);
      n=n/i;
    }
    i=i+1;
  }
  //se i >= raiz de n entao n eh primo
  if(n>1)
    primos.push_back(n);
  return primos;
}

int main(){
  vector <int> primos;
  primos = fatora(120);
  printf("%d",primos[0]);
  for(int i=0;i<primos.size();i++)
    printf("*%d",primos[i]);
  printf("\n");
  system("PAUSE");
} 

SAÍDA 

2*2*2*2*3*5
Pressione qualquer tecla para continuar. . .

Podemos provar que nenhum número composto será adicionado no vetor primos.

Teorema: O algoritmo fatora encontra a decomposição em fatores primos.

Prova: Suponha por absurdo que um número não-primo foi adicionado no vetor primos (vetor que guarda a decomposição em fatores primos). Seja d o menor número não-primo adicionado no vetor primos. Pelo Teorema Fundamental da Aritmética, d também pode ser decomposto em fatores primos. Considere a seguinte decomposição de d = p_1*p_2*..._*p_k onde p_i < d. Logo, existe um primo p_i que foi “pulado” durante o algoritmo uma vez que p_i < d. Absurdo, todos os números menores que d são testados antes de d.  

O pior caso desse algoritmo tem a seguinte complexidade:

A segurança do sistema RSA depende da dificuldade da fatoração de um número n em fatores primos. Vamos imaginar que queiramos quebrar a criptografia de um sistema RSA com uma chave de 1024 bits. Logo, n ~ 2¹⁰²⁴ ~ 2¹⁰⁰⁰~ (2¹⁰)¹⁰⁰ ~ (10³)¹⁰⁰~ 10³⁰⁰. Precisamos testar aproximadamente ~10¹⁵⁰ números. Considerando que seja possível testar 10⁹ fatores por segundo. Para testar todos os valores possíveis ~10¹⁴¹ segundos ~ 10¹³⁰ milênios.  
 

Programa em JavaScript que descobre a decomposição em fatores primos utilizando esse algoritmo

Teste de Primalidade

Na escola, aprendemos que os números primos são especiais: ele tem exatamente dois divisores 1 e n e todo número maior que 1 pode ser representado de maneira única por produto de números primos (chamado de fatores primos). Este resultado é conhecido como Teorema Fundamental da Aritmética. 

Na Grécia Antiga, Euclides já tinha percebido que os primos eram os blocos fundamentais, os átomos, que constituem o edifício dos números inteiros. Vários matemáticos tentaram encontrar uma fórmula para descobrir os primos, mas, até agora, nenhum matemático conseguiu a solução desse antigo quebra-cabeça. 

Euclides (300 AC) formulou e respondeu a seguinte pergunta sobre os primos:

Os números primos são finitos ou infinitos?

Para responder esta pergunta, ele utilizou um método de prova chamada redução ao absurdo.

Prova: Suponha por absurdo que a quantidade de números primos é finita e pode ser listada da seguinte maneira P = {p1,...,pk}. Seja N definido como o produto de todos os primos somando mais 1:

N = p1*...*pk + 1

Nenhum primo pode ser divisor de N, pois para todo pi, N % pi = 1. Como nenhum primo é divisor de N então N não pode ser dividido por nenhum primo. Portanto, N é também um número primo,  mas  N não pertence ao conjunto P. Absurdo!

Até onde eu preciso testar para saber se um número é primo ou não?

Se n não é primo, então o menor divisor de n, d_0, menor ou igual a .

Este resultado também pode provar utilizando a redução por absurdo.

Prova: Suponha por absurdo que n não é primo e o menor divisor de n d₀ maior que  .

d₀ > 

d₀² > n

d₀ > n/d₀

Temos que d₀ é um divisor (e fator primo) de n e podemos escrever n = d₀*k onde k  Z. Temos também n/d₀ é um divisor de n, uma vez que, podemos escrever n = (n/d₀)*d₀  e d₀  Z.

Absurdo, d₀ é o menor divisor de n e existe um outro divisor de n, n/d_0, menor que d_0.

#include <stdio.h>
#include <math.h>
#include <stdlib.h>

/*
funcao ehprimo(n) retorna 0, se n nao eh primo 
                          1, se n eh primo

se i eh divisor de n entao n nao eh primo

se nenhum numero 2<=i<=sqrt(n) eh divisor de n entao n eh primo 
*/
int ehprimo(int n){
  int limite;
  limite = (int)sqrt(n);
  for(int i=2;i<=limite;i++){
    if(n%i==0) return 0;
  }
  return 1;
}

int main(){
  int cont=0;  
  for(int i=2;i<=1000000;i++){
    if(ehprimo(i)) { 
      cont++;
    }

  }
  printf("numeros primos menores que 1000000: %d\n",cont);
  system("PAUSE");

}

Saída
numeros primos menores que 1000000: 78498
Pressione qualquer tecla para continuar. . .

Podemos conferir o resultado aqui:
http://www.wolframalpha.com/input/?i=number+of+prime+from+1+to+1000000


Gauss tentou responder uma pergunta interessante sobre os números primos:
Quantos primos existem até um certo número n?
Gauss era fascinado por um suplemento do seu livro de logaritmos que trazia uma tabela de números primos. Os logaritmos eram previsíveis mas os primos eram completamente aleatórios. Parecia não existir qualquer conexão entre os logaritmos e os primos. Gauss denotou a quantidade de números primos até x e observou essa estranha relação entre e log(x).



A fórmula que ele encontrou é uma boa aproximação para a quantidade de  números de primos até um certo número.

Confira aqui:
Comparação até 100
Comparação até 1000
Comparação até 10000
Comparação até 100000
Comparação até 1000000

Links relacionados:
http://elementosdeteixeira.blogspot.com/2012/02/o-desafio-dos-numeros-primos.html

Equações Diofantinas Lineares

Na matemática, uma equação diofantina linear é uma equação envolvendo soma de variável ou constante, que só podem assumir valores inteiros. As equações diofantinas possuem menos equações do que variáveis e a sua solução envolve descobrir números inteiros que satisfaçam as equações.

Problema 1 Um cachecol custa, na Rússia, 19 rublos, mas o caso é que o comprador só tem notas de 3, e o caixa, só de 5. Nessas condições, será possível pagar a importância da compra, e de que modo?

Podemos traduzir o seguinte problema para a seguinte equação diofantina 3x – 5y = 19 que possui infinitas soluções inteiras e positivas. A compra pode ser feita utilizando 8 notas de 3 rublos e recebendo 1 nota de 5 rublos de troco. Podemos descrever a solução da seguinte maneira:

3x = 19 + 5y 3x = 19 mod 5 3x = 4 mod 5 x = 8 + 5n , n ∈N

5y = 3x – 19 5y = -19 mod 3 5y = 2 mod 3 y = 1 + 3n , n ∈N

Exercício 1 Propõe-se a uma pessoa que multiplique a data do dia do seu nascimento, por 12, e o número que indica o mês correspondente, por 31. Com a soma desses produtos é possível calcula a data de aniversário da dita pessoa. Considere, por exemplo, se a pessoa nasceu em 08 de fevereiro:

8*12 = 96 , 2*31 = 62; 96 + 62 = 158.

Como podemos deduzir a data de nascimento da pessoa?

Problema 2 Quantas quadras de basquete e quantas de vôlei são necessárias para que 80 alunos joguem simultaneamente? E se forem 77? Sabendo que o basquete e vôlei são jogados, respectivamente, por duas equipes com 5 e 6 jogadores em cada uma.

A equação que descreve a solução desse problema é 10x + 12y = 80.

10x = 80 – 12y 10x = 80 mod 12 10x = 8 mod 12 x = 8 e y = 0

12y = 80 – 10x 12y = 0 mod 10 y = 5 , x = 2

A segunda equação 10x + 12y = 77 não tem soluções inteiras.

Problema 3: Para agrupar 13 aviões em filas de 3 ou de 5, quantas filas serão necessárias de cada tipo?

A equação diofantina que traduz o problema é 3x + 5y = 13.

3x = 13 – 5y 3x = 13 mod 5 3x = 3 mod 5 x=1, y=2

5y = 13 – 3x 5y = 13 mod 3 5y = 1 mod 3 y=2, x=1

Esta equação possui uma única solução.

Resultados de divisibilidade:

·         Se um número inteiro d|a, então d|am, para qualquer inteiro m;

·         Se d|a e d|b, então d|a+b

·    (Teorema de Bézout) Se mdc(a,b) = d, então existem inteiros m e n tais que d = am+bn.

O teorema de Bézout nos oferece um método para descobrir se uma equação diofantina tem ou não solução. A equação ax + by = c tem admite solução?

(1) Se c = mdc(a,b), a equação ax + by = c tem solução x = m e y = n do teorema de Bézout.

(2) Se c = dt, onde d = mdc(a,b) e t é um número inteiro, existem inteiros m e n tais que am+bn= d.  Assim,      c = dt = (am+bn)t = a(mt) + (bn)t , a solução será x = mt e y = nt.

Podemos enunciar o seguinte teorema:

Teorema Uma equação diofantina ax+by=c , em que a != 0 e b !=0, admite solução se, e somente se, d = mdc(a,b) | c.

Exemplo

10x + 12y = 77 não tem soluções inteiras.

d = mdc(10,12) = 2 não divide 77

10x+12y = 80 tem soluções inteiras.

d = mdc(10,12) = 2 | 80

Como encontrar os números inteiros m e n tais que d = am + bn, onde d = mdc(a,b)?

Um modo de encontrar esses números m e n é através do algoritmo de Euclides para o cálculo do mdc(a,b).

Se a e b são inteiros, com b > 0, existem q e r, com 0<=r<b tais que a = bq + r (algoritmo da divisão). Supondo que x|a e x|b então x | r = a - bq (combinação linear de a e b).  Da mesma maneira, se x|b e x|r então x| a = bq + r (combinação linear de b e r). Logo, mdc(a,b) se reduz a encontrar mdc(b,r).

Supondo que r_n seja o primeiro resto nulo temos:

mdc(a,b) = mdc(b,r₁) = mdc(r₁,r₂) = … = mdc(r_n-1,r_n) = r_n-1.

a = bq₁ + r₁ (r₁ = a-bq₁)

b = r₁q₂ + r₂ (r₂ = b – (a-bq₁)q₂ = b – aq₂ + bq₁q₂ )

….

r_{n-1 =} r_nq_n (mdc(a,b)=r_n-1 = am + bn , m e n inteiros)

Através de um processo de substituição, podemos encontrar o valor de m e n.

(1) 120 = 23*5 + 5

(2) 23  = 5*4  + 3

(3) 5   = 3*1  + 2

(4) 3   = 2*1  + 1

(5) 2   = 1*2  + 0

(1) 5 = 1*120 - 5*23

(2) 3 = 1*23 - 4*5 Substituindo o 5 temos

    3 = 1*23 - 4*(1*120 - 5*23)

    3 = -4*120 + 21*23

(3) 2 = 1*5 - 1*3 Substituindo o valor de 5 e 3 temos

    2 = 1(1*120 - 5*23) - 1(-4*120 + 21*23)

    2 = 5*120 - 26*23

(4) 1 = 1*3 - 1*2 Novamente substituindo 3 e 2

    1 = 1(-4*120 + 21*23) - 1(5*120 - 26*23)

    1 = -9*120 + 47*23

portanto, m = -9 e n = 47 e temos MDC(120,23) = 120 * ( − 9) + 47 * 23

Algoritmo de Euclides Estendido – Versão 1

int mdc(int a, int b, int *m, int *n) {

int mm, nn, d;

if(b==0) {

*m=1; *n=0;

return a;

}

d = mdc(b, a%b, &xx, &yy);

*m = nn;  *n = mm - a/b*nn;

return d;

 }

Algoritmo de Euclides Versão 2

int eMcd (int a, int b, int *m, int *n){

int x, yAnt, r, aIni, bIni, sr,q;

aIni = a; bIni = b;

x = 1; yAnt = 0;

while (b != 0){

r = a % b;

q=a/b;

a = b;

b = r;

sr = x - yAnt*q;

x = yAnt;

yAnt = sr;

}

*m = x;

*n = (a – x*aIni)/bIni;

return a;

}

Módulo para número negativos

int mod(int a, int n) {

            return (a%n + n)%n;

 }

Teorema1 Seja d = mdc(a,b) e ax’ + by’ = d então a equação ax = c mod b tem como uma solução x₀ = x’ (c/d) é uma solução da equação.

ax₀ = a (x’ (c/d)) mod b

   = ax’ (c/d) mod b

   = d (c/d) mod b

   = c

Teorema2 Suponha ax = b mod c tenha solução e x₀  seja uma solução. Então essa equação tem exatamente d soluções distintas, módulo n, dado por x_i = x₀ + i (b/d)

ax_i   = a (x₀ + i (b/d)) mod b

       = ax₀  +a i/d*b mod b

       = ax₀

Resolvendo o problema:

Problema 1

3x – 5y = 19

b = 5

a = mod(3,5) = 3

c = mod(19,5) = 4

3x = 4 (mod 5)

mdc(3,5) = 1

3x’ + 5y’ = 1

3(2) + 5*(-1) = 1

x₀ = x’ (c/d) = 2*4 = 8

Problema 2

10x + 12y = 80

a = mod(10,12) = 10

c = mod(80,12) = 8

10x + 12y = 8

10x = 8 mod 12

d = mdc(10,12) = 2

10(-1) + 12(1) = 2

x₀  = -1*(8/2) = -4

x₁ = x₀ + 1*(12/2) = -4 + 6 = 2