/* Abril 28. Estruturas algébricas.
Neste file programarmos divisão euclideano dos inteiros de Gauss
(div e div2 - div2 é mais em espirito de Euclides)
e tambem MDC (maximo divisor comum).

As seguintes funções são parte deste programa:
piso, div, div2, mdc

Estou usando só as seguintes 4 funções de PARI/GP:
1. floor(x)	- piso de número real x
   round(x)	- arredondamento de número real x, igual piso(x+1/2)
2. real(z)	- parte real de número complexo z
   imag(z)	- parte imaginária de número complexo z
   norm(z)	- norma de número complexo z 

Também as seguintes funções implicitas (constantes e operações):
1. a = b	- colocar valor b para variavel a
2. a == b	- igual 1 se a e b são iguais, 0 se diferentes
3. a+b		- soma
4. a-b		- diferença
5. a*b		- produto
6. a/b		- racio 
7. I		- raiz de -1, número complexo
8. [1,-1,I,-I]	- vetor de 4 números complexos (cumprimento 4)
   [d,r]	- vetor de cumprimento 2 de d e r
9. #		- cardinalidade (cumprimento de vetor)

Tente ??+ em PARI para ver a seguinte página sobre operações:
	https://pari.math.u-bordeaux.fr/dochtml/html/operators.html
Também "Standard monadic or dyadic operators"

As seguintes funções de programação:
1. return(x)	- saída de função será igual x
2. if(a,b)	- se a é verdade (não igual 0) então fazer b
3. for		- por

Programming in GP: control statements
"?return
return({x=0}): return from current subroutine with result x.

?if
if(a,{seq1},{seq2}): if a is nonzero, seq1 is evaluated, otherwise 
seq2. seq1 and seq2 are optional, and if seq2 is omitted, the preceding 
comma can be omitted also.

?for
for(X=a,b,seq): the sequence is evaluated, X going from a up to b. If b 
is set to +oo, the loop will not stop.
"

*/

\\ Piso (parte inteiro de baixo) de número gaussiano
piso(z) = (floor(real(z))+I*floor(imag(z)));
myround(z) = (round(real(z))+I*round(imag(z)));
arredondamento(z) = piso(z+(1+I)/2);
/* 

div(a,b)

Dividir um número gaussiano (a) doutro (b) com resíduo 
A entrada e um par dos números complexos (a) e (b).
A saída é um vetor [d,r] tal que 
   a = b*d + r 
 e   
   norm(r) < norm(b)   
 por cada   b != 0.

NB. Se usar 'piso' em vez de 'arredondamento' 

*/
div(a,b,inteira=arredondamento) =
{
z = a/b;
d = inteira(z);
r = a-b*d;
return([d,r]);
}

div2(a,b) =
{
if(norm(a)<norm(b),return([0,a]));
v = [1,-1,I,-I];
/* x virá a-b v_n, i.e um dos a-b, a+b, a-b*I, a+b*I */
for(n=1,#v,x=a-b*v[n];if(norm(x)<norm(a),return([v[n],0]+div2(x,b))));
return("Erro: todos os vizinhos têm norma grande!"); /* Este coisa não vai aparecer */ 
}

/* Um algoritmo recursivo euclideano
   para computar máximo divisor comum de (a) e (b),
   dado algoritmo (dividir) de divisão com resto [precisamos só resto]
   Valor default por dividir é div, mas tambem pode usar div2.
 */
mdc(a,b,dividir=div) =
{
if(b==0,return(a));
if(norm(b)>norm(a),return(mdc(b,a)));
return(mdc(b,dividir(a,b)[2]));
}

R1 = mdc(5,3+4*I)
R1 = mdc(3+4*I,5)
R3 = mdc(5,3+4*I,div2)
R4 = mdc(3+4*I,5,div2)
[R2,R3,R4]/R1