
#converts ideal of imaginary quadratic order to quadratic form
   IdealToQuadraticForm:=function(I)
     local O,a,b,c,w,t,n,ba,L,D,d;
     O:=Order(I);
     D:=Discriminant(O);
     d:=Denominator(I);
     if d <> 1 then 
       return IdealToQuadraticForm(d*I);
     fi;
     n:=Norm(I);
     w:=O.2;
     ba:=Basis(I);
     a:=List(ba[1])[1];
     b:=List(ba[2])[1];
     c:=List(ba[2])[2];
     L:=[Coerce(Z,a/c), Coerce(Z,(-2*b/c-Trace(w))),Norm(b+c*w)/(a*c)];
     Assert(L[2]^2-4*L[1]*L[3]=D, " quadratic form has discriminant of order");
     return rec(base:=L,ext1:=O);
   end;

#converts quadratic form to ideal of imaginary quadratic order
   QuadraticFormToIdeal:=function(L)
     local O,w,a,b,I,b2;
     O:=L.ext1;
     w:=O.2;
     L:=Base(L);
     a:=L[1];
     b:=L[2];
     b2:=Coerce(Integers(),(-b-Trace(w))/2) mod a;
     I:=a*O+(b2+w)*O;
     return I;
  end;

#returns ideal J such that I*J = [1]
   IdInvert:=function(I)
     local L,a,b,O,w;
     O:=Order(I);
     if IsPrincipal(I) then
       return 1*O;
     fi; 
     L:=Base(IdealToQuadraticForm(I));
     w:=O.2;
     if not GCD(L) = 1 then
       return Error(I," is not invertible");
     fi;
     a:=L[1];
     b:=L[2];
     return 1*O+((b-Trace(w)+2*w)/(2*a))*O;
   end;


#debugging function
   Assert:= function(b,s)
     if b then return;
     else return Error("Assumption: ",s,"failed!");
     fi;
   end;



#subfunction for IdealToQuadraticForm
   _step3:= function(a,b,c)
     local tmp,_step2;
     _step2:= function(a,b,c)
       local q,r,k;
       k:= 2*a;
       q:= Div(b,k);
       r:= b mod k;
       if r>a then
         r:= r-k;
         q:= q+1;
       fi;
       c:= c-Coerce(Integers(),1/2*(b+r)*q);
       b:=r;
       return _step3(a,b,c);
     end;
     if a>c then
       b:= -b;
       tmp:=a;
       a:=c;
       c:= tmp;
       return _step2(a,b,c);
     elif a=c and b<0 then
       b:=-b;
       return [a,b,c];
     else return [a,b,c];
     fi; 
   end;


#tests if a quadratic form is reduced
   RedTest:= function(L)
     local a,b,c;
     a:= L[1];
     b:= L[2];
     c:= L[3];
     if (Abs(b)<= a) and (a<=c) then
        if (Abs(b) = a) or (a=c) then
           if b>0 then 
              return true;
           else return false;
           fi;
        else return true;
        fi;
     else return false;
     fi;
   end;

#Reduction of a quadratic form
   QuadReduction:= function(R)
     local O,L,a,b,c,I2,ba,phi,s,_step2;
     _step2:= function(a,b,c)
        local q,r,k;
        k:= 2*a;
        q:= Div(b,k);
        r:= b mod k;
        if r>a then
           r:= r-k;
           q:= q+1;
        fi;
        c:= c-Coerce(Integers(),1/2*(b+r)*q);
        b:=r;
        return _step3(a,b,c);
     end;
     O:=R.ext1;
     L:=Base(R);
     a:= Coerce(Integers(), L[1]);
     b:= Coerce(Integers(),L[2]);
     c:= Coerce(Integers(),L[3]);
     if (-a<b and b<=a) then
         L:= _step3(a,b,c);
     else L:= _step2(a,b,c);
     fi;
     return rec(base:=L,ext1:=O);
   end;

#Reduction of an ideal
   IdReduction:=function(I)
     local L;
     L:=IdealToQuadraticForm(I);
     L:=QuadReduction(L);
     return QuadraticFormToIdeal(L);
   end;







