//Copyright (c) 2011, Xavier Caruso and David Lubicz
//All rights reserved.
//
//Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
//
//   * Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
//   * Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
//
//THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//


/*
This package provides tools (including basic operations on free modules) to work over
S_nu = { sum a_i u^i / v_R(a_i) + nu.i>= 0, a_i in Frac R } where R is DVR.
When nu is rational, we set nu = a/b. For convenience, in the following S_nu is denoted
by S.

In the sequel, we will work with two related rings: Snu[1/pi] (which is referred to as Sp)
and the completion of Snu[pi^b/u^a] (referred to as Su).
*/

/*******************************************************************
 *                                                                 *
 *   MISCELLANEOUS FUNCTIONS                                       *
 *   these functions are not meant to be used outside of           *
 *   the library                                                   *
 *   Authors: X. Caruso, D. Lubicz                                 *
 *                                                                 *
 *******************************************************************/

//
//Miscellaneous functions for the arithmetic over S//

//*************************************
// remove the zero coefficients when printing
// an element of S
//*************************************
RemoveZeroesOverS := function (X)
  local l, pow,S,u;

  S := Parent(X);
  u := UniformizingElement(S);
  l,pow := Coefficients (X);
  for i := #l to 1 by -1 do
    if RelativePrecision(l[i]) gt 0 then break; else Prune (~l); end if;
  end for;

  return Parent(X)!l*u^pow;
end function;

//***************************************
//Reduce the coefficients of a series mod p^val
//***************************************
ReduceSeries := function(f,val)
  local coef,S,u,prec,pow;

  S := Parent(f);
  u := UniformizingElement(S);
  prec := Precision(CoefficientRing(S)); 

  coef,pow := Coefficients(f);
  coef := [ ChangePrecision(ChangePrecision(i,val),prec) : i in coef];

  return S!coef*u^pow;
end function;

//************************************
// Lower degree terms of a series
// Internal use
//************************************
LowerTerms := function(f,n)
  local i,P,a,b;

  P := Parent(f);
  a,b := Eltseq(f);
  if n-b ge 1 then
    a := a[[1..Min(n-b,#a)]];
    return P!a*P.1^b;
  end if;
  return 0;
end function;

//************************************
// Higher degree terms of a series
// Internal use
//************************************
HigherTerms := function(f,n)
  local i,P,a,b;

  P := Parent(f);
  a,b := Eltseq(f);
  a := a[[Max(n+1-b,1)..#a]];
  if b-n-1 ge 0 then 
    return P!a*P.1^(b-n);
  end if;
  return P!a;
end function;

/*******************************************************************
 *                                                                 *
 *   ARITHMETICS                                                   *
 *   Authors: X. Caruso, D. Lubicz                                 *
 *                                                                 *
 *******************************************************************/

//*************************************
//Returns true if an element of S is equal to 
//zero
//*************************************
intrinsic IsZeroOverS(f::RngSerPowElt : nu:=0)->BoolElt
{
  IsZeroOverS(f::nu:=0) : RngSerPowElt:FldReElt -> BollElt

  Returns true if f is 0
}
  local coef,i;
  coef := Coefficients(f);
  for i := 1 to #coef do
    if RelativePrecision(coef[i]) gt 0 then return false; end if;
  end for;
  return true;
end intrinsic;

//*************************************
// Return the Gauss valuation
// and the Weierstrass degree
// of the element f in S
// NB: If f = 0, returns 0, -1
//*************************************
intrinsic pValuationDegreeOverS(f::RngSerPowElt : nu:=0)->RngIntElt,RngIntElt
{
pValuationDegreeOverS(f: nu:=0) : RngSerPowElt:FldReElt -> RngIntElt,RngIntElt

Return the Gauss valuation
and the Weierstrass degree
of the element f in S
NB: If f = 0, returns 0, -1
}

  require Type(nu) in {FldRatElt,RngIntElt}:"nu must be in Z or Q";

  coef,pow := Coefficients(f);
  val := 0; deg := 0;
  vadd := pow*nu;
  for i := 1 to #coef do
    if RelativePrecision(coef[i]) eq 0 then continue; end if;
    v := Valuation(coef[i])+vadd;
    if deg le 0 or v lt val then val := v; deg := i; end if;
    vadd := vadd + nu;
  end for;
  return val,deg+pow-1;
end intrinsic;

//*************************************
// Return the Gauss valuation
// of the element f in S
// NB: If f = 0, returns Precision(Parent(f))+1
//*************************************
intrinsic pValuationOverS(f::RngSerPowElt: nu:=0)->RngIntElt
{
pValuationDegreeOverS(f: nu:=0) : RngSerPowElt:FldReElt -> RngIntElt

Return the Gauss valuation
of the element f in S
NB: If f = 0, returns Precision(Parent(f))+1
}
 
  val, deg := pValuationDegreeOverS(f : nu:=nu);

  if deg ge 0 then return val;
  else return Precision(Parent(f))+1;
  end if;
end intrinsic;

//*************************************
// Test if f is distinguished
// and returns its Weiestrass degree
// NB: If f not distinguished returns false, _
//*************************************
intrinsic IsDistinguishedOverS(f::RngSerPowElt: nu:=0)->BoolElt
{
  IsDistinguishedOverS(f:nu:=0)->RngSerPowElt:FldReElt-> BoolElt

Test if f is distinguished
and returns its Weierstrass degree
NB: If f not distinguished returns false, _
}

val, deg := pValuationDegreeOverS(f:nu:=nu);
  if val eq 0 then return true, deg;
  else return false, _;
  end if;
end intrinsic;

//***************************************
//Divide all the elements of a series x by a power
//val of the valuation
//***************************************
DivideSeries := function(f,val)
  local u,prec,i,p,d,res,coef,pow,c;

  S := Parent(f);
  u := UniformizingElement(S);
  p := UniformizingElement(CoefficientRing(S));
  d := p^val;
  pPrecision := Precision(CoefficientRing(S));
  res := S!0;
  coef, pow := Coefficients(f);

  for i := 1 to #coef do
    if Valuation(coef[i]) ge val then
      c := CoefficientRing(S)!(coef[i]/d);
      res := res + S!(u^(pow+i-1)*c);
    else
      ChangePrecision(~res,pow+i-1);
      break i;
    end if;
  end for;
  return res;
end function;

//************************************
// EuclideanWithDegreeWithPrecOverS algorithm returns 
// Q, R such that:
//   A = B*Q + R
//   d is the Weierstrass degree of B
//   we suppose that A and B have infinite u-precision (i.e are
//   polynomials)
//   and ensure that Q and R are known up to piadic-precision piprec
//************************************
intrinsic EuclideanWithPrecWithDegreeOverS(A::RngSerPowElt,B::RngSerPowElt,d::RngIntElt, piprec::RngIntElt : nu:=0)->RngSerPowElt,RngSerPowElt
{
  Q,R := EuclideanWithPrecWithDegreeOverS(A,B,d, piprec : nu:=0):RngSerPowElt,RngSerPowElt,RngIntElt, RngIntElt:FldeElt -> RngSerPowElt,RngSerPowElt

EuclideanWithDegreeWithPrecOverS algorithm returns 
Q, R such that:
A = B*Q + R
d is the Weierstrass degree of B
we suppose that A and B have infinite u-precision (i.e are
polynomials)
and ensure that Q and R are known up to padic-precision piprec
}
  local Blow,Bhiinv,Alow,Ahi,Q,Qadd,prec;

  S := Parent(A);
  prec := Precision(S);
  uprec := d*(piprec+1);

  if d eq 0 then return A/B, S!0; end if;
  ChangePrecision(~B, uprec);
  ChangePrecision(~A, uprec);
  Blow := LowerTerms(B,d);
  Bhiinv := HigherTerms(B,d)^-1;

  Q := S!0;
  while true do
    Alow := LowerTerms(A,d);
    Ahi := HigherTerms(A,d);
    if pValuationOverS(Ahi:nu:=nu) ge piprec then break; end if;
    Qadd := Ahi*Bhiinv;
    Q := Q+Qadd;
    A := Alow - Blow*Qadd;
    
    uprec := uprec -d;
    ChangePrecision(~A, uprec);
    ChangePrecision(~Q, uprec);
  end while;
  
  ChangePrecision(~A,prec);
  ChangePrecision(~Q,prec);
  return Q,A;
end intrinsic;

//************************************
// same thing as before without degree 
//************************************
intrinsic EuclideanWithPrecOverS(A::RngSerPowElt,B::RngSerPowElt,piprec::RngIntElt : nu:=0)->RngSerPowElt,RngSerPowElt
{
Q,R := EuclideanWithPrecWithDegreeOverS(A,B, piprec : nu:=0):RngSerPowElt,RngSerPowElt,RngIntElt:FldeElt -> RngSerPowElt,RngSerPowElt

EuclideanWithDegreeWithPrecOverS algorithm returns 
Q, R such that:
A = B*Q + R
we suppose that A and B have infinite u-precision (i.e are
polynomials)
and ensure that Q and R are known up to padic-precision piprec
}
  local d;
  _, d := pValuationDegreeOverS(B:nu:=nu);
  return EuclideanWithPrecWithDegreeOverS(A,B,d,piprec:nu:=nu); 
end intrinsic;

//************************************
//  same thing as before without precision
//************************************
intrinsic EuclideanOverS(A::RngSerPowElt,B::RngSerPowElt : nu:=0)->RngSerPowElt,RngSerPowElt
{
Q,R := EuclideanWithPrecWithDegreeOverS(A,B, : nu:=0):RngSerPowElt,RngSerPowElt:FldeElt -> RngSerPowElt,RngSerPowElt

EuclideanWithDegreeWithPrecOverS algorithm returns 
Q, R such that:
A = B*Q + R
we suppose that A and B have infinite u-precision (i.e are
polynomials)
and ensure that Q and R are known up to maximal p-adic precision
}
  local piprec;
  piprec := Precision(CoefficientRing(Parent(A)));
  return EuclideanWithPrecOverS(A,B,piprec:nu:=nu); 
end intrinsic;

//************************************
// Weierstrass preparation
// Let f be a distinguished element of S 
// (i.e.  which Gauss valuation is 0)
// with Weierstrass degree d
//write f = UP where U is invertible in S
//and P is a polynomial of degree d 
//************************************
intrinsic WeierstrassOverS(f::RngSerPowElt : nu:=0)->RngSerPowElt, RngSerPowElt
{
  U,P := WeierstrassOverS(f::RngSerPowElt : nu:=0): RngSerPowElt : FldReElt -> RngSerPowElt, RngSerPowElt

Weierstrass preparation
Let f be a distinguished element of S 
(i.e.  which Gauss valuation is 0)
with Weierstrass degree d
write f = UP where U is invertible in S
and P is a polynomial of degree d 
}
  local mbool, d,P,q,r;

  P := Parent(f);
  mbool,d := IsDistinguishedOverS(f:nu:=nu); 

  if mbool eq false then return P!0,P!0; end if;
  q,r := EuclideanOverS(P.1^d,f:nu:=nu);
  return q^(-1),P.1^d - r;
end intrinsic;


//************************************
//let nu=a/b with a,b relatively primes
//let chi such that chi^b=pi
// Extended gcd
// returns d,d',u,u',k,k' such that
// . chi^pia au + chi^pib bk = d with Valuaton(d)=Valuation(b)
// . au' + bk' = d' with Valuation(d') > Valuation(b)
// . uk' - u'k = 1
// Moreover d and d' correspond to two consecutives
// steps in the extended Euclid algorithm
// Here, v is the valuation of b
// and a is assumed to have valuation >= v
// moreover flag = 1 or 2 and flag is the number of
// calls to EnclideanOverS mod 2
//************************************
intrinsic ExtendedGcdWithValuationWithPipowersWithPrecOverS(a::RngSerPowElt,b::RngSerPowElt,v::RngIntElt,pia::FldRatElt,pib::FldRatElt,piprec::RngIntElt : nu:=0)->RngSerPowElt,RngSerPowElt,RngSerPowElt,RngSerPowElt,RngSerPowElt,RngSerPowElt,RngIntElt
{
  d,d1,u,u1,k,k1 := ExtendedGcdWithValuationWithPipowersWithPrecOverS(a,b,v,pia,pib,piprec : nu:=0): RngSerPowElt,RngSerPowElt,RngIntElt,FldRatElt,FldRatElt,RngIntElt : FldReElt-> RngSerPowElt,RngSerPowElt,RngSerPowElt,RngSerPowElt,RngSerPowElt,RngSerPowElt,RngIntElt

 
let nu=a/b with a,b relatively primes
let chi such that chi^b=pi
Extended gcd
returns d,d1,u,u1,k,k1 such that
. chi^pia au + chi^pib bk = d with Valuaton(d)=Valuation(b)
. au1 + bk1 = d1 with Valuation(d1) > Valuation(b)
. uk1 - u1k = 1
Moreover d and d1 correspond to two consecutives
steps in the extended Euclid algorithm
Here, v is the Gauss valuation of b
and a is assumed to have valuation >= v
moreover flag = 1 or 2 and flag is the number of
calls to EuclideanOverS mod 2
}
  local S,vv,q,r,Ma,Mb,Mx,pi,flag;

  S := Parent(a);
  Ma := Matrix(S,1,2,[1,0]);
  Mb := Matrix(S,1,2,[0,1]);
  pi := [pia,pib];
  flag := 2;

  while true do
    vv,d := pValuationDegreeOverS(b:nu:=nu);
    if d lt 0 or vv+pi[flag] gt v then break; end if;
    flag := 3-flag;
    q,r := EuclideanWithPrecOverS(a,b,piprec:nu:=nu);
    a := b; b := -r;
    // We have always Mb [a b]^t = r
    Mx := Ma; Ma := Mb; Mb := -Mx+q*Mb;
  end while;

  return a,b,Ma[1,1],Ma[1,2],Mb[1,1],Mb[1,2],flag;
end intrinsic;

//************************************
//Same as before without pipowers
//without prec
//************************************
intrinsic ExtendedGcdWithValuationOverS(a::RngSerPowElt,b::RngSerPowElt,v::RngIntElt : nu:=0)->RngSerPowElt,RngSerPowElt,RngSerPowElt,RngSerPowElt,RngSerPowElt,RngSerPowElt,RngIntElt

{
d,d1,u,u1,k,k1 := ExtendedGcdWithValuationOverS(a,b,v: nu:=0): RngSerPowElt,RngSerPowElt,RngIntElt: FldReElt-> RngSerPowElt,RngSerPowElt,RngSerPowElt,RngSerPowElt,RngSerPowElt,RngSerPowElt,RngIntElt

 

let nu=a/b with a,b relatively primes
let chi such that chi^b=pi
Extended gcd
returns d,d',u,u',k,k' such that
. chi^pia au + chi^pib bu' = d with Valuaton(d)=Valuation(b)
. ak + bk' = d' with Valuation(d') > Valuation(b)
. uk' - u'k = 1
Moreover d and d' correspond to two consecutives
steps in the extended Euclid algorithm
Here, v is the valuation of b
and a is assumed to have valuation >= v
moreover flag = 1 or 2 and flag is the number of
calls to EnclideanOverS mod 2
}

  local piprec;
  piprec := Precision(CoefficientRing(Parent(a)));
  return ExtendedGcdWithValuationWithPipowersWithPrecOverS(a,b,v,nu,0,0,piprec:nu := nu);
end intrinsic;


//************************************
//Same as before without valuation
//************************************
intrinsic ExtendedGcdWithPipowersWithPrecOverS(a::RngSerPowElt,b::RngSerPowElt,pia::FldRatELt,pib::FldRatElt,piprec::RngIntElt : nu:=0)->RngSerPowElt,RngSerPowElt,RngSerPowElt,RngSerPowElt,RngSerPowElt,RngSerPowElt,RngIntElt

{
   d,d1,u,u1,k,k1 := ExtendedGcdPipowersWithPrecOverS(a,b,pia,pib,piprec : nu:=0): RngSerPowElt,RngSerPowElt,FldRatElt,FldRatElt,RngIntElt : FldReElt-> RngSerPowElt,RngSerPowElt,RngSerPowElt,RngSerPowElt,RngSerPowElt,RngSerPowElt,RngIntElt

 
let nu=a/b with a,b relatively primes
let chi such that chi^b=pi
Extended gcd
returns d,d',u,u',k,k' such that
. chi^pia au + chi^pib bu' = d with Valuaton(d)=Valuation(b)
. ak + bk' = d' with Valuation(d') > Valuation(b)
. uk' - u'k = 1
Moreover d and d' correspond to two consecutives
steps in the extended Euclid algorithm
Here, v is the valuation of b
and a is assumed to have valuation >= v
moreover flag = 1 or 2 and flag is the number of
calls to EnclideanOverS mod 2
}
  local v;
  v,_ := pValuationDegreeOverS(b:nu:= nu); 
  return ExtendedGcdWithValuationWithPipowersWithPrecOverS(a,b,v,pia,pib,piprec:nu := nu);
end intrinsic;

//************************************
//Same as before without pipowers
//************************************
intrinsic ExtendedGcdWithPrecOverS(a::RngSerPowElt,b::RngSerPowElt,piprec::RngIntElt : nu:=0)->RngSerPowElt,RngSerPowElt,RngSerPowElt,RngSerPowElt,RngSerPowElt,RngSerPowElt,RngIntElt

{
d,d1,u,u1,k,k1 := ExtendedGcdWithValuationWithPipowersWithPrecOverS(a,b,v,pia,pib,piprec : nu:=0): RngSerPowElt,RngSerPowElt,RngIntElt,FldRatElt,FldRatElt,RngIntElt : FldReElt-> RngSerPowElt,RngSerPowElt,RngSerPowElt,RngSerPowElt,RngSerPowElt,RngSerPowElt,RngIntElt

let nu=a/b with a,b relatively primes
let chi such that chi^b=pi
Extended gcd
returns d,d',u,u',k,k' such that
. chi^pia au + chi^pib bu' = d with Valuaton(d)=Valuation(b)
. ak + bk' = d' with Valuation(d') > Valuation(b)
. uk' - u'k = 1
Moreover d and d' correspond to two consecutives
steps in the extended Euclid algorithm
Here, v is the valuation of b
and a is assumed to have valuation >= v
moreover flag = 1 or 2 and flag is the number of
calls to EnclideanOverS mod 2
}

  return ExtendedGcdWithValuationWithPipowersWithPrecOverS(a,b,nu,0,0,piprec);
end intrinsic;

//************************************
//Same as before without piprec
//************************************
intrinsic ExtendedGcdOverS(a::RngSerPowElt,b::RngSerPowElt : nu:=0)->RngSerPowElt,RngSerPowElt,RngSerPowElt,RngSerPowElt,RngSerPowElt,RngSerPowElt,RngIntElt

{
d,d1,u,u1,k,k1 := ExtendedGcdWithValuationWithPipowersWithPrecOverS(a,b,v,pia,pib,piprec : nu:=0): RngSerPowElt,RngSerPowElt,RngIntElt,FldRatElt,FldRatElt,RngIntElt : FldReElt-> RngSerPowElt,RngSerPowElt,RngSerPowElt,RngSerPowElt,RngSerPowElt,RngSerPowElt,RngIntElt

let nu=a/b with a,b relatively primes
let chi such that chi^b=pi
Extended gcd
returns d,d',u,u',k,k' such that
. chi^pia au + chi^pib bu' = d with Valuaton(d)=Valuation(b)
. ak + bk' = d' with Valuation(d') > Valuation(b)
. uk' - u'k = 1
Moreover d and d' correspond to two consecutives
steps in the extended Euclid algorithm
Here, v is the valuation of b
and a is assumed to have valuation >= v
moreover flag = 1 or 2 and flag is the number of
calls to EnclideanOverS mod 2
}
  local piprec,v;
  v,_ := pValuationDegreeOverS(b:nu:= nu); 
  piprec := Precision(CoefficientRing(Parent(a)));
  return ExtendedGcdWithValuationWithPipowersWithPrecOverS(a,b,v,RationalField()!0,RationalField()!0,piprec:nu := nu);
end intrinsic;

//************************************
//let nu=a/b with a,b relatively primes
// Extended gcd
// returns u,v,d=pgcd(a,b) such that
//  au + bv = d
//************************************
intrinsic ExtendedGcdWithValuationWithPrecOverSp(a::RngSerPowElt,b::RngSerPowElt,piprec::RngIntElt:nu:=0)->RngSerPowElt, RngSerPowElt, RngSerPowElt, RngSerPowElt, RngSerPowElt, RngSerPowElt

{
  d,d1,u,u1,k,k1 := ExtendedGcdWithValuationWithPrecOverSp(a::RngSerPowElt,b::RngSerPowElt,piprec::RngIntElt:nu:=0) : RngSerPowElt,RngSerPowElt,RngIntElt:FldReElt -> RngSerPowElt, RngSerPowElt, RngSerPowElt, RngSerPowElt, RngSerPowElt, RngSerPowElt

let nu=a/b with a,b relatively primes, this function
  returns d,d1,u,u1,k,k1 :
. chi^pia au + chi^pib bu' = d with Valuaton(d)=Valuation(b)
. ak + bk' = d' with Valuation(d') > Valuation(b)
. uk' - u'k = 1
}
  local S,q,r,Ma,Mb,Mx;

  S := Parent(a);
  Ma := Matrix(S,1,2,[1,0]);
  Mb := Matrix(S,1,2,[0,1]);

  while true do
    if IsZeroOverS(b)  then break; end if;
    q,r := EuclideanWithPrecOverS(a,b,piprec:nu := nu);
    a := b; b := -r;
    Mx := Ma; Ma := Mb; Mb := -Mx+q*Mb;
  end while;

  return a,b,Ma[1,1],Ma[1,2],Mb[1,1],Mb[1,2];
end intrinsic;

//************************************
//same as before without piprec
//************************************
intrinsic ExtendedGcdWithValuationOverSp(a::RngSerPowElt,b::RngSerPowElt:nu:=0)->RngSerPowElt, RngSerPowElt, RngSerPowElt, RngSerPowElt, RngSerPowElt, RngSerPowElt
{
  d,d1,u,u1,k,k1  := ExtendedGcdWithValuationOverSp(a::RngSerPowElt,b::RngSerPowElt:nu:=0) : RngSerPowElt,RngSerPowElt,RngIntElt:FldReElt -> RngSerPowElt, RngSerPowElt, RngSerPowElt, RngSerPowElt, RngSerPowElt, RngSerPowElt

let nu=a/b with a,b relatively primes, this function
returns d,u,v=pgcd(a,b) such that
au + bv = d
}
  local piprec;
  piprec := Precision(CoefficientRing(Parent(a)));
  return ExtendedGcdWithValuationWithPrecOverSp(a,b,piprec:nu := nu);
end intrinsic;
 


//************************************
// same as before without valuation
//************************************
intrinsic ExtendedGcdOverSp(a::RngSerPowElt,b::RngSerPowElt :nu:=0)->RngSerPowElt, RngSerPowElt, RngSerPowElt, RngSerPowElt, RngSerPowElt, RngSerPowElt
{
  d,d1,u,u1,k,k1  := ExtendedGcdOverSp(a::RngSerPowElt,b::RngSerPowElt:nu:=0) : RngSerPowElt,RngSerPowElt:FldReElt -> RngSerPowElt, RngSerPowElt, RngSerPowElt, RngSerPowElt, RngSerPowElt, RngSerPowElt

let nu=a/b with a,b relatively primes, this function
  returns d,d1,u,u1,k,k1 :
. chi^pia au + chi^pib bu' = d with Valuaton(d)=Valuation(b)
. ak + bk' = 0
. uk' - u'k = 1
}
  local v;
  v,_ := pValuationDegreeOverS(b:nu:=nu);
  return ExtendedGcdWithValuationOverSp(a,b:nu:=nu);
end intrinsic;

//************************************
// same as before but returns only the gcd
//************************************
intrinsic GcdOverSp(a::RngSerPowElt,b::RngSerPowElt :nu:=0)->RngSerPowElt
{
d := GcdOverSp(a::RngSerPowElt,b::RngSerPowElt:nu:=0) : RngSerPowElt,RngSerPowElt:FldReElt -> RngSerPowElt

let nu=a/b with a,b relatively primes, this function
  returns d,d1,u,u1,k,k1 :
. chi^pia au + chi^pib bu' = d with Valuaton(d)=Valuation(b)
. ak + bk' = 0
. uk' - u'k = 1
}
  local gcd;

  gcd,_ := ExtendedGcdOverSp(a,b:nu:=nu);
  return gcd; 
end intrinsic;

//*************************************
// pseudo-inverse let s \in W(k)[[u]]
// whose Gauss Valuation is O
// compute t \in W(k)[[u]] and n\in N s.t. s.t=u^n (mod p^pAdicModulo)
// not efficient at all
//*************************************
InvertOverSu := function(f,pAdicModulo,nu)
  local S,upow,p,F,R,val,qdeg,n,uprec,q,list,uval,inv,rem,rem1,pval;

  S := Parent(f);
  upow := UniformizingElement(S);
  p := UniformizingElement(CoefficientRing(S));
  F := ResidueClassField(CoefficientRing(S));
  R := CoefficientRing(S);

  val,qdeg := pValuationDegreeOverS(f);

  if val gt 0 then return 0,-1; end if;

  n:= 0;
  uprec := Precision(S)-qdeg;
  list, uval := Eltseq(f);
  q := [ F!i: i in list | (pValuationOverS(f,nu) eq 0)];
  q := (S![(R!i)*p^(-Floor(nu*i)) : i in q])^(-1);
  inv := q;
  rem := upow^qdeg;
  pval := 0;

  while (not(IsZeroOverS(rem)) and pValuationOverS(rem,nu) le pAdicModulo) do
    rem1 := f*inv-rem;
    val, deg := pValuationDegreeOverS(rem1);
    rem1 := DivideSeries(rem1,val);
    pval := pval + val;
    if deg gt qdeg then
      inv := inv*upow^(deg-qdeg);
      rem := rem*upow^(deg-qdeg);
      rem1 := rem1*upow^(deg-qdeg);
      n := n+deg-qdeg; uprec := uprec - deg+qdeg;
    end if;
    list, _ := Eltseq(f);
    rem1 := S![ R!i: i in list | (pValuationOverS(rem1,nu) eq 0)];
    rem1 := S!Eltseq(rem1); 
    inv := inv + p^pval*S!Eltseq(rem1)*q;
  end while;

  return inv, n, uprec;
end function;

/*******************************************************************
 *                                                                 *
 *   MATRICES over Zp                                              *
 *   Authors: X. Caruso, D. Lubicz                                 *
 *                                                                 *
 *******************************************************************/
LDecompositionOverR := function(M)
  local i,j;
  for  i := 1 to Nrows(M) do
    MultiplyColumn(~M,1/M[i,i],i);
    for j := i+1 to Ncols(M) do
      AddColumn(~M,-M[i,j],i,j);
    end for;
  end for;
  return M;
end function;

Valuation := function(M)
  local i,j;

val := Valuation(M[1,1]);

for i :=1 to Nrows(M) do
  for j := 1 to Ncols(M) do
    if val lt Valuation(M[i,j]) then
    val := Valuation(M[i,j]);
    end if;
  end for;
end for;

return val;
end function;

/*******************************************************************
 *                                                                 *
 *   MATRICES over S                                               *
 *   Authors: X. Caruso, D. Lubicz                                 *
 *                                                                 *
 *******************************************************************/
//*************************************
// remove the zero coefficients when printing
// an element of S
//*************************************
RemoveZeroesOverMatS := function (M)
  local i,j;
  for i := 1 to Nrows(M) do
    for j := 1 to Ncols(M) do
      M[i,j] := RemoveZeroesOverS (M[i,j]);
    end for;
  end for;

  return M;
end function;

//*************************************
//Compute the Gauss valuation of a matrix M
//with coefficients in M
//and returns i,j such that v(M[i,j])=v(M)
//*************************************
pValuationMatOverS := function(M,nu)
  local i,j,val,val1;

  val := pValuationDegreeOverS(M[1,1],nu);
  for i := 1 to Nrows(M) do
    for j := 1 to Ncols(M) do
      val1 := pValuationDegreeOverS(M[i,j],nu);
      if val1 lt val then val := val1; end if;
    end for;
  end for;
  return val;
end function;


//************************************
// Perform simoultaneously the following 
// elementary operations on the columns
// of the matrix M:
//   C_(j1) <- u*C_(j1) + v*C_(j2)
//   C_(j2) <- up*C_(j1) + vp*C_(j2)
///************************************
LinearCombinaisonColumns := procedure (~M,j1,j2,u,v,up,vp) 
  local i, s;
  for i := 1 to Nrows(M) do
    s := M[i,j1];
    M[i,j1] := u*M[i,j1] + v*M[i,j2];
    M[i,j2] := up*s + vp*M[i,j2];
  end for;
end procedure;

//************************************
// Perform simoultaneously the following 
// elementary operations on the rows
// of the matrix M:
//   R_(j1) <- u*R_(j1) + v*R_(j2)
//   R_(j2) <- up*R_(j1) + vp*R_(j2)
//************************************
LinearCombinaisonRows := procedure (~M,j1,j2,u,v,up,vp) 
  local i, s;
  for i := 1 to Ncols(M) do
    s := M[j1,i];
    M[j1,i] := u*M[j1,i] + v*M[j2,i];
    M[j2,i] := up*s + vp*M[j2,i];
  end for;
end procedure;

//*************************************
//Test whether a matrix with coefficient in S
//is zero
//*************************************
intrinsic IsZeroMatOverS(M::Mtrx) -> BoolElt
{
IsZeroMatOverS : AltMatElt -> BoolElt
Test if a Matrix M with coefficients in S is the zero matrix
}
  local i,j;
  for i := 1 to Nrows(M) do
    for j := 1 to Ncols(M) do
      if not(IsZeroOverS(M[i,j])) then return false; end if;
    end for;
  end for;
  return true;
end intrinsic;

//***************************************
//Compute the p-valuation of column i starting
//from row j to row k and returns the first position
//of an element with a smallest valuation
//***************************************
pColValuationWithPipowersOverS := function(M,Pipowers,i,j,k,nu)
  local m, l, val, val1;

  val := pValuationDegreeOverS(M[j,i],nu)+Pipowers[1,j];
  l:=j;
  for m := j+1 to k do
    val1 := pValuationDegreeOverS(M[m,i],nu)+Pipowers[1,m];
    if val1 lt val then
      val := val1; l := m;
    end if;
  end for;
  return val, l;
end function;

//***************************************
//Compute the p-valuation of column i starting
//from row j to row k and returns the first position
//of an element with a smallest valuation
//***************************************
pColValuationOverS := function(M,i,j,k)
  local m, l, val, val1;

  val := pValuationDegreeOverS(M[j,i]);
  l:=j;
  for m := j+1 to k do
    val1 := pValuationDegreeOverS(M[m,i]);
    if val1 lt val then
      val := val1; l := m;
    end if;
  end for;
  return val, l;
end function;



//***************************************
//Compute the p-valuation of row i starting
//from col j to col k and returns the first position
//of an element with a smallest valuation
//***************************************
pRowValuationOverS := function(M,Pipowers,i,j,k,nu)
  local m, l, val, val1;

  val := pValuationDegreeOverS(M[i,j],nu)+Pipowers[1,j];
  l:=j;
  for m := j+1 to k do
    val1 := pValuationDegreeOverS(M[i,m],nu)+Pipowers[1,m];
    if val1 lt val then
      val := val1; l := m;
    end if;
  end for;
  return val, l;
end function;

//**************************************
//Compute the row echelon form of a matrix over Sp
//***************************************
intrinsic RowEchelonOverSp(M::Mtrx,boolND::BoolElt:nu:=0) -> Mtrx,Mtrx,RngIntElt
{
  M1,ND,pAdicModulo:= RowEchelonOverSp(M,boolND:nu:=0) : Mtrx, BoolElt:FldReElt -> Mtrx,Mtrx,RngIntElt

Compute the row echelon form of a matrix over Sp
Compute a matrix M1 in row echelon form such that there exists an
invertible matrix ND such that M1=M ND
If boolND is equel to true it also returns the matrix ND
pAdicModulo is the maximum of the Gauss valuation of the non null elements in the diagonal of M1
}
  local M1,ND,S,pPrecision,p,i,j,vii,vij,u,v,up,vp,inc,pAdicModulo;

  pAdicModulo := 0;
  M1 := M;
  S := Parent(M1[1,1]);
  pPrecision := Precision(CoefficientRing(S));
  p := UniformizingElement(CoefficientRing(S));

  if boolND then ND := IdentityMatrix(S, Ncols(M1)); else ND := 0; end if;

  inc := 1;
  for i := 1 to Nrows(M1) do
    j := inc;
    //find some non zero M1[i,j] 
    while (j le Ncols(M1)) and IsZeroOverS(M1[i,j]) do j := j+1; end while; // Euh, il se passe quoi s'il n'y a que des 0 ? // ben on passe a la ligne suivante
    if j le Ncols(M1) then
      SwapColumns (~M1, inc, j);
      if boolND then SwapColumns (~ND, inc, j); end if;
      //cancel all M1[i,j] for j>inc
      j := inc+1;
      while (j le Ncols(M1)) do
        if not(IsZeroOverS(M1[i,j])) then
	_,_,u,v,up,vp := ExtendedGcdOverS(M1[i,inc], M1[i,j]:nu:=nu);
          LinearCombinaisonColumns (~M1, inc, j, u, v, up, vp);
          if boolND then LinearCombinaisonColumns (~ND, inc, j, u, v, up, vp); end if;
        end if;
        j := j + 1;
      end while;
      pAdicModulo := Max(pAdicModulo, pValuationOverS(M1[i,inc]:nu:=nu));
      inc := inc+1;
    end if;
  end for;
  return M1,ND,pAdicModulo;
end intrinsic;

//**************************************
//Compute the Hermite form H of a matrix M with respect to p
//returns the transformation matrix such that
//H = M ND 
//if boolND is true and pAdicModulo
//***************************************
intrinsic HermiteOverSp(M::Mtrx,boolND::BoolElt:nu:=0) -> Mtrx,Mtrx,RngIntElt
{
  M1,ND,pAdicModulo:= HermiteOverSp(M,boolND:nu:=0) : Mtrx, BoolElt:FldReElt -> Mtrx,Mtrx,RngIntElt

Compute the Hermite form of a matrix over Sp
Compute a matrix M1 in Hermite form such that there exists an
invertible matrix ND such that M1=M ND
If boolND is equel to true it also returns the matrix ND
pAdicModulo is the maximum of the Gauss valuation of the non null elements in the diagonal of M1
}

  local ND, val, M1,S,pPrecision,p,inc,i,j,val1,g,f,q,r;

  M1 := M;
  S := Parent(M1[1,1]);
  pPrecision := Precision(CoefficientRing(S));
  p := UniformizingElement(CoefficientRing(S));

  //first compute the RowEchelon form
  M1,ND,pAdicModulo := RowEchelonOverSp(M,boolND:nu:=nu);
  val1 := pValuationOverS(M1[1,1]:nu:=nu);
  M1[1,1] := DivideSeries(M1[1,1],val1);
  inc := 2;
  //then compute the Hermite form
  for i := 2 to Ncols(M1) do
    while IsZero(M1[inc,i]) do
        inc := inc +1;
    end while;
    val := pValuationOverS(M1[inc,i]:nu:=nu);
    for j := 1 to i-1 do
     val1 := pValuationOverS(M1[inc,j]:nu:=nu);
     if val1 lt  val then
       MultiplyColumn(M,p^(val-pValuationOverS(M1[inc,j]:nu:=nu)),j);
       val1 := val;
       if boolND then
	 MultiplyRow(ND,p^(val-pValuationOverS(M1[inc,j]:nu:=nu)),j);
       end if;
     end if;
     g := DivideSeries(M1[inc,j],val1);
     f := DivideSeries(M1[inc,i],val);
     q,r := EuclideanOverS(g,f:nu:=0);
     LinearCombinaisonColumns(~M1,j,i,1,-q*p^(val1-val),0,1);
     if boolND then
      LinearCombinaisonColumns(~ND,j,i,1,-q*p^(val1-val),0,1); 
     end if;
    end for;
    if pValuationOverS(M1[inc,i]:nu:=nu) gt pAdicModulo then
		    pAdicModulo := pValuationOverS(M1[inc,i]:nu:=nu);
    end if;
  end for;

  return M1,ND,pAdicModulo;
end intrinsic;

//**************************************
//Compute the Smith form S of a matrix M with respect to p
//if boolN returns the transformation matrices 
//NG and ND such that
//S = NG M ND
//***************************************
intrinsic SmithOverSp(M::Mtrx,boolN::BoolElt:nu:=0) -> Mtrx,Mtrx,RngIntElt
{
  M1,NG,ND,pAdicModulo:= SmithOverSp(M,boolN:nu:=0) : Mtrx, BoolElt:FldReElt -> Mtrx,Mtrx,RngIntElt

Compute the Hermite form of a matrix over Sp
Compute a matrix M1 in Hermite form such that there exists an
invertible matrices NG ND such that M1=NG M ND
If boolND is equel to true it also returns the matrices NG,ND
pAdicModulo is the maximum of the Gauss valuation of the non null elements in the diagonal of M1
}
  local M1,NG,ND,S,pPrecision,p,i,j,vii,vij,u,v,up,vp,inc,pAdicModulo,NotZeroRow;

  pAdicModulo := 0;
  M1 := M;
  S := Parent(M1[1,1]);
  pPrecision := Precision(CoefficientRing(S));
  p := UniformizingElement(CoefficientRing(S));

  if boolN then 
    ND := IdentityMatrix(S, Ncols(M1)); 
    NG := IdentityMatrix(S, Nrows(M1)); 
  else 
    ND := 0; 
    NG := 0;
  end if;

  inc := 1;
  i := 1;
  while i le Nrows(M1) do
  NotZeroRow := true;
    while (NotZeroRow) do
      NotZeroRow := false;
      j := inc;
      //find some non zero M1[i,j] 
      while (j le Ncols(M1)) and IsZeroOverS(M1[i,j]) do j := j+1; end while; 
      if j le Ncols(M1) then
        SwapColumns (~M1, inc, j);
        if boolN then SwapColumns (~ND, inc, j); end if;
        //cancel all M1[i,j] for j>inc
        j := inc+1;
        while (j le Ncols(M1)) do
          if not(IsZeroOverS(M1[i,j])) then
            NotZeroRow := true;
	    _,_,u,v,up,vp := ExtendedGcdOverSp(M1[i,inc], M1[i,j]:nu:=nu);
            LinearCombinaisonColumns (~M1, inc, j, u, v, up, vp);
            if boolN then LinearCombinaisonColumns (~ND, inc, j, u, v, up, vp); end if;
          end if;
          j := j + 1;
        end while;
      end if;
      //cancel all M1[j,inc] for j>i
      j := i+1;
      while (j le Nrows(M1)) do
        if not(IsZeroOverS(M1[j,inc])) then
	   _,_,u,v,up,vp := ExtendedGcdOverS(M1[i,inc], M1[j,inc]:nu:=nu);
           LinearCombinaisonRows (~M1, i, j, u, v, up, vp);
           if boolN then LinearCombinaisonRows (~NG, i, j, u, v, up, vp); end if;
        end if;
        j := j + 1;
      end while;  
    end while;
    pAdicModulo := Max(pAdicModulo, pValuationOverS(M1[i,inc]:nu:=nu));
    if not(IsZero(M1[i,inc])) then
      inc := inc+1;
    end if;
     i:= i+1;
  end while;
  return M1,NG,ND,pAdicModulo;
end intrinsic;

//*************************************
// Compute the relations of a matrix over Sp
// that is given a matrix M compute R such that
// M.R =0
//*************************************
intrinsic RelationsOverSp(M::Mtrx:nu:=0) -> Mtrx,RngIntElt
{
  M1,pAdicModulo:= RelationsOverSp(M:nu:=0) : Mtrx:FldReElt -> Mtrx,RngIntElt

Compute the relations matrix over Sp of the matrix M that is given a matrix M compute R of maximal rang such that M.R=0
pAdicModulo is the maximum of the Gauss valuation of the non null elements in the diagonal of M once put in row echolon form.
}

  local M1,ND1,R,i,first,pAdicModulo;

  M1 := M;
  M1,ND1,pAdicModulo := RowEchelonOverSp(M1,true:nu:=nu);
  R := Matrix(Parent(M1[1,1]), Ncols(M1), 0, []);

  for i := 1 to Ncols(M1) do
    if (IsZeroMatOverS(ColumnSubmatrix(M1,i,1))) then
      R := HorizontalJoin(R,ColumnSubmatrix(ND1,i,1));
    end if;
  end for;

  return R, pAdicModulo;
end intrinsic;


//**************************************
//Compute the row echelon form of a matrix with coefficients
//in Su
//***************************************
intrinsic RowEchelonOverSu(M::Mtrx,boolND::BoolElt,pAdicModulo::RingIntElt:nu:=0) -> Mtrx,Mtrx
{
  M1,ND:= RowEchelonOverSu(M,boolND,pAdicModulo:nu:=0) : Mtrx, BoolElt:FldReElt -> Mtrx,Mtrx

Compute the row echelon form of a matrix over Su. More precisely, this function
computes a matrix M1 in row echelon form such that there exists an
invertible matrix ND such that M1=M ND modulo p^pAdicModulo
If boolND is equel to true it also returns the matrix ND.
}

  local M1,ND,S,uPrecision,p,i,j,vii,vij,d,u,v,qii,qij,minval,col,m,inc;

  M1 := M;
  S := Parent(M1[1,1]);
  uPrecision := Precision(S);
  upow := UniformizingElement(S);
  p := UniformizingElement(CoefficientRing(S));
  pPrecision := Precision(CoefficientRing(S));

  if boolND then ND := IdentityMatrix(S, Ncols(M1)); else ND := 0; end if;

  inc := 1;
  for i := 1 to Nrows(M1) do
    minval := pValuationOverS(M1[i,inc]:nu:=nu);
    col := inc;
    for j := inc to Ncols(M1) do
      if pValuationOverS(M1[i,j]:nu:=nu) lt minval then
	minval := pValuationOverS(M1[i,j]:nu:=nu);
        col := j;
      end if;
    end for;
    if minval lt pPrecision then
      SwapColumns (~M1, inc, col);
      if boolND then SwapColumns (~ND, inc, col); end if;
      M1ii := DivideSeries(M1[i,inc], minval);
      inv, n := InvertOverSu(M1ii,pAdicModulo); // NDXavier: je pense qu'on peut mettre pAdicModulo-minval
      //cancel all M1[i,j] for j>i
      for j := inc+1 to Ncols(M1) do
        if not(IsZeroOverS(M1[i,j])) then
          vij := Valuation(M1[i,j]);
          if vij lt n then
            MultiplyColumn (~M1, upow^(n-vij),j);
            if boolND then MultiplyColumn (~ND, upow^(n-vij),j); end if;
            vij := n;
          end if;
          c := -upow^(vij-n)*inv*(Parent(M1[1,1])!Eltseq(M1[i,j]));
          AddColumn (~M1, c, i, j); 
          if boolND then AddColumn (~ND, c, i, j); end if;
        end if;
      end for;	
      MultiplyColumn(~M1,inv,inc);
      if boolND then MultiplyColumn(~ND,inv,inc); end if;
      inc := inc+1;
    end if;
  end for;

  return M1,ND;
end intrinsic;

//**************************************
//Compute the Hermite form of a matrix with respect to u
//with precision pAdicModulo
//***************************************
intrinsic HermiteOverSu(M::Mtrx,boolND::BoolElt,pAdicModulo::RingIntElt:nu:=0) -> Mtrx,Mtrx
{
  M1,ND:= HermiteOverSu(M,boolND,pAdicModulo:nu:=0) : Mtrx, BoolElt:FldReElt -> Mtrx,Mtrx

Compute the row echelon form of a matrix over Su. More precisely, this function
computes a matrix M1 in Hermite form such that there exists an
invertible matrix ND such that M1=M ND modulo p^pAdicModulo
If boolND is equel to true it also returns the matrix ND.
}
  local S,upow,p,pPrecision,M1,i,j,M1ii,inv,v,uval,uval1,minval,series;

  S := Parent(M[1,1]);
  upow := UniformizingElement(S);
  p := UniformizingElement(CoefficientRing(S));
  pPrecision := Precision(CoefficientRing(S));

  //first compute the RowEchelon form
  M1,ND := RowEchelonOverSu(M,boolND,pAdicModulo:nu:=nu);
  inc := 1;
  //then compute the Hermite form
  for i := 1 to Ncols(M1) do
    while IsZero(M1[inc,i]) do
       inc := inc +1;
    end while;
    minval := pValuationOverS(M1[inc,i]:nu:=nu);

    //now reduce columns to the right of i
    j := 1;
    uval := Valuation(M1[inc,i]);
    while j lt i do
      uval1 := Valuation(M1[inc,j]);
      if uval1 lt uval then
	MultiplyColumn(~M1,upow^(uval-uval1),j);
	if boolND then 
	MultiplyColumn(~ND,upow^(uval-uval1),j);
	end if;
	uval1 := uval;
      end if;
        series := DivideSeries(M1[inc,j]-ReduceSeries(M1[inc,j],minval),minval);
        series := S!Coefficients(series)*upow^(uval1-uval);
        LinearCombinaisonColumns(~M1,j,i,1,-series,0,1);
        if boolND then LinearCombinaisonColumns(~ND,j,i,1,-series,0,1);
	end if;
      j := j+1;
    end while;
  end for; 
  
  return M1,ND;
end intrinsic;

//**************************************
//Compute the Smith form of a matrix with respect to u
//with precision pAdicModulo
//if boolN returns the transformation matrices
//NG and ND such that
//S = NG M ND
//***************************************
intrinsic SmithOverSu(M::Mtrx,boolN::BoolElt,pAdicModulo::RingIntElt:nu:=0) -> Mtrx,Mtrx,Mtrx
{
  M1,ND:= SmithOverSu(M,boolN,pAdicModulo:nu:=0) : Mtrx, BoolElt:FldReElt -> Mtrx,Mtrx,Mtrx

Compute the Smith form of a matrix over Su. More precisely, this function
computes a matrix M1 in Hermite form such that there exists
invertible matrices NG, ND such that M1=NG M ND modulo p^pAdicModulo
If boolN is equel to true it also returns the matrices NG,ND.
}


  local M1,ND,NG,S,uPrecision,p,i,j,vii,vij,d,u,v,qii,qij,minval,col,m,inc;

  M1 := M;
  S := Parent(M1[1,1]);
  uPrecision := Precision(S);
  upow := UniformizingElement(S);
  p := UniformizingElement(CoefficientRing(S));
  pPrecision := Precision(CoefficientRing(S));

  if boolN then 
    NG := IdentityMatrix(S, Nrows(M1)); 
    ND := IdentityMatrix(S, Ncols(M1)); 
  else 
    NG := 0;
    ND := 0; 
  end if;

  inc := 1;
  cont_flag := true;
  while (cont_flag) do
  //look for the element of M1[i,j] with the smallest p-valuation
    minval := pValuationOverS(M1[inc,inc]:nu:=nu); col := inc; row := inc;
    for i := inc to Nrows(M1) do
      for j :=inc to Ncols(M1) do 
	if pValuationOverS(M1[i,j]:nu:=nu) lt minval then
	  minval := pValuationOverS(M1[i,j]:nu:=nu);
	  col := j;
	  row := i;
	end if;
      end for;
    end for;

    if minval eq pPrecision then
      cont_flag := false;
    else
      SwapColumns(~M1,inc,col);
      SwapRows(~M1,inc,row);
      if boolN then 
      SwapRows(~NG,inc,row);
      SwapColumns (~ND,inc,col); 
      end if;
      M1ii := DivideSeries(M1[inc,inc], minval);
      inv, n := InvertOverSu(M1ii,pAdicModulo);
      MultiplyColumn(~M1,inv,inc);
      if boolN then MultiplyColumn(~ND,inv,inc); end if;
      //now the pivot is M[inc,inc]
      //cancel all M[inc,i] for i > inc
      for j := inc+1 to Ncols(M1) do
        if not(IsZeroOverS(M1[inc,j])) then
          vij := Valuation(M1[inc,j]);
          if vij lt n then
            MultiplyColumn (~M1, upow^(n-vij),j);
            if boolN then MultiplyColumn(~ND, upow^(n-vij),j); end if;
            vij := n;
          end if;
          c := -upow^(vij-n)*inv*(Parent(M1[1,1])!Eltseq(M1[inc,j]));
          AddColumn (~M1, c, inc, j); 
          if boolN then AddColumn(~ND, c, inc, j); end if;	  
        end if;
      end for;
      //cancel all M[i,inc] for i > inc
      for j := inc+1 to Nrows(M1) do
        if not(IsZeroOverS(M1[j,inc])) then
          vij := Valuation(M1[j,inc]);
          if vij lt n then
            MultiplyRow (~M1, upow^(n-vij),j);
            if boolN then MultiplyRow(~NG, upow^(n-vij),j); end if;
            vij := n;
          end if;
          c := -upow^(vij-n)*inv*(Parent(M1[1,1])!Eltseq(M1[j,inc]));
          AddRow (~M1, c, inc, j); 
          if boolN then AddRow(~NG, c, inc, j); end if;	  
        end if;
      end for;
      inc := inc+1;
      if (inc gt Nrows(M1)) or (inc gt Ncols(M1)) then
	cont_flag := false;
      end if;
    end if;
  end while;

  return M1,NG,ND;
end intrinsic;

//*************************************
// Compute the relations of a Matrix over Su
//*************************************
intrinsic RelationsOverSu(M::Mtrx,pAdicModulo::RngIntElt:nu:=0) -> Mtrx
{
  M1:= RelationsOverSp(M,pAdicModulo:nu:=0) : Mtrx,RingIntElt:FldReElt -> Mtrx

Compute the relations matrix over Su of the matrix M that is given a matrix M compute R of maximal rang such that M.R=0 modulo p^pAdicModulo
}

  local M1,ND1,R,i;

  M1 := M;
  M1,ND1 := RowEchelonOverSu(M1,true,pAdicModulo:nu:=nu);
  R := Matrix(Parent(M1[1,1]), Ncols(M1), 0, []);

  for i := 1 to Ncols(M1) do
    if (pValuationMatOverS(ColumnSubmatrix(M1,i,1)) ge pAdicModulo) then
      R := HorizontalJoin(R,ColumnSubmatrix(ND1,i,1));
    end if;
  end for;

  return R;
end intrinsic;


/*******************************************************************
 *                                                                 *
 *   FREE MODULES                                                  *
 *   Authors: X. Caruso, D. Lubicz                                 *
 *                                                                 *
 *******************************************************************/
/* A module is given as a submodule of S^d (with S = Snu, Sp or Su). 
It is represented by the matrix M of a generating family. When S = Snu,
it can also be given as a couple (Mp,Mu) where Mp (resp. Mu) is a matrix
of a generating family of M[1/pi] (resp. M[pi^b/u^a]^) (NB : M = Mp \cap Mu)
*/

/*******************************************************************
 *                                                                 *
 *   MODULES OVER Sp                                               *
 *                                                                 *
 *******************************************************************/

//*************************************
//Compute the reduced determinant of Sp a module
//in HNF representation
//In the case that the module is full rank and T1,...,Tn 
//are diagonal terms of the HNF
//Returns \prod (T_i/p^(pValuation(T_i))
//*************************************
DetHNFOverSp := function(M,nu)
  local result,i,j,Sp;

  Sp := Parent(M[1,1]);
  result := DivideSeries(M[1,1],pValuationOverS(M[1,1],nu));

  j := 1;
  for i := 2 to Ncols(M) do
    while (j le Nrows(M)) and IsZeroOverS(M[j,i]) do 
      j := j+1;
    end while;
    if j gt Nrows(M) then
      break i;
    else
      result := result*DivideSeries(M[j,i],pValuationOverS(M[j,i],nu));
    end if;
  end for;

  return result;
end function;


//*************************************
//Compute the sum of two free modules over Sp
//*************************************
intrinsic SumOverSp(M1::Mtrx,M2::Mtrx:nu:=0) -> Mtrx
{
M := SumOverSp(M1,M2:nu:=0) : Mtrx,Mtrx:FldReElt -> Mtrx

Returns the matrix M representing the sum of the two free Sp-modules represented by the matrices M1 and M2. 
}
  local S,M,ND,Col;
  S := Parent(M1[1,1]);
  M := HorizontalJoin(M1,M2);
  M := RowEchelonOverSp(M,false:nu:=nu);
  M1 := Matrix(S,Nrows(M),0,[]);
  for i := 1 to Ncols(M) do
    Col := ColumnSubmatrix(M,i,1);
    if not(IsZeroMatOverS(ColumnSubmatrix(M,i,1))) then
      M1 := HorizontalJoin(M1,Col);
    end if;
  end for;
  return M1;
end intrinsic;

//*************************************
//Compute the intersection of two free modules over Sp
//*************************************
intrinsic InterOverSp(M1::Mtrx,M2::Mtrx:nu:=0) -> Mtrx
{
M := InterOverSp(M1,M2:nu:=0) : Mtrx,Mtrx:FldReElt -> Mtrx

Returns the matrix M representing the intersection of the two free Sp-modules represented by the matrices M1 and M2. 
}
  local R;
  R := RelationsOverSp(HorizontalJoin(M1,M2):nu:=nu);
  return M1*ExtractBlock(R,1,1,Ncols(M1),Ncols(R));
end intrinsic;

//*************************************
//Returns a basis of the kernel of a morphism 
//of modules over Sp given by a matrix M
//*************************************
intrinsic KernelOverSp(M::Mtrx:nu:=0) -> Mtrx
{
M := KernelOverSp(M:nu:=0) : Mtrx:FldReElt  -> Mtrx

Returns the matrix M whose column vectors is a basis of the kernel of a morphism of modules over Sp given by the matrix M.
}
  local M1,ND,X,i;

  M1,ND := RowEchelonOverSp(M,true:nu:=nu);
  K := ZeroMatrix(CoefficientRing(M),Ncols(M),0);

  i :=1;
  while (i le Ncols(M1)) do
    if IsZeroMatOverS(ColumnSubmatrix(M1,i,1)) then
      InsertBlock(K, ColumnSubmatrix(ND,i,1), 0, 1);
    end if;
    i := i +1;
  end while;

  return K;
end intrinsic;

//*************************************
//Returns a basis of the image of a morphism 
//of modules over Sp given by a matrix M
//*************************************
intrinsic ImageOverSp(M::Mtrx,boolND::BoolElt:nu:=0) -> Mtrx,Mtrx
{
  M1,ND := ImageOverSp(M,boolND,nu) : Mtrx,BoolElt: -> Mtrx,Mtrx

Returns the matrix M1 whose column vectors is a basis of the image of
  the morphism of modules over Sp given by the matrix M. If
  boolND=true returns the matrix ND such that the matrix M1 is a
  submatrix of M ND.  }

  local M1,i,ND;
  M1,ND := RowEchelonOverSp(M,boolND:nu:=nu);
  i :=1;
  while (i le Ncols(M1)) do
    if IsZeroMatOverS(ColumnSubmatrix(M1,i,1)) then
      RemoveColumn(~M1,i);
      if boolND then RemoveColumn(~ND,i); end if;
    else
      i := i +1;
    end if;
  end while;
  return M1,ND;
end intrinsic;

//*************************************
// Test if Y belongs to the module generated by M
// Returns (false, 0, 0) if Y does not belong to M
// Otherwise, returns
// - (true, X, val) where M*X = p^val*Y if boolX = true
// - (true, 0, 0) if boolX = false
//*************************************

// First, we assume that M is p-row-echelonned.
_IsInOverSp := function (M, Y, boolX,nu)
  local S, pPrecision, p, i, j, inc, val, X, v, q, r,vq;

  S := Parent(M[1,1]);
  p := UniformizingElement (CoefficientRing(S));
  pPrecision := Precision (CoefficientRing(S));
  if boolX then X := ZeroMatrix(S, Ncols(M), 1); else; X := 0; end if;
  inc := 1; val := 0;
  for i := 1 to Ncols(M) do
    if inc gt Nrows(M) or IsZeroOverS (M[i,inc]) then
      if (not IsZeroOverS (Y[i,1])) then
        return false,0,0;
      end if;
    else
      v := pValuationOverS(M[i,inc],nu);
      q,r := ExtendedGcdOverSp (Y[i,1], DivideSeries(M[i,inc], v),nu);
      if IsZeroOverS(r) then
        vq := pValuationOverS(q,nu);
        if v ge vq then
          val := val + v-vq;
          for j := 1 to Nrows(M) do Y[j,1] := Y[j,1] * p^(v-vq); end for;
          if boolX then
            for j := 1 to i-1 do X[j,1] := X[j,1] * p^(v-vq); end for;
          end if;
          v := vq;
        end if;
        if boolX then
          X[i,1] := DivideSeries(q,v);
          Y := Y - X[i,1]*ColumnSubmatrix(M,inc,1);
        else 
          Y := Y - DivideSeries(q,v)*ColumnSubmatrix(M,inc,1);
        end if;
      else
        return false,0,0;
      end if;
      inc := inc + 1;
    end if;
  end for;

  return true,X,val;
end function;

IsInOverSp := function (M, Y, boolX,nu)
  local Mp, ND, bool, X, val;

  Mp,ND := ImageOverSp(M, boolX,nu);
  bool,X,val := _IsInOverSp(Mp, Y, boolX,nu);
  if (not bool) then return false,0,0; end if;
  if boolX then return true, ND*X, val; end if;
  return true, 0, 0; 
end function;

//*************************************
// Test if M1 is included in M2 as modules over Sp
// if boolP returns also (P,v) such that
// M2*P = p^v M1
//*************************************
intrinsic IsIncludedOverSp(M1::Mtrx,M2::Matrx,boolP::BoolElt:nu:=0) -> BoolElt, Mtrx, RngIntElt
{
  Bool,P,v := IsIncludedOverSp(M1,M2,boolP:nu:=0) Mtrx,Mtrx,BollElt:FldReElt -> BoolElt, Mtrx, RngIntElt

Test if the Sp-module represented M1 is included in the Sp-module
represented by M2.  If Bool=true it also returns (P,v) such that M2*P = p^v M1 if it exists.
}
  local M2p,i,ND,v,X,p,P,S;

  S := Parent(M1[1,1]);
  p := UniformizingElement(CoefficientRing(S));
  M2p,ND := ImageOverSp(M2, boolP:nu:=nu);
 
  v := 0; 
  if boolP then
  P := Matrix(S,Ncols(M2p),0,[]);
  else
  P := 0;
  end if;

  for i := 1 to Ncols(M1) do
    bool,X,val := _IsInOverSp (M2p, ColumnSubmatrix (M1, i, 1), boolP,nu);
    if (bool) then
      if (boolP) then
        if (val gt v) then
          P := p^(val-v)*P;
          val := v;
	else 
          X := p^(v-val)*X;
	end if;
        P := HorizontalJoin(P,X);
      end if;
    else
      return false,0,0;
    end if;
  end for;

  return true,ND*P,v;
end intrinsic;

//*************************************
// Test if M1 is equal to M2 as modules over Sp
// if boolP returns also (P,v) such that
// M2*P = p^v M1
//*************************************
intrinsic IsEqualOverSp(M1::Mtrx,M2::Matrx,boolP::BoolElt:nu:=0) -> BoolElt, Mtrx, RngIntElt
{
  Bool,P,v := IsEqualOverSp(M1,M2,boolP:nu:=0) Mtrx,Mtrx,BollElt:FldReElt -> BoolElt, Mtrx, RngIntElt

Test if the Sp-module represented M1 is equal to the Sp-module
represented by M2.  If Bool=true it also returns (P,v) such that M2*P = p^v if it exists.
}

  local mybool,P,v;
  mybool, P,v := IsIncludedOverSp(M1,M2,boolP:nu:=nu);
  mybool := mybool and IsIncludedOverSp(M1,M2,false);
  return mybool, P, v;
end intrinsic;

//*************************************
// Returns the direct sum of M1 and M2
// modules over Sp
//*************************************
intrinsic DirectSumOverSp(M1::Mtrx,M2::Mtrx) -> Mtrx
{
M := DirectSumOverSp(M1,M2) Mtrx,Mtrx -> Mtrx

Compute the module represented by M equal to the direct sum of the modules represented respectively by M1 and M2.
}
  return HorizontalJoin(M1,M2);
end intrinsic;

/*******************************************************************
 *                                                                 *
 *   MODULES OVER Su                                               *
 *                                                                 *
 *******************************************************************/


//*************************************
//Compute the reduced determinant of Su a module
//in HNF representation
//In the case that the module is full rank and T1,...,Tn 
//are diagonal terms of the HNF
//Returns \prod (T_i/u^(uValuation(T_i))
//*************************************
DetHNFOverSu := function(M)
  local result,i,j,Su;

  Su := Parent(M[1,1]);
  result := Su!Eltseq(M[1,1]);

  j := 1;
  for i := 2 to Ncols(M) do
    while (j le Nrows(M)) and IsZeroOverS(M[j,i]) do 
      j := j+1;
    end while;
    if j gt Nrows(M) then
      break i;
    else
      result := result*(Su!Eltseq(M[j,i]));
    end if;
  end for;

  return result;
end function;


//*************************************
//Compute a basis of the kernel of a morphism
//over Su given by the matrix M
//*************************************
intrinsic KernelOverSu(M::Mtrx:nu:=0) -> Mtrx
{
M := KernelOverSu(M:nu:=0) : Mtrx:FldReElt  -> Mtrx

Returns the matrix M whose column vectors is a basis of the kernel of a morphism of modules over Su given by the matrix M.
}

  local M1,ND,X,i;

  M1,ND := RowEchelonOverSu(M,true:nu:=nu);
  K := ZeroMatrix(CoefficientRing(M),Ncols(M),0);

  i :=1;
  while (i le Ncols(M1)) do
    if IsZeroMatOverS(ColumnSubmatrix(M1,i,1)) then
      InsertBlock(K, ColumnSubmatrix(ND,i,1), 0, 1);
    end if;
    i := i +1;
  end while;

  return K;
end intrinsic;

//*************************************
//Compute a basis of the image of 
//a morphism over Su given by the matrix M
//*************************************
intrinsic ImageOverSu(M::Mtrex,boolND::BoolElt:nu:=0) -> Mtrx,Mtrx
{
  M1,ND := ImageOverSp(M,boolND,nu) : Mtrx,BoolElt: -> Mtrx,Mtrx

Returns the matrix M1 whose column vectors is a basis of the image of
  the morphism of modules over Su given by the matrix M. If
  boolND=true returns the matrix ND such that the matrix M1 is a
  submatrix of M ND.  }


  local M1,i,ND;
  M1,ND := RowEchelonOverSu(M,boolND:nu:=nu);
  i :=1;
  while (i le Ncols(M1)) do
    if IsZeroMatOverS(ColumnSubmatrix(M1,i,1)) then
      RemoveColumn(~M1,i);
      if boolND then
        RemoveColumn(~ND,i);
      end if;
    else
      i := i +1;
    end if;
  end while;
  return M1,ND;
end intrinsic;

//*************************************
//Compute the intersection of two free modules over Su
//*************************************
intrinsic InterOverSu(M1::Mtrx,M2::Mtrx:nu:=0) -> Mtrx
{
M := InterOverSu(M1,M2:nu:=0) : Mtrx,Mtrx:FldReElt -> Mtrx

Returns the matrix M representing the intersection of the two free Su-modules represented by the matrices M1 and M2. 
}

  local R;
  R := RelationsOverSu(HorizontalJoin(M1,M2));
  return M1*ExtractBlock(R,1,1,Ncols(M1),Ncols(R));
end intrinsic;

//*************************************
//Compute the sum of two free modules over Su
//*************************************
intrinsic SumOverSu(M1::Mtrx,M2::Mtrx:nu:=0) -> Mtrx
{
M := SumOverSu(M1,M2:nu:=0) : Mtrx,Mtrx:FldReElt -> Mtrx

Returns the matrix M representing the sum of the two free Su-modules represented by the matrices M1 and M2. 
}

  local S,M,ND,Col;
  S := Parent(M1[1,1]);
  M := HorizontalJoin(M1,M2);
  M := RowEchelonOverSu(M,false:nu:=nu);
  M1 := Matrix(S,Nrows(M),0,[]);
  for i := 1 to Ncols(M) do
    Col := ColumnSubmatrix(M,i,1);
    if not(IsZeroMatOverS(ColumnSubmatrix(M,i,1))) then
      M1 := HorizontalJoin(M1,Col);
    end if;
  end for;
  return M1;
end intrinsic;

//*************************************
//Compute the sum of two free modules over Su
//*************************************
intrinsic DirectSumOverSu(M1::Mtrx,M2::Mtrx) -> Mtrx
{
M := DirectSumOverSu(M1,M2) Mtrx,Mtrx -> Mtrx

Compute the module represented by M equal to the direct sum of the modules represented respectively by M1 and M2.
}

  return HorizontalJoin(M1,M2);
end intrinsic;

//*************************************
// Test if Y belongs to the module generated by M[1/u] mod pAdicModulo
// Returns (false, 0, 0) if Y does not belong to M[1/u] mod
// pAdicModulo
// Otherwise, returns
// - (true, X, val) where M*X = u^val*Y if boolX = true
// - (true, 0, 0) if boolX = false
//*************************************

// First, we assume that M is u-row-echelonned.
_IsInOverSu := function(M, Y, boolX,pAdicModulo,nu)
local inc,S,p,X,i,uval,uval1;

  S := Parent(M[1,1]);
  upow := UniformizingElement(Parent(S));
  p := UniformizingElement(CoefficientRing(Parent(S)));

  if boolX then
    X := ZeroMatrix(S, Ncols(M), 1);
  else
    X := 0;
  end if;

  inc := 1; n :=0;
  for i := 1 to Ncols(M) do
    while(IsZero(M[inc,i])) do
      if not(IsZeroOverS(Y[inc,1])) then
	return false,0,0;
      end if;
      inc := inc +1;
    end while;

    if inc gt Nrows(M) then 
      break i;
    end if;
    pval := pValuationOverS(M[inc,i],nu);
    if pval gt pValuationOverS(Y[inc,1],nu) then
      return false, 0,0;
    else
      uval := Valuation(M[inc,i]);
      uval1 := Valuation(Y[inc,1]);
      if uval gt uval1 then
	Y := upow^(uval-uval1);
	n := n +uval-uval1;
	uval1 := uval;
      end if;
      Y := Y - upow^(uval1-uval)*S!Coefficients(DivideSeries(Y[inc,1],pval))*ColumnSubmatrix(M,i,i+1);
      if boolX then
        X[i,1] := upow^(uval1-uval)*S!Coefficients(DivideSeries(Y[inc,1],pval));
      end if; 
    end if;
  end for;

return true,X,n;
end function;

IsInOverSu := function (M, Y, boolX)
  local Mp, ND, bool, X, val;

  Mp,ND := ImageOverSu(M, boolX);
  bool,X,val := _IsInOverSu(Mp, Y, boolX);
  if (not bool) then return false,0,0; end if;
  if boolX then return true, ND*X, val; end if;
  return true, 0, 0; 
end function;

//*************************************
// Test if M1 is included in M2
// if boolP returns also (P,v) such that
// M2*P = u^v M1
//*************************************
intrinsic IsIncludedOverSu(M1::Mtrx,M2::Mtrx,boolP::BoolElt:nu:=0)-> BoolElt, Mtrx, RngIntElt

{
  Bool,P,v := IsIncludedOverSu(M1,M2,boolP:nu:=0) Mtrx,Mtrx,BollElt:FldReElt -> BoolElt, Mtrx, RngIntElt

Test if the Su-module represented M1 is included in the Su-module
represented by M2.  If Bool=true it also returns (P,v) such that M2*P = u^v M1 if it exists.
}

  local M2u,i,ND,v,X,p,P,S;

  S := Parent(M1[1,1]);
  upow := UniformizingElement(S);
  M2u,ND := ImageOverSu(M2, boolP:nu:=nu);
 
  v := 0; 
  if boolP then
  P := Matrix(S,Ncols(M2u),0,[]);
  else
  P := 0;
  end if;

  for i := 1 to Ncols(M1) do
    bool,X,val := _IsInOverSu (M2u, ColumnSubmatrix (M1, i, 1), boolP,nu);
    if (bool) then
      if (boolP) then
        if (val gt v) then
          P := upow^(val-v)*P;
          val := v;
	else 
          X := upow^(v-val)*X;
	end if;
        P := HorizontalJoin(P,X);
      end if;
    else
      return false,0,0;
    end if;
  end for;

  return true,ND*P,v;
end intrinsic;

//*************************************
// Test if M1 is equal to M2
// if boolP returns also (P,v) such that
// M2*P = p^v M1
//*************************************
intrinsic IsEqualOverSu(M1::Mtrx,M2::Mtrx, boolP::BoolElt:nu:=0)-> BoolElt, Mtrx, RngIntElt
{
  Bool,P,v := IsEqualOverSu(M1,M2,boolP:nu:=0) Mtrx,Mtrx,BollElt:FldReElt -> BoolElt, Mtrx, RngIntElt

Test if the Su-module represented M1 is equal to the Su-module
represented by M2.  If Bool=true it also returns (P,v) such that M2*P = u^v if it exists.
}

  local mybool,P,v;
  mybool, P,v := IsIncludedOverSu(M1,M2,boolP:nu:=nu);
  mybool := mybool and IsIncludedOverSu(M1,M2,false);
  return mybool, P, v;
end intrinsic;

/*******************************************************************
 *                                                                 *
 *   MODULES OVER S                                                *
 *                                                                 *
 *******************************************************************/

//*************************************
// Compute the free module associated to a module
// given par a matrix M
// The relations of M is given by R
//*************************************
FreeWithROverS := function(M,R)
  local S,pPrecision,p,i,j,k,inc,P;
  local u,v,up,vp;

  if (Ncols(R) eq 0) then return M; end if;

  S := Parent(M[1,1]);
  pPrecision := Precision(CoefficientRing(S));
  p := UniformizingElement(CoefficientRing(S));

  P := IdentityMatrix(S, Ncols(M));
  inc := 1;
  for i := 1 to Ncols(R) do
    //check if the column contains non zero elements
    val,j := pColValuationOverS(R,i,inc,Nrows(R));
    if not(IsZeroMatOverS(ExtractBlockRange(R, inc, i, Nrows(R), i+1))) then
      // first we prepare the row so that it contains a distinguished element 
      SwapRows (~R, inc, j); SwapColumns (~M, inc, j); SwapRows(~P, inc,j); 
      // cancel all R[j,i] for j>inc
      while true do
        vij,j := pColValuationOverS(R,i,inc+1,Nrows(R));
        if  IsZeroMatOverS(ExtractBlockRange(R, inc+1, i, Nrows(R), i+1))then break; end if;
        vij := vij - val;
        if vij gt 0 then
          // Perform OR1 operation
          for k := 1 to Nrows(R) do
            if k ne inc then R[k,i] := R[k,i]/p^vij; end if;
          end for;
          for k := 1 to Ncols(R) do
            if k ne i then R[inc,k] := R[inc,k]*p^vij; end if;
          end for;
          MultiplyColumn(~M,1/p^vij,inc);
	  MultiplyRow(~P, p^vij,inc); 
        end if;
	
        _,_,u,v,up,vp := ExtendedGcdWithValuationOverS(R[i,inc],R[i,j],val);
        LinearCombinaisonRows (~R, inc, j, u, v, up, vp);
        LinearCombinaisonColumns (~M, inc, j, vp, -up, -v, u);
	LinearCombinaisonRows(~P, inc,j,u,v,up,vp);
      end while;
      inc := inc+1;
    end if;
  end for;

  //extract the free part of M
  inc := 1;
  for i := 1 to Nrows(R) do
    if IsZeroMatOverS(ExtractBlockRange(R,i,1,i+1,Ncols(R))) and (not(IsZeroMatOverS(ExtractBlockRange(M,1,i,Nrows(M),i+1)))) then
      inc := inc+1;
    else
      RemoveColumn (~M, inc);
      RemoveRow (~P, inc);
    end if;
  end for;

  return M,P;
end function;

intrinsic FreeOverS(M::Mtrx) -> Mtrx,Mtrx
{
  M,P := FreeOverS(M1) Mtrx -> Mtrx,Mtrx 

Compute a matrix M representing the free module associated to a module
given par a matrix M1 as well as P such that M is a submatrix of M1*P.
}
  return FreeWithROverS(M,RelationsOverSp(M));
end intrinsic;

//*************************************
//Compute a basis over S of an S module
//in Su Sp representation
//*************************************
intrinsic BasisOverS(MS::Mtrx:nu:=0) -> Mtrx
{
MF := BasisOverS(MS:nu:=0) Mtrx: FldReElt -> Mtrx

Returns a matrix giving a basis over S of an S module in Su Sp representation.
}

  local M,MF,MSu,MSp;

  MSu := MS[1];
  MSp := MS[2];

  D1 := DetHNFOverSp(MSp,nu);
  D2 := DetHNFOverSu(MSu); 
  MSu := D1*MSu;
  MSp := D2*MSp;

  M := HorizontalJoin(MSu,MSp);
  R := RelationsOverSp(M:nu:=nu);

  MF := FreeWithROverS(M,R:nu:=nu);

  return MF;
end intrinsic;
