Ci-dessous un exemple de comment écrire un premier p (qui est 1 modulo
4) comme une somme de deux carrés, en utilisant Maple. Le logiciel Maple
n'est pas disponible gratuitement, mais est installé sur les machines
des travaux pratiques.

Le lecteur constatera que l'utilisation semble peut-être facile, mais
n'est pas du tout mathématiquement rigoureuse. Par exemple, les objets
ne sont pas typés: on ne sait pas ce que sont les choses que l'on
manipule. Cela fait un grand bordel. Malheureusement, les logiciels
utilisés pour l'épreuve de modélisation de l'Agrégation sont Maple et
Mupad (qui lui est gratuit, mais pas `libre' dans le sens que la
source n'est pas disponible). 

Les commentaires dans les lignes qui suivent commencent par ';', et
ont étés ajoutées après le calcul. Le `prompt' de Maple est `>'. 



> restart;

; Malheureusement, ceci est souvent nécessaire. 

> p:=nextprime(10^10);

; Ceci calcule ne premier nombre premier après 10^10. 

                           p := 10000000019

; Notons que cela marche tres vite. Ce p est -1 mod 4. On va prendre des
; p plus grands.  

> p:=nextprime(10^100);

  p := 10000000000000000000000000000000000000000000000000000000000\
        000000000000000000000000000000000000000267

; Cela prend plus de temps, mais c'est encore raisonnable. D'ailleurs,
; maple ne donne pas de preuve que ce p est premier (voir la commande
; isprime).  De toute façon, ce p est -1 mod 4.

> p:=nextprime(p);

  p := 10000000000000000000000000000000000000000000000000000000000\
        000000000000000000000000000000000000000949

; Maintenant p est 1 mod 4. Calculons la partie impaire de p-1. 

> m:=(p-1)/4;

  m := 25000000000000000000000000000000000000000000000000000000000\
        00000000000000000000000000000000000000237

; On a donc p-1 = 2^2*m avec m impair.  Prenons un a dans F_p au hasard,
; et calculons a^m.

> a:=2; b:=Power(a,m) mod p;

                                a := 2


  b := 54205298266425948739380045800247485111179100381642542689455\
        30069756410577665995572552521420196783014

> b1:=b^2 mod p;

  b1 := 1000000000000000000000000000000000000000000000000000000000\
        0000000000000000000000000000000000000000948

; Notons que b1 =-1 dans F_p.  On a donc réussi en un seul coup: b est
; une racine de -1 dans F_p.

> with(GaussInt):
> z:=GIgcd(p,b-I);

  z := 7766881905507050845172598218029833369440123277895

         + 99697921470138519447541656418848509184628524016382 I

; Les coordonnées sont u et v tels que u^2+v^2=p.  Voyons comment
; extraire u et v de z.

> nops(z);

                                  2

> u:=op(1,z);

        u := 7766881905507050845172598218029833369440123277895

> v:=op(2,z);

      v := 99697921470138519447541656418848509184628524016382 I

> nops(v);

                                  2

> w:=op(1,v);

       w := 99697921470138519447541656418848509184628524016382

> u^2+w^2;

  1000000000000000000000000000000000000000000000000000000000000000\
        0000000000000000000000000000000000949

> u^2+w^2-p;

                                  0

; Il existe aussi d'autres logiciels de calcul scientifique. Par
; exemple: mupad, pari.