res += string[current_index:]
return res;
-def _get_equation_poly(poly, dic):
- raise DevelopementError("_get_equation_poly");
+def _get_equation_poly(poly, dic):
+ '''
+ Method that computes the differential equation that the evaluation
+ of a polynomial over some functions satisfies.
+
+ This procedure has 5 main steps:
+ (1) Compute a vector space $W$ generated by a vector of generators $(g1,...,gk)$
+ (2) Compute a derivation matrix $D$ for $W$ w.r.t. the generators
+ (3) Compute a initial vector $v_0$ that represents poly in $W$.
+ (4) Compute the ansatz system $S$
+ (5) Solve the ansatz system
+
+ The steps (1), (2), (3) can be done simultaneously in the method _get_derivation_matrix(poly,dic).
+ The step (4) use the inductive formula $v_{i+1} = Dv_i + \partial v_i$. Then $S = (v0 | v1 | ...)$.
+ The step (5) will be done using Bareiss Algorithm.
+ '''
+ ### Steps (1)-(3)
+ D, v0 = _get_derivation_matrix(poly, dic);
+
+ ### Step 4
+ R = D[0][0].parent();
+ if(sage.rings.fraction_field.is_FractionField(R)):
+ R = R.base();
+ F = R.fraction_field();
+
+ v = [v0];
+ for i in range(D.nrows()):
+ v += [D*v[-1] + vector(F, [el.derivative() for el in v[-1]])];
+
+ S = Matrix(F, v).transpose();
+
+ ### Step 5
+ ## TODO We rely on the system solver of the operator of the functions we are
+ ## plug-in to the polynomial.
+ solution = dic.keys()[0].equation._get_element_nullspace(S);
+
+ ## We create the final operator
+ return dic.keys()[0].parent().element([el for el in solution]).equation;
+
+def _get_derivation_matrix(poly, dic, simpl = False):
+ '''
+ Method that computes the derivation matrix of a vector space that
+ contains the vector space generated by the evaluation of poly
+ using the dictionary dic and a vector that represent poly into
+ that vector space.
+
+
+ This method has a recursive structure:
+ - If for any key in dic, dic[key][i] appears in poly for i > 0,
+ then we simplify the polynomial so it only contains the variables
+ dic[key][0]. We return the matrix obtained with that polyomial and
+ the associated vector.
+ - If poly has multiple terms, we return the direct sum of the matrices
+ for each term and the direct sum of the vectors.
+ - If poly is a monomial with multiple variables, we return the Kronecker
+ sum of the matrices of the result and the tensor product of the vectors.
+ - If poly is a monomial with just one variable, we compute the
+ Kronecker sum of the companion matrix associated with the function
+ related with that variable degree times and return the vector (1,0,...,0).
+
+ WARNING:
+ - The argument simpl must never be called with True but within this method
+ - The very base case can be improve using the fact that we are computing
+ the power of an element.
+ '''
+ ## Simplification step
+ #if(not simpl):
+ # non_simple = sum([dic[key][1:] for key in dic], []);
+ # if(any(var in poly.variables() for var in non_simple)):
+ # to_eval = {}
+ # for key in dic:
+ # current = dic[key];
+ # for i in range(1,len(current)):
+ # to_eval[str(current[i])] = current[0];
+ # D, v = _get_derivation_matrix(poly(**to_eval, dic, simpl=True);
+ # ## TODO
+ # raise DevelopementError("_get_derivation_matrix");
+
+ ## From this point on, we can assume only simplified polynomials appears
+ ## Case non-monomial
+ if(not poly.is_monomial()):
+ coeffs = poly.coefficients();
+ sols = [_get_derivation_matrix(mon, dic, simpl=True) for mon in poly.monomials()];
+ D = Matrix.block_diagonal(*[sol[1] for sol in sols]);
+ v = vector(D[0][0].parent(), sum([[coeffs[i]*el for el in sols[i][2]] for i in range(len(sols))], []);
+
+ return D,v;
+ ## More than one variable
+ elif(len(poly.variables()) > 1):
+ var = poly.variables()[0];
+ deg = poly.degree(var);
+ oth = poly.parent()(poly/(var**deg));
+
+ D1,v1 = _get_derivation_matrix(var**deg, dic);
+ D2,v2 = _get_derivation_matrix(oth, dic);
+
+ D = ajpastor.misc.Matrix.kronecker_sum(D1, D2);
+ v = vector(v1.parent().base(), sum(([el for el in row] for row in v1.tensor_product(v2)), []));
+
+ return D,v;
+ ## A power of a variable (using Faa di Bruno's formula)
+ elif(len(poly.variables()) == 1):
+ var = poly.variables()[0];
+ deg = poly.degree();
+ ## Getting the function and the offset for the initial values
+ func = None;
+ for key in dic.keys():
+ if(var in dic[key]):
+ func = key;
+ offset = dic[key].index(var);
+ break;
+ ## Checking the variable can be easily represent
+ if(offset > func.getOrder()):
+ raise DevelopementError("_get_derivation_matrix");
+ ## TODO Go on from here
+
+
+ raise DevelopementError("_get_derivation_matrix");
+
+## TODO Saved function from the console
+def get_sequences(a,b):
+ M = [[tensor_power(get_vector_prime(n+1),i+1) for i in range(a)] for n in range(b)]
+ return [[[list(M[i][j]).index(k^(j+1))+1 for k in [1] + [el for el in primes_first_n(i)]] for j in range(5)] for i in range(5)];
+
__CACHED_INIT_POLY = {};
def _get_initial_poly(poly, dic, m):
if(var in dic[key]):
func = key;
offset = dic[key].index(var);
+ break;
init = lambda n : func.getInitialValue(offset+n);
projects / ajpastor / diff_defined_functions.git / commitdiff