errore := function (arg)

## Riporta un messaggio di errore o informazione, senza ritornare
## Messaggio di default molto utile

local TESTO;

if Length(arg) > 0
then
TESTO := arg[1];
else
TESTO := "Very useful error message\n";
fi;
Print(TESTO, "\n");
return;

end;

errore("");
errore("giMCD (a, b) calcola verbosamente il MCD dei due interi di Gauss");
errore("a e b. In realtà funziona anche in qualsiasi dominio euclideo, tipo Z");
errore("");

giMCD := function (a, b)

local aorig, borig,
r, q, 
c11, c12, c21, c22, 
d1, d2;

c11 := 1; c12 := 0;
c21 := 0; c22 := 1;

aorig := a; borig := b;

while b <> 0 do
r := EuclideanRemainder(a, b);
q := EuclideanQuotient (a, b);
# Qui devo tener conto anche del quoziente

Print(a, " = (", b, ") * (", q, ") + (", r, ")\n");

a := b;
b := r;

d1:= c11 - c21 * q;
d2:= c12 - c22 * q;

c11:= c21; c12 := c22;
c21:= d1;  c22 := d2;

od;

Print("Il massimo comun divisore fra ", aorig, " e ", borig, 
	  " è ", a, ", e si ha:\n");

Print(a, " = (", aorig, ") * (", c11, ") + (", borig, ") * (", c12, ")\n");
# Si spiega il risultato

return [a,c11, c12];
#Qui restituisco sia il MCD a, che i coefficienti della combinazione lineare

end;

realPart := function(a)

return (a + ComplexConjugate(a))/2;

end;

imPart := function(a)

return E(4)*(-a + ComplexConjugate(a))/2;

end;

errore("");
errore("sqrtMinusOne(p) calcola una radice quadrata di -1 modulo p");
errore("Viene controllato che p sia primo, e congruo a 1 mod 4");
errore("");

sqrtMinusOne := function (p)

local x, tmp;

if not IsPrimeInt(p) then
errore("Sorry, argument is not prime");
return;
fi;
if p mod 4 > 1 then 
errore("Sorry, argument is prime but not 1 mod 4");
return; 
fi;

for x in Primes do

if PowerMod(x, (p-1)/2, p) > 1 then 
  tmp := PowerMod(x, (p-1)/4, p); 
  return [Minimum([tmp, p - tmp]),x];
else
tmp := PowerMod(x, (p-1)/4, p);
if tmp > 1 then tmp := -1; fi;
Print(x, "^(p-1)/4 is ", tmp, " mod p = ", p, ", not good\n");
fi;

od;

return;

end;

errore("");
errore("sumOfSquares(p) scrive il primo p come somma di due quadrati");
errore("verifica che p sia primo, e congruo a 1 mod 4");
errore("");

sumOfSquares := function (p)

local a, b, c, ra, ia;

if not IsPrimeInt(p) then
errore("Sorry, argument is not prime");
return;
fi;
if p mod 4 > 1 then 
errore("Sorry, argument is prime but not 1 mod 4");
return; 
fi;

c := sqrtMinusOne(p);
b := c[1]; c := c[2];

Print("Square root of -1 mod ", p, " is ", b, ",\n");
Print("Obtained as +- ", c, "^(", (p-1)/4, ")\n");

a:= giMCD(p, b + E(4))[1];

ra := realPart(a); ia := imPart(a);

Print(p, " = ", ra, "^2 + ", ia, "^2\n");

return [ra, ia];

end;

errore("");
errore("NextPrimeIntSame(p) trova il prossimo primo che sia nella");
errore("stessa classe di congruenza modulo 4 del numero dispari p");
errore("");

NextPrimeIntSame := function(p)
local pr;
if p mod 2 = 0 then
errore("p deve essere dispari");
return;
fi;
pr := p;
repeat pr := NextPrimeInt(pr);
until (p mod 4) = (pr mod 4);
return pr;
end;