Other deprecated files (i.e., never used) were deleted.
Files adapted: bareiss.py, matrix.py, ring_w_sequence.py, verbose.py
Files deleted: euclidean.py, hermite.py, nullspace_wrapper.py,
storj_villard.py, timing.py
-# This file was *autogenerated* from the file ./bareiss.sage
-from sage.all_cmdline import * # import sage library
+from sage.all import * # import sage library
-_sage_const_2 = Integer(2); _sage_const_1 = Integer(1); _sage_const_0 = Integer(0)
from ajpastor.misc.matrix import swap_rows;
from ajpastor.misc.matrix import swap_cols;
This class receives several parameters:
- `parent`: the polynomial ring where the operations will be done. A fraction field is allowed if the base domain is a polynomial ring.
- `M`: the matrix we want to treat with the Algorithm
- - `method`: a Boolean function
+ - `method`: a Boolean function such that I = {g in parent | method(g) == True}
This algorithm also works in a quotient of the polynomial ring (K[X]/I) and the method provided must be a method to check if some polynomial is in the ideal I.
'''
## Checking the matrix input
self.base_matrix = Matrix(self.parent(), M);
- self.change_of_columns = Permutation(range(_sage_const_1 ,self.base_matrix.ncols()));
+ self.change_of_columns = Permutation(range(1 ,self.base_matrix.ncols()));
## Storing the checking method
self.__in_ideal = method;
#################################################
#@wverbose
def echelon_form(self):
+ '''
+ Method to compute the right echelon form of the matrix self.base_matrix.
+
+ This method store the result for future callings in a private variable.
+ '''
if(self.__echelon is None):
self.__compute_echelon();
return self.__echelon;
- @wverbose
+ #@wverbose
def __compute_echelon(self):
- sverbose("Computing the echelon form of a matrix of size %dx%d" %(self.base_matrix.nrows(), self.base_matrix.ncols()));
+ '''
+ Method to really compute the right echelon form of self.base_matrix. This method
+ is private and it is not recommended to be called by the user.
+ '''
+ #sverbose("Computing the echelon form of a matrix of size %dx%d" %(self.base_matrix.nrows(), self.base_matrix.ncols()));
self.__actions = [("base")];
self.__steps = [self.base_matrix];
## If we have initial relations, we perform the first reduction
if(self.__have_ideal(self.__gb)):
- sverbose("Initial reduction of the matrix"); ## Verbose
- self.__steps += [self.__simplify_matrix(self.__steps[-_sage_const_1 ])];
+ #sverbose("Initial reduction of the matrix"); ## Verbose
+ self.__steps += [self.__simplify_matrix(self.__steps[-1 ])];
self.__actions += [("f_reduce", self.__gb)];
- tr = self.__steps[-_sage_const_1 ].nrows(); tc = self.__steps[-_sage_const_1 ].ncols();
- cr = _sage_const_0 ; cc = _sage_const_0 ; i = -_sage_const_1 ;
+ tr = self.__steps[-1 ].nrows(); tc = self.__steps[-1 ].ncols();
+ cr = 0 ; cc = 0 ; i = -1 ;
- sverbose("Starting the iterations");
- sverbose.start_iteration(min(tr,tc)+_sage_const_1 , True, False);
+ #sverbose("Starting the iterations");
+ #sverbose.start_iteration(min(tr,tc)+1 , True, False);
while(i < min(tr,tc)):
- i = i + _sage_const_1 ;
+ i = i + 1 ;
- sverbose.increase_depth();
+ #sverbose.increase_depth();
- pivot = self.__choose_pivot(self.__steps[-_sage_const_1 ], i,i);
+ pivot = self.__choose_pivot(self.__steps[-1 ], i,i);
- cr,cc = pivot[_sage_const_0 ];
- new_rels = pivot[_sage_const_2 ];
+ cr,cc = pivot[0 ];
+ new_rels = pivot[2 ];
## If there are new relations, we simplify the matrix
- if(len(new_rels) > _sage_const_0 ):
- sverbose("New relations found looking for a pivot");
- sverbose("\t%s" %new_rels);
+ if(len(new_rels) > 0 ):
+ #sverbose("New relations found looking for a pivot");
+ #sverbose("\t%s" %new_rels);
self.__gb = ideal(self.parent(), tuple(new_rels) + tuple(self.__gb)).groebner_basis();
- self.__steps += [self.__simplify_matrix(self.__steps[-_sage_const_1 ])];
+ self.__steps += [self.__simplify_matrix(self.__steps[-1 ])];
self.__actions += [("reduce", new_rels)];
- sverbose("New reductions applied");
+ #sverbose("New reductions applied");
## If no pivot was found, we finish
- if(cr == -_sage_const_1 or cc == -_sage_const_1 ):
- sverbose("No pivot found.");
- sverbose.decrease_depth();
- sverbose.next_iteration();
+ if(cr == -1 or cc == -1 ):
+ #sverbose("No pivot found.");
+ #sverbose.decrease_depth();
+ #sverbose.next_iteration();
break;
## If we have a new pivot, we moved everythin so it is in the proper position
- swap_actions = pivot[_sage_const_1 ];
+ swap_actions = pivot[1 ];
for action in swap_actions:
- if(action[_sage_const_0 ] == "sw_r"):
- self.__steps += [swap_rows(self.__steps[-_sage_const_1 ], action[_sage_const_1 ],action[_sage_const_2 ])];
- elif(action[_sage_const_0 ] == "sw_c"):
- self.__steps += [swap_cols(self.__steps[-_sage_const_1 ], action[_sage_const_1 ],action[_sage_const_2 ])];
- self.change_of_columns = Permutation((action[_sage_const_1 ]+_sage_const_1 ,action[_sage_const_2 ]+_sage_const_1 ))*self.change_of_columns;
+ if(action[0 ] == "sw_r"):
+ self.__steps += [swap_rows(self.__steps[-1 ], action[1 ],action[2 ])];
+ elif(action[0 ] == "sw_c"):
+ self.__steps += [swap_cols(self.__steps[-1 ], action[1 ],action[2 ])];
+ self.change_of_columns = Permutation((action[1 ]+1 ,action[2 ]+1 ))*self.change_of_columns;
self.__actions += [action];
## One the pivot is in position, we proceed with the Bareiss elimination
- M = self.__steps[-_sage_const_1 ];
+ M = self.__steps[-1 ];
## We save the new matrix and go on
- self.__steps += [self.__bareiss(self.__steps[-_sage_const_1 ], i)];
+ self.__steps += [self.__bareiss(self.__steps[-1 ], i)];
self.__actions += [("bareiss", i)];
- sverbose.decrease_depth();
- sverbose.next_iteration();
+ #sverbose.decrease_depth();
+ #sverbose.next_iteration();
- sverbose("Performing gcd simplification");
- gcds = [gcd(row) for row in self.__steps[-_sage_const_1 ]];
+ #sverbose("Performing gcd simplification");
+ gcds = [gcd(row) for row in self.__steps[-1 ]];
for i in range(len(gcds)):
- if(gcds[i] == _sage_const_0 ):
- gcds[i] = _sage_const_1 ;
- self.__steps += [Matrix(self.parent(), [[el/gcds[i] for el in self.__steps[-_sage_const_1 ][i]] for i in range(self.__steps[-_sage_const_1 ].nrows())])];
+ if(gcds[i] == 0 ):
+ gcds[i] = 1 ;
+ self.__steps += [Matrix(self.parent(), [[el/gcds[i] for el in self.__steps[-1 ][i]] for i in range(self.__steps[-1 ].nrows())])];
self.__actions += [("gcd_simpl")];
- self.__echelon = self.__steps[-_sage_const_1 ];
- sverbose("Finished computation of echelon form");
+ self.__echelon = self.__steps[-1 ];
+ #sverbose("Finished computation of echelon form");
return;
@derived_property
def rank(self):
+ '''
+ Method to compute the property "rank" of self.base_matrix. This method
+ uses the right echelon form of the matrix and check how many
+ rows are all zeros.
+ '''
for i in range(self.base_matrix.nrows()):
- if(any((el != _sage_const_0 ) for el in self.echelon_form()[-(i+_sage_const_1 )])):
+ if(any((el != 0 ) for el in self.echelon_form()[-(i+1 )])):
return self.base_matrix.nrows()-i;
def relations(self):
+ '''
+ Method to get the realtions between the variables found in the matrix.
+
+ A relation is a polynomial p in self.parent() such that self.is_in_ideal(p) is True.
+ We return a Groebner basis of all the relations found while computing the echelon
+ form of the matrix.
+
+ WARNING: It is important to remark that this method does not return the basis
+ of the ideal that self.is_in_ideal() defines. This is because only found relations
+ are taking into account at this step.
+ '''
self.echelon_form();
return self.__gb;
@cached_method
def right_kernel_matrix(self):
+ '''
+ Method to compute the right kernel matrix of self.base_matrix.
+
+ This method computes the echelon form of the matrix and then starts a
+ simple computation to clean denominator and obtain a basis formed with
+ elements of the integral domain self.parent().
+ '''
sol_dimension = self.base_matrix.ncols()-self.rank;
M = self.echelon_form();
- if(sol_dimension > _sage_const_0 ):
+ if(sol_dimension > 0 ):
## Compute the product of the nice diagonal
A = self.__get_lcm([M[i][i] for i in range(self.rank)]);
to_mult = [-A/M[i][i] for i in range(self.rank)];
- ker_basis = [vector(self.parent(), [to_mult[j]*M[j][i+self.rank] for j in range(self.rank)] + [_sage_const_0 for j in range(i)] + [A] + [_sage_const_0 for j in range(i+_sage_const_1 , sol_dimension)]) for i in range(sol_dimension)];
+ ker_basis = [vector(self.parent(), [to_mult[j]*M[j][i+self.rank] for j in range(self.rank)] + [0 for j in range(i)] + [A] + [0 for j in range(i+1 , sol_dimension)]) for i in range(sol_dimension)];
ch = self.change_of_columns;
## If there were a change of columns (i.e. ch is not the trivial permutation)
## we swap the result
- if(Permutation(range(_sage_const_1 ,self.base_matrix.ncols())) != ch):
+ if(Permutation(range(1 ,self.base_matrix.ncols())) != ch):
ch = ch.inverse();
- rows = [[_sage_const_0 for i in range(M.ncols())] for j in range(len(ker_basis))];
+ rows = [[0 for i in range(M.ncols())] for j in range(len(ker_basis))];
for j in range(M.ncols()):
- new_col = ch(j+_sage_const_1 )-_sage_const_1 ;
+ new_col = ch(j+1 )-1 ;
for i in range(len(ker_basis)):
rows[i][new_col] = ker_basis[i][j];
return ker_basis;
else:
- return [vector(self.parent(), [_sage_const_0 for i in range(M.ncols())])];
+ return [vector(self.parent(), [0 for i in range(M.ncols())])];
#################################################
### Other cached methods
#################################################
@cached_method
def is_in_ideal(self, p):
+ '''
+ Method that defines the ideal over this algorithm is working on.
+
+ This method is a boolean method (returns True or False) in such a way that
+ I = {g in self.parent() : self.is_in_ideal(g) == True} is an ideal
+ of self.parent().
+ '''
p = self.parent()(p);
try:
#################################################
### Private methods for Bareiss Algorithm
#################################################
- @wverbose
+ #@wverbose
def __choose_pivot(self, M, ir, ic):
'''
This method computes the next pivot element for the algorithm and returns the information to prepare the matrix for the next step.
The ir and ic are parameters to begin the search from the position (ir,ic)
+
+ A valid pivot is any element on the matrix with position strictly greater than (ir,ic)
+ (i.e., both coordinates must be greater or equal to those) and must not be in the
+ ideal defined by self.is_in_ideal().
+
+ INPUT:
+ - M: the current matrix we are choosing the pivot.
+ - ir: initial index for the rows.
+ - ic: initial index fot the columns.
+ OUTPUT:
+ This method returns a tuple of 4 elements:
+ - A pair (fr,fc) for the position of the chosen pivot.
+ (-1,-1) means no pivot was found.
+ - A tuple of actions to transform M in order to put
+ the pivot in the appropriate position.
+ - A set of new relations found during this pivot chosing.
'''
- sverbose("Looking for a new pivot from position (%d,%d)" %(ir,ic));
+ #sverbose("Looking for a new pivot from position (%d,%d)" %(ir,ic));
relations = set();
actions = [];
end = False;
- fc = -_sage_const_1 ; fr = -_sage_const_1 ;
+ fc = -1 ; fr = -1 ;
## Checking rows
for cc in range(ic, M.ncols()):
## Checking columns
for cr in range(ir, M.nrows()):
- sverbose("Checking if position (%d,%d) is zero (%s)" %(cr,cc, M[cr][cc]));
+ #sverbose("Checking if position (%d,%d) is zero (%s)" %(cr,cc, M[cr][cc]));
to_check = M[cr][cc];
- if(len(relations) > _sage_const_0 ):
- sverbose("Reducing the selected element with the found relations");
+ if(len(relations) > 0 ):
+ #sverbose("Reducing the selected element with the found relations");
to_check = M[cr][cc].reduce(relations);
- sverbose("New element to check: %s" %(to_check));
+ #sverbose("New element to check: %s" %(to_check));
if(not to_check == self.parent().zero()):
if(not self.is_in_ideal(to_check)):
- sverbose("Non-zero element found. Pivot found");
+ #sverbose("Non-zero element found. Pivot found");
end = True;
fr = cr;
break;
else:
- sverbose("Non trivial zero found. New relation: %s" %(to_check));
+ #sverbose("Non trivial zero found. New relation: %s" %(to_check));
relations.add(to_check);
if(end):
fc = cc;
break;
- if(fc != -_sage_const_1 and fc != ic):
+ if(fc != -1 and fc != ic):
actions += [("sw_c", ic, cc)];
- for i in range(_sage_const_1 ,(min(cc-ic, M.ncols()-cc))):
+ for i in range(1 ,(min(cc-ic, M.ncols()-cc))):
actions += [("sw_c", ic+i, M.ncols()-i)];
- if(fr != -_sage_const_1 and fr != ir):
+ if(fr != -1 and fr != ir):
actions += [("sw_r", ir, cr)];
- sverbose("Finished search of pivot: (%s,%s)" %(fr,fc));
+ #sverbose("Finished search of pivot: (%s,%s)" %(fr,fc));
return ((fr,fc), tuple(actions), relations);
def __bareiss(self, M, i):
+ '''
+ Method that applies a Bareiss reduction to the matrix M on the position (i,i).
+ This method assumes that M[i][i] is a valid pivot for this reduction.
+
+ The output matrix will preserve the structure of M for positions (j,k) with
+ j,k < i and will create zeros on the i-th column.
+
+ A final reduction with gcds for each column is performed to keep the degrees
+ of the polynomials as small as possible.
+
+ These computations are not in-place, so the resulting matrix needs more
+ memory allocation to be computed.
+
+ INPUT:
+ - M: matrix where we want perform the Bareiss reduction.
+ - i: position where we have a pivot.
+ '''
rows = [];
+ ## Reduction of the rows before i
for j in range(i):
- if(M[j][i] != _sage_const_0 ):
+ if(M[j][i] != 0 ):
rows += [[M[j][k]*M[i][i]-M[j][i]*M[i][k] for k in range(M.ncols())]];
else:
rows += [[M[j][k] for k in range(M.ncols())]];
+
+ ## The i-th row remains the same as before
rows += [[M[i][k] for k in range(M.ncols())]];
- for j in range(i+_sage_const_1 , M.nrows()):
- if(M[j][i] != _sage_const_0 ):
+
+ ## Reduction of the rows after i
+ for j in range(i+1 , M.nrows()):
+ if(M[j][i] != 0 ):
rows += [[M[j][k]*M[i][i]-M[j][i]*M[i][k] for k in range(M.ncols())]];
else:
rows += [[M[j][k] for k in range(M.ncols())]];
+ ## GCD simplification of the rows
try:
gcds = [gcd(row) for row in rows];
rows = [[el/gcd[i] for el in rows[i]] for i in range(len(rows))];
return Matrix(self.parent(), rows);
def __have_ideal(self, basis):
- if(len(basis) == _sage_const_1 ):
- return basis[_sage_const_0 ] != _sage_const_0 ;
- return len(basis) > _sage_const_0 ;
+ '''
+ Auxiliar method to know if some relation have been already found.
+
+ This method returns True if the current Groebner basis we have computed
+ gives some non-trivial ideal.
+ '''
+ if(len(basis) == 1 ):
+ return basis[0 ] != 0 ;
+ return len(basis) > 0 ;
def __simplify_matrix(self, M):
+ '''
+ Method to simplify a polynomial matrix using the relations found.
+
+ INPUT:
+ - M: the matrix to simplify.
+ OUTPUT:
+ A matrix N such that N[i][j] is the reduction via the current
+ Groebner basis (self.relations()) of M[i][j].
+ '''
rows = [[el for el in row] for row in M];
if(self.__have_ideal(self.__gb)):
rows = [[el.reduce(self.__gb) for el in row] for row in rows];
return Matrix(self.parent(), rows);
def __get_lcm(self, input):
+ '''
+ Auxiliar method for computing the lcm of a sequence.
+ '''
try:
return lcm(input);
except AttributeError:
try:
## Relying on gcd
p = self.__parent;
- res = p(_sage_const_1 );
+ res = p(1 );
for el in input:
res = p((res*el)/gcd(res,el));
return res;
+++ /dev/null
-
-# This file was *autogenerated* from the file ./euclidean.sage
-from sage.all_cmdline import * # import sage library
-
-_sage_const_1 = Integer(1); _sage_const_0 = Integer(0)
-
-####################################################################################################
-####################################################################################################
-###
-### Euclidean module
-###
-### ------------------------------------------------------------------------------------------------
-###
-### This file contains the code for some euclidean operations over for polynomials.
-###
-### ------------------------------------------------------------------------------------------------
-###
-### Version: 0.0
-### Date of begining: 29-03-2017
-####################################################################################################
-####################################################################################################
-
-## Basic polynomial -- Private functions
-def __coefficients(PolR, p, var):
- p = PolR(p);
- if(not var in PolR.gens()):
- raise ValueError("The variable given (%s) is not in the ring %s"%(var, PolR));
-
- var = PolR(var);
-
- try:
- return [p.coefficient({var:i}) for i in range(p.degree(var)+_sage_const_1 )];
- except TypeError:
- return p.coefficients(sparse = False);
-
-def __lc(PolR, p, var):
- return __coefficients(PolR, p, var)[-_sage_const_1 ];
-
-def __degree(PolR, p, var):
- p = PolR(p);
- if(not var in PolR.gens()):
- raise ValueError("The variable given (%s) is not in the ring %s"%(var, PolR));
-
- var = PolR(var);
- try:
- return p.degree(var);
- except TypeError:
- return p.degree();
-
-## Polynomial Division
-def PolyDivide(PolR, A,B):
- A = PolR(A); B = PolR(B);
- x = PolR.gens()[-_sage_const_1 ];
-
- Q = PolR.zero(); R = A;
- d = __degree(PolR, R, x) - __degree(PolR, B, x);
- while(R != PolR.zero() and d >= _sage_const_0 ):
- T = (__lc(PolR, R, x)/__lc(PolR, B, x))*x**d;
- Q += T;
- R -= B*T;
- d = __degree(PolR, R, x) - __degree(PolR, B, x);
-
- return (Q,R);
-
-def PolyPseudoDivide(PolR, A, B):
- A = PolR(A); B = PolR(B);
- x = PolR.gens()[-_sage_const_1 ];
-
- b = __lc(PolR, B, x); N = __degree(PolR, A, x) - __degree(PolR, B, x)+_sage_const_1 ;
- Q = PolR.zero(); R = A;
- d = __degree(PolR, R, x) - __degree(PolR, B, x);
- while(R != PolR.zero() and d >= _sage_const_0 ):
- T = __lc(PolR, R, x)*x**d;
- N -= _sage_const_1 ;
- Q = b*Q + T;
- R = b*R - T*B;
- d = __degree(PolR, R, x) - __degree(PolR, B, x);
-
- return (b**N*Q, b**N*R);
-
-## Euclidean Algorithms
-def Euclidean(PolR, a, b):
- a = PolR(a); b = PolR(b);
- while(b != PolR.zero()):
- (q,r) = PolyDivide(PolR,a,b);
- a = b;
- b = r;
- return a;
-
-def ExtendedEuclidean_slow(PolR, a, b):
- a = PolR(a); b = PolR(b);
- a1 = PolR.one(); a2 = PolR.zero(); b1 = PolR.zero(); b2 = PolR.one();
-
- while(b != PolR.zero()):
- (q,r) = PolyDivide(PolR,a,b);
- a = b; b = r;
- r1 = a1-q*b1; r2 = a2-q*b2;
- a1 = b1; a2 = b2; b1 = r1; b2 = r2;
-
- return (a1,a2,a);
-
-def HalfExtendedEuclidean(PolR, a, b):
- a = PolR(a); b = PolR(b);
- a1 = PolR.one(); b1 = PolR.zero();
-
- while(b != PolR.zero()):
- (q,r) = PolyDivide(PolR,a,b);
- a = b; b = r;
- r1 = a1-q*b1;
- a1 = b1; b1 = r1;
-
- return (a1,a);
-
-def ExtendedEuclidean(PolR, a, b):
- a = PolR(a); b = PolR(b);
- (s,g) = HalfExtendedEuclidean(PolR, a,b);
- (t,r) = PolyDivide(PolR, g-s*a,b);
-
- return (s,t,g);
-
-def ExtendedEuclidean_diophantine_slow(PolR, a,b,c):
- a = PolR(a); b = PolR(b); c = PolR(c);
- x = PolR.gens()[-_sage_const_1 ];
-
- (s,t,g) = ExtendedEuclidean(PolR, a,b);
- (q,r) = PolyDivide(PolR, c,g);
- if(r != PolR.zero()):
- raise ValueError("(%s) is not in the ideal generated by (%s) and (%s)" %(c,a,b));
- s = q*s; t = t*q;
-
- if(s != PolR.zero() and (__degree(PolR, s, x) >= __degree(PolR, b, x))):
- (q,r) = PolyDivide(PolR, s,b);
- s = r;
- t = t+q*a;
- return (s,t);
-
-def HalfExtendedEuclidean_diophantine(PolR, a,b,c):
- a = PolR(a); b = PolR(b); c = PolR(c);
- x = PolR.gens()[-_sage_const_1 ];
-
- (s,g) = HalfExtendedEuclidean(PolR, a,b);
- (q,r) = PolyDivide(PolR, c,g);
- if(r != PolR.zero()):
- raise ValueError("(%s) is not in the ideal generated by (%s) and (%s)" %(c,a,b));
- s = q*s;
-
- if(s != PolR.zero() and (__degree(PolR, s, x) >= __degree(PolR, b, x))):
- (q,r) = PolyDivide(PolR, s,b);
- s = r;
- return s;
-
-def ExtendedEuclidean_diophantine(PolR, a,b,c):
- a = PolR(a); b = PolR(b); c = PolR(c);
- s = HalfExtendedEuclidean_diophantine(PolR, a,b,c);
- (t,r) = PolyDivide(PolR, c-s*a, b);
- return (s,t);
-
-## Partial Fractions
-def PartialFraction_simple(PolR, a, d):
- a = PolR(a); d = [PolR(el) for el in d];
- if(len(d) == _sage_const_0 ):
- return [a];
-
- d1 = prod(d); d2 = d1/d[_sage_const_0 ];
-
- (q,r) = PolyDivide(PolR, a,d1);
- if(len(d) == _sage_const_1 ):
- return [q,r];
-
- (a1,t) = ExtendedEuclidean_diophantine(PolR, d2, d[_sage_const_0 ],r);
- aux = PartialFraction_simple(PolR, t,d[_sage_const_1 :])
-
- return [aux[_sage_const_0 ]+q, a1] + aux[_sage_const_1 :];
-
-def PartialFraction(PolR, a, d, e = None):
- a = PolR(a); d = [PolR(el) for el in d];
- if(e is None):
- e = [_sage_const_1 for i in range(len(d))];
- else:
- for i in range(len(e)):
- e[i] = int(e[i]);
- if(e[i] < _sage_const_0 ):
- raise ValueError("No negative exponents allowed in the PartialFraction method. The %d-th exponent is negative." %i);
-
- if(len(d) != len(e)):
- raise TypeError("PartialFraction need the same amount of exponents as denominators");
-
- dens = [d[i]**e[i] for i in range(len(d))];
- pf = PartialFraction_simple(PolR, a, dens);
-
- res = [pf[_sage_const_0 ]];
- for i in range(len(e)):
- part_res = [];
- for j in range(e[i],_sage_const_0 ,-_sage_const_1 ):
- (q,aux) = PolyDivide(PolR, pf[i+_sage_const_1 ], d[i]);
- pf[i+_sage_const_1 ] = q;
- part_res += [aux];
- part_res.reverse();
- res[_sage_const_0 ] = res[_sage_const_0 ] + pf[i+_sage_const_1 ];
- res += part_res;
-
- return res;
-
-## Primitive part and content of polynomials
-def content(PolR, p):
- p = PolR(p);
- x = PolR.gens()[-_sage_const_1 ];
-
- coeffs = __coefficients(PolR, p, x);
- return gcd(coeffs);
-
-def pp_only(PolR, p):
- p = PolR(p);
- return p/content(PolR, p);
-
-def pp(PolR, p):
- c = content(PolR, p);
- return (c, p/c);
-
-## Squarefree factorization
-def Squarefree_Musser(PolR, A):
- A = PolR(A);
- x = PolR.gens()[-_sage_const_1 ];
-
- (c,S) = pp(PolR, A);
-
- S_minus = gcd(S, S.derivative(x));
- S_star = PolR(S/S_minus);
- res = []
- while(__degree(PolR, S_minus, x) > _sage_const_0 ):
- Y = gcd(S_star, S_minus);
- res += [PolR(S_star/Y)];
- S_star = Y;
- S_minus = PolR(S_minus/Y);
- res += [S_star];
-
- res[_sage_const_0 ] = res[_sage_const_0 ]*(c*S_minus);
- return res;
-
-def Squarefree_Yun(PolR, A):
- A = PolR(A);
- x = PolR.gens()[-_sage_const_1 ];
-
- (c,S) = pp(PolR, A);
- S_prime = S.derivative(x);
- S_minus = gcd(S,S_prime);
- S_star = PolR(S/S_minus);
- Y = PolR(S_prime/S_minus);
- res = [];
- Z = Y - S_star.derivative(x);
- while(Z != _sage_const_0 ):
- res += [gcd(S_star, Z)];
- S_star = PolR(S_star/res[-_sage_const_1 ]);
- Y = PolR(Z/res[-_sage_const_1 ]);
- Z = Y - S_star.derivative(x);
-
- res += [S_star];
- res[_sage_const_0 ] = res[_sage_const_0 ]*c;
-
- return res;
-
-
+++ /dev/null
-
-# This file was *autogenerated* from the file ./hermite.sage
-from sage.all_cmdline import * # import sage library
-
-_sage_const_1 = Integer(1); _sage_const_0 = Integer(0)
-from ajpastor.misc.euclidean import *;
-from ajpastor.misc.verbose import *;
-from ajpastor.misc.matrix import swap_cols as swap;
-
-####################################################################################
-###
-### Two Stages for Hermite Normal Form
-###
-####################################################################################
-@wverbose
-def stage1(A,i):
- parent = A.parent().base();
- sverbose("Initiating the loop...");
- sverbose.start_iteration(i,False,False);
- sverbose.increase_depth();
- for j in range(i):
- (p,q,r) = ExtendedEuclidean(parent, A[j][j],A[j][i]);
-
- sverbose("Computed the EEA: (%s)(%s) + (%s)(%s) = %s" %(p,A[j][j],q,A[j][i],r));
-
- new_rows = [[el for el in row] for row in A];
- p2 = parent(-A[j][i]/r); q2 = parent(A[j][j]/r);
-
- sverbose("Column operation matrix:\n%s" %(Matrix(parent,[[p,p2],[q,q2]])));
- for k in range(A.nrows()):
- new_rows[k][j] = p*A[k][j] + q*A[k][i];
- new_rows[k][i] = p2*A[k][j] + q2*A[k][i];
- coeff = new_rows[j][j].leading_coefficient();
- coeff2 = new_rows[j+_sage_const_1 ][i].leading_coefficient();
- for k in range(A.nrows()):
- new_rows[k][j] = new_rows[k][j]/coeff;
- new_rows[k][i] = new_rows[k][i]/coeff2;
- A = Matrix(parent,new_rows);
- sverbose("New matrix computed:\n%s" %A);
- sverbose.next_iteration();
- sverbose.decrease_depth();
- sverbose("Finished the loop");
- return A;
-
-@wverbose
-def stage2(Ao, A, i):
- parent = Ao.parent().base();
- m = i+_sage_const_1 ;
- S = Matrix(parent,[[Ao[i][j] for j in range(m)] for i in range(m)]);
- D = Matrix(parent,[([_sage_const_0 for j in range(l)] + [A[l][l]] + [_sage_const_0 for j in range(l+_sage_const_1 ,m)]) for l in range(m)]);
- K = S * BackslashOperator() * D;
- new_rows = [[el for el in row] + [_sage_const_0 for j in range(m,Ao.ncols())] for row in K];
- for j in range(m,Ao.ncols()):
- new_rows += [[_sage_const_0 for k in range(j)] + [_sage_const_1 ] + [_sage_const_0 for k in range(j+_sage_const_1 ,Ao.ncols())]];
- K = Matrix(K.parent().base(), new_rows);
- return Ao*K;
-
-####################################################################################
-###
-### Rearrangement methods
-###
-####################################################################################
-@wverbose_w_data
-def index_n_from_m(n,m,subset):
- '''
- Let S(n,m) be the subset of the parts of {1,...,m} with n elements. It is well known that |S(n,m)| = binomial(m,n). Let A(n,m) = {1,2,...,binomial(m,m)}. Then there are bijections between S(n,m) and A(n,m).
-
- Consider f a bijection that orders the elements of S(n,m) in a lexicographic order. This method implements such map.
-
- INPUT:
- - n: size fo the subsets we are considering
- - m: number of elements we can take elements
- - subset: an element of S(n,m).
- '''
- if(not (isinstance(subset,list) or not isinstance(subset,set))):
- raise TypeError("The argument subset must be a set or a list of elements");
-
- subset = list(set(subset)); subset.sort();
-
- if(len(subset) != n):
- raise TypeError("The argument subset must have %d elements, but it has %d" %(n,len(subset)));
-
- for el in subset:
- if((not (isinstance(el,Integer) or type(el)==int)) or (el <= _sage_const_0 ) or (el > m)):
- raise ValueError("The elements of subset must be integers from 1 to %d" %m);
-
- if(n == _sage_const_0 ):
- return _sage_const_0 ;
- elif(n == _sage_const_1 ):
- return subset[_sage_const_0 ];
-
- return sum(binomial(m-j,n-_sage_const_1 ) for j in range(_sage_const_1 ,subset[_sage_const_0 ])) + index_n_from_m(n-_sage_const_1 ,m-subset[_sage_const_0 ],[subset[i]-subset[_sage_const_0 ] for i in range(_sage_const_1 ,len(subset))]);
-
-
-@wverbose_w_data
-def subset_n_from_m(n,m,i):
- '''
- Let S(n,m) be the subset of the parts of {1,...,m} with n elements. It is well known that |S(n,m)| = binomial(m,n). Let A(n,m) = {1,2,...,binomial(m,m)}. Then there are bijections between S(n,m) and A(n,m).
-
- Consider f a bijection that orders the elements of S(n,m) in a lexicographic order (see method index_n_from_m). This method inverts such method.
-
- INPUT:
- - n: size fo the subsets we are considering
- - m: number of elements we can take elements
- - i: index of the subset we want to rebuild.
- '''
- if(n < _sage_const_1 or m < n):
- raise TypeError("Imposible compute this: not valid configuration");
- if(n == _sage_const_1 ):
- sverbose("Base case (n = 1). Returning [%d]" %i);
- return [i];
- bin = binomial(m-_sage_const_1 ,n-_sage_const_1 )
- current_sum = _sage_const_0 ;
- index = _sage_const_1 ;
- while(current_sum + bin < i):
- index += _sage_const_1 ;
- current_sum += bin;
- bin = binomial(m-index,n-_sage_const_1 );
-
- sverbose("Index found: %d" %index);
-
- output = [index] + [el+index for el in subset_n_from_m(n-_sage_const_1 ,m-index,i-current_sum)];
- sverbose("Output: %s" %output);
-
- return [index] + [el+index for el in subset_n_from_m(n-_sage_const_1 ,m-index,i-current_sum)];
-
-@wverbose_w_data
-def permutation_from_m(subset,m):
- '''
- Given a subset of {1,2,...,m} on the form [a_1<...<a_n], this method computes a permutation p in S_m such that for all i = 1,...,n
- - p(a_i) = i.
- '''
- perm = [el for el in subset];
- for j in range(_sage_const_1 ,m+_sage_const_1 ):
- if(not j in perm):
- perm += [j];
- return Permutation(perm).inverse();
-
-@wverbose_w_data
-def get_rearrange_from_minor(nrows,ncols,size,i):
- '''
- Given the size of a matrix, a size of a minor that is not zero and the index of such minor we compute the permutations that move rows and columns in order to have that minor in the left top corner of the matrix.
- '''
- bin = binomial(ncols,size);
- for_rows = subset_n_from_m(size,nrows,ceil(i/bin));
-
- to_cols = i%bin;
- if(to_cols == _sage_const_0 ):
- to_cols = bin;
- for_cols = subset_n_from_m(size,ncols,to_cols);
-
- return (permutation_from_m(for_rows,nrows),permutation_from_m(for_cols,ncols));
-
-@wverbose
-def get_rearrange(M):
- '''
- Given a matrix M, this method computes the maximal non-zero minor (so, the rank of the matrix) and return the permutations needed to put that minor in the left top corner of the matrix.
- '''
- if(M.nrows() > M.ncols()):
- (o,j) = get_rearrange(M.transpose());
- return ((o[_sage_const_1 ],o[_sage_const_0 ]),j);
-
- done = False;
- j = M.nrows()+_sage_const_1 ;
- while((not done) and (j > _sage_const_1 )):
- j -= _sage_const_1 ;
- i = _sage_const_1 ;
- minors = M.minors(j);
- for el in minors:
- if(el != _sage_const_0 ):
- done = True;
- break;
- i += _sage_const_1 ;
-
- return (get_rearrange_from_minor(M.nrows(),M.ncols(),j,i), j);
-
-@wverbose
-def rearrange(M):
- '''
- Method that returns the same matrix M after changing rows and columns in order to have the rank defining minor in the top left corner of the matrix.
- '''
- (orders, size) = get_rearrange(M);
-
- return Matrix(M.parent().base(), [[M[orders[_sage_const_0 ](i+_sage_const_1 )-_sage_const_1 ][orders[_sage_const_1 ](j+_sage_const_1 )-_sage_const_1 ] for j in range(M.ncols())] for i in range(M.nrows())]), size;
-
-
-@wverbose
-def deep_rearrange(M):
- '''
- Method that returns the same matrix M after changing rows and columns in order to have all the possible principal minors non zero.
- '''
- M,rank = rearrange(M);
-
- nrows = M.nrows(); ncols = M.ncols();
- for i in range(nrows):
- minors = M.minors(i+_sage_const_1 );
- j = i;
- for j in range(i,ncols):
- if(minors[j-i] != _sage_const_0 ):
- break;
- if(j != i):
- M = swap(M,i+_sage_const_1 ,j+_sage_const_1 );
-
- return M;
-
-####################################################################################
-###
-### Hermite Normal Form Method
-###
-####################################################################################
-@wverbose
-def HNF_square(A, check_argument = True):
- '''
- Method that computes the Hermite Normal Form of a non-singular square matrix.
-
- We can disable the checking for non-singularity changing the argument 'check_argument'.
- '''
- if(A.ncols() != A.nrows()):
- raise TypeError("This method requires a square matrix");
-
- if(check_argument and A.determinant() == _sage_const_0 ):
- raise ValueError("This method requires a non-singular matrix");
-
- sverbose("Matrix before rearrange:\n%s" %A);
- A = deep_rearrange(A);
- sverbose("Matrix after rearrange:\n%s" %A);
-
- current = Matrix(A.parent().base(), [[el for el in row] for row in A]);
- sverbose.start_iteration(A.nrows()-_sage_const_1 , False, False);
- for i in range(_sage_const_1 ,A.nrows()):
- current = stage2(A,stage1(current,i),i);
- sverbose.next_iteration();
-
- return current;
-
-@wverbose
-def HNF(A):
- if(A.ncols() < A.nrows()):
- sverbose("Case with more rows than columns. Using the transpose case");
- return HNF(A.transpose()).transpose();
-
- parent = A.parent().base();
-
- A,rank = rearrange(A);
- sverbose("Matrix rearrange. The rank computed is %d" %rank);
- new_rows = [[el for el in A[i]] for i in range(rank)];
- for i in range(rank,A.ncols()):
- new_rows += [[kronecker_delta(j,i) for j in range(A.ncols())]];
-
- B = HNF_square(Matrix(parent,new_rows));
- sverbose("Computed the square Hermite Normal Form, we return the first %d rows" %rank);
-
- return Matrix(parent,[[el for el in B[i]] for i in range(rank)] + [[_sage_const_0 for j in range(B.ncols())] for i in range(rank,A.nrows())]);
-
-####################################################################################
-###
-### Package exports
-###
-####################################################################################
-#__all__ = ["HNF"];
-
-# This file was *autogenerated* from the file ./matrix.sage
-from sage.all_cmdline import * # import sage library
-
-_sage_const_1 = Integer(1); _sage_const_0 = Integer(0)
+from sage.all import * # import sage library
from sage.structure.element import is_Matrix;
The way to extend that simple elements is marked by argument `constant_or_identity`: if it is False, then we put a constant matrix fullfilled with such element. Otherwise, we build a diagonal matrix with that element in the diagonal (i.e. it is the identity matrix times the element).
'''
## we check that rows is a list of lists of the same size
- d = len(rows[_sage_const_0 ]);
- for i in range(_sage_const_1 , len(rows)):
+ d = len(rows[0 ]);
+ for i in range(1 , len(rows)):
if(d != len(rows[i])):
raise ValueError("Te rows provided can not be seen as a matrix");
## We check the sizes
rows_hights = [__check_row(row, parent) for row in rows];
- cols_widths = [__check_col([rows[i][j] for i in range(len(rows))], parent) for j in range(len(rows[_sage_const_0 ]))];
+ cols_widths = [__check_col([rows[i][j] for i in range(len(rows))], parent) for j in range(len(rows[0 ]))];
rows_with_matrix = [];
for i in range(len(rows)):
row_with_matrix = [];
- for j in range(len(rows[_sage_const_0 ])):
+ for j in range(len(rows[0 ])):
if(rows[i][j] in parent):
if(constant_or_identity):
row_with_matrix += [Matrix(parent, [[rows[i][j]*kronecker_delta(k,l) for l in range(cols_widths[j])] for k in range(rows_hights[i])])];
- If args has one element and it it a list, we build a diagonal matrix using such elements
- If args has one element and it is not a list, we build a diagonal matrix using such element repeated `size` times.
'''
- if(len(args) > _sage_const_1 ):
+ if(len(args) > 1 ):
return diagonal_matrix(parent, args);
- if(len(args) == _sage_const_1 and (not (type(args[_sage_const_0 ]) == list))):
- return diagonal_matrix(parent, [args[_sage_const_0 ] for i in range(kwds.get("size", _sage_const_1 ))]);
+ if(len(args) == 1 and (not (type(args[0 ]) == list))):
+ return diagonal_matrix(parent, [args[0 ] for i in range(kwds.get("size", 1 ))]);
- list_of_elements = args[_sage_const_0 ];
+ list_of_elements = args[0 ];
d = len(list_of_elements);
rows = [];
for i in range(d):
- rows += [[_sage_const_0 for j in range(i)] + [list_of_elements[i]] + [_sage_const_0 for j in range(i+_sage_const_1 ,d)]];
+ rows += [[0 for j in range(i)] + [list_of_elements[i]] + [0 for j in range(i+1 ,d)]];
return block_matrix(parent, rows);
def vandermonde_matrix(parent, *args, **kwds):
- if(len(args) > _sage_const_1 ):
+ if(len(args) > 1 ):
return vandermonde_matrix(parent, args, **kwds);
else:
- casted = [parent(el) for el in args[_sage_const_0 ]];
+ casted = [parent(el) for el in args[0 ]];
try:
ncols = kwds["size"];
return Matrix(parent, rows);
def random_matrix(parent, *args, **kwds):
- if(len(args) > _sage_const_1 ):
+ if(len(args) > 1 ):
return random_matrix(parent, args, **kwds);
else:
- casted = [int(el) for el in args[_sage_const_0 ]];
+ casted = [int(el) for el in args[0 ]];
try:
ncols = kwds["size"];
ncols = kwds["ncols"];
nrows = kwds["nrows"];
except KeyError:
- if(len(casted) > _sage_const_1 ):
- nrows = casted[_sage_const_0 ];
- ncols = casted[_sage_const_1 ];
+ if(len(casted) > 1 ):
+ nrows = casted[0 ];
+ ncols = casted[1 ];
else:
- nrows = casted[_sage_const_0 ];
- ncols = casted[_sage_const_1 ];
+ nrows = casted[0 ];
+ ncols = casted[1 ];
try:
min = kwds["min"];
except KeyError:
- min = -sys.maxint-_sage_const_1 ;
+ min = -sys.maxint-1 ;
try:
max = kwds["max"];
except KeyError:
elif(isinstance(X, list)):
try:
is_matrix = True;
- d = len(X[_sage_const_0 ]);
- i = _sage_const_1 ;
+ d = len(X[0 ]);
+ i = 1 ;
while(is_matrix and i < len(X)):
is_matrix = (is_matrix and d == len(X[i]));
- i += _sage_const_1 ;
+ i += 1 ;
if(is_matrix):
return True;
elif(not (el in parent)):
raise ValueError("The row has not the proper format -- non-castable element");
if(hight is None):
- return _sage_const_1 ;
+ return 1 ;
return hight;
def __check_col(col,parent):
raise ValueError("The col has not the proper format -- non-castable element");
if(width is None):
- return _sage_const_1 ;
+ return 1 ;
return width;
def __ncols(X):
try:
return X.ncols();
except AttributeError:
- return len(X[_sage_const_0 ]);
+ return len(X[0 ]);
except Exception:
raise TypeError("Impossible to compute the number of columns for element %s" %X);
def turn_matrix(M, vertical=False):
if(vertical):
- new_rows = [[M[i][j] for j in range(M.ncols())] for i in range(-_sage_const_1 ,-M.nrows()-_sage_const_1 ,-_sage_const_1 )];
+ new_rows = [[M[i][j] for j in range(M.ncols())] for i in range(-1 ,-M.nrows()-1 ,-1 )];
else:
- new_rows = [[M[i][j] for j in range(-_sage_const_1 ,-M.ncols()-_sage_const_1 ,-_sage_const_1 )] for i in range(M.nrows())];
+ new_rows = [[M[i][j] for j in range(-1 ,-M.ncols()-1 ,-1 )] for i in range(M.nrows())];
return Matrix(M.parent().base(), new_rows);
def matrix_of_dMovement(M,v,D, cols):
res = [v];
- for j in range(_sage_const_1 ,cols):
- res += [__dm(M,res[-_sage_const_1 ],D)];
+ for j in range(1 ,cols):
+ res += [__dm(M,res[-1 ],D)];
return Matrix(M.parent().base(), res).transpose();
+++ /dev/null
-
-# This file was *autogenerated* from the file ./nullspace_wrapper.sage
-from sage.all_cmdline import * # import sage library
-
-_sage_const_1 = Integer(1); _sage_const_0 = Integer(0)
-from .matrix import swap_rows;
-from .matrix import swap_cols;
-
-from .cached_property import derived_property;
-
-from .verbose import *;
-
-from ore_algebra.nullspace import *;
-
-from sage.rings.polynomial.polynomial_ring import is_PolynomialRing as isUniPolynomial;
-from sage.rings.polynomial.multi_polynomial_ring import is_MPolynomialRing as isMPolynomial;
-
-class ElementInIdealFound(ValueError):
- # Used for returning an exception with an element that is zero but is not the zero polynomial
- def __init__(self, element):
- self.__element = element;
- print "ElementInIdealFound raised"
-
- def __str__(self):
- return "A zero found: %s" %self.__element;
-
- def element(self):
- return self.__element;
-
-class NullspaceAlgorithm(object):
- ### Initialization method
- def __init__(self, parent, M, method=None, relations = []):
- '''
- This class represent the application of the Bareiss Algorithm over a matrix with polynomial coefficients.
- This class receives several parameters:
- - `parent`: the polynomial ring where the operations will be done. A fraction field is allowed if the base domain is a polynomial ring.
- - `M`: the matrix we want to treat with the Algorithm
- - `method`: a Boolean function
-
- This algorithm also works in a quotient of the polynomial ring (K[X]/I) and the method provided must be a method to check if some polynomial is in the ideal I.
- '''
- ## Checking the parent parameter
- if(parent.is_field()):
- parent = parent.base();
- if(not (isUniPolynomial(parent) or isMPolynomial(parent))):
- raise TypeError("The parent for this algorithm must be a polynomial ring.\n\t Got: %s" %parent);
- if((not (parent.base() is QQ)) and (not(parent.base() is ZZ))):
- raise TypeError("The parent for these algorithms must be based on rational numbers\n\t Got: %s" %parent);
- self.__parent = PolynomialRing(ZZ,parent.gens());
-
- ## Checking the matrix input
- self.base_matrix = Matrix(parent, M);
- ## If the parent are the Rational numbers we clean denominators
- ## Otherwise, the parent are the Integers, so we can go on.
- if(parent.base() is QQ):
- ## Cleaning rational denominators
- new_rows = [];
- for row in M:
- aux_list = [];
- for el in row:
- aux_list += [coeff.denominator() for coeff in el.coefficients()];
- den_lcm = lcm(aux_list);
- new_rows += [den_lcm*row];
- M = new_rows;
- self.__matrix = Matrix(self.parent(), M);
- else:
- self.__matrix = self.base_matrix;
-
- ## Storing the checking method
- self.__in_ideal = method;
-
- ## Cached elements
- self.__gb = ideal(self.parent(), relations).groebner_basis();
-
- ### Getter methods
- def parent(self):
- '''
- Method to have the same interface that the majority of objects of SAGE. Returns the polynomial ring where the matrix have its coefficients.
- '''
- return self.__parent;
-
- #################################################
- ### Linear algebra methods
- #################################################
- @derived_property
- def rank(self):
- return self.base_matrix.nrows()-self.right_kernel_matrix().ncols();
-
- def relations(self):
- self.right_kernel_matrix();
- return self.__gb;
-
- @cached_method
- def right_kernel_matrix(self):
- solver = self.__decide_solver();
-
- finished = False;
- output = None;
- while(not finished):
- try:
- output = solver(self.__matrix);
- finished = True;
- except ElementInIdealFound as e:
- if(e.element().reduce(self.__gb) == _sage_const_0 ):
- raise RuntimeError("Same zero detected twice. Probable infinite loop");
- self.__gb = ideal(self.parent(), self.__gb+[e.element()]).groebner_basis();
- self.__matrix = self.__simplify_matrix(self.__matrix);
- ## Putting the output in matrix form
- rows = [[_sage_const_0 for i in range(len(output))] for j in range(self.__matrix.ncols())];
- for i in range(len(output)):
- for j in range(self.__matrix.ncols()):
- rows[j][i] = output[i][j];
- return Matrix(self.parent(), rows);
-
- #################################################
- ### Other cached methods
- #################################################
- @cached_method
- def is_in_ideal(self, p):
- p = self.parent()(p);
-
- try:
- return self.__in_ideal(p) is True;
- except Exception:
- return False;
-
- #################################################
- ### Private methods for the Wrapper
- #################################################
- @wverbose
- def __pivot(self,M, r, n, c, m, zero):
- # We delegate to the method _pivot of ore_algebra.nullspace
- from ore_algebra.nullspace import _pivot;
- (i,j) = _pivot(M,r,n,c,m,zero);
- sverbose("A pivot found at position (%d,%d): %s" %(i,j,repr(M[i][j])));
- if(self.is_in_ideal(M[i][j])):
- sverbose("The pivot is zero");
- raise ElementInIdealFound(M[i][j]);
-
- return (i,j);
-
- def __decide_solver(self):
- return kronecker(gauss(lambda M,r,n,c,m,zero : self.__pivot(M,r,n,c,m,zero)));
-
- def __have_ideal(self, basis):
- if(len(basis) == _sage_const_1 ):
- return basis[_sage_const_0 ] != _sage_const_0 ;
- return len(basis) > _sage_const_0 ;
-
- def __simplify_matrix(self, M):
- rows = [[el for el in row] for row in M];
- if(self.__have_ideal(self.__gb)):
- rows = [[el.reduce(self.__gb) for el in row] for row in rows];
-
- return Matrix(self.parent(), rows);
-
-
-
-# This file was *autogenerated* from the file ./ring_w_sequence.sage
-from sage.all_cmdline import * # import sage library
-
-_sage_const_1 = Integer(1); _sage_const_0 = Integer(0)
+from sage.all import * # import sage library
########################################################
-###
from sage.rings.ring import Ring;
def getSequenceElement(self, el, n):
self_gen = 'x';
try:
- self_gen = self.base().gens()[-_sage_const_1 ];
+ self_gen = self.base().gens()[-1 ];
except:
pass;
- if(self_gen == _sage_const_1 ):
- if(n > _sage_const_0 ):
- return _sage_const_0 ;
- elif(n==_sage_const_0 ):
+ if(self_gen == 1 ):
+ if(n > 0 ):
+ return 0 ;
+ elif(n==0 ):
return el;
else:
raise ValueError("Impossible get negative element of the sequence");
res = derivative(res,self_gen);
except AttributeError:
## Not derivative possible. Considering a constant
- if(n == _sage_const_0 ):
+ if(n == 0 ):
return res;
- return _sage_const_0 ;
+ return 0 ;
try:
- return res(**{repr(self_gen):_sage_const_0 })/factorial(n);
+ return res(**{repr(self_gen):0 })/factorial(n);
except TypeError:
## Not callable element. Returning the element without evaluation
return res/factorial(n);
+++ /dev/null
-
-# This file was *autogenerated* from the file ./storj_villard.sage
-from sage.all_cmdline import * # import sage library
-
-_sage_const_100 = Integer(100); _sage_const_2 = Integer(2); _sage_const_1 = Integer(1); _sage_const_0 = Integer(0)
-from .cached_property import derived_property;
-
-from sage.rings.polynomial.polynomial_ring import is_PolynomialRing as isUniPolynomial;
-
-class StorjVillardAlgorithm(object):
- ### Initialization method
- def __init__(self, parent, M):
- '''
- This class represent the application of the Storjohann-Villard algorithm (ISSAC 2005 paper) over a matrix with polynomial coefficients.
- This class receives several parameters:
- - `parent`: the polynomial ring where the operations will be done. A fraction field is allowed if the base domain is a polynomial ring.
- - `M`: the matrix we want to treat with the Algorithm
- '''
- ## Checking the parent parameter
- if(parent.is_field()):
- parent = parent.base();
- if(not (isUniPolynomial(parent))):
- raise TypeError("The parent for this algorithm must be a univariate polynomial ring.\n\t Got: %s" %parent);
- self.__parent = parent;
-
- self.__c_rank = None;
- self.__c_nulls = None;
-
- ## Checking the matrix input
- self.base_matrix = Matrix(self.parent(), M);
-
- ### Getter methods
- def parent(self):
- '''
- Method to have the same interface that the majority of objects of SAGE. Returns the polynomial ring where the matrix have its coefficients.
- '''
- return self.__parent;
-
- #################################################
- ### Linear algebra methods
- #################################################
- @derived_property
- def rank(self):
- if(self.__c_rank is None):
- return self.__nullspace()[_sage_const_0 ];
- return self.__c_rank;
-
- @cached_method
- def right_kernel_matrix(self):
- return self.symmetric_problem.left_kernel_matrix();
-
- @cached_method
- def left_kernel_matrix(self):
- if(self.__c_nulls is None):
- return self.__nullspace()[_sage_const_1 ];
- return self.__c_nulls;
-
- #################################################
- ### Other cached methods
- #################################################
- @derived_property
- def symmetric_problem(self):
- return StorjVillardAlgorithm(self.parent(), M.transpose());
-
- #################################################
- ### Private methods for Storjohann-Villard Algorithm
- #################################################
- def __null_min_vector(self, M, thres):
- '''
- This method receives a matrix M and a degree threshold and computes some vectors in the
- left null-space of M.
-
- This requires that M has more rows than columns and a maximal rank.
-
- This operation use some random matrices.
- '''
- if(M.nrows() <= M.ncols()):
- raise TypeError("Invalid matrix: more columns than rows");
- n = M.ncols();
- p = M.nrows()-n;
- d = max([max([el.degree() for el in row]) for row in M]);
-
- ## M = QM for a random matrix Q
- Q = random_matrix(self.parent(), n+p, min=-_sage_const_100 ,max=_sage_const_100 );
- M = Q*M;
- ## M = M(x+x0) for a random x0
- x0 = self.__rand_pol(_sage_const_0 );
- M = Matrix(self.parent(), [[el(x=x+x0) for el in row] for row in M]);
-
- ## Splitting the matrix into a square matrix and remaining rows
- A = Matrix(self.parent(), [[M[i][j] for j in range(n)] for i in range(n)]);
- B = Matrix(self.parent(), [[M[i][j] for j in range(n)] for i in range(n, n+p)]);
-
- ## Checking rank of M
- if(A.determinant() == _sage_const_0 ):
- raise ValueError("Failure in the algorithm. Probably rank(M) lower than require");
-
- ## Computing some degree bound
- nu = thres+d+int(ceil(n*d/p));
-
- ## Computing a expansion of BA^-1 mod x^nu
- H = raise NotImplementedError("Do not know how to implement this");
-
- ## Compressing the matrix H
- Hp = H*Matrix(self.parent(), [[__rand_pol(d) for i in range(p)] for j in range(n)]);
- t = [d-_sage_const_1 for i in range(p)] + [_sage_const_0 for i in range(p)];
-
- ## Compute a sigma-basis with respect to t for [-Ip Hp^T]^T of order nu
- L = raise NotImplementedError("Do not know how to implement this");
-
- ## Computing kappa = #{rows of L such t-degree is at most thres}
- ordered_rows = {};
- kappa = _sage_const_0
- for row in L:
- t_degree = max([row[i] - t[i] for i in range(len(t))]);
- if(t_degree <= thres):
- kappa += _sage_const_1 ;
-
- if(ordered_rows.get(t_degree) is None):
- ordered_rows[t_degree] = [];
- ordered_rows[t_degree] = ordered_rows[t_degree] + [row];
- S = [];
- for i in range(thres+_sage_const_1 ):
- S += ordered_rows[i];
-
- ## Computing the partial reconstructions Ni
- N = [el*block_matrix(self.parent(),H,_sage_const_1 ) for i in el];
- N = [Matrix(self.parent(), [[el%((self.parent().gens()[_sage_const_0 ])**(thres+_sage_const_1 )) for el in row] for row in Ni]) for Ni in N];
- N = [Matrix(self.parent(), [[el(x=x-x0) for el in row] for row in Ni]) for Ni in N];
- N = [Ni*Q for Ni in N];
- l = [Ni*M].count(_sage_const_0 );
-
- if(l != kappa):
- raise ValueError("Failure in the algorithm. Error certifiating kappa");
-
- N = Matrix(self.parent(), N);
-
- lc_N = Matrix(self.parent().base(), [[el.leading_coefficient() for el in row] for row in N]);
- if(lc_N.rank() != max(lc_N.nrows(), lc_N.ncols())):
- raise ValueError("Failure in the algorithm. No minimality got");
-
- return kappa, N;
-
- def __nullspace2n(self, M):
- '''
- This method compute some `small` linearly independent polynomial vectors in the nullspace of a matrix M.
-
- Such matrix M must have more rows than columns, but no more than twice.
-
- Moreover, M must have maximal rank.
- '''
- if(M.nrows() <= M.ncols() or M.nrows() > _sage_const_2 *M.ncols()):
- raise TypeError("Invalid matrix: invalid amount of rows with respect to columns");
-
- n = M.ncols();
- q = M.nrows()-n;
- d = max([max([el.degree() for el in row]) for row in M]);
-
- ## M = QM for a random matrix Q
- Q = random_matrix(self.parent(), n+q, min=-_sage_const_100 ,max=_sage_const_100 );
- M = Q*M;
-
- ## Splitting the matrix into a square matrix and remaining rows
- A = Matrix(self.parent(), [[M[i][j](x=self.__rand_pol(_sage_const_0 )) for j in range(n)] for i in range(n)]);
-
- if(A.determinant() == _sage_const_0 ):
- raise ValueError("Failure in the algorithm. Non singular upper part of matrix");
-
- I = set();
- p = q;
-
- res = [];
-
- while(len(I) < q):
- ## Getting indices left
- index_to_check = set([n+i for i in range(_sage_const_1 ,q+_sage_const_1 )])-I;
-
- ## Computing a degree bound
- delta = _sage_const_2 *n*d/p;
-
- ## Building a matrix C
- rows = [[_sage_const_0 for j in range(n+q)] for i in range(n+p)];
- for i in range(n):
- rows[i][i] = _sage_const_1 ;
- for j in range(p):
- rows[n+j][index_to_check[j]] = _sage_const_1 ;
- C = Matrix(self.parent().base(), rows);
-
- ## Using C and M to use the null_min_vector method
- M1 = C*M;
- kappa, D = self.__null_min_vector(M1, delta);
- N = Matrix(self.parent(),[Di*C for Di in D]);
-
- ## Getting new indices
- J = raise NotImplementedError("Do not know how to implement this");
-
- ## Updating loop variables
- I = I + set(J);
- p = p - kappa;
-
- ## Nullspace update
- res = block_matrix(self.parent(), res.transpose(), N.transpose()).transpose();
-
- return res*Q;
-
- def __nullspace(self):
- '''
- This method privately compute the left kernel of the matrix and its rank and stores the result in cached variables of the structure.
-
- This method returns both the rank and the minimal basis.
- '''
- m = self.base_matrix.nrows();
- n = self.base_matrix.ncols();
- d = max([max([el.degree() for el in row]) for row in M]);
-
- ## Computing a guess r0 for rank
- r0 = raise NotImplementedError("Do not know how to implement this");
- M1 = raise NotImplementedError("Do not know how to implement this");
-
- ## If the guess is too big, we can return an empty nullspace
- if(m == r0):
- self.__c_rank = m;
- self.__c_nulls = [];
-
- return self.__c_rank, self.__c_nulls;
-
- ## Randomly ensure that the top r0xr0 submatrix of M1 is non-singular
- raise NotImplementedError("Do not know how to implement this");
-
- s = int(ceil((m-_sage_const_2 *r0)/r0));
-
- ## Getting the
- N = [];
- for k in range(s):
- N += [row for row in self.__null_min_vector(raise NotImplementedError("Do not know how to implement this"), d)[_sage_const_1 ]];
- N += [row for row in self.__nullspace2n(raise NotImplementedError("Do not know how to implement this"))];
-
- N = Matrix(self.parent(), N);
- if(N*M != _sage_const_0 ):
- raise ValueError("Failure in the algorithm. No kernel computed");
-
- self.__c_rank = r0;
- self.__c_nulls = N;
-
- return self.__c_rank, self.__c_nulls;
-
-
- def __rand_pol(self, degree):
- return self.parent()([randint(-_sage_const_100 ,_sage_const_100 ) for i in range(randint(_sage_const_0 ,degree)+_sage_const_1 )]);
-
+++ /dev/null
-
-# This file was *autogenerated* from the file ./timing.sage
-from sage.all_cmdline import * # import sage library
-
-_sage_const_3 = Integer(3); _sage_const_2 = Integer(2); _sage_const_1 = Integer(1); _sage_const_0 = Integer(0); _sage_const_60 = Integer(60); _sage_const_24 = Integer(24)
-
-
-##########################################
-###
-### TIME OPERATIONS
-###
-##########################################
-from timeit import default_timer as timer;
-
-__GLOBAL_TIMER = None;
-__GENERAL_TIMER = None;
-def start_timer(reset = False):
- global __GLOBAL_TIMER;
- global __GENERAL_TIMER;
-
- if(reset):
- __GENERAL_TIMER = timer();
-
- __GLOBAL_TIMER = timer();
-
-def finish_timer(text = "an action"):
- global __GLOBAL_TIMER;
- t = timer() - __GLOBAL_TIMER;
- __print_timer(t, "Perfomed %s: it took" %(text));
-
-def get_general_timing():
- global __GENERAL_TIMER;
- t = timer() - __GENERAL_TIMER;
- __print_timer(t, "General timing:");
-
-def __print_timer(T, text):
- END = "secs.";
- if(T > _sage_const_60 ):
- END = "mins.";
- T = float(T/_sage_const_60 );
-
- if(T > _sage_const_60 ):
- END = "h.";
- T = float(T/_sage_const_60 );
-
- if(T > _sage_const_24 ):
- END = "days"
- T = float(T/_sage_const_24 );
- print "%s %s %s" %(text, T, END);
-
-##########################################
-###
-### COMPARISON BETWEEN METHODS
-###
-##########################################
-def __compare_methods(arguments, *args):
- if(len(args) == _sage_const_1 and (type(args[_sage_const_0 ]) == list)):
- args = args[_sage_const_0 ];
-
- res = [arguments];
- for method in args:
- try:
- init = timer();
- method(*arguments);
- time = timer()-init;
- res += [time];
- except Exception as e:
- print e
- res += [-_sage_const_1 ];
- return res;
-
-def __compare_for_entries(entries, *args):
- if(len(args) == _sage_const_0 and type(args[_sage_const_0 ]) == list):
- args = args[_sage_const_0 ];
- try:
- res = [["Entry\Alg"] + [el.__name__ for el in args]];
- for arguments in entries:
- res += [__compare_methods(arguments, *args)];
-
- return res;
- except AttributeError:
- raise ValueError("The args must be a list of methods");
-
-def __compute_global_stats(results):
- try:
- res = [["Stat\Alg"] + [el for el in results[_sage_const_0 ][_sage_const_1 :]]];
- reverse = [[results[i][j] for i in range(len(results))] for j in range(len(results[_sage_const_0 ]))];
- stats = [];
- for i in range(_sage_const_1 , len(results[_sage_const_0 ])):
- stats += [__compute_stats(reverse[i][_sage_const_1 :])];
-
- res += [["Mean"] + [stats[j][_sage_const_0 ] for j in range(len(stats))]];
- res += [["Min"] + [stats[j][_sage_const_1 ] for j in range(len(stats))]];
- res += [["Max"] + [stats[j][_sage_const_2 ] for j in range(len(stats))]];
- res += [["N. Err"] + [stats[j][_sage_const_3 ] for j in range(len(stats))]];
-
- return res;
- except (TypeError, IndexError):
- raise TypeError("The input for compute stats must be a list of lists with at least one row");
-
-def __compute_stats(results):
- cMin = float('inf');
- cMax = float('-inf');
- cSum = _sage_const_0 ;
-
- for el in results:
- cSum += el;
- cMin = min(cMin, el);
- cMax = max(cMax, el);
-
- return [cSum/len(results), cMin, cMax, results.count(-_sage_const_1 )];
-
-def full_comparison(entries, *args):
- if(len(args) == _sage_const_1 and (type(args[_sage_const_0 ]) == list)):
- args = args[_sage_const_0 ];
-
- return __compute_global_stats(__compare_for_entries(entries, *args));
-
-# This file was *autogenerated* from the file ./verbose.sage
-from sage.all_cmdline import * # import sage library
+from sage.all import * # import sage library
-_sage_const_100 = Integer(100); _sage_const_1 = Integer(1); _sage_const_0 = Integer(0)
####################################################################
###
-### verbose.sage
+### verbose.py
###
### Code for perform an easy conditional printing.
###
### - To change the configuration of the sverbose, call the public methods of the class
### - Added a decorator to easily set up a function to verbose
###
+### Version 1.1 (18-06-2018)
+### - Move file from .sage to .py
+###
####################################################################
import sys;
class __verbose_session(object):
def __init__(self):
self.__time = True;
- self.__current_depth = _sage_const_0 ;
- self.__level = _sage_const_1 ;
+ self.__current_depth = 0 ;
+ self.__level = 1 ;
self.__default_beg = "";
self.__begs = ["** ","## ","-- ","%% ","&& ","@@ "];
self.__end = "\n";
- self.__deep_wrapper = _sage_const_0 ;
+ self.__deep_wrapper = 0 ;
- self.__n_iterations = _sage_const_0 ;
+ self.__n_iterations = 0 ;
self.__perc = False;
self.__inline = True;
- self.__current_iteration = _sage_const_0 ;
+ self.__current_iteration = 0 ;
## Basic call method
def __call__(self,string, end=None):
print "|",str(time.now()), "|",
for i in range(self.__current_depth):
print "\t",
- if(self.__current_depth > _sage_const_0 ):
- print self.__begs[(self.__current_depth-_sage_const_1 )%len(self.__begs)],
+ if(self.__current_depth > 0 ):
+ print self.__begs[(self.__current_depth-1 )%len(self.__begs)],
else:
print self.__default_beg,
## Changing depth methods
def reset_depth(self):
- self.__current_depth = _sage_const_0 ;
+ self.__current_depth = 0 ;
def change_depth(self,amount):
- self.__current_depth = max(_sage_const_0 , self.__current_depth+amount);
+ self.__current_depth = max(0 , self.__current_depth+amount);
- def increase_depth(self,amount = _sage_const_1 ):
- if(amount < _sage_const_0 ):
+ def increase_depth(self,amount = 1 ):
+ if(amount < 0 ):
raise ValueError("To increase the tab depth it is needed a positive number, not %d" %amount);
self.change_depth(amount);
- def decrease_depth(self,amount = _sage_const_1 ):
- if(amount < _sage_const_0 ):
+ def decrease_depth(self,amount = 1 ):
+ if(amount < 0 ):
raise ValueError("To decrease the tab depth it is needed a positive number, not %d" %amount);
self.change_depth(-amount);
self.__n_iterations = iterations;
self.__perc = percentage;
self.__inline = in_line;
- self.__current_iteration = _sage_const_1 ;
+ self.__current_iteration = 1 ;
def next_iteration(self):
msg = "Performed %d-th iteration " %self.__current_iteration;
- if(self.__perc and self.__n_iterations > _sage_const_0 ):
- msg += "(%.2f %%)" %(self.__current_iteration*_sage_const_100 /self.__n_iterations);
+ if(self.__perc and self.__n_iterations > 0 ):
+ msg += "(%.2f %%)" %(self.__current_iteration*100 /self.__n_iterations);
else:
msg += "(%d/%d)" %(self.__current_iteration, self.__n_iterations);
end = None;
self(msg, end);
- self.__current_iteration += _sage_const_1 ;
+ self.__current_iteration += 1 ;
## Conditional definition of the verbose variables
sverbose = __verbose_session();
sverbose.increase_depth(sverbose.get_deep_wrapper());
sverbose("Calling method %s" %(func.__name__));
sverbose("Arguments:");
- if(len(args) > _sage_const_0 ):
+ if(len(args) > 0 ):
sverbose("%s" %(str(args)));
- if(len(kwds)> _sage_const_0 ):
+ if(len(kwds)> 0 ):
sverbose("%s" %(str(kwds)));
state = sverbose.get_iteration_state();
sverbose.increase_depth();
try:
set_verbose(get_verbose());
except Exception:
- __VERBOSE_LEVEL = _sage_const_0 ;
+ __VERBOSE_LEVEL = 0 ;
def set_verbose(level):
global __VERBOSE_LEVEL;
- __VERBOSE_LEVEL = max(_sage_const_0 , level);
- if(__VERBOSE_LEVEL >= _sage_const_1 ):
+ __VERBOSE_LEVEL = max(0 , level);
+ if(__VERBOSE_LEVEL >= 1 ):
print "-- Changed verbose level to %d" %(__VERBOSE_LEVEL);
def get_verbose():
return __VERBOSE_LEVEL;
__all__ += ["set_verbose", "get_verbose"];
-
-
-## Getters and setters of the verbose level
-#def set_verbose(level):
-# global __VERBOSE_LEVEL;
-#
-# __VERBOSE_LEVEL = max(0, level);
-# verbose("-- Changed verbose level to %d" %(__VERBOSE_LEVEL), 1);
-#
-#def get_verbose():
-# global __VERBOSE_LEVEL;
-# return __VERBOSE_LEVEL;
-
-## Main printing function
-#def verbose(text, level = 0, end = '\n', depth = 0, beg = '', time = False):
-# global __VERBOSE_LEVEL;
-#
-# if(level <= __VERBOSE_LEVEL):
-# if(time):
-# from datetime import datetime as time;
-# print str(time.now()), " --",
-# for i in range(depth):
-# print "\t",
-# if(len(beg) > 0):
-# print beg, " ",
-# print text,
-# print end,
-# sys.stdout.flush();
-
-## Other useful methods
-#__sv_active = False;
-#__sv_time = False;
-#__sv_current_depth = 0;
-#__sv_current_level = 0;
-#__sv_beg = "";
-#__sv_end = "\n";
-
-## Making smart verbose
-#def sverbose(text, end=None):
-# global __sv_active;
-# global __sv_current_depth;
-# global __sv_current_level;
-# global __sv_beg;
-# global __sv_time;
-# global __sv_end;
-# if(end is None):
-# tend = __sv_end;
-# else:
-# tend = end;
-#
-# if(__sv_active):
-# verbose(text, __sv_current_level, end=tend, depth=__sv_current_depth, beg=__sv_beg, time=__sv_time);
-
-## Starting and finishing sessions
-#def start_smart_verbose(level, time=True):
-# global __sv_active;
-# global __sv_current_depth;
-# global __sv_current_level;
-# global __sv_time;
-#
-# if(__sv_active):
-# verbose("A smart verbose session was already initiated. Reseting session...", 1, beg = "--");
-
-# __sv_active = True;
-# __sv_time = time;
-# __sv_current_depth = 0;
-# __sv_current_level = level;
-#
-#def end_smart_verbose():
-# global __sv_active;
-#
-# __sv_active = False;
-
-## Changing depth in smart session
-#def change_depth(amount):
-# global __sv_active;
-# global __sv_current_depth;
-#
-# if(__sv_active):
-# __sv_current_depth = max(0, __sv_current_depth+amount);
-#
-#def increase_depth(amount = 1):
-# if(amount <= 0):
-# raise ValueError("To increase the tab depth it is needed a positive number, not %d" %amount);
-# change_depth(amount);
-#
-#def decrease_depth(amount = 1):
-# if(amount <= 0):
-# raise ValueError("To decrease the tab depth it is needed a positive number, not %d" %amount);
-# change_depth(-amount);
-#
-#def change_beginning(beg = ""):
-# global __sv_beg;
-# __sv_beg = beg;
-
-## Printing iterations in smart session
-#__sv_iteration = False;
-#__sv_n_iterations = 0;
-#__sv_perc = False;
-#__sv_inline = True;
-#__sv_current_iteration = 0;
-#def start_iteration(iterations, percentage = False, in_line = True):
-# global __sv_iteration;
-# global __sv_n_iterations;
-# global __sv_perc;
-# global __sv_inline;
-# global __sv_current_iteration;
-#
-# __sv_iteration = True;
-# __sv_n_iterations = iterations;
-# __sv_perc = percentage;
-# __sv_inline = in_line;
-# __sv_current_iteration = 1;
-#
-#def next_iteration():
-# global __sv_iteration;
-# global __sv_n_iterations;
-# global __sv_perc;
-# global __sv_inline;
-# global __sv_current_iteration;
-#
-# if(__sv_iteration):
-# msg = "Performed %d-th iteration " %__sv_current_iteration;
-#
-# if(__sv_perc and __sv_n_iterations > 0):
-# msg += "(%.2f %%)" %(__sv_current_iteration*100/__sv_n_iterations);
-# else:
-# msg += "(%d/%d)" %(__sv_current_iteration, __sv_n_iterations);
-#
-# if(__sv_inline):
-# end = "\r";
-# else:
-# end = None;
-#
-# sverbose(msg, end);
-# __sv_current_iteration += 1;
-#
-#def end_iteration():
-# global __sv_iteration;
-# global __sv_inline;
-# global __sv_current_level;
-# global __VERBOSE_LEVEL;
-#
-# if(__sv_iteration):
-# __sv_iteration = False;
-# if(__sv_inline and __sv_current_level <= __VERBOSE_LEVEL):
-# print "";
-
-
projects / ajpastor / diff_defined_functions.git / commitdiff