import random

#################-- General --###########################################################

## This code implements an algorithm taken from the paper "An elementary proof that
## rationally isometric quaratic forms are isometric". Its imput are two hermitian forms
## over a valuation ring, and an isometry between them which is defined over the fraction
## field (or division ring), and its output an isometry which is defined over the
## valuation ring. Both hermitian forms and isometries are represented as matrices.

## The file features four classes:
##   1. QQ -- implementing rational numbers,
##   2. QQi -- implementing elements of the field Q[i] with i=sqrt(-1),
##   3. Quat -- implementing elements of a quaternion algebra (-1,-3) over Q (the
##      values -1 and -3 can be changed during running time or manually).
##   4. rational_iso_algo -- implements the algorithm above, and several related methods
##      such as generating random input for the algorithm.

## Valuation rings in QQ, QQi and Quat are not implemented. To "caputure" a valuation
## ring in these claases, one uses an additive valuation. Several such valuations are
## implemented below. 

## For convenience, within the file are also included the following particular instances
## of rational_iso_alo:
##   1. rational_iso_algo_QQ(p) -- the algorithm class over QQ, with the p-adic valuation
##   2. rational_iso_algo_QQi(p) -- the alogrithm class over QQi, with the complex
##      conjugation involution, and the p-adic valuation (p is a prime number congruent
##      3 mod 4)
##   3. rational_iso_algo_Quat3 -- the algorithm class over Quat, with the *orthogonal*
##      involution (a+bi+cj+dij)-->(a+bi+cj-dij), and the unique extension of the
##      3-adic valuation. (In general one can take any prime for which Quat does not
##      ramify.)

## Here is a sample code:

def test_rational_iso_algo_instance(alg):
    A,B,U = alg.random_input(4,[-2,-1,1,2]) # worst case complexity
    print ("A = ")
    alg.mat_print(A)
    print ("B = ")
    alg.mat_print(B)
    print ("U = ")
    alg.mat_print(U)
    print ("UAU* == B? ...", alg.test_isometry(A,B,U))
    print ("valuation(U) = ")
    alg.mat_print([[alg.valuation(x) for x in r] for r in U])
    print ("Finding isometry V over valuation ring...")
    V = alg.find_isometry(A,B,U)
    print ("Done.")
    print ("V = ")
    alg.mat_print(V)
    print ("VAV* == B? ...", alg.test_isometry(A,B,V))
    print ("valaution(V) = ")
    alg.mat_print([[alg.valuation(x) for x in r] for r in V])

def test_rational_iso_algo():
    print ("----- Testing rational_iso_algo_QQ(3) -----")
    alg = rational_iso_algo_QQ(3)
    test_rational_iso_algo_instance(alg)
    print()

    print ("----- Testing rational_iso_algo_QQi(3) -----")
    alg = rational_iso_algo_QQi(3)
    test_rational_iso_algo_instance(alg)
    print()

    print ("----- Testing rational_iso_algo_Quat3 -----")
    alg = rational_iso_algo_Quat3
    test_rational_iso_algo_instance(alg)
    print()

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

#################-- QQ class: implements rational numbers --#############################

class QQ(object):
    """
    Implements the rational numbers.
    Disclaimer: I know that there are built in implementations.
    """
    def __init__(self, nom=0, denom=1):
        if denom==0:
            raise "Zero denominator!"
        if denom < 0:
            denom = -denom
            nom = -nom
        self.nom = nom
        self.denom = denom
        self._red()

    def _red(self):
        "Reduces nominator and denominator by the greatest common divisor."
        a = self.nom
        b = self.denom
        while a%b != 0:
            a, b = b, a%b
        self.nom = self.nom//b
        self.denom = self.denom//b

    def __str__(self):
        return str(self.nom)+"/"+str(self.denom)

    def __repr__(self):
        return str(self)

    def __float__(self):
        return float(self.nom)/float(self.denom)

    def __add__(self,other):
        return QQ(self.nom*other.denom+self.denom*other.nom,self.denom*other.denom)

    def __sub__(self,other):
        return QQ(self.nom*other.denom-self.denom*other.nom,self.denom*other.denom)

    def __mul__(self,other):
        return QQ(self.nom*other.nom,self.denom*other.denom)

    def __truediv__(self,other):
        return QQ(self.nom*other.denom,self.denom*other.nom)

    def __neg__(self):
        return QQ(-self.nom,self.denom)

    def __ne__(self,other):
        return (self.nom!=other.nom)or(self.denom!=other.denom)

    def __eq__(self,other):
        return (self.nom==other.nom)and(self.denom==other.denom)

    def inverse(self):
        return QQ(self.denom,self.nom)


def QQ_padic_val(x,p):
    "Returns the p-adic valuation of x."
    i=0
    nom = x.nom
    denom = x.denom
    if nom==0:
        return 100000000 # infinity
    while nom%p==0:
        i+=1
        nom//=p
    while denom%p==0:
        i-=1
        denom//=p
    return i

def QQ_element_factory_for_padic_val(v,p):
    "Returns a QQ instance with p-adic valuation v."
    if v>=0:
        return QQ(p**v)
    else:
        return QQ(1,p**(-v))

def rand_QQ(_range=10):
    return QQ(random.randrange(-_range,_range+1),random.randrange(1,_range+1))


###################-- QQi class: implement the field Q[i] --####################

class QQi(object):
    """
    Implements the field Q[i] where Q is the rationals and i=sqrt(-1).
    Disclaimer: I know that there are built in implementations.
    """
    def __init__(self, real=0, im=0):
        if type(real) is QQ:
            self.real = real
        else:
            self.real = QQ(real)
        if type(im) is QQ:
            self.im = im
        else:
            self.im = QQ(im)

    def __str__(self):
        return "("+str(self.real)+")+("+str(self.im)+")i"

    def __repr__(self):
        return str(self)

    def __add__(self,other):
        return QQi(self.real+other.real, self.im+other.im)

    def __sub__(self,other):
        return QQi(self.real-other.real, self.im-other.im)

    def __mul__(self,other):
        return QQi(self.real*other.real - self.im*other.im, self.real*other.im + self.im*other.real)

    def __truediv__(self,other):
        return self * other.inverse()

    def __neg__(self):
        return QQi(-self.real,-self.im)

    def __ne__(self,other):
        return (self.real!=other.real)or(self.im!=other.im)

    def __eq__(self,other):
        return (self.real==other.real)and(self.im==other.im)

    def inverse(self):
        a = self.real*self.real+self.im*self.im
        return QQi(self.real/a,-self.im/a)

    def conjugate(self):
        return QQi(self.real,-self.im)


def QQi_padic_val(x,p):
    "Returns the p-adic valuation of x. p must be a prime congruent to 3 mod 4"
    return min(QQ_padic_val(x.real,p),QQ_padic_val(x.im,p))

def QQi_element_factory_for_padic_val(v,p):
    "Returns a QQi instance with p-adic valuation v. p must be a prime congruent to 3 mod 4"
    if v>=0:
        return QQi(p**v)
    else:
        return QQi(QQ(1,p**(-v)))

def rand_QQi(_range=2):
    return QQi(rand_QQ(_range),rand_QQ(_range))


#################-- Quat: Implementing the a quaternion algebra over QQ --################

# the algebra implemented here is (-1,-3). One can change ALL instances at once by
# setting the values of Quat.i_squared and Quat.j_squared

class Quat(object):
    """
    Implements the quaternion algebra (i_squared,j_squared) over the rational numbers.
    This the 4-dimensional algebra spanned over Q by 1,i,j,ij and subject to the relations:
    i*i = i_squared (default: -1)
    j*j = j_squared (default: -3)
    i*j = -j*i
    The values of i_squared and j_squared can be changed, causing all instances of Quat
    to "switch" to the new quaternion algebra.
    """
    
    i_squared=QQ(-1) # the square of i
    j_squared=QQ(-3) # the square of j
    
    def __init__(self, a=0,b=0,c=0,d=0):
        "represents a+b*i+c*j+d*ij"
        if type(a) is QQ:
            self.a=a
        else:
            self.a = QQ(a)
        if type(b) is QQ:
            self.b=b
        else:
            self.b = QQ(b)
        if type(c) is QQ:
            self.c=c
        else:
            self.c = QQ(c)
        if type(d) is QQ:
            self.d=d
        else:
            self.d = QQ(d)

    def __str__(self):
        return "("+str(self.a)+")+("+str(self.b)+")i+("+str(self.c)+")j+("+str(self.d)+")ij"

    def __repr__(self):
        return str(self)

    def __add__(self,other):
        return Quat(self.a+other.a,self.b+other.b,self.c+other.c,self.d+other.d)

    def __sub__(self,other):
        return Quat(self.a-other.a,self.b-other.b,self.c-other.c,self.d-other.d)

    def __mul__(self,other):
        return Quat(self.a*other.a + self.b*other.b*Quat.i_squared + \
                    self.c*other.c*Quat.j_squared - self.d*other.d*Quat.i_squared*Quat.j_squared, \
                    self.a*other.b + self.b*other.a - self.c*other.d*Quat.j_squared + self.d*other.c*Quat.j_squared, \
                    self.a*other.c + self.c*other.a + self.b*other.d*Quat.i_squared - self.d*other.b*Quat.i_squared, \
                    self.a*other.d + self.d*other.a + self.b*other.c - self.c*other.b)

    def __neg__(self):
        return Quat(-self.a,-self.b,-self.c,-self.d)

    def __ne__(self,other):
        return (self.a!=other.a)or(self.b!=other.b)or(self.c!=other.c)or(self.d!=other.d)

    def __eq__(self,other):
        return (self.a==other.a)and(self.b==other.b)and(self.c==other.c)and(self.d==other.d)

    def norm(self):
        return self.a*self.a - self.b*self.b*Quat.i_squared - self.c*self.c*Quat.j_squared + \
               self.d*self.d*Quat.i_squared*Quat.j_squared

    def inverse(self):
        N = self.norm()
        return Quat(self.a/N,-self.b/N,-self.c/N,-self.d/N)

    def conjugate(self):
        "Returns a-b*i-c*j-d*ij."
        return Quat(self.a,-self.b,-self.c,-self.d)

    def orth_inv(self):
        "Returns a+b*i+c*j-d*ij."
        return Quat(self.a,self.b,self.c,-self.d)

    def __pow__(self,n):
        if n==0:
            return Quat(1)
        if n>=0:
            tmp = self
            res = Quat(1)
            while n>0:
                if n%2==1:
                    res = res*tmp
                tmp = tmp*tmp
                n = n//2
            return res
        else:
            return (self.__pow__(-n)).inverse()
                
def Quat_padic_val(x,p):
    """
    Returns the p-adic valuation of x.
    Here, the p-adic valuation stands for the unique extension of the padic valuation over Q to the
    Quaternions. It exists if and only if the quaternion algebra (Quat.i_squared,Quat.j_squared) ramifies at p,
    in which case it is given by the p-adic valuation of the norm.
    """
    return QQ_padic_val(x.norm(),p)

def Quat_element_factory_for_padic_val(v,element_of_norm_p):
    """
    Returns an instance of Quat with p-adic valuation v.
    element_of_norm_p must be a quaternion elemenet of norm p, and the quaterion algebra must ramify at p.
    To produce elements which are invariant some involution, element_of_norm_p must be invariant under the
    that involution.
    """
    return element_of_norm_p**v

def rand_Quat(_range=2):
    return Quat(rand_QQ(_range),rand_QQ(_range),rand_QQ(_range),rand_QQ(_range))



##################--- Rational Isomorphism algorithm ---########################################

class rational_iso_algo(object):
    """
    The class implement an algorithm derived from the paper "An elementary proof that
    rationally isometric quaratic forms are isometric". Its imput are two hermitian forms
    over a valuation ring, and an isometry between them which is defined over the fraction
    field (or division ring), and its output an isometry which is defined over the
    valuation ring. Both hermitian forms and isometries are represented as matrices
    (i.e. lists of lists).
    
    If the hermitian forms are represented by matrices A and B, then an isometry is a matrix
    U satisying UAU*=B where U* stands for the "star involution of U". (The star involution
    is applying the field involution to all entries of the matrix and then taking transpose.)

    The main methods of the class are:
        1. find_isometry(A,B,U) --- implements the algorithm
        2. test_isometry(A,B,U) --- tests whether U is an isometry from A to B
        3. random_input() --- returns random valid input A,B,U for the algorithm
        4. mat_print() --- print a matrix
        5. KAK_decomp(U) --- returns the KAK-decomposition of a matrix U

    The class's members include:
        1. field --- the class implementing the base field (or division ring)
        2. valuation --- an additive valuation on the field (which defines a valuation ring
           in the field)
        3. involution --- an involution of the field.
        4. sym_element_factory --- a function taking a valuation and returning an element
           of the field with given valuation and which is invariant under the involution.
           This is needed only if the involution is not the identity map.
        5. random_element_factory --- produces random elements of the field.
           This is needed only for random_input(), the algorithm implemented by the class
           is deterministic.
    The precise requirements from the memebers of the class appear in the documentation
    of the constructor.
    """

    def __init__(self, \
                 field=QQ, \
                 valuation=None, \
                 involution=None, \
                 sym_element_factory=None, \
                 random_element_factory=None
                 ):
        """
        Arguments:
            1. field (madatory, default: QQ) --- a class implementing a field (or division ring).
            2. valaution (mandatory) --- an additive valauation on "field"
            3. involution (default: None, meaning identity) --- an involution on "field"
            4. sym_element_factory (optional if involution is identity) --- a function
               taking an argument in the range of "valuation" and returning an element of
               "field" which is invaraint under "involution" and whose valaution is the argument given. 
            5. random_element_factory (optional) --- a function returning a random instance of "field"

        Requirements:
            [1] The class "field" must support the following methods:
                arithmetic operation (+,-,*),
                equting elements (==,!=),
                constructor taking an integer,
                a function inverse() returning the inverse of an object.
            [2] The input and return value of "involution" is a single object of type "field".
            [3] The input of "valuation" is a single object of type "field".
                The type returned by "valuation" must support the following methods:
                arithmetic operations (+,-),
                comparison of elements (<,>,==,!=,<=,>=)
            [4] The input of "sym_element_factory" is an object of the type returned by
                "valuation" and its output is an object of type "field".
            [5] "rand_element_factory" has no input and it returns an instance of "field".
        """
        self.field = field
        self.valuation = valuation

        if involution is None:
            self.involution = lambda x:x
        else:
            self.involution = involution

        self.sym_element_factory = sym_element_factory
        self.random_element_factory = random_element_factory

        self._one = field(1)
        self._zero = field(0)
        self._zero_val = self.valuation(field(1))

    

    ##################--- Operations on matrices ---###########################################

    def _mat_dims(self,A):
        "Returns the dimensions of a matrix."
        if len(A)==0:
            return 0,0
        return len(A), len(A[0])

    def _mat_trans(self,A):
        "Returns the transpose of a matrix."
        n,m = self._mat_dims(A)
        return [[A[i][j] for i in range(n)] for j in range(m)]

    def _mat_star(self,A):
        """
        Returns the star of the matrix, namely, the transpose of the matrix after "involution"
        is applied to all entries.
        """
        n,m = self._mat_dims(A)
        return [[self.involution(A[i][j]) for i in range(n)] for j in range(m)]

    def _mat_sub(self,A,B):
        if self._mat_dims(A)!=self._mat_dims(B):
            raise "Bad dimensions"
        return [[x-y for x,y in zip(r,s)] for r,s in zip(A,B)]

    def _mat_neg(self,A):
        return [[-x for x in r] for r in A]

    def _mat_mul(self,A,B):
        if A == []:
            return []
        nA,mA = self._mat_dims(A)
        nB,mB = self._mat_dims(B)
        if mA != nB:
            raise "Bad dimensions!"
        Bt = self._mat_trans(B)
        return [[sum([x*y for (x,y) in zip(rA,rBt)],self._zero) for rBt in Bt] for rA in A]

    def _mat_lscalar_mul(self,A,a):
        return [[a*x for x in r] for r in A]

    def _mat_rscalar_mul(self,A,a):
        return [[x*a for x in r] for r in A]

    def _mat_id(self,n):
        "Returns an identitiy matrix of size n."
        A=[[self._zero for i in range(n)] for j in range(n)]
        for i in range(n):
            A[i][i]=self._one
        return A

    def _mat_inv(self,A):
        "Returns the inverse of A, or None if A is not invertbile."
        n,m = self._mat_dims(A)
        if n != m:
            return None
        one = self._one
        zero = self._zero
        E = [(r + (n*[zero])) for r in A]
        for i in range(n):
            E[i][n+i]=one
        for i in range(n):
            for j in range(i,n):
                if E[j][i] != zero:
                    break
            if E[j][i]==zero:
                return None
            self._row_mul(E,j,E[j][i].inverse())
            for k in range(n):
                if k != j:
                    self._row_op(E,k,j,E[k][i])
            if i != j:
                E[i], E[j] = E[j], E[i]
        return [r[n:] for r in E]

    def _mat_conj(self,U,A):
        "Returns UAU*."
        return self._mat_mul(self._mat_mul(U,A),self._mat_star(U))

    def _mat_clone(self,A):
        "Returns a copy of A."
        return [[x for x in r] for r in A]

    ##################--- Row and column operations ---###########################################

    def _row_op(self,A,i,j,a):
        A[i] = [x-a*y for x,y in zip(A[i],A[j])]

    def _col_op(self,A,i,j,a):
        for r in A:
            r[i] = r[i]-r[j]*a

    def _row_swap(self,A,i,j):
        A[i],A[j] = A[j],A[i]

    def _col_swap(self,A,i,j):
        for r in A:
            r[i],r[j] = r[j],r[i]
            
    def _row_mul(self,A,i,a):
        A[i] = [a*x for x in A[i]]

    def _col_mul(self,A,i,a):
        for r in A:
            r[i] = r[i]*a

    ##################--- Random objects ---###########################################

    def _rand_R_element(self):
        "Returns a random element of R, the valuation ring. Based on random_element_factory."
        a = self.random_element_factory()
        if self.valuation(a) < self._zero_val:
            a = a.inverse()
        return a

    def _rand_R_nonzero_element(self):
        "Returns a random element of R, the valuation ring. Based on random_element_factory."
        a = self.random_element_factory()
        if self.valuation(a) < self._zero_val:
            a = a.inverse()
        if a == self._zero:
            a = self._one
        return a

    def _rand_R_inv_element(self):
        "Returns a random invertible element of R, the valuation ring."
        a = self.random_element_factory()
        if a == self._zero:
            a = self._one
        return a.inverse()*self.sym_element_factory(self.valuation(a))

    def _rand_R_mat(self,n,m):
        "Returns a random matrix over R, the valuation ring."
        return [[self._rand_R_element() for i in range(m)] for j in range(n)]

    def _rand_R_sym_mat(self,n):
        "Returns a random star-symmetric matrix of R, the valuation ring."
        A = self._rand_R_mat(n,n)
        return [[A[i][j]+self.involution(A[j][i]) for i in range(n)] for j in range(n)]

    def _rand_R_inv_mat(self,n):
        "Returns a random invertible matrix of R, the valuation ring."
        # Disclaimer: I know that there are other ways to do this ...
        
        flag = True
        while flag:
            A = self._rand_R_mat(n,n)
            A_inv = self._mat_inv(A)
            flag = (A_inv is None)
            if not flag:
                for x in sum(A_inv,[]):
                    if self.valuation(x) < self._zero_val:
                        flag = True
                        break
        return A

    def _rand_R_inv_sym_mat(self,n):
        "Returns a random star-symmetric invertible matrix of R, the valuation ring."
        # Disclaimer: I know that there are better ways to do this ...
        
        flag = True
        while flag:
            A = self._rand_R_sym_mat(n)
            A_inv = self._mat_inv(A)
            flag = (A_inv is None)
            if not flag:
                for x in sum(A_inv,[]):
                    if self.valuation(x) < self._zero_val:
                        flag = True
                        break
        return A


    ####################--- Public Methods ---##############################################

    def mat_print(self,A):
        "Prints a matrix A."
        for r in A:
            print(r)

    def random_input(self,n,KAK_digonal_valuations=None):
        """
        Generates a random input of size n x n for the algorithm. Namely, returns three matrices
        A,B,U such that A,B are symmetric and invertible over R, the valuation ring,
        and UAU*=B (U may not be defined over R).
        The optional variable "KAK_digonal_valuations" specifies the valuations of the diagonal
        etries of the "A part" of the KAK decomposition of the isometry constructed. This list
        must satisfy the condition that for every entry v, the frequency of v and -v
        in the list is the same. If "KAK_digonal_valuations" consists of distinct values,
        then the algorithm reaches its worst case complexity. On the other hand, a list of
        zeroes will cause U to be defined over R, thus making the algorithm's complexity minimal.
        """
        d_vals = []
        
        if KAK_digonal_valuations is None:
            L = sorted([self.valuation(self._rand_R_nonzero_element()) for i in range(n//2)])
            if n%2==0:
                d_vals = [-x for x in L]+L
            else:
                d_vals = [-x for x in L[::-1]]+[self._zero_val]+L
        else:
            d_vals = sorted(KAK_digonal_valuations)
            if d_vals != [-x for x in d_vals[::-1]]:
                raise "Impossible diagonal valuations!"

        # producing blocks sizes
        block_sizes = []

        i=0
        while i<n:
            j=i
            while d_vals[i]==d_vals[j]:
                i+=1
                if i==n:
                    break
            block_sizes.append((i-j,d_vals[j]))

        # producing blocks
        t = len(block_sizes)
        blocks = [[None for i in range(t)] for j in range(t)]
        for i in range(t):
            s,v = block_sizes[i]

            if v == self._zero_val:
                blocks[i][i] = self._rand_R_inv_sym_mat(s)
            else:
                blocks[i][i] = self._rand_R_sym_mat(s)
                if v < self._zero_val:
                    a = self.sym_element_factory(-v)
                    blocks[i][i] = self._mat_lscalar_mul(blocks[i][i],a)
                    blocks[i][i] = self._mat_rscalar_mul(blocks[i][i],a)

            for j in range(i+1,t):
                s1,v1 = block_sizes[j]
                if v+v1 == self._zero_val:
                    blocks[i][j] = self._rand_R_inv_mat(s)
                else: 
                    blocks[i][j] = self._rand_R_mat(s,s1)
                    if v+v1 < self._zero_val:
                        blocks[i][j] = self._mat_lscalar_mul(blocks[i][j],self.sym_element_factory(-v-v1))
                blocks[j][i] = self._mat_star(blocks[i][j])

        # glueing blocks
        
        A = []
        for i in range(t):
            s,v = block_sizes[i]
            r = blocks[i]
            for j in range(s):
                A.append(sum([x[j] for x in r],[]))

        # generating the isometry
        U = [[self._zero for i in range(n)] for j in range(n)]

        for i in range(n):
            U[i][i] = self.sym_element_factory(d_vals[i])
            
        # Adding "noise"
        P = self._rand_R_inv_mat(n)
        Q = self._rand_R_inv_mat(n)

        A = self._mat_conj(self._mat_inv(P),A)
        U = self._mat_mul(U,P)
        U = self._mat_mul(Q,U)

        B = self._mat_conj(U,A)

        return A,B,U

    def KAK_decomp(self,A):
        """
        Returns the KAK decomposition of a matrix A over "field".
        Namely, the function returns three matrices K1, D, K2 such that
        K1 and K2 are defined and invertbile over R, the valuation ring,
        D is digonal (by not necessarily defined over R), and A = K1 x D x K2.
        """
        n,m=self._mat_dims(A)
        if n!=m:
            raise "Bad dimensions!"
        K1 = self._mat_id(n)
        K2 = self._mat_id(n)
        A = self._mat_clone(A)

        for k in range(n,1,-1):
            
            i1=0
            j1=0
            v1 = self.valuation(A[0][0])
            for i in range(k):
                for j in range(k):
                    v = self.valuation(A[i][j])
                    if v < v1:
                        i1,j1,v1=i,j,v
            
            self._row_swap(A,k-1,i1)
            self._col_swap(K1,k-1,i1)
            self._col_swap(A,k-1,j1)
            self._row_swap(K2,k-1,j1)
            a=A[k-1][k-1].inverse()

            for i in range(k-1):
                b = A[i][k-1]*a
                self._row_op(A,i,k-1,b)
                self._col_op(K1,k-1,i,-b)
                b = a*A[k-1][i]
                self._col_op(A,i,k-1,b)
                self._row_op(K2,k-1,i,-b)

        return K1,A,K2

    def test_isometry(self,A,B,U):
        """
        Returns True if UAU* == B
        """
        return self._mat_conj(U,A) == B

    def find_isometry(self,A,B,U):
        """
        Implementing the algorithm of the class.
        
        Let R denote the valuation ring of "field".
        The input consists of two matrices A,B over R, which are star-symmetric and invertible
        over R, and a matrix U over "field" such that UAU*==B.
        The function returns a matrix V defined over R and satisfying VAV*=B
        """

        if not self.test_isometry(A,B,U):
            raise "Forms not isometric!"

        n,m = self._mat_dims(U)

        # Throughout the algorithm, we will conjugate A and B. The transitions prerformed
        # are recorded in the variables P and Q. Namely, PAP* and Q^{-1}BQ^{-1}* will give us the
        # original value of A and B.
        
        K1,U,K2 = self.KAK_decomp(U)
        P = K2
        A = self._mat_conj(K2,A)
        Q = K1
        B = self._mat_conj(self._mat_inv(K1),B)

        # d_val are the valuation diagonal etries of U (which is now diagonal)
        # v is the maximal valution in d_val
        
        d_val = [self.valuation(U[i][i]) for i in range(n)]
        v = max(d_val)
        
        while v > self._zero_val: 

            # rearranging U,A,B to make d_val be of the form [v,v,v,..., -v,-v,-v,..., <all other valuations>]
            # s records the number of times v (and -v) occur in d_vals
            # t = 2*s (if no exception is raised)
            
            s=0
            for i in range(n):
                
                if d_val[i] == v:
                    if i != s:
                        d_val[s],d_val[i] = d_val[i],d_val[s]
                        U[s][s],U[i][i] = U[i][i],U[s][s]
                        self._row_swap(A,i,s)
                        self._col_swap(A,i,s)
                        self._row_swap(B,i,s)
                        self._col_swap(B,i,s)
                        self._row_swap(P,i,s)
                        self._col_swap(Q,i,s)
                    s+=1
                
            t=s
            for i in range(s,n):
                if d_val[i] == -v:
                    if i != t:
                        d_val[t],d_val[i] = d_val[i],d_val[t]
                        U[t][t],U[i][i] = U[i][i],U[t][t]
                        self._row_swap(A,i,t)
                        self._col_swap(A,i,t)
                        self._row_swap(B,i,t)
                        self._col_swap(B,i,t)
                        self._row_swap(P,i,t)
                        self._col_swap(Q,i,t)
                    t+=1

            if t!=s+s:
                raise "One A,B is not invertibe over the valuation ring!"

            # v2 is the second largest valuation in d_vals
            # p is an element of "field" with valaution v-v2
            
            v2 = max(d_val[t:]+[self._zero_val])
            for i in range(t,n):
                if d_val[i]==v2:
                    break
            p = U[0][0]*U[i][i].inverse()
            if (p != self.involution(p)) or not (self.sym_element_factory is None):
                p = self.sym_element_factory(v-v2)
            p_inv = p.inverse()

            # for the code to follow, it is good to keep in mind the following block matrices
            #     [A11      A12         A13   ]  
            # A = [... (p x A22 x p) (p x A23)]  
            #     [...      ...         A33   ] 

            A11 = [r[:s] for r in A[:s]]
            A12 = [r[s:t] for r in A[:s]]
            A13 = [r[t:] for r in A[:s]]
            A22 = self._mat_lscalar_mul(self._mat_rscalar_mul([r[s:t] for r in A[s:t]],p_inv),p_inv)
            A23 = self._mat_lscalar_mul([r[t:] for r in A[s:t]],p_inv)
            A33 = [r[t:] for r in A[t:]]
            A33_inv = self._mat_inv(A33)

            # We now construct three matrice P1,P2,P3 with which we shall conjugate A.
            # The effect of this conjugation will be equivalent to having A conjugated by
            # diag(p,p,p,..., p_inv,p_inv,p_inv,...., 1,1,1,...). The outer loop will then repeat itself.
            # See the paper "An elementary proof that rationally isometric quadratic forms are isometric"
            # for further explantions about the computations to follow.

            # for brevity:
            _mul2 = lambda x,y: self._mat_mul(x,y)
            _mul3 = lambda x,y,z: self._mat_mul(self._mat_mul(x,y),z)
            _sub  = lambda x,y: self._mat_sub(x,y)
            _star = lambda x: self._mat_star(x)
            _rsm  = lambda x,a: self._mat_rscalar_mul(x,a)
            _lsm  = lambda x,a: self._mat_lscalar_mul(x,a)

            X = A11
            Y = A12
            Z = A22

            if t!=n:
                X = _sub(A11, _mul3(A13,A33_inv,_star(A13)))
                Y = _sub(A12, _mul3(A13,A33_inv,_rsm(_star(A23),p)))
                Z = _sub(A22, _mul3(A23,A33_inv,_star(A23)))

            Y_inv = self._mat_inv(Y)

            V13 = self._mat_neg(_mul3(Y_inv,A13,A33_inv))
            V23 = _mul2(_lsm(A23,-p),A33_inv)

            W13 = _mul2(_lsm(A13,p),A33_inv)
            W23 = _mul2(A23,A33_inv)
            W11 = _lsm(_rsm(Y,p_inv),p)

            Id1 = self._mat_id(s)
            Id2 = self._mat_id(n-t)
            z = [self._zero for i in range(s)]

            P1 = [Y_inv[i] + z      + V13[i] for i in range(s)] + \
                 [z        + Id1[i] + V23[i] for i in range(s)] + \
                 [z        + z      + Id2[i] for i in range(n-t)]

            P3 = [W11[i] + z      + W13[i] for i in range(s)] + \
                 [z      + Id1[i] + W23[i] for i in range(s)] + \
                 [z      + z      + Id2[i] for i in range(n-t)]

            X = _mul3(Y_inv,X,_star(Y_inv))
            pId = self._mat_lscalar_mul(Id1,p)
            pId_inv = self._mat_lscalar_mul(Id1,p_inv)

            # W = 2(2-pX-pZ)^{-1}
            W = self._mat_inv(_sub(_sub(Id1, \
                                        _lsm(X,p*self.field(2).inverse())), \
                                   _lsm(Z,p*self.field(2).inverse())))
            
            # U11 = p - (1-pX)Wp = p - (1-pX)2(2-pX-pZ)^{-1}p
            U11 = _sub(pId,_rsm(_mul2(_sub(Id1,_lsm(X,p)),W),p))

            # U12 = (1-pX)W = (1-pX)2(2-pX-pZ)^{-1} 
            U12 = _mul2(_sub(Id1, _lsm(X,p)), W)

            # U21 = (p^{-1} - Z)Wp = (p^{-1} - Z)2(2-pX-pZ)^{-1}p
            U21 = _rsm(_mul2(_sub(pId_inv,Z),W),p)

            # U22 = p^{-1} - (p^{-1}-Z)W = p^{-1} - (p^{-1}-Z)2(2-pX-pZ)^{-1}
            U22 = _sub(pId_inv,_mul2(_sub(pId_inv,Z),W))
            
            P2 = [U11[i] + U12[i] + [self._zero]*(n-t) for i in range(s)] + \
                 [U21[i] + U22[i] + [self._zero]*(n-t) for i in range(s)] + \
                 [z + z + Id2[i] for i in range(n-t)]


            # Updating P

            P = _mul2(P1,P)
            P = _mul2(P2,P)
            P = _mul2(P3,P)

            # Updating A, U, and d_val

            u = self.valuation(p)
            v = v-u

            for i in range(s):
                self._row_mul(A,i,p)
                self._row_mul(A,i+s,p_inv)
                self._col_mul(A,i,p)
                self._col_mul(A,i+s,p_inv)
                U[i][i] = U[i][i]*p_inv
                U[i+s][i+s] = U[i+s][i+s]*p
                d_val[i] = d_val[i] - u
                d_val[i+s] = d_val[i+s] + u

                                
        if v < self._zero_val: 
            raise "One A,B is not invertibe over the valuation ring or not star-symmetric!" 

        return self._mat_mul(self._mat_mul(Q,U),P)


################-- Particular instances of rational_iso_algo --#########################

def small_rand_QQ():
    return [QQ(0),QQ(1),QQ(-1)][random.randrange(3)]

def small_rand_QQi():
    return [QQi(0),QQi(0),QQi(0), \
            QQi(1),QQi(-1),QQi(0,1),QQi(0,-1)][random.randrange(7)]

def small_rand_Quat():
    return [Quat(0),Quat(0),Quat(0), \
            Quat(0),Quat(0),Quat(0),Quat(0), \
            Quat(1),Quat(-1),Quat(0,1),Quat(0,-1), \
            Quat(0,0,1),Quat(0,0,-1),Quat(0,0,0,1),Quat(0,0,0,-1)][random.randrange(15)]

def rational_iso_algo_QQ(p):
    """
    The rational_iso_algo class over QQ, with the p-adic valuation.
    p must be a prime number.
    """
    return rational_iso_algo(QQ, \
                             lambda x:QQ_padic_val(x,p), \
                             None, \
                             lambda v:QQ_element_factory_for_padic_val(v,p), \
                             small_rand_QQ)                                          

def rational_iso_algo_QQi(p):
    """
    The rational_iso_algo class over QQi, with the with the complex conjugation involution
    and the p-dic valuation.
    p must be a prime number congruent to 3 mod 4.
    """
    return rational_iso_algo(QQi, \
                             lambda x:QQi_padic_val(x,p), \
                             lambda x:x.conjugate(), \
                             lambda v:QQi_element_factory_for_padic_val(v,p), \
                             small_rand_QQi)                                          


##  rational_iso_algo_Quat3 implement rational_iso_also over Quat, with the *orthogonal*
##  involution (a+bi+cj+dij)-->(a+bi+cj-dij), and the unique extension of the
##  3-adic valuation. The variables Quat.i_squared and Quat.j_squared must be set to
##  their default values, -1 and -3 (resp.), before using this instance.

rational_iso_algo_Quat3 = rational_iso_algo(Quat, \
                                            lambda x:Quat_padic_val(x,3), \
                                            lambda x:x.orth_inv(), \
                                            lambda v:Quat_element_factory_for_padic_val(v,Quat(0,0,1,0)), \
                                            small_rand_Quat)