############################################
############################################
# Gesucht ist das Primideal p, das entsteht, wenn ein Primideal P aus der
# Maximalordnung mit o geschnitten wird:

intersection := function(ORDER, P)
	local Bo, BP, n, M, N, H, i, Bp;

	n := ORDER.n;
	Bo := TransformationMatrix(ORDER.o,ORDER.O);
	BP := BasisMatrix(P);

# 1.) Setze die Basismatrizen von o und P zu M uebereinander:
	M := VerticalJoin(Matrix(ORDER.Ri,Bo),Matrix(ORDER.Ri,BP));
# 2.) Bestimme eine Matrix mit maximalen Rang N, so dass N*M = 0:
	N := Matrix(ORDER.Ri,NullspaceMatrix(M));

# 3.) Streiche die letzten n Spalten von N weg:
	for i in [1..n] do
		N := RemoveColumn(N, n+1);
	od;

# 4.) Bp ist die Basismatrix von o geschnitten P:
	return Ideal(ORDER.o, N);
end;

##########################################################################

# Fuer Testbeispiele in kash3, fuer die allgemein der Schnitt zweier
# Ideale ueber der selben Ordnung benoetigt wird:

InstallMethod(
rec(name := "Intersection",
    kind := "FUNCTION",
    sin :=  [[elt-ids^int/ord^num,"A"],[elt-ids^int/ord^num,"B"]],
    sou :=  [[elt-ids^int/ord^num,"C"]],
    short:= "The intersection 'C' of the ideals 'A' and 'B'. Both ideals must be in the same order.",
    author := ["Carolin Just"],
    ex := ["o := EquationOrder(X^2+5);\nA :=256*o;\nB := 6*o;\nIntersection(A,B);"]
),
function(a,b)
	local A, B, n, o, M, N, i, AB;

	A := BasisMatrix(a);
	B := BasisMatrix(b);
	n := NumberOfRows(A);
	o := Order(a);

    if o <> Order(b) then
      Error("Intersection: Both ideals must be in the same order");
    fi;
# 1.) Setze die Basismatrizen von a und b zu M uebereinander:
	M := VerticalJoin(Matrix(Z,A),Matrix(Z,B));
# 2.) Bestimme eine Matrix mit maximalen Rang N, so dass N*M = 0:
	N := Matrix(Z,NullspaceMatrix(M));

# 3.) Streiche die letzten n Spalten von N weg:
	for i in [1..n] do
		N := RemoveColumn(N, n+1);
	od;

# 4.) AB ist die Basismatrix von a geschnitten b:
	AB := N*A;
	return Ideal(o, AB);
end
);

##########################################################################

# Gesucht ist das Ideal a, das entsteht, wenn ein Ideal A aus der
# mit einer Unterordnung o geschnitten wird:


InstallMethod(
rec(name := "Intersection",
    kind := "FUNCTION",
    sin :=  [[ord^num,"o"],[elt-ids^int/ord^num,"A"]],
    sou :=  [[elt-ids^int/ord^num,"a"]],
    short:= "The intersection 'a' of the ideal 'A' and the order 'o'. The order 'o' must be a suborder of the order of 'A'",
    author := ["Carolin Just"],
    ex := ["o := EquationOrder(X^2-5);\nO := MaximalOrder(o);\nA :=256*O;\na := Intersection(o,A);"]
),
function(o, A)
	local Ri, Bo, BA, n, M, N, H, i;

        if not FieldOfFractions(o)=FieldOfFractions(Order(A)) then
          Error("Intersection: the order of A and the order o are not compatible");
        fi;

        n := Degree(o);
	Bo := TransformationMatrix(o,Order(A));

        if not IsUnit(Bo.ext1) then
          Error("Intersection: the order o must be a suborder of the order of A");
        fi;

        BA := BasisMatrix(A);
	Ri := CoefficientRing(o);


# 1.) Setze die Basismatrizen von o und A zu M uebereinander:
	M := VerticalJoin(Matrix(Ri,Bo),Matrix(Ri,BA));
# 2.) Bestimme eine Matrix mit maximalen Rang N, so dass N*M = 0:
	N := Matrix(Ri,NullspaceMatrix(M));

# 3.) Streiche die letzten n Spalten von N weg:
	for i in [1..n] do
		N := RemoveColumn(N, n+1);
	od;

# 4.) Ba ist die Basismatrix von o geschnitten A:
	return Ideal(o, N);
end
);


##########################################################################
##########################################################################
# Gesucht sind die Primideale p, die a enthalten.

prime_ideals_of_a := function(ORDER, a)
	local amax, prid_amax, prid_a, P, p; 

# Berechne aO in der MaxOrd O:
	amax := Ideal(ORDER.O, Generators(a));

# Zerlege aO in Primideale:
	prid_amax := Factorization(amax);
	prid_a := [];

# Speichere Primideale nicht doppelt ab: 
	for P in prid_amax do
		p := intersection(ORDER, P[1]);
		if not (p in prid_a) then
			Add_(prid_a,p); 
		fi;
	od;

# Gib die Liste der Primideale p zurueck:
	return prid_a;
end;

###############################################################################

# Bestimmung des Exponenten m mit p^mOp < aOp
# 	Test: Ist ein Ideal einer Ordnung in einer Primidealpotenz der selben 
# 	Ordnung enthalten?

mSearch := function(o,a,p)
	local ch_p, my, apmy, m, s, ende, pm, gm, pm_in_apmy, alpha;

	ch_p := IsPrimePower(Index(o,p)).ext1;
	my := Valuation(Index(o,a),ch_p).base;
	apmy := a+p^my;
	s := 0;
	ende := Ilog(2,my)+1;

	while s<=ende do
		m := Minimum(2^s,my);
		pm := p^m;
		gm := Generators(pm);
		pm_in_apmy := true;
		for alpha in gm do
			if not(alpha in apmy) then
				pm_in_apmy := false;
			fi;
		od;
		if pm_in_apmy then
			s := ende+1;
		else
			s := s+1;
		fi;
	od;
	return m;
end;

################################################################################
################################################################################


# Diese Funktion erzeugt ein record, das alle wichtigen Infos, die nicht doppelt
# berechnet werden sollen, enthaelt

infos_zu_p := function(ORDER,a,p,m)

	local 	o, n, Ri, Frac, omega_vec, p_k, P_k, P_k_inv, P_invQ, tau_multi, Quot_k, 
		Relationen_k_l, ha, ende, s, k, pk, Pk, S, Q_lk, Q_lk_inv, Bas_vec, Bas_m, 
		i, Bas_apm, iota;

	o := ORDER.o;
	n := ORDER.n;
	Ri := ORDER.Ri;
	Frac := ORDER.Frac;

# Initialisierung der record-Elemente
	omega_vec := Matrix(o,n,1,Basis(o));
	p_k := [];	# Liste der Ideale p^k
	P_k := []; 	# Liste der TrafoMatrizen Pk von p^k
	P_k_inv := [];	# Liste der Inversen Pk^-1
	P_invQ := [];	# Liste der Matrizen Pk^-1*Q_lk. Diese werden fuer den
			# diskreten Logarithmus in (1+p^k)/(1+p^l) gebraucht.
	tau_multi := [];# Liste der Erzeuger (mult. Gruppe)
	Quot_k := [];	# Liste der Restklassenringe o/p^k
	Relationen_k_l := [];	# Liste der Relationenmatrizen von (1+p^k)/(1+p^l)
	ha := [];	# Liste der Erzeuger von (1+a+p^m)/(1+p^m)
	ende := Ilog(2,m);

# Fuellen der record-Elemente
	for s in [0..ende+1] do
		k := Minimum(2^s,m);
# Das Ideal b^k, seine Trafomatrix und deren Inverses:
		pk := p^k;	# p^k = p^(2^s)
   # Achtung: TrafoMatrizen werden bei Kash3 immer schon in HNF angegeben!  
		Pk := TransformationMatrix(pk);
		p_k[s+1] := pk;
		P_k[s+1] := Pk.base;
		P_k_inv[s+1] := Matrix(Frac,Pk)^-1;

# Relationen von (1+b^k)/(1+b^l). Ausnahmsweise eingetragen an der Stelle fuer l (!!)
# Umgeaendert zu einer SNF
		if (s<>0) then 
			S := SmithForm(Matrix(ORDER.Ri, Pk * P_k_inv[s]));
				Relationen_k_l[s+1] := S.base;
			Q_lk := S.ext2;		# Zur Anpassung der Erzeugervektoren tau
			Q_lk_inv := Matrix(Ri,Matrix(Frac,Q_lk)^(-1));
				P_invQ[s] := P_k_inv[s]*Q_lk;	# fuer diskreten Log. 
								# in (1+p^k)/(1+p^l)

# Neue Konstruktion der Basen:
# Die additiven Basisvektoren - an die SNF von Pl*Pk^(-1) angepasst
			Bas_vec := Matrix(o,Q_lk_inv*P_k[s])*omega_vec;

# Die multiplikativen Erzeuger von (1+p^k)/(1+p^l). Als Vektor abgespeichert:
#			Bas_m := Bas_vec;
			for i in [1..n] do
#				SetEntry(Bas_m,i,1, Bas_m[i][1]+1);
				SetEntry(Bas_vec,i,1, Bas_vec[i][1]+1);
			od;
#			tau_multi[s] := Bas_m;
			tau_multi[s] := Bas_vec;
		fi;

# Die Restklassenringe o/p^k
		Quot_k[s+1] := Quotient(o, pk);

# Erzeuger von (1+(a+p^m)^k)/(1+(a+p^m)^l)
# Brauche ich fuer die Berechnung der Matrix A in unitGroup
		if k<>m then
			Bas_apm := Basis((a+p^m)^k);	
				# Liste der Basiselemente von (a+p^m)^k
			for iota in [1..n] do
				ha[s*n+iota] := Bas_apm[iota]+1;     
					# = Basiselemente +1
			od;
		fi;	
	od;

	return 	    rec(p_k:=p_k, 
			P_k:=P_k,
		   	P_k_inv:=P_k_inv,
			P_invQ := P_invQ,
		   	Relationen_k_l:=Relationen_k_l, 
		   	tau_multi:=tau_multi, 
			Quot_k:= Quot_k, 
		   	ha:=ha);
end;

#############################################################################################
#############################################################################################


##################################################################################
# Der diskrete Logarithmus fuer (1+b)/(1+b^m)

# Es gibt 1 Hilfsfunktion:
#     Die Fkt. fuer den diskreten Log. in (1+b^k)/(1+b^l) gibt einen Vektor c
#     mit Koeffizienten fuer die Darstellung eines Elementes in (1+b^k)/(1+b^l)
#     aus.
#
# Die Hauptfunktion gibt eine Sequenz a aus, die alle Koeffizienten enhaelt,
# so dass [eta] aus (1+b)/(1+b^m) darstellbar ist als [a*g], g sind die Erzeuger
# von (1+b)/(1+b^m)

###################################################################################
# Hilfsfunktion:

# s zeigt auf k!!
# Berechnet wird der diskrete Logarithmus in (1+b^k)/(1+b^l):

discrete_logarithm_bk_bl := function(ORDER,IDEAL,s,eta)
	local beta;

# 1.) Bestimme die Koeffizienten von eta-1 bzgl. omega. Sie werden fuer (2.) gebraucht:
	beta := ElementToSequence(Coerce(ORDER.o,eta-1));

# 2.) Bestimme die Koffizienten von eta-1 bzgl. tau_k und gib sie aus:
	return Matrix(Z,Matrix(ORDER.Frac, 1, ORDER.n, beta) * IDEAL.P_invQ[s]);
		# Achtung, hier steht Z - zu beachten fuer R = F_q[t]!!! 
end;

#######################################################################################

diskreter_logarithmus:= function(ORDER,IDEAL,eta,m)

	local k,coeff,s,Quot_m,l,c,g;

	k := 1;
	coeff := [];
	s := 1;	# s zeigt auf die Stelle fuer b^k,
		# s+1 zeigt auf die Stelle fuer b^l
	Quot_m := IDEAL.Quot_k[Ilog(2,m)+2];	# = o/p^m
	eta := Coerce(Quot_m,eta);

	if k=m then	# Fuer den Sonderfall m=1
		coeff := ElementToSequence(discrete_logarithm_bk_bl(ORDER,IDEAL,s,eta));
	fi;

	while k<>m do
		l := Minimum(2*k,m);
		c := discrete_logarithm_bk_bl(ORDER,IDEAL,s,eta);
		g := Transpose(Matrix(Quot_m, IDEAL.tau_multi[s]));	
			# die Eintraege der Erzeuger nach o/p^m schieben, damit PowerProduct 
			# schneller laeuft
		eta := eta*(PowerProduct(g,c)^(-1));
		coeff := Concatenation(coeff,ElementToSequence(c));
		k := l;
		s := s+1;
	od;

	return coeff;
end;


################################################################################
################################################################################


##############################################################################
# Mg gehoeren zu (1+b)/(1+b^k).
# Nh gehoeren zu (1+b^k)/(1+b^l).
# s zeigt auf k!!

exact_sequence := function(ORDER, IDEAL, s, Mg, Nh)

	local M, N, g, g_hoch_M_Seq, g_hoch_M, k, r, P_seq, i, P, c, Rel;

	M := Mg[1];
	r := NumberOfRows(M);
	N := Nh[1];

# 1.) Den Vektor M*g berechnen: - Bis jetzt funktioniert PowerProduct nicht, wenn M eine Matrix ist,
#				  deswegen erstmal noch so umstaendlich...
	g := Matrix(IDEAL.Quot_k[s+1],Transpose(Mg[2]));
	g_hoch_M_Seq := [];
	for k in [1..r] do
		g_hoch_M_Seq[k] := PowerProduct(g,M[k]);
	od;
	g_hoch_M := Matrix(ORDER.o,r,1,g_hoch_M_Seq);

# 2.) Fuer jeden Eintrag von g_hoch_M den diskreten Logarithmus in (1+b^k)/(1+b^l)
#     loesen und in einer Liste abspeichern:
	P_seq := [];
	for i in [1..r] do
		P_seq[i] := 
		discrete_logarithm_bk_bl(ORDER,IDEAL,s,g_hoch_M[i][1]);
	od;
# 3.) Aus der Liste die Matrix P bauen:
	P := VerticalJoin(P_seq);

# 4.) M und N diagonal anordnen, danach -P einfuegen:
	c := NumberOfColumns(M);
	Rel := InsertBlock(DiagonalJoin(M,N),Matrix(Z,-P),1,c+1);

# 5.) Ausgegeben werden die neue Relationenmatrix M und der neue Erzeugervektor g:
	return [HermiteForm(Rel).base, VerticalJoin(Mg[2],Nh[2])];
end;

##############################################################################
# s zeigt auf l!!!

computing_1_plus_p_nach_1_plus_p_m := function(ORDER, IDEAL,m)
	local ende, Mg, k, s, l, Nh;

		ende := Ilog(2,m);

# 1.) Berechne Erzeuger und Relationen zu (1+p)/(1+p^2):
		Mg := [IDEAL.Relationen_k_l[2], IDEAL.tau_multi[1]];

# 2.) Von (1+p)/(1+p^k) und (1+p^k)/(1+p^l) nach (1+p)/(1+p^l). 
		k := 2;
		s := 2;
		while k<m do
			l := Minimum(2*k,m);
			Nh := [IDEAL.Relationen_k_l[s+1], IDEAL.tau_multi[s]];
			Mg := exact_sequence(ORDER, IDEAL, s, Mg, Nh);
			k := l;
			s := s+1;
		od;

	return Mg;
end;

#############################################################################
#############################################################################
# Ein Testbeispiel:

# Hier ein Beispiel fuer eine Nicht-Maximalordnung und Nicht-Gleichungsordnung:
# (laeuft auch durch!!)

#oe := EquationOrder(X^2-5);
#o := Order(oe,Matrix(Z,2,2,[2,0,0,6]),2);
#n := 2;
#Frac := Q;
#Ri := Z;

#o := MaximalOrder(X^20 -5*X^19 +3*X^18 -7*X^17 +12*X^16 -3*X^15 +7*X^14 +17*X^13-6*X^12 
#-9*X^11 +5*X^10 +11*X^9 +3*X^8 -10*X^7 -18*X^6 
#+10*X^5 +17*X^4 +16*X^3 -14*X^2 +3*X -15);
#Frac := Q;
#Ri := Z;
#n := 20;
#b := Ideal(o,7,1+o.2); 
#Gens := Generators(b);
#eta := 789*Gens[1] + 456233*Gens[2] +1;
#m := 15;

#o := MaximalOrder(X^14-15*X^13+14*X^12-5*X^11-17*X^10 +17*X^9 +4*X^8 +8*X^7 +16*X^6 
#-17*X^5 +8*X^4 +12*X^3 -4*X^2 -16*X -10);
#Frac := Q;
#Ri := Z;
#n := 14;
#b := Ideal(o,5,o.2);  # laut kash3 ebenfalls ganz o??

#b := Ideal(o,o.2); 
#eta := 45629074*o.2 +1;
#m := 17;

#o := EquationOrder(X^6+11*X^5+20*X^4+3*X^3-18*X^2+7*X+11);
#O := MaximalOrder(X^6+11*X^5+20*X^4+3*X^3-18*X^2+7*X+11);
#Frac := Q;
#Ri := Z;
#n := 6;
#b := Ideal(o,47,[40,20,1,0,0,0]); 
#Gens := Generators(b);
#eta := 789*Gens[1] + 456233*Gens[2] +1;
#m := 40;

# b := Ideal(o,3,[1,0,0,1,0,0]); Das ist bereits die ganze Ordnung!!

#o := MaximalOrder(X^3+23*X^2+88*X+51);
#Frac := Q;
#Ri := Z;
#n := 3;
#p := 7*o;
#b := Ideal(o,29,[18,1,0]);
#Gens := Generators(b);
#eta := 789*Gens[1] + 456233*Gens[2] +1;
#m := 31;

#o := MaximalOrder(X^4+14*X^3+8*X^2+4*X+11);
#Frac := Q;
#Ri := Z;
#n := 4;
#b := Ideal(o,2,[1,1,0,0]);
#Gens := Generators(b);
#eta := Element(o,[4,1,0,0]);
#eta := 789*Gens[1] + 456233*Gens[2] +1;

#o := EquationOrder(X^2-5);
#m := 5;
#b := 5*o;
#b3 := b^3;
#b5 := b^5;

#############################################################################
################################################################################
# Das ist die Hauptfunktion, die alle bisher bereit gestellten Techniken 
# zusammenfuehrt. Ihr Ablauf entspricht dem Ablauf des Pseudocodes auf S.9

#############################################################################
# Hilfsfunktion: Chinesischer Restsatz fuer unseren speziellen Fall
# Input:  Ordnung o
#	  Liste der Ideale a+p^m,
#	  den i-ten Erzeugervektor g,
#	  den Zeiger i
#	  die Anzahl s der Primideale, die a enthalten
# Output: angepasster Erzeugervektor h

chin_rem := function(o, L_apm, g, i, s) 
	local b, j, apm, len_g, h, k;

# 1.) Multipliziere die (a+p_j^m_j), j <> i
	b := 1;
	for j in [1..i-1] do
		b := b*L_apm[j];
	od;
	for j in [i+1..s] do
		b := b*L_apm[j];
	od;

# 2.) Loese fuer jeden Eintrag von g den Chinesischen Restsatz
	apm := L_apm[i];
	len_g := Length(g);
	h := [];
	for k in [1..len_g] do
		Add_(h, ChineseRemainderTheorem(apm,b,Coerce(o,g[k][1]), Coerce(o,1)));
	od;

# 3.) Gib h als Vektor zurueck
	return Matrix(len_g,1,h);
end;  

###############################################################################
# Hauptfunktion:
# Input:	Restklassenring o/a
# Output:	Einheitengruppe (o/a)* einschliesslich der Abbildungen zwischen
#		(o/a)* und (o/a) 

unitGroup := function(oa)

	local 	a, o, ORDER, omega, Ri, PrIds, s, m_List, i, P_IDEAL, Rel, Erz, L_apm,
		o_by_p, p, m, K, q, UK, Kphi, UKphi, Zeta, vp, t, Xi_mod_pm, Xi,
		RES_FIELD, Mg, M, g, len_ha, A_entry, j, A, MA, H, Erz_fuer_psi,
		Rel_fuer_psi, S, Rel_fuer_phi, Q_smith, Erz_fuer_phi, zeiger, nrow_rel,
		G, App,	U, phi, psi, disc_log;

# 0.) Grundsaetzlich brauchen wir:
	a := Modulus(oa);
	o := Order(a);

# 1.) Falls a=o, gib Fehlermeldung aus:
	if IsUnit(Determinant(BasisMatrix(a))) then
		Error("ideal = order. This quotient ring has no unit group.");
	fi;

# 2.) Konstruiere ein record mit den wichtigen Angaben zur Ordnung o:
	ORDER := rec();
		omega := Basis(o);
		Ri := CoefficientRing(o);
			ORDER.o := o;
			ORDER.omega := omega;
			ORDER.n := Length(omega);
			ORDER.Ri := Ri;
			ORDER.Frac := FieldOfFractions(Ri);
			ORDER.O := MaximalOrder(o);

# 3.a) Bestimme alle Primideale p, die a enthalten 
	PrIds := prime_ideals_of_a(ORDER,a);
# 3.b) Bestimme die Anzahl s der Primideale
	s := Length(PrIds);
# 3.c) Bestimme alle Exponenten m
	m_List := [];
	for i in [1..s] do
		m_List[i] := mSearch(o,a,PrIds[i]);
	od;
# 3.d) Sammle zu jedem Primideal alle noetigen Informationen
	P_IDEAL:= [];
	for i in [1..s] do
		P_IDEAL[i] := infos_zu_p(ORDER,a,PrIds[i],m_List[i]);
	od;

# 4) Berechne fuer jedes Primideal p die Einheitengruppe von (O_p/aO_p):
	Rel := [];	# Liste aller Relationen der O_p_i/aO_p_i
	Erz := [];	# Liste aller Erzeuger der O_p_i/aO_p_i
	L_apm := [];	# Liste aller a+p_i^(m_i)
	o_by_p := [];	# Liste der Abbildungen fuer die (o/p_i)*

	for i in [1..s] do
		p := PrIds[i];	# Primideal p
		m := m_List[i];	# Exponent m zu p
		L_apm[i] := a+p^m;

# 4.1.) Berechne die multiplikative Gruppe von (O/p):
	# Koerper K := (o/p) und dessen Groesse q. 
	# Dabei gibt (q-1) Relation des Erzeugers von (o/p)* an 
		K := ResidueClassField(p);
		q := Size(Base(K));

	# Einheitengruppe UK := (o/p)*
		UK := UnitGroup(K);

	# Fuer den diskreten Logarithmus in (o/p)*
		Kphi := K.ext1;
		UKphi := UK.ext1;
		
	# Repraesentant von (o/p)* in o
		Zeta := Preimage(UKphi(UK.1),Kphi); 

# 4.1.a) Falls m=1, sind wir mit der Berechnung von (o/p)* fertig
		if m=1 then
			Erz[i] := Matrix([[Zeta]]);	
				# = Erzeugervektor
			Rel[i] := Matrix([[q-1]]);
				# = Relationenmatrix
			RES_FIELD := rec();
				RES_FIELD.Kphi := Kphi;
				RES_FIELD.UKphi := UKphi;
			o_by_p[i] := RES_FIELD;
				# = Record, das die Abbildungen zwischen (o/p)* und o enthaelt.

# 4.1.b) Fuer alle anderen m berechnen wir (1+p)/(1+a+p^m)
		else
	# Dafuer Hensellifting von Zeta:
			vp := Valuation(Coerce(o,q),p);	
			t := 0;
			while (Minimum(q^t,1+t*vp)<m) do
				t := t+1;
			od;
			Xi_mod_pm := Coerce(Quotient(o,p^m),Zeta)^(q^t);  

	# Der Repraesentant von Xi_mod_pm als Matrix
			Xi := Matrix([[Preimage(Xi_mod_pm,Quotient(o,p^m).ext1)]]);

	# Die Abbildungen fuer (o/p)* in einem record sammeln:
			RES_FIELD := rec();
				RES_FIELD.Kphi := Kphi;
				RES_FIELD.UKphi := UKphi;
			o_by_p[i] := RES_FIELD;

# 4.2.) Berechne (1+p)/(1+p^m):
			Mg := computing_1_plus_p_nach_1_plus_p_m(ORDER,P_IDEAL[i],m);
			M := Mg[1];
			g := Mg[2];

# 4.3.) Berechne (1+a+p^m)/(1+p^m):
# a) Erzeuger fuer (1+(a+p^m))/(1+(a+p^m)^2), ..., (1+(a+p^m)^(2^ende))/(1+(a+p^m)^m) 
#     sind unter IDEAL.ha bereits abgespeichert.
# b) Berechne fuer jeden einzelnen Erzeuger den diskreten Logarithmus bzgl. der
#     Erzeuger g von (1+p)/(1+p^m)
# c) Setze die erhaltenen Sequenzen zur Matrix A zusammen:

			len_ha := Length(P_IDEAL[i].ha);
			if (len_ha>0) then
				A_entry := [];
				for j in [1..len_ha] do
					A_entry[j] := Matrix(Ri, 1, len_ha, 
					diskreter_logarithmus(ORDER,P_IDEAL[i],P_IDEAL[i].ha[j],m));
				od;
				A := VerticalJoin(A_entry);

# 4.4.) Setze alle Relationenmatrizen zusammen:
				# Setze M und A uebereinander =: MA, gib MA in HNF zurueck:
				MA := HermiteForm(VerticalJoin(M,A)).base;

				# Die untere Haelfte von MA ist Nullmatrix, sie wird abgeschnitten:
				MA := ExtractBlock(MA, 1, 1, len_ha, len_ha);

			else
				MA := M;
			fi;

			# Setze (q-1) und MA zu Diagonalmatrix zusammen:
			Rel[i] := DiagonalJoin(Matrix([[q-1]]), MA);

# 4.5.) Setze die Erzeuger Xi von (o/p)* und g zusammen: 
			Erz[i] := VerticalJoin(Xi,g);
		fi;

	od;

## Hier muessen vier Listen herauskommen:
#  a) Erz: Liste aller Erzeuger, also zu jedem p_i
#  b) Rel: Liste aller Relationen
#  c) L_apm : Liste der a+p^m
#  d) o_by_p: Liste mit Abbildungen zu (o/p_i)*


# 5) Jetzt noch den Chinesischen Restsatz anwenden:
# (Dafuer habe ich oben die Liste L_apm der Ideale a+p^m angelegt)

	if s=1 then
		H := Erz;
	else
		H := [];
		for i in [1..s] do
			H[i] := chin_rem(o,L_apm,Erz[i],i,s);
		od;
	fi;

# 6.) Zusammensetzen von Erzeugern und Relationen
#	Fuer psi brauche ich genau den zusammengesetzten Vektor aus H
#	Fuer phi brauche ich die Basis, die entsteht, wenn die SNF angewandt
#	wird. - Nicht die Erzeuger

# 6.1) Daten zu psi
	Erz_fuer_psi := VerticalJoin(H);
	Rel_fuer_psi := DiagonalJoin(Rel);

# 6.2) Daten zu phi
	S := SmithForm(Rel_fuer_psi);
	Rel_fuer_phi := S.base;
	Q_smith := Matrix(Z,Matrix(Q,S.ext2)^(-1));
	Erz_fuer_phi := [];
	for i in [1..NumberOfRows(Erz_fuer_psi)] do
		Erz_fuer_phi[i] := 
			PowerProduct(Transpose(Matrix(oa,Erz_fuer_psi)), Q_smith[i]);
	od;
	zeiger := 1;
	while (Rel_fuer_phi[zeiger][zeiger] = 1) do
		Remove_(Erz_fuer_phi,1);
		zeiger := zeiger+1;
	od;


# 7.) Konstruktion der Einheitengruppe U := (o/a)*

	nrow_rel := NumberOfRows(Rel_fuer_phi);
	G := FreeAbelianGroup(nrow_rel);
	App := Apply([1..nrow_rel], x-> Coerce(G,List(Rel_fuer_phi[x])));
	U := Quotient(G,App);


# 8.) Diskreter Logarithmus:

# Input:	Ein Element u aus U := (o/a)*
# Output:	Eine "Map" in kash. Das ist: 
#		eine Abbildung phi mit: 
#		Urbildbereich (o/a)*,
#		Bildbereich (o/a),
#		Umkehrabbildung psi, die in kash durch "Preimage" aufgerufen wird 

# 8.1.) phi: Von (o/a)* nach o/a
	# Input:	Element u aus (o/a)*
	# Output:	Repraesentant in o/a (?)

	# Ich brauche hierfuer den per SNF veraenderten Erzeugervektor Erz_fuer_phi,
	# am besten bereits nach (o/a) geschoben und als Sequenz. 

	phi := function(u)
		local c;

		c := ElementToSequence(u);
		return PowerProduct(Erz_fuer_phi,c);		 
	end;

# 8.2.) psi: Von o/a nach (o/a)* 
	# Input:	Repraesentant alpha aus o
	# Output:	Koeffizientenvektor von alpha bzgl. der Erzeuger von (o/a)*

	# Ich brauche hierfuer den unveraenderten Erzeugervektor Erz_fuer_psi

	psi := function(alpha)
		local dl, i, a0, a1, dl_smith;

	# I.) Wenn [alpha] kein Element der Einheitengruppe ist, wirf Fehlermeldung
		if (not IsUnit(alpha)) then
			Error("Element is no unit.");
		fi;

	# II.) Setze den diskreten Logarithmus aus den diskreten Logarithmen der opi/aopi
	#     zusammen:
		dl := [];
		for i in [1..s] do
			# diskreter Logarithmus von alpha in (o/p)*
				a0 := ElementToSequence(Preimage(o_by_p[i].Kphi(Coerce(o,alpha)),o_by_p[i].UKphi));
				alpha := alpha*(Coerce(oa,H[i][1][1])^(-a0[1]));	# so ungefaehr...
				dl := Concatenation(dl,a0);
			# diskreter Logarithmus von alpha in (1+p)/(1+a+p^m)				
				g := Remove(ElementToSequence(Matrix(oa,H[i])),1);
					# (1.Eintrag von H[i] wegstreichen liefert Erzeugervektor von (1+p_i)/(1+p_i^m_i))
					# (H[i] nach oa schieben verkuerzt Rechenaufwand beim Potenzieren)
				if Length(g)>0 then
					a1 := diskreter_logarithmus(ORDER,P_IDEAL[i],alpha,m_List[i]);
					alpha := alpha*PowerProduct(g,a1);
					dl := Concatenation(dl,a1);
				fi;
		od;

	# III.) Passe dl an die Basis (Erz_fuer_psi) an:
		dl_smith := ElementToSequence(Matrix(Z,1,Length(dl),dl)*S.ext2);
		zeiger := 1;
		while (Rel_fuer_phi[zeiger][zeiger] = 1) do
			Remove_(dl_smith,1);
			zeiger := zeiger+1;
		od;

	# IV.) Wandle dl_smith in kanonische Repraesentanten um und gib zurueck:
#		return ElementToSequence(Coerce(U,dl_smith));
		return Coerce(U,dl_smith);
	end;

# 8.3.) disc_log gibt phi als Umkehrabbildung von psi aus.
	disc_log := Map(U,oa,phi,psi);

# 9.) Rueckgabe eines records, das als Basis die Einheitengruppe U = (o/a)* und als 
#     Erweiterung die Map disc_log enthaelt
	return rec(base := U.base, ext1 := disc_log); 
end;


#InstallMethod(
#rec(name:="UnitGroup",
#    kind:="FUNCTION",
#    sin :=[[res^num,"oa"]],
#    sou :=[[grp^abl,"U"],[map(grp^abl,res^num),"phi"]],
#    short := "carolin weiss wie es geht",
#    ex := [],
#    see := [],
#    author :=["Carolin Just"]
#    ),
#unitGroup);

_UnitGroup := UnitGroup;
UnitGroup := function(arg)
  if Is(Type(arg[1]),res^num) then
    if not IsMaximal(Order(Modulus(arg[1]))) then
      return unitGroup(arg[1]);
    fi;
  fi;
  return _CallFunction(_UnitGroup,arg);
end;


