#
#--> Gaussian(R)
# Creates the domain of Gaussians over the ring R e.g. Z, Q
# Author MBM: 89
#
macro(real=1,imag=2);
Gaussian := proc() local F,G; option remember;

	F := args[1];
	if hasCategory(F,Field) then G := Field()
	elif hasCategory(F,Ring) then G := Ring()
	else ERROR(`1st argument must be a Ring`)
	fi;

	defOperation( Conj,  G &-> G, G );
	defOperation( Make,  [F,F] &-> G, G );
	defOperations( {Re,Im}, G &-> F, G );

	# Rep is a Record(real:F,imag:F) implemented as Maple List
	G[DomainName] := Gaussian;
	G[Characteristic] := F[Characteristic];
	G['_i'] := [F[0],F[1]];
	G[Make] := proc(x,y) [x,y] end;
	G[Re] := proc(x) x[real] end;
	G[Im] := proc(x) x[imag] end;

	G[Output] := subs('R' = F, proc(x)
		R[Output](x[real]) + R[Output](x[imag]) * _i
	   end);
	G[Random] := subs('R' = F, proc() [R[Random](),R[Random]()] end);
	G[Input] := subs('R' = F,proc(x) local i,r,y;
		y := subs(_i=(-1)^(1/2),(-1)^(1/2)=_i,expand(x));
		r := traperror(coeff(y,_i,0));
		if r = lasterror then RETURN(FAIL) fi;
		r := R[Input](r);
		if r = FAIL then RETURN(FAIL) fi;
		i := traperror(coeff(y,_i,1));
		if i = lasterror then RETURN(FAIL) fi;
		i := R[Input](i);
		if i = FAIL then RETURN(FAIL) fi;
		[r,i]
	    end);


	if hasProperty(F,CannonicallyUniquelyRepresented) then
		addProperty(G,CannonicallyUniquelyRepresented);
		G[`=`] := <evalb(x=y)>
	else	G[`=`] := subs('R' = F, proc(x,y)
			R[`=`](x[real],y[real]) and R[`=`](x[imag],y[imag])
		    end)
	fi;

	G[0] := [F[0],F[0]];
	G[1] := [F[1],F[0]];

	# Inline arithmetic operations for these special cases
	if isDomain(F,Integer) then
		RETURN( GaussianInteger(op(G)) )
	elif isDomain(F,Rational) or isDomain(F,Floats) then
		RETURN( MapleGaussianArithmetic(op(G)) )
	fi;

	G[`+`] := subs('R' = F, proc(x,y)
		[R[`+`](x[real],y[real]),R[`+`](x[imag],y[imag])]
	    end);
	G[`-`] := subs('R' = F, proc(x,y)
		if nargs = 1 then [R[`-`](x[real]),R[`-`](x[imag])]
		else [R[`-`](x[real],y[real]),R[`-`](x[imag],y[imag])] fi
	    end);
	G[`*`] := subs('R' = F, proc(x,y)
		if type(x,integer) then
		    [R[`*`](x,y[real]),R[`*`](x,y[imag])]
		else
		    [R[`-`](R[`*`](x[real],y[real]),R[`*`](x[imag],y[imag])),
		     R[`+`](R[`*`](x[real],y[imag]),R[`*`](x[imag],y[real]))]
		fi
	    end);
	if hasCategory(G,Field) then
		G[EuclideanNorm] := subs('R' = F, proc(x)
		    R[`+`](R[`^`](x[real],2),R[`^`](x[imag],2))
		end);
		G[Inv] := subs('R' = F, 'D' = G, proc(x) local d;
		    if type(x,integer) then [R[Inv](x),0] else
			d := D[EuclideanNorm](x);
			[R[`/`](x[real],d),R[`-`](R[`/`](x[imag],d))]
		    fi
	 	end);
	fi;
	G[Conj] := subs('R' = F, proc(x) [x[real],R[`-`](x[imag])] end);
	op(G)
end:

MapleGaussianArithmetic := proc() local G;
	# Inline Code +, -, *, Inv, EuclideanNorm, conj
	G := args[1];
	G[`+`] := proc(x,y) [x[real]+y[real],x[imag]+y[imag]] end;
	G[`-`] := proc(x,y)
		if nargs = 1 then [-x[real],-x[imag]]
		else [x[real]-y[real],x[imag]-y[imag]] fi
	    end;
	G[`*`] := proc(x,y)
		if type(x,integer) then [x*y[real],x*y[imag]] else
			[x[real]*y[real]-x[imag]*y[imag],
		 	x[real]*y[imag]+x[imag]*y[real]]
		fi
	    end;
	if hasCategory(G,Field) then
		G[EuclideanNorm] := proc(x) [x[real]^2,x[imag]^2] end;
		G[Inv] := proc(x) local d;
		    if type(x,integer) then [1/x,0] else
			d := x[real]^2+x[imag]^2;
			[x[real]/d,-x[imag]/d]
		    fi
	 	end
	fi;
	G[Output] := proc(x) x[real] + x[imag] * _i end;
	G[Conj] := proc(x) [x[real],-x[imag]] end;
	op(G)
end:

# Special Implementation for Gaussian Integers
GaussianInteger := proc() local G; option remember;

	G := EuclideanDomain(args);
	G := MapleGaussianArithmetic(G);

	# add Euclidean Domain operations
	G[Unit] := proc(x)
		if x[imag] > 0 and x[real] <= 0 then [0,1]
		elif x[real] < 0 and x[imag] <= 0 then [-1,0]
		elif x[imag] < 0 and x[real] >= 0 then [0,-1]
		else [1,0]
		fi
	    end;
	G[Inv] := proc(x)
		if x = [1,0] then [1,0]
		elif x = [-1,0] then [-1,0]
		elif x = [0,1] then [0,-1]
		elif x = [0,-1] then [0,1]
		elif x = [0,0] then ERROR(`division by zero`)
		else FAIL
		fi
	    end;
	G[EuclideanNorm] := proc(x) x[real]^2+x[imag]^2 end;
	G[Rem] := proc(a,b,q) local nb,cr,ci,qr,qi,r;
		nb := b[real]^2+b[imag]^2;
		cr := a[real]*b[real]+a[imag]*b[imag];
		ci := a[imag]*b[real]-a[real]*b[imag];
		qr := iquo(cr,nb,'r');
		if cr >= 0 then if r > iquo(cr,2) then qr := qr + 1 fi fi;
		if cr < 0 then if r < iquo(cr,2) then qr := qr - 1 fi fi;
		qi := iquo(ci,nb,'r');
		if ci >= 0 then if r > iquo(ci,2) then qi := qi + 1 fi fi;
		if ci < 0 then if r < iquo(ci,2) then qi := qi - 1 fi fi;
		if nargs = 3 then q := [qr,qi] fi;
		[a[real]-b[real]*qr+b[imag]*qi,a[imag]-b[imag]*qr-b[real]*qi]
	    end;
	G[Prime] := subs('F' = G, proc(x) local y;
		y := F[Normal](x);
		if y[imag] = 0 then irem(y[real],4) = 3 and isprime(y[real])
		else isprime(F[EuclideanNorm](y))
		fi
	    end);
	op(G)
end:

save `G.m`;
quit
