Ce code permet de résoudre l'équation x^2 + y^2 = p, p premier. Très pratique en cryptographie ;o)
Source / Exemple :
/* CORNACCHIA Algorithm
*
- GOAL: given d and p prime, find (x,y) such that x^2 + y^2 = d * p
*
- in this implementation, d = 1. This code does not accept d <> 1 !!!
- see "A Course in Computational Algebraic Number Theory" by Henri Cohen
- coded by malik@hammoutene.com, 2004,the 5th of August
*
#pragma comment(lib, "gmp.lib")
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include "gmp.h"
#define _WIN32_WINNT 0x0400
#include <windows.h>
#include <wincrypt.h>
int genPrime(){ // GENERATION DE NOMBRE PREMIER
mpz_t t,t2, prime_n;
int d2;
d2=getrand();
mpz_init(t);
mpz_set_ui(t,d2);
mpz_mul_ui(t,t,4);
mpz_add_ui(t,t,1);
if (mpz_probab_prime_p(t,20)==0){
while (mpz_probab_prime_p(t,20)==0){
d2=getrand();
mpz_set_ui(t,d2);
mpz_mul_ui(t,t,4);
mpz_add_ui(t,t,1);
}
}
mpz_init(t2);
mpz_set(t2,t);
mpz_set(t,t2); // petite finte pour le return
mpz_init(prime_n);
mpz_set(prime_n,t);
mpz_clear(t);
mpz_clear(t2);
return (prime_n);
}
int getrand() // GENERATION DE NOMBRE ALEATOIRE
{
HCRYPTPROV hProv = 0;
int rnd;
CryptAcquireContext(&hProv, 0, 0, PROV_RSA_FULL, CRYPT_VERIFYCONTEXT);
CryptGenRandom(hProv, sizeof(rnd), (BYTE*)&rnd);
CryptReleaseContext(hProv, 0);
return abs(rnd/1000); // pour être dans un ordre honnête (~ 2^20 et plus)
}
int power(int val, int pow)
{ int ret_val = 1;
int i;
for(i = 0; i < pow; i++)
ret_val *= val;
return(ret_val);
}
//************* SHANKS-TONELLI Algorithm
int shanks(mpz_t a2, mpz_t prime_number){
int random_n;
mpz_t a,a_prime,b,x,y,z,q,t,temp,temp2;
mpz_t rand_n;
int e = 0;
int compteur=1;
int r,e2,m,mp,r2;
mpz_init(a);
mpz_sub_ui(a,prime_number,1);
mpz_init(q);
mpz_init(a_prime);
mpz_mod_ui(q,a,2);
while(mpz_cmp_ui(q,0) == 0)
{
mpz_div_ui(a_prime, a, 2);
e = e+1;
mpz_set(a,a_prime);
mpz_mod_ui(q,a,2);
}
e2 = power(2,e);
mpz_div_ui(q, prime_number, e2);
// 1 Find generator
random_n = getrand();
mpz_init(rand_n);
mpz_set_ui(rand_n,random_n);
if (mpz_legendre(rand_n, prime_number) != -1){
while (mpz_legendre(rand_n, prime_number) != -1){
random_n = getrand();
mpz_set_ui(rand_n,random_n);
}
}
mpz_init(z);
mpz_powm(z,rand_n,q,prime_number);
// 2 Initialize
mpz_init(y);
mpz_set(y,z);
r = e;
mpz_init(temp);
mpz_sub_ui(temp, q, 1);
mpz_div_ui(temp, temp, 2);
mpz_init(x);
mpz_powm(x,a2,temp,prime_number);
mpz_pow_ui(temp,x,2);
mpz_mul(temp, temp, a2);
mpz_init(b);
mpz_mod(b,temp,prime_number);
mpz_mul(temp, a2, x);
mpz_mod(x,temp,prime_number);
// 3 Find exponent
mpz_mod(temp, b, prime_number);
mpz_init(t);
if (mpz_cmp_ui(temp,1) != 0){
mpz_init(temp2);
while(mpz_cmp_ui(temp,1) != 0){
m=1;
mp=power(2,m);
mpz_powm_ui(temp2,b,mp,prime_number);
while(mpz_cmp_ui(temp2,1) !=0)
{
m++;
mp=power(2,m);
mpz_powm_ui(temp2,b,mp,prime_number);
}
if (m==r){
printf("a is not a quadratic residue mod p...\n");
exit(0);}
// 4 Reduce exponent
r2 = r-m-1;
r2 = power(2,r2);
mpz_powm_ui(t,y,r2,prime_number);
mpz_powm_ui(y,t,2,prime_number);
r=m;
mpz_mul(x,x,t);
mpz_mod(x,x,prime_number);
mpz_mul(b,b,y);
mpz_mod(b,b,prime_number);
mpz_mod(temp, b, prime_number);
}
}
mpz_clear(a);
mpz_clear(a2);
mpz_clear(a_prime);
mpz_clear(b);
mpz_clear(rand_n);
mpz_clear(y);
mpz_clear(z);
mpz_clear(q);
mpz_clear(t);
mpz_clear(temp);
mpz_clear(temp2);
return (x);
}
struct cornAlgo{
mpz_t x_sol;
mpz_t y_sol;
};
struct cornAlgo cornacchia(mpz_t prime){
struct cornAlgo resultats;
int d = 1;
int k2;
mpz_t temp,temp3,a2,a,x01,b,l,c,r;
mpz_t prime_number;
mpz_init(prime_number);
mpz_set(prime_number,prime);
//************* CORNACCHIA ALGORITHM
// 1 Test if residue
mpz_init(temp);
mpz_set_ui(temp,d);
mpz_neg(temp,temp);
k2= mpz_legendre(temp,prime_number);
mpz_clear(temp);
if (k2 == -1) {
printf("the equation has no solution\n");
exit (0);}
// 2 Compute square root
mpz_init(a2);
mpz_set_ui(a2,d);
mpz_neg(a2,a2);
mpz_init(x01);
mpz_set(x01, shanks(a2,prime_number));
if (mpz_cmp(x01,prime_number)>0){
mpz_neg(x01,x01);
}
mpz_init(temp3);
mpz_div_ui(temp3, prime_number,2);
mpz_add_ui(temp3, temp3, 1);
if (mpz_cmp(x01,temp3)>=0){
while (!(mpz_cmp(x01,temp3)>=0)){
mpz_add(x01, x01, prime_number);
}
}
mpz_init(a);
mpz_set(a,prime_number);
mpz_init(b);
mpz_set(b,x01);
mpz_init(l);
mpz_sqrt(l,prime_number);
// 3 Euclidean algorithm
mpz_init(r);
if (mpz_cmp(b,l)>0){
while ((mpz_cmp(b,l)>0)){
mpz_mod(r,a,b);
mpz_set(a,b);
mpz_set(b,r);
}
}
// 4 Test solution: ici, comme d = 1, les tests ne sont pas utiles
mpz_pow_ui(temp3,b,2);
mpz_neg(temp3,temp3);
mpz_add(temp3,temp3,prime_number);
mpz_init(c);
mpz_set(c, temp3);
mpz_sqrt(c,c);
mpz_init(resultats.x_sol);
mpz_set(resultats.x_sol,b);
mpz_init(resultats.y_sol);
mpz_set(resultats.y_sol,c);
mpz_clear(temp3);
mpz_clear(a);
mpz_clear(x01);
mpz_clear(l);
mpz_clear(r);
mpz_clear(prime_number);
return resultats;
}
//************* MAIN
main(int argc, char *argv[])
{
struct cornAlgo resultatsXY;
mpz_t prime_number;
// génération d'un nombre premier tel que p = 1 + 4k
mpz_init(prime_number);
mpz_set(prime_number,genPrime());
// exécution de l'algorithme de Cornacchia
resultatsXY = cornacchia(prime_number);
// affichage du résultat
mpz_out_str (stdout, 10, prime_number); printf(" = ");
mpz_out_str (stdout, 10, resultatsXY.x_sol); printf("^2 + ");
mpz_out_str (stdout, 10, resultatsXY.y_sol); printf("^2\n");
}
Conclusion :
N'oubliez pas d'installer GMP pour faire fonctionner ce code (j'ai mis une marche à suivre dans le forum).
Merci de me tenir au courant en cas d'amélioration (par exemple, si quelqu'un traite le cas d>1)
Vous n'êtes pas encore membre ?
inscrivez-vous, c'est gratuit et ça prend moins d'une minute !
Les membres obtiennent plus de réponses que les utilisateurs anonymes.
Le fait d'être membre vous permet d'avoir un suivi détaillé de vos demandes et codes sources.
Le fait d'être membre vous permet d'avoir des options supplémentaires.