Vendredi 17 nov. 2023

Entiers de Gauss

Prétexte pour revisiter C++ que je ne pratique plus depuis quelques années et réviser un peu d'algèbre que je suis en train d'oublier. Un entier de Gauss est un nombre complexe dont la partie réelle et la partie imaginaire sont des entiers relatifs.

Introduction

Cette partie est largement inspirée par le papier de Vincent Lefèvre et le cours de P. Caldero.

Définition

L'ensemble des entiers de Gauss est l'ensemble $\mathbb{Z}[i]=\{a+bi\,|\,(a,b)\in\mathbb{Z}^2\}$.

Anneau commutatif

$\mathbb{Z}[i]$ est un anneau commutatif unitaire. (Très facile à établir : $\mathbb{Z}[i]\subset\mathbb{C}$)

Norme

La norme algébrique dans $\mathbb{Z}[i]$ est définie par $N(a+bi)=a^2+b^2$.

Propriété $N(zz')=N(z)N(z')$. (Preuve : dans ℂ, $|zz'|^2=|z|^2+|z'|^2, ...)

Inversibles

Les éléments inversibles de $\mathbb{Z}[i]$ sont $1, -1, i, -i$.

Preuve : $N(z.z^{-1})=1~\Rightarrow~N(z).N(z^{-1})=1$, ...

Irréductibles

Si $x\in\mathbb{Z}[i]$ et $p=N(x)$ est premier, alors $x$ est irréductible.

Soit $x=uv$, et $p=N(x)=N(u)N(v)$ entier naturel premier, alors $N(u)$ ou $N(v)$ est 1, et du coup $u$ ou $v$ est inversible.

Si $x$ n'est ni réel pur ni imaginaire pur, la réciproque est vraie.

Remarque : $x$ irréductible signifie $x=uv~\Rightarrow~u$ ou $v$ inversible.

Division euclidienne

Il s'agit de trouver $q$ et $r$ tels que $x=y.q+r$ et $N(r)\ltN(y)$. En fait on peut faire mieux et faire en sorte que $N(r)\le\dfrac{N(y)}{2}$. Malgrè cette inéquation le couple $(q,r)$ n'est pas unique. On peut cependant faire que ce choix soit plus fin, en faisant l'arrondi le plus proche(round()).

Dans ℂ on calcule $\dfracxy=\alpha+\betai$ où $\alpha$ et $\beta$ ne sont pas nécessairement entiers. On prend $a$ et $b$ entiers relatifs tels que $|a-\alpha|\le\dfrac{1}{2}$ et $|b-\beta|\le\dfrac{1}{2}$. D'où $\dfrac{N(r)}{N(y)}\le\dfrac{1}{4}+\dfrac{1}{4}=\dfrac{1}{2}$. (Comme me disaient mes élèves : ça se voit sur la figure)

ℤ[i] est un anneau principal

Soit $I$ un idéal de $\mathbb{Z}[i]$. On peut choisir un élément $m\ne0$ de $I$ de norme minimale.

$m\mathbb{Z}[i]\subsetI$ est claire.

Montrons que $I\subsetm\mathbb{Z}[i]$ : Pour tout $x\inI$ on considère le reste $r$ de la division de $x$ par $m$. Alors $r\inI$ et sa norme est inférieure à celle de $m$. Donc $N(r)=0$ et donc $r=0$. D'où $x$ est multiple de $m$ et donc $x\inm\mathbb{Z}[i]$.

La classe Gi

Squelette

Comme c'est pour jouer, j'ai défini une struct plutôt qu'une class, ainsi tout est public par défaut. L'identifiant de type Gi est un entier de Gauss (Gaussian integer).

#include <string>
#include <iostream>
using namespace std;
struct Gi
{
   int re, im; // partie réelle, partie imaginaire
   Gi (int a = 0, int b = 0) : re(a), im(b) {}; // constructeur par défaut
};
// affichage
ostream& operator << (ostream& os, Gi z) {
   os << z.re << showpos << z.im << "i";
   return os;
}
// addition (a + bi) + (a' + b'i) = a + a' + (b + b')i
Gi operator + (const Gi& x, const Gi& y)
{
   return Gi(x.re + y.re, x.im + y.im);
}
// main
int main(int argc, char * argv[]) 
{
   Gi x(3, 2), y(-1, 2);
   cout << "x = " << x << endl;
   cout << "y = " << y << endl;
   cout << "x + y = " << x + y << endl;
}

Résultat :

x = 3+2i
y = -1+2i
x + y = +2+4i

Erreurs

#define NAI 2147483647

NAI est un alias de INT_MAX définie dans /usr/include/limits.h qui correspond à la valeur maximale qu'un int peut contenir. C'est une sorte de not a number.

const Gi naGi(NAI,NAI);

naGi (not a Gi) sera le retour d'une fonction qui renvoie un Gi lorsqu'elle rencontre une erreur de dépassement de capacité ou de division par zéro, etc. C'est là aussi une sorte de n'est pas un entier de Gauss.

Cosmétique

Je me suis employé à singer les entrées-sorties de Giac en console, en particulier pour l'utiliser pour contrôler mes calculs. Voici la fonction membre Str() et l'opérateur << : 

ostream& operator << (ostream& os, const Gi& g)
{
   if(g.isnan()) return os << "nan";
   if(g == Gi(0, 0)) return os << "0";
   if(g.re) os << g.re;
   if(!g.im) return os;
   if(g.im == 1) return g.re ? os << "+i" : os << "i";
   if(g.im == -1) return os << "-i";
   if(g.re != 0 && g.im > 0) os << "+";
   return os << g.im << "i";
}
string Gi::str() const
{
   ostringstream os;
   os << *this;
   return os.str();
}

Et pour faire bonne mesure j'ai ajouté un constructeur qui prend comme paramètre une string, qui utilise la même syntaxe.

La fonction stoGi reconnaît des chaînes comme : "2", "i", "-i", "3i", "-3i", "4*i", "1+2i". Elle retourne un Gi ou naGi en cas d'erreur de syntaxe.

Le constructeur Gi(string) utilise la fonction stoGi pour instancier un Gi.

Gi stoGi(string s)
{
   if(s.empty()) return naGi;
   if(s == "+i" || s == "i") return Gi(0, 1); // or return i;
   if(s == "-i") return Gi(0, -1); // or return -i;
   istringstream is(s);
   int v;
   string tail;
   Gi next;
   if(is >> v) {
      char c = is.peek();
      if(c == 'i') return Gi(0, v);
      else if(c == '+' || c == '-') {
         is >> tail;
         return Gi(v, 0) += stoGi(tail);
      } else return Gi(v, 0);
   }
   else return naGi;
}
Gi::Gi(string s)
{
   Gi g = stoGi(s);
   re = g.re;
   im = g.im;
}

opérateur de cast en complex

L'opérateur de cast en complex<int> a pour vocation d'utiliser la librairie numérique complex pour vérifier les calculs.

   operator complex<int>() const { return complex<int>(re, im); }

Remarque : Normalement le patron complex accepte les types float, double et long double. J'ai tenté complex<int> et gcc n'a pas bronché.

Classe Gi

Voici l'état de la classe Gi avec les ajouts expliqués ci-dessus.

Je l'ai complétée par la surcharge des opérateurs - et * et deux fonctions membres :

/* Gi.cc 
:w | !g++ % -c -Wall
(C) 2023 by marnout à free pt fr  Released under the GPL
*/
#ifndef MRD_GINT_H
#define MRD_GINT_H 1
#include <string>
#include <iostream>
#include <complex>
#define NAI 2147483647
using namespace std;
struct Gi
{
   int re, im;
   Gi(int a = 0, int b = 0) : re(a), im(b) {};
   Gi(string s); 
   bool isnan() const { return re == NAI || im == NAI; }
   int N() const { return isnan() ? NAI : re*re + im*im; }
   Gi conj() const;
   string str() const; 
   Gi& operator +=(const Gi& g) ;
   bool operator == (const Gi& g) const { return re == g.re && im == g.im; }
   operator complex<int>() const { return complex<int>(re, im); }
};
const Gi naGi(NAI,NAI);
const Gi i(0, 1);
Gi
Gi::conj() const  { 
   return isnan() ? naGi :  Gi(re, -im); 
}
Gi
operator -(const Gi& g)
{
   if(g.isnan()) return naGi;
   return Gi(-g.re, -g.im);
}
Gi
operator + (const Gi& g1, const Gi& g2)
{
   if(g1.isnan() || g2.isnan()) return naGi;
   return Gi(g1.re + g2.re, g1.im + g2.im);
}
Gi
operator - (const Gi& g1, const Gi& g2)
{
   if(g1.isnan() || g2.isnan()) return naGi;
   return Gi(g1.re - g2.re, g1.im - g2.im);
}
Gi
operator * (const Gi& g1, const Gi& g2)
{
   if(g1.isnan() || g2.isnan()) return naGi;
   return Gi(g1.re*g2.re - g1.im*g2.im , g1.re*g2.im + g1.im*g2.re);
}
Gi&
Gi::operator += (const Gi& g2)
{
   if(isnan() || g2.isnan()) { re = NAI; im = NAI; }
   else { re += g2.re; im += g2.im; }
   return *this;
}
Gi 
operator / (const Gi& g1, const Gi& g2)
{   // !naGi && !0
   if(g1.isnan() || g2.isnan() || (g2.re == 0 && g2.im == 0)) return naGi;
   Gi num = g1*(g2.conj());
   int denom = g2.N();
   return Gi(round((float)num.re/denom), round((float)num.im/denom));
}
Gi 
operator % (const Gi& g1, const Gi& g2)
{
   return g1 - g2*(g1/g2);
}
ostream&
operator << (ostream& os, const Gi& g)
{
   if(g.isnan()) return os << "nan";
   if(g == Gi(0, 0)) return os << "0";
   if(g.re) os << g.re;
   if(!g.im) return os;
   if(g.im == 1) return g.re ? os << "+i" : os << "i";
   if(g.im == -1) return os << "-i";
   if(g.re != 0 && g.im > 0) os << "+";
   return os << g.im << "i";
}
string
Gi::str() const
{
   ostringstream os;
   os << *this;
   return os.str();
}
Gi stoGi(string s)
{
   if(s.empty()) return naGi;
   if(s == "+i" || s == "i") return Gi(0, 1); // or return i;
   if(s == "-i") return Gi(0, -1); // or return -i;
   istringstream is(s);
   int v;
   string tail;
   Gi next;
   if(is >> v) {
      char c = is.peek();
      if(c == 'i') return Gi(0, v);
      else if(c == '+' || c == '-') {
         is >> tail;
         return Gi(v, 0) += stoGi(tail);
      } else return Gi(v, 0);
   }
   else return naGi;
}
Gi::Gi(string s)
{
   Gi g = stoGi(s);
   re = g.re;
   im = g.im;
}
Gi 
gcd ( Gi a, Gi b )
{
  Gi c;
  while (a.re || a.im) {
     c = a; 
     a = b%a;  
     b = c;
  }
  return b;
}
#endif

Division euclidienne

Ce qui fait l'intérêt des entiers de Gauss et qui justifie qu'on puisse les appeler «entiers» c'est la division euclidienne. Par définition $q$ et $r$ sont respectivement le quotient et le reste de la division euclidienne de $a$ par $b$ si $a=bq+r$ et $N(r)<N(b)$, où N est la norme algébrique. La norme est aux entiers de Gauss ce que la valeur absolue est aux entiers relatifs.

La méthode que j'ai utilisée consiste à calculer le quotient $\dfrac{a}{b}=\dfrac{a\bar{b}}{b\bar{b}}=\dfrac{b\barb}{N(b)}$ et à choisir l'entier de Gauss le plus proche que me donne la fonction round.

Cette division n'est possible que si $b\ne0$. 

Gi operator / (const Gi& g1, const Gi& g2)
{	// !naGi && !0
   if(g1.isnan() || g2.isnan() || (g2.re == 0 && g2.im == 0)) return naGi;
   Gi num = g1*(g2.conj());
   int denom = g2.N();
   return Gi(round((float)num.re/denom), round((float)num.im/denom));
}
Gi operator % (const Gi& g1, const Gi& g2)
{
   return g1 - g2*(g1/g2);
}

Avec ces notations l'opérateur / donne le quotient et \% le reste.

Utilisation de la classe

Pour tester la classe Gi, j'ai introduit une fonction alea() qui exploite la toute nouvelle fonction randint, qui n'a encore que le statut experimental et qui permet de de délivrer un entier au hasard entre les bornes qui lui sont fournies. (Cf. les lignes 7 et 12 du programme ci-dessous).

Pour vérifier les opérateurs + , - et *   j'utilise l'opérateur de cast de Gi en complex : affichage cmplx: à la suite du résultat.

Enfin, toujours pour vérifier les résultats des opérateurs / et % j'utilise l'excellente librairie de Giac-xcas, lignes 8, 9, 11 et de 33 à 40. Ceci me permet les affichages Giac.iquo: et Giac.irem.

Remarque il faut compiler avec la ligne : g++ main.cc -o main -lgiac -lgmp.

Voici maintenant un programme qui utilise la classe Gi :

/* gi.cc
:w | !g++ % -o .%< -lgiac -lgmp -Wfatal-errors -Wall
!./.%<
(C) 2023 by marnout à free pt fr  Released under the GPL
*/
#include "Gi.cc"
#include <experimental/random>
#include <giac/config.h>
#include <giac/giac.h>
using namespace std;
using namespace giac;
int alea() { return experimental::randint(-8, 9); }
string giac_div(const string s, const string s1, const string s2);
int main(int argc, char* argv[])
{
   Gi g1, g2;
   if(argc > 2) { 
      g1 = Gi(argv[1]); g2 = Gi(argv[2]); 
   } else {
      g1 = Gi(alea(), alea()); g2 = Gi(alea(), alea()); 
   }
   cout << "Entiers de Gauss\n";
   cout << "g1 := " << g1 << " ; g2 := " << g2 << endl;
   cout << "(" << g1 << ") + (" << g2 << ") = " << g1 + g2;
   cout <<"\tcmplx: " << complex<int>(g1 + g2) << endl;
   cout << "(" << g1 << ") - (" << g2 << ") = " << g1 - g2;
   cout <<"\tcmplx: " << complex<int>(g1 - g2) << endl;
   cout << "(" << g1 << ")(" << g2 << ") = " << g1 * g2;
   cout <<"\tcmplx: " << complex<int>(g1 * g2) << endl;
   cout << "(" << g1 << ")/(" << g2 << ") = " << g1 / g2;
   cout << "\tGiac.iquo: " << giac_div("iquo", g1.str(), g2.str()) << endl;
   cout << "(" << g1 << ")%(" << g2 << ") = " << g1 % g2;
   cout << "\tGiac.irem: " << giac_div("irem", g1.str(), g2.str()) << endl;
   cout << "gcd(" << g1 << ", " << g2 << ") = " << gcd(g1, g2);
   cout << "\tGiac.igcd: " << giac_div("igcd", g1.str(), g2.str()) << endl;
}
string
giac_div(const string s, const string s1, const string s2)
{
   context ct;
   gen e(s + "(" + s1 + ", " + s2 + ")", &ct);
   e = eval(e, 1, &ct);
   return print(e, &ct);
}

Et un exemple de sortie

Appuyez sur ENTRÉE ou tapez une commande pour continuer
Entiers de Gauss
g1 := -2-2i ; g2 := 7
(-2-2i) + (7) = 5-2i	cmplx: (5,-2)
(-2-2i) - (7) = -9-2i	cmplx: (-9,-2)
(-2-2i)(7) = -14-14i	cmplx: (-14,-14)
(-2-2i)/(7) = 0	Giac.iquo: 0
(-2-2i)%(7) = -2-2i	Giac.irem: -2-2*i
gcd(-2-2i, 7) = -1	Giac.igcd: 1

Au fil de mes tests, je suis tombé sur la sortie suivante :

Entiers de Gauss
g1 := 4+4i ; g2 := 5
(4+4i) + (5) = 9+4i	cmplx: (9,4)
(4+4i) - (5) = -1+4i	cmplx: (-1,4)
(4+4i)(5) = 20+20i	cmplx: (20,20)
(4+4i)/(5) = 1+i	Giac.iquo: 0
(4+4i)%(5) = -1-i	Giac.irem: 4+4*i
gcd(4+4i, 5) = 1	Giac.igcd: 1

Les résultats pour la division euclidienne ne concordent pas. Mystère.

P.S. J'ai déjà traité le pgcd. Il me reste au moins Bezout. Il y a longtemps que je n'ai pas programmé du Euclide-Bezout, la dernière fois c'était dans les années 90, et c'était en pascal.

Entiers quadratiques

Définition

Nous dirons qu'un nombre réel $x$ est un entier quadratique s'il existe trois entiers relatifs $a$, $b$ et $n$ tels que $x=a+b\sqrt{n}$.

Pour $n\in\mathbb{Z}$ fixé, $\mathbb{Z}(n)=\{ a+b\sqrt{n}\,|\,(a,b)\in\mathbb{Z}^2\}$ est l'ensemble des entiers quadratiques. Ce n'est pas la définition académique.

Moyennant certaines restrictions sur $n$, les entiers quadratiques peuvent avoir des propriétés communes avec les entiers de Gauss. Je ne connais pas ces restrictions (sans facteur carré ?). Toujours est-il, voici deux sources C++ qui mettent en œuvre ces entiers.

La classe Z_n

/* Z_n.cc
:w | !g++ -c % -Wall
Entiers quadratiques, entiers de Gauss
(C) 2023 by marnout à free pt fr  Released under the GPL
*/
#ifndef Z_N_CC
#define Z_N_CC
#ifndef N2
#define N2 2
#endif
#include <string>
#include <iostream>
#include <sstream>
#define NAI 2147483647
using namespace std;
// Entiers quadratiques Z_n : {a + b√n | a ∈ ℤ, b ∈ ℤ}
// n ∈ ℤ, n != 1 , ∀p ∈ ℕ, p² ∤ n
class Z_n
{
   static const int n2 = N2; // n²
public:
   int a, b;
   Z_n(int p = 0, int q = 0) : a(p), b(q) {} // constructeur
   Z_n(string);
   string str(); // pour affichage
   bool badZ_n() const { return a == NAI || b == NAI; }
   // norme N = a² + b² non √(a² + b²) !
   int N() const { return badZ_n() ? NAI : a*a + b*b; } 
   friend Z_n operator + (const Z_n u, const Z_n v);
   friend Z_n operator - (const Z_n u, const Z_n v); 
   friend Z_n operator * (const Z_n u, const Z_n v);
};
Z_n naZ_n(NAI, NAI);
// str : a ± b√n 
string
Z_n::str()
{
   if(badZ_n()) return "nan";
   ostringstream sstr;
   sstr << a;
   if(b >= 0) sstr << "+";
   sstr << b; 
   sstr << "√" << Z_n::n2;
   return sstr.str();
}
Z_n operator + (const Z_n u, const Z_n v) 
{
   if(u.badZ_n() || v.badZ_n()) return naZ_n;
   return Z_n(u.a + v.a, u.b + v.b);
}
Z_n operator - (const Z_n u, const Z_n v) 
{
   if(u.badZ_n() || v.badZ_n()) return naZ_n;
   return Z_n(u.a - v.a, u.b - v.b);
}
Z_n operator * (const Z_n u, const Z_n v)
{
   if(u.badZ_n() || v.badZ_n()) return naZ_n;
   return Z_n(u.a*v.a + Z_n::n2*u.b*v.b, u.a*v.b + u.b*v.a);
}
Z_n::Z_n(string s)
{
   istringstream is(s);
   string tail;
   is >> a >> b >> tail;
   if(s.empty() || tail != "√"+to_string(n2)) a = b = NAI;
}
#endif

Programme utilisateur de Z_n.cc

/* z_n.cc
:w | !g++ % -o .%< -Wfatal-errors -Wall
!./.%<
(C) 2023 by marnout à free pt fr  Released under the GPL
*/
#include "Z_n.cc"
#include <experimental/random>
using namespace std;
int alea() { return experimental::randint(-4, 5); }
int main(int argc, char * argv[]) 
{
	Z_n x, y;
	if(argc > 2) { 
		x = Z_n(argv[1]); 
		y = Z_n(argv[2]); 
	} else {
		x = Z_n(alea(), alea()); y = Z_n(alea(), alea()); 
	}
	cout << "Entiers quadratiques\n";
	cout << "x := " << x.str() << endl;
	cout << "y := " << y.str() << endl;
	cout << "x + y = " << (x + y).str() << endl;
	cout << "x - y = " << (x - y ).str() << endl;
	cout << "xy = " << (x * y).str() << endl;
}


Réalisé avec Qlam - LGPL