From 76b4ce2a9d466bb4491ddeafca1a6c8dd3b405a9 Mon Sep 17 00:00:00 2001 From: Antonio Jimenez Pastor Date: Mon, 25 Jun 2018 20:39:43 +0200 Subject: [PATCH] The initial values for the polynomial computation is finished. It is based in a recursive structure using the generalized Leibniz rule for decomposing the products in the monomials, and Bruno di Faa's formula for the power of variables. --- ajpastor/dd_functions/ddFunction.py | 91 +++++++++++++++++++++++++++++++------ 1 file changed, 76 insertions(+), 15 deletions(-) diff --git a/ajpastor/dd_functions/ddFunction.py b/ajpastor/dd_functions/ddFunction.py index 97f4235..14e13db 100644 --- a/ajpastor/dd_functions/ddFunction.py +++ b/ajpastor/dd_functions/ddFunction.py @@ -816,10 +816,7 @@ class ParametrizedDDRing(DDRing): def ngens(self): return len(self.parameters()); - - #def variables(self, as_symbolic=False): - # return self.base_ddRing().variables(as_symbolic); - + # Override to_depth method from DDRing def to_depth(self, dest): return ParametrizedDDRing(self.base_ddRing().to_depth(dest), self.parameters(True)); @@ -1981,7 +1978,8 @@ class DDFunction (IntegralDomainElement): ''' Hash method for DDFunctions. ''' - return int("%d%d%d" %(self.getOrder(),self.equation.get_jp_fo(),self.size())); + return sum(self.getInitialValueList(20)); + #return int("%d%d%d" %(self.getOrder(),self.equation.get_jp_fo(),self.size())); #return sum([hash(coeff) for coeff in self.getOperator().getCoefficients()]); def __getitem__(self, key): @@ -2272,9 +2270,9 @@ def polynomial_computation(poly, functions): if(len(dics) == 1): if(hasattr(f.equation, "annihilator_of_polynomial")): ## Reordering the variables - parent_poly = PolynomialRing(chains[0][0].parent().original_ring(), dics[chains[0][0]]); - poly = parent_poly(poly); - newOperator = f.equation.annihilator_of_polynomial(poly); + parent_poly2 = PolynomialRing(chains[0][0].parent().original_ring(), dics[chains[0][0]]); + poly2 = parent_poly2(poly); + newOperator = f.equation.annihilator_of_polynomial(poly2); if(newOperator is None): newOperator = _get_equation_poly(poly, dics); @@ -2282,6 +2280,7 @@ def polynomial_computation(poly, functions): needed_initial = newOperator.get_jp_fo()+1 ; ### Getting as many initial values as posible until the new order + newInit = [_get_initial_poly(poly, dics, i) for i in range(needed_initial)]; ## If some error happens, we delete all results after it if(None in newInit): @@ -2294,7 +2293,7 @@ def polynomial_computation(poly, functions): ## Building the element result = parent.element(newOperator, newInit, check_init=False, name=newName); - result.change_built("polynomial", (poly, {str(gens[i]) : functions[i]})); + result.change_built("polynomial", (poly, {str(gens[i]) : functions[i] for i in range(len(gens))})); return result; @@ -2342,17 +2341,79 @@ def _m_replace(string, to_replace, escape=("(",")")): def _get_equation_poly(poly, dic): raise DevelopementError("_get_equation_poly"); -def _get_initial_poly(poly, dic, i): +__CACHED_INIT_POLY = {}; +def _get_initial_poly(poly, dic, m): + ''' + Method that computes the ith initial value of the polynomial + given when evaluated using the dictionary dic. + ''' + ## Constant case + if(poly.is_constant()): + #print "Calling (%s +++ %s +++ %d)" %(poly, dic, m); + if(m > 0): + return 0; + return poly.parent().base()(poly); + + ## Non-constant polynomial + ## Checking the Cache + global __CACHED_INIT_POLY; + if(not poly in __CACHED_INIT_POLY): + __CACHED_INIT_POLY[poly] = []; + + found = False; + for el in __CACHED_INIT_POLY[poly]: + if(el[0] == dic): + found = True; + if(m in el[1]): + #print "Calling (%s +++ %s +++ %d) -- CACHED" %(poly, dic, m); + return el[1][m]; + break; + if(not found): + __CACHED_INIT_POLY[poly] += [(dic, {})]; + #print "Calling (%s +++ %s +++ %d)" %(poly, dic, m); + + result = None; ## General case (poly is not a monomial) if(not (poly.is_monomial())): - c = poly.coefficients(); m = poly.monomials(); - return sum(c[i]*_get_equation_poly(m[i],dic,i) for i in range(len(m))); + c = poly.coefficients(); mo = poly.monomials(); + result = sum(c[j]*_get_initial_poly(mo[j],dic,m) for j in range(len(mo))); ## More than one variable - if(len(poly.variables()) > 1): - raise DevelopementError("_get_initial_poly"); + elif(len(poly.variables()) > 1): + var = poly.variables()[0]; + deg = poly.degree(var); + oth = poly.parent()(poly/(var**deg)); - raise DevelopementError("_get_initial_poly"); + result = sum(binomial(m,k)*_get_initial_poly(var**deg, dic, k)*_get_initial_poly(oth, dic, m-k) for k in range(m+1)); + + ## 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); + + init = lambda n : func.getInitialValue(offset+n); + + ## Computing the final result + if(m > 0): + result = sum(falling_factorial(deg,k)*(init(0)**(deg-k))*bell_polynomial(m,k)(*[init(j+1) for j in range(m-k+1)]) for k in range(min(deg,m)+1)); + else: + result = init(0)**deg; + + if(not(result is None)): + ## Saving the result + for el in __CACHED_INIT_POLY[poly]: + if(el[0] == dic): + el[1][m] = result; + break; + return result; + else: + raise ValueError("An impossible point was reached"); ################################################################################################### ### COMMAND CONVERSION OF DD_FUNCTIONS -- 2.1.4