From 36ca13b9ec9c184040e11b6b1360d106cfe0f986 Mon Sep 17 00:00:00 2001 From: Antonio Jimenez Pastor Date: Mon, 25 Jun 2018 16:58:57 +0200 Subject: [PATCH] Modified the method to decide parents for the example functions Started to add the method to compute polynomials directly (Not finished) --- ajpastor/dd_functions/ddExamples.py | 19 +- ajpastor/dd_functions/ddFunction.py | 662 +++++++++++++++++++++--------------- ajpastor/operator/oreOperator.py | 9 + 3 files changed, 415 insertions(+), 275 deletions(-) diff --git a/ajpastor/dd_functions/ddExamples.py b/ajpastor/dd_functions/ddExamples.py index 808cebd..69d6cdc 100644 --- a/ajpastor/dd_functions/ddExamples.py +++ b/ajpastor/dd_functions/ddExamples.py @@ -2,6 +2,10 @@ # This file was *autogenerated* from the file ./ddExamples.sage from sage.all_cmdline import * # import sage library + +from sage.rings.polynomial.polynomial_ring import is_PolynomialRing; +from sage.rings.polynomial.multi_polynomial_ring import is_MPolynomialRing; + _sage_const_3 = Integer(3); _sage_const_2 = Integer(2); _sage_const_1 = Integer(1); _sage_const_0 = Integer(0); _sage_const_6 = Integer(6) from ajpastor.dd_functions.ddFunction import *; @@ -566,12 +570,9 @@ def DAlgebraic(polynomial, init=[], dR=None): mon *= y; ## Building the derivation matrix of <1,y,y^2,...> - print rows M = Matrix(F, rows).transpose(); - print M ## Creating the vector representing y y_vector = vector(F, [_sage_const_0 ,_sage_const_1 ] + [_sage_const_0 ]*(polynomial.degree()-_sage_const_2 )); - print y_vector; ## Building ans solving the system to_solve = move(M, y_vector, lambda p : p.derivative(), M.ncols()+_sage_const_1 ); v = to_solve.right_kernel_matrix()[_sage_const_0 ]; @@ -624,20 +625,26 @@ def Federschwinger(a,b,c,k,init=(_sage_const_0 ,_sage_const_0 )): ### ################################################################################## ################################################################################## -def __decide_parent(input, parent = None, depth = _sage_const_1 ): +def __decide_parent(input, parent = None, depth = 1): if(parent is None): R = input.parent(); if(isinstance(R, sage.symbolic.ring.SymbolicRing)): parameters = set([str(el) for el in input.variables()])-set(['x']); - if(len(parameters) > _sage_const_0 ): + if(len(parameters) > 0 ): parent = ParametrizedDDRing(DDRing(DFinite, depth=depth), parameters); else: parent = DDRing(PolynomialRing(QQ,x), depth=depth); + elif(is_MPolynomialRing(R) or is_PolynomialRing(R)): + parameters = [str(gen) for gen in R.gens()[1:]]; + if(len(parameters) > 0): + parent = ParametrizedDDRing(DDRing(PolynomialRing(QQ,R.gens()[0]), depth=depth), parameters); + else: + parent = DDRing(PolynomialRing(QQ,R.gens()[0]), depth = depth); else: try: parent = DDRing(R, depth = depth); except Exception: - raise TypeError("The object provided is not in a valid Parent"); + raise TypeError("The object provided is not in a valid Parent", e); return parent.base()(input), parent; diff --git a/ajpastor/dd_functions/ddFunction.py b/ajpastor/dd_functions/ddFunction.py index 8642901..97f4235 100644 --- a/ajpastor/dd_functions/ddFunction.py +++ b/ajpastor/dd_functions/ddFunction.py @@ -1,17 +1,19 @@ -# This file was *autogenerated* from the file ./ddFunction.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_5 = Integer(5); _sage_const_100 = Integer(100); _sage_const_12 = Integer(12); _sage_const_11 = Integer(11); _sage_const_1en10 = RealNumber('1e-10') +_sage_const_1en10 = RealNumber('1e-10'); # Python imports import warnings; #SAGE imports from sage.rings.ring import IntegralDomain; -from sage.rings.polynomial.polynomial_ring import is_PolynomialRing as isPolynomial; +from sage.rings.polynomial.polynomial_element import is_Polynomial; +from sage.rings.polynomial.multi_polynomial import is_MPolynomial; +from sage.rings.polynomial.polynomial_ring import is_PolynomialRing; from sage.structure.element import IntegralDomainElement; from sage.categories.integral_domains import IntegralDomains; +from sage.categories.pushout import pushout; from sage.categories.pushout import ConstructionFunctor; #ajpastor imports @@ -26,6 +28,12 @@ from ajpastor.misc.ring_w_sequence import Ring_w_Sequence; from ajpastor.misc.ring_w_sequence import Wrap_w_Sequence_Ring; from ajpastor.misc.ring_w_sequence import sequence; +##################################################### +### Definition of some exceptions +##################################################### +class DevelopementError(NotImplementedError): + def __init__(self, name): + NotImplementedError.__init__(self, "Method %s still under construction. Please be patient and wait until it is finished." %name); ##################################################### ### Definition of the particular warnings we are interested to raise @@ -51,7 +59,7 @@ def _aux_derivation(p): R = p.parent(); return derivative(p,R(x)); except Exception: - return _sage_const_0 ; + return 0 ; ##################################################### ### Class for DD-Rings @@ -60,7 +68,7 @@ class DDRing (Ring_w_Sequence, IntegralDomain): Element = None; _Default_Base = PolynomialRing(QQ,x); - _Default_Depth = _sage_const_1 ; + _Default_Depth = 1 ; _Default_Base_Field = None; _Default_Invertibility = None; _Default_Derivation = None; @@ -93,14 +101,14 @@ class DDRing (Ring_w_Sequence, IntegralDomain): from sage.categories.pushout import AlgebraicExtensionFunctor as algebraic; const_S = S.construction(); const_R = R.construction(); - if(not(isinstance(const_S[_sage_const_0 ],algebraic) and isinstance(const_R[_sage_const_0 ], algebraic))): + if(not(isinstance(const_S[0 ],algebraic) and isinstance(const_R[0 ], algebraic))): return R==S; - if(const_S[_sage_const_1 ] != const_R[_sage_const_1 ]): + if(const_S[1 ] != const_R[1 ]): return False; - polys_R = const_R[_sage_const_0 ].polys; polys_S = const_S[_sage_const_0 ].polys; - names_R = const_R[_sage_const_0 ].names; names_S = const_S[_sage_const_0 ].names; + polys_R = const_R[0 ].polys; polys_S = const_S[0 ].polys; + names_R = const_R[0 ].names; names_S = const_S[0 ].names; if(len(polys_R) != len(polys_S)): return False; @@ -127,12 +135,12 @@ class DDRing (Ring_w_Sequence, IntegralDomain): res = {'algebraic' : [], 'polynomial': [], 'other' : []}; - while((not(const is None)) and (not (_sage_const_1 in gens))): - if(isinstance(current, DDRing) or isinstance(const[_sage_const_0 ], polynomial) or isinstance(const[_sage_const_0 ], multi_polynomial)): + while((not(const is None)) and (not (1 in gens))): + if(isinstance(current, DDRing) or isinstance(const[0 ], polynomial) or isinstance(const[0 ], multi_polynomial)): res['polynomial'] += list(current.gens()); - elif(isinstance(const[_sage_const_0 ], algebraic)): - for i in range(len(const[_sage_const_0 ].polys)): - name = const[_sage_const_0 ].names[i]; + elif(isinstance(const[0 ], algebraic)): + for i in range(len(const[0 ].polys)): + name = const[0 ].names[i]; found = None; for gen in current.gens(): if(current(name) == gen): @@ -140,23 +148,23 @@ class DDRing (Ring_w_Sequence, IntegralDomain): break; if(found is None): raise TypeError("Impossible error: no generator for a name!!"); - res['algebraic'] += [(found, const[_sage_const_0 ].polys[i], const[_sage_const_1 ])]; - elif(not isinstance(const[_sage_const_0 ], fraction)): + res['algebraic'] += [(found, const[0 ].polys[i], const[1 ])]; + elif(not isinstance(const[0 ], fraction)): res['other'] += list(current.gens()); - current = const[_sage_const_1 ]; + current = const[1 ]; const = current.construction(); ## Cleaning "1" from the result ## Put everything as tuples - if(_sage_const_1 in res['algebraic']): + if(1 in res['algebraic']): raise TypeError("1 is a generator of an algebraic extension"); res['algebraic'] = tuple(res['algebraic']); - while(_sage_const_1 in res['polynomial']): - res['polynomial'].remove(_sage_const_1 ); + while(1 in res['polynomial']): + res['polynomial'].remove(1 ); res['polynomial'] = tuple(set(res['polynomial'])); - while(_sage_const_1 in res['other']): - res['other'].remove(_sage_const_1 ); + while(1 in res['other']): + res['other'].remove(1 ); res['other'] = tuple(set(res['other'])); return res, current; @@ -169,26 +177,26 @@ class DDRing (Ring_w_Sequence, IntegralDomain): final_args = []; ## We first compile the parameters the user send - current = _sage_const_0 ; - if(len(args) != _sage_const_1 or args[_sage_const_0 ] != None): + current = 0 ; + if(len(args) != 1 or args[0 ] != None): for i in range(len(args)): - final_args += [[DDRing._Default_Parameters[i][_sage_const_0 ], args[i]]]; + final_args += [[DDRing._Default_Parameters[i][0 ], args[i]]]; current = len(args); for i in range(current, len(DDRing._Default_Parameters)): - final_args += [[DDRing._Default_Parameters[i][_sage_const_0 ], kwds.get(DDRing._Default_Parameters[i][_sage_const_0 ],DDRing._Default_Parameters[i][_sage_const_1 ])]]; + final_args += [[DDRing._Default_Parameters[i][0 ], kwds.get(DDRing._Default_Parameters[i][0 ],DDRing._Default_Parameters[i][1 ])]]; ## Managing the depth over DDRings - if(isinstance(final_args[_sage_const_0 ][_sage_const_1 ], DDRing)): - final_args[_sage_const_1 ][_sage_const_1 ] = final_args[_sage_const_1 ][_sage_const_1 ]+final_args[_sage_const_0 ][_sage_const_1 ].depth(); - final_args[_sage_const_0 ][_sage_const_1 ] = final_args[_sage_const_0 ][_sage_const_1 ]._DDRing__original; + if(isinstance(final_args[0 ][1 ], DDRing)): + final_args[1 ][1 ] = final_args[1 ][1 ]+final_args[0 ][1 ].depth(); + final_args[0 ][1 ] = final_args[0 ][1 ]._DDRing__original; ## Casting to tuple (so is hashable) - to_hash = tuple(tuple(el) for el in final_args[:_sage_const_2 ]); + to_hash = tuple(tuple(el) for el in final_args[:2 ]); final_args = tuple(tuple(el) for el in final_args); - if(final_args[_sage_const_1 ][_sage_const_1 ] == _sage_const_0 ): - return final_args[_sage_const_0 ][_sage_const_1 ]; + if(final_args[1 ][1 ] == 0 ): + return final_args[0 ][1 ]; if(not cls in DDRing._CACHED_DD_RINGS): DDRing._CACHED_DD_RINGS[cls] = {}; @@ -211,8 +219,8 @@ class DDRing (Ring_w_Sequence, IntegralDomain): ## Other attributes initialization self.__variables = None; - if(depth > _sage_const_1 ): - DDRing.__init__(self,DDRing(base, depth-_sage_const_1 , base_field, invertibility, derivation, default_operator), _sage_const_1 , base_field, invertibility, derivation, default_operator); + if(depth > 1 ): + DDRing.__init__(self,DDRing(base, depth-1 , base_field, invertibility, derivation, default_operator), 1 , base_field, invertibility, derivation, default_operator); else: ### We call the builders of the superclasses Ring_w_Sequence.__init__(self, base, method = lambda p,n : self(p).getSequenceElement(n)); @@ -234,7 +242,7 @@ class DDRing (Ring_w_Sequence, IntegralDomain): # If the base ring is already a DDRing, we take the correspond ring from base if isinstance(base, DDRing): self.base_field = base.base_field; - self.__depth = base.__depth + _sage_const_1 ; + self.__depth = base.__depth + 1 ; self.__original = base.__original; # Otherwise, we set this simplest ring as our current coefficient ring else: @@ -242,8 +250,8 @@ class DDRing (Ring_w_Sequence, IntegralDomain): from sage.categories.pushout import PolynomialFunctor, FractionField; current = base; const = base.construction(); - while((not (const is None)) and (isinstance(const[_sage_const_0 ], PolynomialFunctor) or isinstance(const[_sage_const_0 ],FractionField))): - current = const[_sage_const_1 ]; + while((not (const is None)) and (isinstance(const[0 ], PolynomialFunctor) or isinstance(const[0 ],FractionField))): + current = const[1 ]; const = current.construction(); if(not current.is_field()): @@ -254,29 +262,31 @@ class DDRing (Ring_w_Sequence, IntegralDomain): else: self.base_field = base_field; - self.__depth = _sage_const_1 ; + self.__depth = 1 ; self.__original = base; ### Saving the invertibility criteria if(invertibility is None): try: - self_var = self.variables(True,False)[_sage_const_0 ]; - self.base_inversionCriteria = lambda p : p(**{self_var : _sage_const_0 }) == _sage_const_0 ; + self_var = self.variables(True,False)[0 ]; + self.base_inversionCriteria = lambda p : p(**{self_var : 0 }) == 0 ; except IndexError: - self.base_inversionCriteria = lambda p : self.base()(p)==_sage_const_0 ; + self.base_inversionCriteria = lambda p : self.base()(p)==0 ; else: self.base_inversionCriteria = invertibility; ### Saving the base derivation if(derivation is None): try: - self_var = self.variables(True,False)[_sage_const_0 ]; + self_var = self.variables(True,False)[0 ]; self.base_derivation = lambda p : p.derivative(self.base()(self_var)); except IndexError: - self.base_derivation = lambda p : _sage_const_0 ; + self.base_derivation = lambda p : 0 ; else: self.base_derivation = derivation; + ### Setting the default "get_recurence" method + self.__get_recurrence = None; ### Setting new conversions @@ -319,7 +329,7 @@ class DDRing (Ring_w_Sequence, IntegralDomain): return None; ## Using the 'other' generators - if(len(gens_S['other']) > _sage_const_0 ): + if(len(gens_S['other']) > 0 ): return self.base()._coerce_map_from_(S); ## Comparing the algebraic construction @@ -327,7 +337,7 @@ class DDRing (Ring_w_Sequence, IntegralDomain): looking = gens_S['algebraic'][i]; found = False; for el in gens_self['algebraic']: - if(el[_sage_const_1 ] == looking[_sage_const_1 ] and str(el[_sage_const_0 ]) == str(looking[_sage_const_0 ])): + if(el[1 ] == looking[1 ] and str(el[0 ]) == str(looking[0 ])): found = True; break; if(not found): @@ -370,7 +380,6 @@ class DDRing (Ring_w_Sequence, IntegralDomain): - TypeError will be raised if a problem with the algebraic extension is found. - ''' - from sage.categories.pushout import pushout; from sage.categories.integral_domains import IntegralDomains; if(isinstance(S, sage.symbolic.ring.SymbolicRing)): @@ -383,7 +392,7 @@ class DDRing (Ring_w_Sequence, IntegralDomain): except RuntimeError: return None; - if(len(gens_S['other']) > _sage_const_0 or len(gens_self['other']) > _sage_const_0 ): + if(len(gens_S['other']) > 0 or len(gens_self['other']) > 0 ): raise TypeError("Impossible to compute a pushout: generators not recognized found.\n\t- %s" %(list(gens_S['other']) + list(gens_self['other']))); ########################################## @@ -404,11 +413,11 @@ class DDRing (Ring_w_Sequence, IntegralDomain): F = F.fraction_field(); ## Extending the field with algebraic extensions - polys = {str(el[_sage_const_0 ]):el[_sage_const_1 ] for el in gens_self['algebraic']} + polys = {str(el[0 ]):el[1 ] for el in gens_self['algebraic']} for el in gens_S['algebraic']: - if(polys.get(str(el[_sage_const_0 ]), el[_sage_const_1 ]) != el[_sage_const_1 ]): - raise TypeError("Incompatible names in algebraic extensions:\n\t- %s\n\t- %s" %(el,(el[_sage_const_0 ],polys[str(el[_sage_const_0 ])]))); - polys[str(el[_sage_const_0 ])] = el[_sage_const_1 ]; + if(polys.get(str(el[0 ]), el[1 ]) != el[1 ]): + raise TypeError("Incompatible names in algebraic extensions:\n\t- %s\n\t- %s" %(el,(el[0 ],polys[str(el[0 ])]))); + polys[str(el[0 ])] = el[1 ]; sorted_names = sorted(polys); sorted_polys = [polys[el] for el in sorted_names]; @@ -432,10 +441,10 @@ class DDRing (Ring_w_Sequence, IntegralDomain): ########################################## ## Building the final DDRing ########################################## - if(len(variables) > _sage_const_0 ): + if(len(variables) > 0 ): F = PolynomialRing(F,[str(el) for el in variables]); R = DDRing(F, depth = depth); - if(len(params) > _sage_const_0 ): + if(len(params) > 0 ): R = ParametrizedDDRing(R, params); return R; @@ -446,20 +455,20 @@ class DDRing (Ring_w_Sequence, IntegralDomain): def _element_constructor_(self, *args, **kwds): el_class = self.element_class; - if(len(args) < _sage_const_1 ): + if(len(args) < 1 ): raise ValueError("Impossible to build a lazy element without arguments"); - if(args[_sage_const_0 ] is self): - X = args[_sage_const_1 ]; + if(args[0 ] is self): + X = args[1 ]; else: - X = args[_sage_const_0 ]; + X = args[0 ]; try: if(isinstance(X, DDFunction)): if(X.parent() is self): return X; else: - return self.element([coeff for coeff in X.equation.getCoefficients()], X.getInitialValueList(X.equation.get_jp_fo()+_sage_const_1 ), name=X._DDFunction__name); + return self.element([coeff for coeff in X.equation.getCoefficients()], X.getInitialValueList(X.equation.get_jp_fo()+1 ), name=X._DDFunction__name); else: try: try: @@ -475,10 +484,10 @@ class DDRing (Ring_w_Sequence, IntegralDomain): name = X._DDFunction__name; except: pass; - return el.change_init_values([sequence(X,i)*factorial(i) for i in range(el.equation.get_jp_fo()+_sage_const_1 )], name=name); + return el.change_init_values([sequence(X,i)*factorial(i) for i in range(el.equation.get_jp_fo()+1 )], name=name); except AttributeError: try: - return self.element([_sage_const_1 ],[], self.base()(X), name=str(X)); + return self.element([1 ],[], self.base()(X), name=str(X)); except Exception: return self(str(element)); except TypeError as e: @@ -488,7 +497,7 @@ class DDRing (Ring_w_Sequence, IntegralDomain): return (); def ngens(self): - return _sage_const_0 ; + return 0 ; def construction(self): return (DDRingFunctor(self.depth(), self.base_field), self.__original); @@ -567,7 +576,7 @@ class DDRing (Ring_w_Sequence, IntegralDomain): return self.element(coeffs,init_values); ## Otherwise, we need to know the initial value condition equation = self.element(coeffs).equation; - warnings.warn("Not-simple random element not implemented. Returning zero", DDFunctionWarning, stacklevel=_sage_const_2); + warnings.warn("Not-simple random element not implemented. Returning zero", DDFunctionWarning, stacklevel=2); return self.zero(); @@ -615,7 +624,7 @@ class DDRing (Ring_w_Sequence, IntegralDomain): ''' return self.base_inversionCriteria(x); - def element(self,coefficients=[],init=[],inhomogeneous=_sage_const_0 , check_init=True, name=None): + def element(self,coefficients=[],init=[],inhomogeneous=0 , check_init=True, name=None): ''' Method to create a DDFunction contained in self with some coefficients, inhomogeneous term and initial values. This method just call the builder of DDFunction just puting the base argument as self. ''' @@ -627,23 +636,27 @@ class DDRing (Ring_w_Sequence, IntegralDomain): element = self(element); rx = X; - self_var = str(self.variables(True)[_sage_const_0 ]); + self_var = str(self.variables(True)[0 ]); if((rx is None) and (self_var in input)): rx = input.pop(self_var); if not (rx is None): - if(rx == _sage_const_0 ): - return element.getInitialValue(_sage_const_0 ); + if(rx == 0 ): + return element.getInitialValue(0 ); elif(rx in self.base_field): - return element.numerical_solution(rx,**input)[_sage_const_0 ]; + return element.numerical_solution(rx,**input)[0 ]; try: return element.compose(rx); except ValueError: pass; - elif(len(input) == _sage_const_0 ): + elif(len(input) == 0 ): return element; raise NotImplementedError("Not implemented evaluation of an element of this ring (%s) with the parameters %s and %s" %(self,rx,input)); + def get_recurrence(self, *args, **kwds): + if(self.__get_recurrence is None): + raise NotImplementedError("Recurrence method not implemented for %s" %self); + return self.__get_recurrence(*args, **kwds); def variables(self, as_symbolic=False, fill = True): if(self.__variables is None): @@ -655,11 +668,11 @@ class DDRing (Ring_w_Sequence, IntegralDomain): self.__variables = tuple(set(self.__variables)); if(as_symbolic): - if(len(self.__variables) == _sage_const_0 and fill): + if(len(self.__variables) == 0 and fill): return tuple([var(DDRing._Default_variable)]); return self.__variables; else: - if(len(self.__variables) == _sage_const_0 and fill): + if(len(self.__variables) == 0 and fill): return tuple([self.base()(var(DDRing._Default_variable))]); return tuple(self.base()(el) for el in self.__variables); @@ -684,9 +697,9 @@ class ParametrizedDDRing(DDRing): @staticmethod def __classcall__(cls, *args, **kwds): ## In order to call the __classcall__ of DDRing we treat the arguments received - base_ddRing = args[_sage_const_0 ]; - if(len(args) > _sage_const_1 ): - parameters = args[_sage_const_1 ]; + base_ddRing = args[0 ]; + if(len(args) > 1 ): + parameters = args[1 ]; else: parameters = kwds.get('parameters',None); names = kwds.get('names',None); @@ -704,7 +717,7 @@ class ParametrizedDDRing(DDRing): parameters = [parameters]; else: parameters = list(parameters); - if(len(parameters) == _sage_const_0 ): + if(len(parameters) == 0 ): raise TypeError("The list of variables can not be empty"); for i in range(len(parameters)): @@ -723,21 +736,21 @@ class ParametrizedDDRing(DDRing): base_field = base_ddRing.base_field; constructions = [base_ddRing.construction()]; # Here is the DD-Ring operator - while(constructions[-_sage_const_1 ][_sage_const_1 ] != base_field): - constructions += [constructions[-_sage_const_1 ][_sage_const_1 ].construction()]; + while(constructions[-1 ][1 ] != base_field): + constructions += [constructions[-1 ][1 ].construction()]; new_basic_field = PolynomialRing(base_field, parameters).fraction_field(); current = new_basic_field; - for i in range(_sage_const_1 ,len(constructions)): - current = constructions[-i][_sage_const_0 ](current); + for i in range(1 ,len(constructions)): + current = constructions[-i][0 ](current); - if(constructions[_sage_const_0 ][_sage_const_0 ].depth() > _sage_const_1 ): - base_ring = ParametrizedDDRing(DDRing(base_ddRing.original_ring(),depth=constructions[_sage_const_0 ][_sage_const_0 ].depth()-_sage_const_1 ), parameters); - ring = DDRing.__classcall__(cls, base_ring, _sage_const_1 , base_field = new_basic_field, default_operator = base_ddRing.default_operator); + if(constructions[0 ][0 ].depth() > 1 ): + base_ring = ParametrizedDDRing(DDRing(base_ddRing.original_ring(),depth=constructions[0 ][0 ].depth()-1 ), parameters); + ring = DDRing.__classcall__(cls, base_ring, 1 , base_field = new_basic_field, default_operator = base_ddRing.default_operator); Ring_w_Sequence.__init__(ring, base_ring, method = lambda p,n : ring(p).getSequenceElement(n)); IntegralDomain.__init__(ring, base_ring, category=IntegralDomains()); else: - ring = DDRing.__classcall__(cls, current, depth = constructions[_sage_const_0 ][_sage_const_0 ].depth(), base_field = new_basic_field, default_operator = base_ddRing.default_operator); + ring = DDRing.__classcall__(cls, current, depth = constructions[0 ][0 ].depth(), base_field = new_basic_field, default_operator = base_ddRing.default_operator); ring.__init__(base_ddRing, parameters); return ring; @@ -778,12 +791,12 @@ class ParametrizedDDRing(DDRing): return self.__baseDDRing; def _repr_(self): - res = "(%s" %str(self.parameters()[_sage_const_0 ]); - for i in range(_sage_const_1 ,len(self.parameters())): + res = "(%s" %str(self.parameters()[0 ]); + for i in range(1 ,len(self.parameters())): res += ", " + str(self.parameters()[i]); res += ")"; - if(len(self.parameters()) > _sage_const_1 ): + if(len(self.parameters()) > 1 ): return "%s with parameters %s" %(self.base_ddRing(),res); else: return "%s with parameter %s" %(self.base_ddRing(),res); @@ -814,16 +827,16 @@ class ParametrizedDDRing(DDRing): # Override eval method from DDRing def eval(self, element, X=None, **input): rx = X; - self_var = str(self.variables(True)[_sage_const_0 ]); + self_var = str(self.variables(True)[0 ]); if(X is None and self_var in input): rx = input.pop(self_var); ### If X is not None, then we evaluate the variable of the ring if(not (rx is None)): - if(len(input) > _sage_const_0 ): + if(len(input) > 0 ): return self.eval(element, **input)(**{self_var:rx}); else: return super(ParametrizedDDRing, self).eval(element, rx); - elif(len(input) > _sage_const_0 ): + elif(len(input) > 0 ): ### Otherwise, we evaluate the parameters element = self(element); @@ -837,14 +850,14 @@ class ParametrizedDDRing(DDRing): if(self_var in destiny_parameters): raise TypeError("Parameters can only be evaluated to constants.\n\t- Got: %s" %(input)); - if(len(destiny_parameters) == _sage_const_0 ): + if(len(destiny_parameters) == 0 ): destiny_ring = self.base_ddRing(); else: destiny_ring = ParametrizedDDRing(self.base_ddRing(), tuple(destiny_parameters)); new_equation = destiny_ring.element([el(**input) for el in element.equation.getCoefficients()]).equation; - new_init = [el(**input) for el in element.getInitialValueList(new_equation.get_jp_fo()+_sage_const_1 )]; + new_init = [el(**input) for el in element.getInitialValueList(new_equation.get_jp_fo()+1 )]; new_name = None; if(element._DDFunction__name is not None): new_name = din_m_replace(element._DDFunction__name, {key: str(input[key]) for key in input}, True); @@ -859,7 +872,7 @@ class DDFunction (IntegralDomainElement): ##################################### ### Init and Interface methods ##################################### - def __init__(self, parent, input = _sage_const_0 , init_values = [], inhomogeneous = _sage_const_0 , check_init = True, name = None): + def __init__(self, parent, input = 0 , init_values = [], inhomogeneous = 0 , check_init = True, name = None): ''' Method to initialize a DD-Function insade the DD-Ring given in 'parent' @@ -879,32 +892,32 @@ class DDFunction (IntegralDomainElement): ## init_values -> init ## inhomogeneous -> inhom zero = False; - if((type(input) == list and len(input) == _sage_const_0 ) or input == _sage_const_0 ): + if((type(input) == list and len(input) == 0 ) or input == 0 ): zero = True; - elif(type(input) == list and all(el == _sage_const_0 for el in input)): + elif(type(input) == list and all(el == 0 for el in input)): zero = True; - elif((type(input) == list and len(input) == _sage_const_1 ) or (isinstance(input, Operator) and input.getOrder() == _sage_const_0 )): - zero = (inhomogeneous == _sage_const_0 ); + elif((type(input) == list and len(input) == 1 ) or (isinstance(input, Operator) and input.getOrder() == 0 )): + zero = (inhomogeneous == 0 ); if(zero): ### The equation is zero, we need inhomogeneous equals to zero - if(not inhomogeneous == _sage_const_0 ): + if(not inhomogeneous == 0 ): raise ValueError("Incompatible equation (%s) and inhomogeneous term (%s)" %(input, inhomogeneous)); else: ### Also, all the initial values must be zero for el in init_values: - if (not el == _sage_const_0 ): + if (not el == 0 ): raise ValueError("Incompatible equation (%s) and initial values (%s)" %(input, init_values)); - init = [_sage_const_0 ]; - equation = [_sage_const_0 ,_sage_const_1 ]; - inhom = _sage_const_0 ; + init = [0 ]; + equation = [0 ,1 ]; + inhom = 0 ; else: equation = input; init = [el for el in init_values]; inhom = inhomogeneous; ### Cached elements - self.__pows = {_sage_const_0 :_sage_const_1 , _sage_const_1 :self}; # Powers-cache + self.__pows = {0 :1 , 1 :self}; # Powers-cache self.__derivative = None; # The derivative of a function ### Assigning the differential operator @@ -920,7 +933,7 @@ class DDFunction (IntegralDomainElement): # We cast the inhomogeneous term to an element of self.parent() # If that is not zero, then we compute the new operator and initial values needed # to define the equation. - if(inhom != _sage_const_0 ): + if(inhom != 0 ): ## Getting the basic elements inhom = self.parent()(inhom); field = parent.base_field; @@ -930,11 +943,11 @@ class DDFunction (IntegralDomainElement): ## Get the number of elements to check ## If init is too big, we use all the elements for a better checking - n_init = new_eq.get_jp_fo()+_sage_const_1 ; + n_init = new_eq.get_jp_fo()+1 ; to_check = max(n_init, len(init)); ## Getting the matrix of the current equation - M = Matrix(field, self.equation.get_recursion_matrix(to_check-_sage_const_1 -self.equation.forward_order)); + M = Matrix(field, self.equation.get_recursion_matrix(to_check-1 -self.equation.forward_order)); v = vector(field, inhom.getSequenceList(M.nrows())); ## Solving the system MX = v @@ -957,17 +970,17 @@ class DDFunction (IntegralDomainElement): ## After creating the original operator, we check we can not extract an "x" factor - coeff_gcd = _sage_const_1 ; - if(isPolynomial(self.parent().base())): + coeff_gcd = 1 ; + if(is_PolynomialRing(self.parent().base())): l = []; for el in self.equation.getCoefficients(): l += el.coefficients(x); coeff_gcd = gcd(l); - if(coeff_gcd == _sage_const_0 ): - coeff_gcd = _sage_const_1 ; + if(coeff_gcd == 0 ): + coeff_gcd = 1 ; g = coeff_gcd*gcd(self.equation.getCoefficients()); - if(g != _sage_const_1 and g != _sage_const_0 ): + if(g != 1 and g != 0 ): coeffs = [coeff/g for coeff in self.equation.getCoefficients()]; self.equation = self.__buildOperator(coeffs); @@ -977,8 +990,8 @@ class DDFunction (IntegralDomainElement): self.__calculatedSequence = {}; if(len(init) > self.equation.get_jp_fo()): initSequence = [init[i]/factorial(i) for i in range(len(init))]; - M = self.equation.get_recursion_matrix(len(initSequence)-self.equation.forward_order-_sage_const_1 ); - if(M*vector(initSequence) == _sage_const_0 ): + M = self.equation.get_recursion_matrix(len(initSequence)-self.equation.forward_order-1 ); + if(M*vector(initSequence) == 0 ): for i in range(len(initSequence)): self.__calculatedSequence[i] = initSequence[i]; else: @@ -1032,11 +1045,11 @@ class DDFunction (IntegralDomainElement): @derived_property def ps_order(self): if(self.is_null): - return -_sage_const_1 ; + return -1 ; else: - i = _sage_const_0 ; - while(self.getInitialValue(i) == _sage_const_0 ): - i += _sage_const_1 ; + i = 0 ; + while(self.getInitialValue(i) == 0 ): + i += 1 ; return i; @@ -1054,23 +1067,23 @@ class DDFunction (IntegralDomainElement): else: ## Otherwise, we try to get as usual the value d = self.getOrder(); - i = max(n-d,_sage_const_0 ); + i = max(n-d,0 ); rec = self.equation.get_recursion_row(i); - while(rec[n] == _sage_const_0 and i <= self.equation.jp_value): - i += _sage_const_1 ; + while(rec[n] == 0 and i <= self.equation.jp_value): + i += 1 ; rec = self.equation.get_recursion_row(i); - if(rec[n] == _sage_const_0 ): + if(rec[n] == 0 ): raise TypeError("Impossible to compute recursively the required value"); ## Checking that we only need previous elements - for i in range(n+_sage_const_1 , len(rec)): - if(not (rec[i] == _sage_const_0 )): + for i in range(n+1 , len(rec)): + if(not (rec[i] == 0 )): raise TypeError("Impossible to compute recursively the required value"); ## We do this operation in a loop to avoid computing initial values ## if they are not needed res = self.parent().base_field.zero(); for i in range(n): - if(not (rec[i] == _sage_const_0 )): + if(not (rec[i] == 0 )): res -= rec[i]*(self.getSequenceElement(i)); self.__calculatedSequence[n] = self.parent().base_field(res/rec[n]); @@ -1082,7 +1095,7 @@ class DDFunction (IntegralDomainElement): Extra method that returns the list of the first `n` initial coefficients of the power-serie expansion of the function. If it is not possible to get so many values, the method DO NOT return an error. It returns just a truncated list. ''' ### We check the argument is correct - if n < _sage_const_0 : + if n < 0 : raise ValueError("The argument must be a non-negative integer"); res = []; @@ -1107,7 +1120,7 @@ class DDFunction (IntegralDomainElement): Extra method that returns the list of the first `n` initial values for the function. If it is not possible to get so many values, the method DO NOT return an error. It returns just a truncated list. ''' ### We check the argument is correct - if n < _sage_const_0 : + if n < 0 : raise ValueError("The argument must be a non-negative integer"); res = []; @@ -1126,7 +1139,7 @@ class DDFunction (IntegralDomainElement): else: for coeff in self.equation.getCoefficients(): try: - if(coeff != _sage_const_0 ): + if(coeff != 0 ): to_sum += [coeff.degree()] except Exception: pass; @@ -1157,6 +1170,9 @@ class DDFunction (IntegralDomainElement): def change_name(self, new_name): self.__name = new_name; + + def has_name(self): + return not(self.__name is None); ##################################### ### Arithmetic methods @@ -1177,7 +1193,7 @@ class DDFunction (IntegralDomainElement): ''' Method that compute and return a DD-Function `f` such `f*self == 1`, i.e. this method computes the multiplicative inverse of `self`. ''' - if(self.getInitialValue(_sage_const_0 ) == _sage_const_0 ): + if(self.getInitialValue(0 ) == 0 ): raise ValueError("Can not invert a function with initial value 0 --> That is not a power serie"); coeffs = self.getOperator().getCoefficients(); @@ -1186,13 +1202,13 @@ class DDFunction (IntegralDomainElement): newName = None; if(not(self.__name is None)): newName = DinamicString("1/(_1)",self.__name); - if(self.getOrder() == _sage_const_0 ): + if(self.getOrder() == 0 ): raise ZeroDivisionError("Impossible to invert the zero function"); - elif(self.getOrder() == _sage_const_1 ): - return self.parent().element([-coeffs[_sage_const_0 ],coeffs[_sage_const_1 ]], [_sage_const_1 /self.getInitialValue(_sage_const_0 )], check_init=False, name=newName); + elif(self.getOrder() == 1 ): + return self.parent().element([-coeffs[0 ],coeffs[1 ]], [1 /self.getInitialValue(0 )], check_init=False, name=newName); else: newDDRing = DDRing(self.parent()); - return newDDRing.element([self.derivative(), self], [_sage_const_1 /self.getInitialValue(_sage_const_0 )], check_init=False, name=newName); + return newDDRing.element([self.derivative(), self], [1 /self.getInitialValue(0 )], check_init=False, name=newName); def add(self, other): ''' @@ -1214,7 +1230,7 @@ class DDFunction (IntegralDomainElement): ### Computing the new name newName = None; if(not(self.__name is None) and (not(other.__name is None))): - if(other.__name[_sage_const_0 ] == '-'): + if(other.__name[0 ] == '-'): newName = DinamicString("_1_2", [self.__name, other.__name]); else: newName = DinamicString("_1+_2", [self.__name, other.__name]); @@ -1225,21 +1241,21 @@ class DDFunction (IntegralDomainElement): result = None; if(self.is_constant and other.is_constant): parent = self.parent(); - newOperator = [_sage_const_0 ,_sage_const_1 ]; - newInit = [self(_sage_const_0 )+other(_sage_const_0 )]; + newOperator = [0 ,1 ]; + newInit = [self(0 )+other(0 )]; result = parent.element(newOperator, newInit, check_init = False, name=newName); elif(other.is_constant): parent = self.parent(); - newOperator = parent.element(self.equation, inhomogeneous=other(_sage_const_0 )*self.equation.getCoefficient(_sage_const_0 )).equation; - newInit = [self(_sage_const_0 )+other(_sage_const_0 )] + [self.getInitialValue(i) for i in range(_sage_const_1 ,newOperator.get_jp_fo()+_sage_const_1 )]; + newOperator = parent.element(self.equation, inhomogeneous=other(0 )*self.equation.getCoefficient(0 )).equation; + newInit = [self(0 )+other(0 )] + [self.getInitialValue(i) for i in range(1 ,newOperator.get_jp_fo()+1 )]; result = parent.element(newOperator, newInit, check_init = False, name=newName); - result.change_built("polynomial", (PolynomialRing(self.parent().base_field,'x1')("x1+%s" %other(_sage_const_0 )), {'x1':self})); + result.change_built("polynomial", (PolynomialRing(self.parent().base_field,'x1')("x1+%s" %other(0 )), {'x1':self})); elif(self.is_constant): parent = other.parent(); - newOperator = parent.element(other.equation, inhomogeneous=self(_sage_const_0 )*other.equation.getCoefficient(_sage_const_0 )).equation; - newInit = [self(_sage_const_0 )+other(_sage_const_0 )] + [other.getInitialValue(i) for i in range(_sage_const_1 ,newOperator.get_jp_fo()+_sage_const_1 )]; + newOperator = parent.element(other.equation, inhomogeneous=self(0 )*other.equation.getCoefficient(0 )).equation; + newInit = [self(0 )+other(0 )] + [other.getInitialValue(i) for i in range(1 ,newOperator.get_jp_fo()+1 )]; result = parent.element(newOperator, newInit, check_init = False, name=newName) - result.change_built("polynomial", (PolynomialRing(self.parent().base_field,'x1')("x1+%s" %self(_sage_const_0 )), {'x1':other})); + result.change_built("polynomial", (PolynomialRing(self.parent().base_field,'x1')("x1+%s" %self(0 )), {'x1':other})); return result; ## We build the new operator @@ -1249,7 +1265,7 @@ class DDFunction (IntegralDomainElement): newOperator = self.equation.add_solution(other.equation); ### Getting the needed initial values for the solution - needed_initial = newOperator.get_jp_fo()+_sage_const_1 ; + needed_initial = newOperator.get_jp_fo()+1 ; ### Getting as many initial values as posible until the new order op1Init = self.getInitialValueList(needed_initial); @@ -1286,21 +1302,21 @@ class DDFunction (IntegralDomainElement): result = None; if(self.is_constant and other.is_constant): parent = self.parent(); - newOperator = [_sage_const_0 ,_sage_const_1 ]; - newInit = [self(_sage_const_0 )-other(_sage_const_0 )]; + newOperator = [0 ,1 ]; + newInit = [self(0 )-other(0 )]; result = parent.element(newOperator, newInit, check_init = False, name=newName); elif(other.is_constant): parent = self.parent(); - newOperator = parent.element(self.equation, inhomogeneous=other(_sage_const_0 )*self.equation.getCoefficient(_sage_const_0 )).equation; - newInit = [self(_sage_const_0 )-other(_sage_const_0 )] + [self.getInitialValue(i) for i in range(_sage_const_1 ,newOperator.get_jp_fo()+_sage_const_1 )]; + newOperator = parent.element(self.equation, inhomogeneous=other(0 )*self.equation.getCoefficient(0 )).equation; + newInit = [self(0 )-other(0 )] + [self.getInitialValue(i) for i in range(1 ,newOperator.get_jp_fo()+1 )]; result = parent.element(newOperator, newInit, check_init = False, name=newName); - result.change_built("polynomial", (PolynomialRing(self.parent().base_field,'x1')("x1-%s" %other(_sage_const_0 )), {'x1':self})); + result.change_built("polynomial", (PolynomialRing(self.parent().base_field,'x1')("x1-%s" %other(0 )), {'x1':self})); elif(self.is_constant): parent = other.parent(); - newOperator = parent.element(other.equation, inhomogeneous=self(_sage_const_0 )*other.equation.getCoefficient(_sage_const_0 )).equation; - newInit = [self(_sage_const_0 )-other(_sage_const_0 )] + [-other.getInitialValue(i) for i in range(_sage_const_1 ,newOperator.get_jp_fo()+_sage_const_1 )]; + newOperator = parent.element(other.equation, inhomogeneous=self(0 )*other.equation.getCoefficient(0 )).equation; + newInit = [self(0 )-other(0 )] + [-other.getInitialValue(i) for i in range(1 ,newOperator.get_jp_fo()+1 )]; result = parent.element(newOperator, newInit, check_init = False, name=newName) - result.change_built("polynomial", (PolynomialRing(self.parent().base_field,'x1')("-x1+%s" %self(_sage_const_0 )), {'x1':other})); + result.change_built("polynomial", (PolynomialRing(self.parent().base_field,'x1')("-x1+%s" %self(0 )), {'x1':other})); return result; ## We build the new operator @@ -1310,7 +1326,7 @@ class DDFunction (IntegralDomainElement): newOperator = self.equation.add_solution(other.equation); ### Getting the needed initial values for the solution - needed_initial = newOperator.get_jp_fo()+_sage_const_1 ; + needed_initial = newOperator.get_jp_fo()+1 ; ### Getting as many initial values as posible until the new order op1Init = self.getInitialValueList(needed_initial); @@ -1334,28 +1350,28 @@ class DDFunction (IntegralDomainElement): ''' ### We check some simplifications: if one of the functions is zero, then we can return directly 0 if ((other.is_null) or (self.is_null)): - return _sage_const_0 ; + return 0 ; if(self.is_one): return other; elif(other.is_one): return self; elif(self.is_constant and other.is_constant): - return self.getInitialValue(_sage_const_0 )*other.getInitialValue(_sage_const_0 ); + return self.getInitialValue(0 )*other.getInitialValue(0 ); elif(self.is_constant): - return other.scalar(self.getInitialValue(_sage_const_0 )); + return other.scalar(self.getInitialValue(0 )); elif(other.is_constant): - return self.scalar(other.getInitialValue(_sage_const_0 )); + return self.scalar(other.getInitialValue(0 )); ### We build the new operator newOperator = self.equation.mult_solution(other.equation); ### Getting the needed initial values for the solution - needed_initial = newOperator.get_jp_fo()+_sage_const_1 ; + needed_initial = newOperator.get_jp_fo()+1 ; ### Getting as many initial values as posible until the new order op1Init = self.getInitialValueList(needed_initial); op2Init = other.getInitialValueList(needed_initial); - newInit = [sum([binomial(i,j)*op1Init[j] * op2Init[i-j] for j in range(i+_sage_const_1 )]) for i in range(min(len(op1Init),len(op2Init)))]; + newInit = [sum([binomial(i,j)*op1Init[j] * op2Init[i-j] for j in range(i+1 )]) for i in range(min(len(op1Init),len(op2Init)))]; ### Computing the new name newName = None; @@ -1381,26 +1397,26 @@ class DDFunction (IntegralDomainElement): - ``self.scalar(k) == self.mult(DDFunction.getConstantFunction(k))``. ''' ### We first check if the scalar is zero, because the return is direct. - if(r == _sage_const_0 ): + if(r == 0 ): return self.parent().zero(); ### Otherwise, we check that the element is a constant or not because the algorithm is ### much easier r = self.parent().base()(r); - if(self.parent().base_derivation(r) == _sage_const_0 ): - if(r == _sage_const_1 ): + if(self.parent().base_derivation(r) == 0 ): + if(r == 1 ): return self; - n = self.equation.get_jp_fo()+_sage_const_1 ; + n = self.equation.get_jp_fo()+1 ; coeffs = self.getOperator().getCoefficients(); init = self.getInitialValueList(n); if(isinstance(r, DDFunction)): - r = r.getInitialValue(_sage_const_0 ); + r = r.getInitialValue(0 ); ### Computing the new name newName = None; if(not(self.__name is None)): - if(r == -_sage_const_1 ): + if(r == -1 ): newName = DinamicString("-(_1)" ,self.__name); else: newName = DinamicString("(_1)*(_2)", [repr(r), self.__name]); @@ -1414,20 +1430,20 @@ class DDFunction (IntegralDomainElement): @derived_property def zero_extraction(self): if(self == self.parent().zero()): - return (_sage_const_0 ,self); + return (0 ,self); else: - n = _sage_const_0 ; - while(self.getInitialValue(n) == _sage_const_0 ): - n = n+_sage_const_1 ; + n = 0 ; + while(self.getInitialValue(n) == 0 ): + n = n+1 ; - X = self.parent().variables()[_sage_const_0 ]; - if(n == _sage_const_0 ): - return (_sage_const_0 ,self); + X = self.parent().variables()[0 ]; + if(n == 0 ): + return (0 ,self); else: d = self.getOrder(); coeffs = self.getOperator().getCoefficients(); - newEquation = self.parent().element([sum([coeffs[j+l]*(binomial(j+l, j)*falling_factorial(n,l)*X**(n-l)) for l in range(min(n,d-j)+_sage_const_1 )]) for j in range(d+_sage_const_1 )], check_init = False).equation; - newInit = [factorial(i)*self.getSequenceElement(i+n) for i in range(newEquation.get_jp_fo()+_sage_const_1 )]; + newEquation = self.parent().element([sum([coeffs[j+l]*(binomial(j+l, j)*falling_factorial(n,l)*X**(n-l)) for l in range(min(n,d-j)+1 )]) for j in range(d+1 )], check_init = False).equation; + newInit = [factorial(i)*self.getSequenceElement(i+n) for i in range(newEquation.get_jp_fo()+1 )]; ### Computing the new name newName = None; @@ -1442,13 +1458,13 @@ class DDFunction (IntegralDomainElement): if(self.is_null): return self.parent().zero(); if(other.is_constant): - return self.scalar(_sage_const_1 /other.getInitialValue(_sage_const_0 )); + return self.scalar(1 /other.getInitialValue(0 )); s_ze = self.zero_extraction; o_ze = other.zero_extraction; - if(s_ze[_sage_const_0 ] >= o_ze[_sage_const_0 ]): - result = self.parent()(x)**(s_ze[_sage_const_0 ]-o_ze[_sage_const_0 ])*(s_ze[_sage_const_1 ]*o_ze[_sage_const_1 ].inverse); + if(s_ze[0 ] >= o_ze[0 ]): + result = self.parent()(x)**(s_ze[0 ]-o_ze[0 ])*(s_ze[1 ]*o_ze[1 ].inverse); ### Computing the new name newName=None; @@ -1458,7 +1474,7 @@ class DDFunction (IntegralDomainElement): result.__name = newName; return result; else: - raise ValueError("Impossible perform a division if the x factor of the denominator (%d) is greater than in the numerator (%d)"%(o_ze[_sage_const_0 ],s_ze[_sage_const_0 ])); + raise ValueError("Impossible perform a division if the x factor of the denominator (%d) is greater than in the numerator (%d)"%(o_ze[0 ],s_ze[0 ])); def gcd(self, other): if(other in self.parent()): @@ -1471,12 +1487,12 @@ class DDFunction (IntegralDomainElement): elif(other.is_null): return self; - X = self.parent().variables()[_sage_const_0 ]; + X = self.parent().variables()[0 ]; se = self.zero_extraction; oe = other.zero_extraction; - return self.parent()(X**min(se[_sage_const_0 ],oe[_sage_const_0 ])*gcd(se[_sage_const_1 ].getInitialValue(_sage_const_0 ),oe[_sage_const_1 ].getInitialValue(_sage_const_0 ))); + return self.parent()(X**min(se[0 ],oe[0 ])*gcd(se[1 ].getInitialValue(0 ),oe[1 ].getInitialValue(0 ))); ##################################### ### Differential methods @@ -1491,19 +1507,19 @@ class DDFunction (IntegralDomainElement): if(self.__derivative is None): if(self.is_constant): ### Special case: is a constant - self.__derivative = self.parent()(_sage_const_0 ); + self.__derivative = self.parent()(0 ); else: ### We get the new operator newOperator = self.equation.derivative_solution(); ### We get the required initial values (depending of the order of the next derivative) - initList = self.getInitialValueList(newOperator.get_jp_fo()+_sage_const_2 ); - newInit = [initList[i+_sage_const_1 ] for i in range(min(len(initList)-_sage_const_1 ,newOperator.get_jp_fo()+_sage_const_1 ))]; + initList = self.getInitialValueList(newOperator.get_jp_fo()+2 ); + newInit = [initList[i+1 ] for i in range(min(len(initList)-1 ,newOperator.get_jp_fo()+1 ))]; ### Computing the new name newName = None; if(not(self.__name is None)): - if(self.__name[-_sage_const_1 ] == "'"): + if(self.__name[-1 ] == "'"): newName = DinamicString("_1'", self.__name); else: newName = DinamicString("(_1)'", self.__name); @@ -1515,7 +1531,7 @@ class DDFunction (IntegralDomainElement): return self.__derivative; - def integrate(self, constant=_sage_const_0 ): + def integrate(self, constant=0 ): ''' Method to get a DDFunction `g` that satisfies `D(g) = self` and `g(0) = constant`. @@ -1527,12 +1543,12 @@ class DDFunction (IntegralDomainElement): ### We get the initial values for the integral just adding at the begining of the list the constant value # If not enough initial values can be given, we create the integral with so many initial values we have - newInit = [self.parent().base_field(constant)] + self.getInitialValueList(newOperator.get_jp_fo()+_sage_const_1 ); + newInit = [self.parent().base_field(constant)] + self.getInitialValueList(newOperator.get_jp_fo()+1 ); ### Computing the new name newName = None; if(not(self.__name is None)): - if(constant == _sage_const_0 ): + if(constant == 0 ): newName = DinamicString("int(_1)", self.__name); else: newName = DinamicString("int(_1) + _2", [self.__name, repr(constant)]); @@ -1574,19 +1590,19 @@ class DDFunction (IntegralDomainElement): ## If we have the basic function we return the same element ## Also, if self is a constant, the composition is again the constant - self_var = self.parent().variables(True)[_sage_const_0 ]; + self_var = self.parent().variables(True)[0 ]; if(self_var == other or self.is_constant): return self; ## Checking that 'other' at zero is zero value = None; try: - value = other(**{str(self_var):_sage_const_0 }); + value = other(**{str(self_var):0 }); except Exception as e: - warnings.warn("Be careful my friend!! Evaluation at zero were not possible!!\nElement %s" %other, DDFunctionWarning, stacklevel=_sage_const_2 ); + warnings.warn("Be careful my friend!! Evaluation at zero were not possible!!\nElement %s" %other, DDFunctionWarning, stacklevel=2 ); raise e; - if(value != _sage_const_0 ): - raise ValueError("Can not compose with a power series with non-zero constant term. Obtained: %s" %(other(**{str(self_var):_sage_const_0 }))); + if(value != 0 ): + raise ValueError("Can not compose with a power series with non-zero constant term. Obtained: %s" %(other(**{str(self_var):0 }))); ###################################### ## Computing the destiny ring @@ -1605,8 +1621,8 @@ class DDFunction (IntegralDomainElement): else: destiny_ring = push.to_depth(sp.depth()); - if(destiny_ring.depth() >= _sage_const_3 ): - warnings.warn("You are going too deep (depth=%d). The abyss of hell is closer than this function... This may not finish." %destiny_ring.depth(), TooDeepWarning, stacklevel=_sage_const_2 ); + if(destiny_ring.depth() >= 3 ): + warnings.warn("You are going too deep (depth=%d). The abyss of hell is closer than this function... This may not finish." %destiny_ring.depth(), TooDeepWarning, stacklevel=2 ); ###################################### ## Computing the new operator @@ -1620,16 +1636,16 @@ class DDFunction (IntegralDomainElement): ###################################### ## Computing the new initial values ###################################### - required = new_equation.get_jp_fo()+_sage_const_1 ; + required = new_equation.get_jp_fo()+1 ; ## Getting as many initial values as we can and need init_f = self.getInitialValueList(required); init_g = None; try: init_g = g.getInitialValueList(required); except AttributeError: - init_g = [_sage_const_0 ] + [factorial(n)*new_equation.base().sequence(g,n) for n in range(_sage_const_1 ,required)]; + init_g = [0 ] + [factorial(n)*new_equation.base().sequence(g,n) for n in range(1 ,required)]; ## Computing the new initial values - new_init = [init_f[_sage_const_0 ]]+[sum([init_f[j]*bell_polynomial(i,j)(*init_g[_sage_const_1 :i-j+_sage_const_2 ]) for j in range(_sage_const_1 ,i+_sage_const_1 )]) for i in range(_sage_const_1 ,min(len(init_f), len(init_g), required))]; ## See Faa di Bruno's formula + new_init = [init_f[0 ]]+[sum([init_f[j]*bell_polynomial(i,j)(*init_g[1 :i-j+2 ]) for j in range(1 ,i+1 )]) for i in range(1 ,min(len(init_f), len(init_g), required))]; ## See Faa di Bruno's formula ###################################### @@ -1654,7 +1670,7 @@ class DDFunction (IntegralDomainElement): The method used is comparing the first initial values of the function and check they are zero. If some initial value can not be computed, then False is returned. ''' - return self.is_fully_defined and all(self.getInitialValue(i) == _sage_const_0 for i in range(self.equation.get_jp_fo()+_sage_const_1 )); + return self.is_fully_defined and all(self.getInitialValue(i) == 0 for i in range(self.equation.get_jp_fo()+1 )); @derived_property def is_one(self): @@ -1663,7 +1679,7 @@ class DDFunction (IntegralDomainElement): The metod used is checking that the element is a constant and its first initial value is one. ''' - return (self.is_constant and self.getInitialValue(_sage_const_0 ) == _sage_const_1 ); + return (self.is_constant and self.getInitialValue(0 ) == 1 ); @derived_property def is_constant(self): @@ -1676,7 +1692,7 @@ class DDFunction (IntegralDomainElement): ## Test zero and if not, check if a constant can be solution try: if(not self.is_null): - return self.equation[_sage_const_0 ] == _sage_const_0 and all(self.getInitialValue(i) == _sage_const_0 for i in range(_sage_const_1 , self.equation.get_jp_fo()+_sage_const_1 )); + return self.equation[0 ] == 0 and all(self.getInitialValue(i) == 0 for i in range(1 , self.equation.get_jp_fo()+1 )); return True; except TypeError: return False; @@ -1700,7 +1716,7 @@ class DDFunction (IntegralDomainElement): This means that, given some initial values and the differential equation, the solution of such problem is unique (True) or not (False) ''' - max_init = self.equation.get_jp_fo()+_sage_const_1 ; + max_init = self.equation.get_jp_fo()+1 ; return len(self.getInitialValueList(max_init)) == max_init; def equals(self,other): ### TO REVIEW @@ -1730,12 +1746,12 @@ class DDFunction (IntegralDomainElement): ## Special case with constants if(other.is_constant): - return self.is_constant and self(_sage_const_0 ) == other(_sage_const_0 ); + return self.is_constant and self(0 ) == other(0 ); elif(self.is_constant): return False; if(not (self.is_fully_defined and other.is_fully_defined)): - return (self.equation == other.equation) and (self.getInitialValueList(self.equation.get_jp_fo()+_sage_const_1 ) == other.getInitialValueList(other.equation.get_jp_fo()+_sage_const_1 )); + return (self.equation == other.equation) and (self.getInitialValueList(self.equation.get_jp_fo()+1 ) == other.getInitialValueList(other.equation.get_jp_fo()+1 )); ### THREE OPTIONS ## 1. Compare with the JP-Values of each function @@ -1757,11 +1773,11 @@ class DDFunction (IntegralDomainElement): ### If an error occur, then other can not be substracted from self, so they can not be equals return False; - def numerical_solution(self, value, delta=_sage_const_1en10 , max_depth=_sage_const_100 ): + def numerical_solution(self, value, delta=_sage_const_1en10 , max_depth=100 ): try: ## We try to delegate the computation to the operator (in case of OreOperators) - if(self.equation.getCoefficients()[-_sage_const_1 ](_sage_const_0 ) != _sage_const_0 ): - sol = self.equation.operator.numerical_solution(self.getSequenceList(self.getOrder()), [_sage_const_0 ,value],delta); + if(self.equation.getCoefficients()[-1 ](0 ) != 0 ): + sol = self.equation.operator.numerical_solution(self.getSequenceList(self.getOrder()), [0 ,value],delta); return sol.center(), sol.diameter(); else: return self.__basic_numerical_solution(value, delta, max_depth); @@ -1771,20 +1787,20 @@ class DDFunction (IntegralDomainElement): return self.__basic_numerical_solution(value, delta, max_depth); def __basic_numerical_solution(self, value, delta,max_depth): - res = _sage_const_0 ; - to_mult = _sage_const_1 ; - to_sum = _sage_const_0 ; - step = _sage_const_0 ; - find = _sage_const_0 ; - while(find < _sage_const_5 and step < max_depth): + res = 0 ; + to_mult = 1 ; + to_sum = 0 ; + step = 0 ; + find = 0 ; + while(find < 5 and step < max_depth): to_sum = self.getSequenceElement(step)*to_mult; res += to_sum; - if(self.getSequenceElement(step) != _sage_const_0 and abs(to_sum) < delta): - find += _sage_const_1 ; - elif(self.getSequenceElement(step) != _sage_const_0 ): - find = _sage_const_0 ; + if(self.getSequenceElement(step) != 0 and abs(to_sum) < delta): + find += 1 ; + elif(self.getSequenceElement(step) != 0 ): + find = 0 ; - step += _sage_const_1 ; + step += 1 ; to_mult *= value; return float(res),float(abs(to_sum)); @@ -1810,7 +1826,7 @@ class DDFunction (IntegralDomainElement): if(not isinstance(other, DDFunction)): other = self.parent()(other); - m = self.equation.get_jp_fo()+other.equation.get_jp_fo()+_sage_const_1 ; + m = self.equation.get_jp_fo()+other.equation.get_jp_fo()+1 ; return self.getInitialValueList(m) == other.getInitialValueList(m); except Exception: return False @@ -1909,9 +1925,9 @@ class DDFunction (IntegralDomainElement): self.__pows[other] = f; elif(other in ZZ): ## Trying integer power other = int(other); - if(other >= _sage_const_0 ): - a = other >> _sage_const_1 ; - b = a+(other&_sage_const_1 ); + if(other >= 0 ): + a = other >> 1 ; + b = a+(other&1 ); self.__pows[other] = self.__pow__(a)*self.__pow__(b); newName = None; if(not(self.__name is None)): @@ -1928,14 +1944,14 @@ class DDFunction (IntegralDomainElement): try: newDDRing = DDRing(self.parent()); other = self.parent().base_ring()(other); - self.__pows[other] = newDDRing.element([(-other)*f.derivative(),f], [el**other for el in f.getInitialValueList(_sage_const_1 )], check_init=False); + self.__pows[other] = newDDRing.element([(-other)*f.derivative(),f], [el**other for el in f.getInitialValueList(1 )], check_init=False); newName = None; if(not(self.__name is None)): newName = DinamicString("(_1)^%s" %(other), self.__name); self.__pows[other].__name = newName; except TypeError: - raise TypeError("Impossible to compute (%s)^(%s) within the basic field %s" %(f.getInitialValue(_sage_const_0 ), other, f.parent().base_ring())); + raise TypeError("Impossible to compute (%s)^(%s) within the basic field %s" %(f.getInitialValue(0 ), other, f.parent().base_ring())); except ValueError: raise NotImplementedError("Powering to an element of %s not implemented" %(other.parent())); return self.__pows[other]; @@ -1945,7 +1961,7 @@ class DDFunction (IntegralDomainElement): def __eq__(self,other): if (other is self): return True; - if((type(other) == sage.rings.integer.Integer or type(other) == int) and other == _sage_const_0 ): + if((type(other) == sage.rings.integer.Integer or type(other) == int) and other == 0 ): return self.is_null; return self.equals(other); @@ -1955,8 +1971,8 @@ class DDFunction (IntegralDomainElement): Method that return the value of the function in the point given. As the function as a power-serie may not converge, this method only works if the argument is 0. Further implementation can be done for DDFunctions that we know about the convergence. ''' - if ((not(X is None)) and X == _sage_const_0 ): - return self.getInitialValue(_sage_const_0 ); + if ((not(X is None)) and X == 0 ): + return self.getInitialValue(0 ); return self.parent().eval(self, X, **input); @@ -1984,7 +2000,7 @@ class DDFunction (IntegralDomainElement): ''' Missing method for DDFuntions. ''' - return _sage_const_0 ; + return 0 ; def __str__(self, detail=True): ''' @@ -1997,24 +2013,24 @@ class DDFunction (IntegralDomainElement): res = "%s %s in %s:\n" %(self.__critical_numbers__(),repr(self),self.parent()); else: res = "%s" %(repr(self)); - if(res[_sage_const_0 ] != '('): + if(res[0 ] != '('): res += " in %s" %(self.parent()); res += ":\n"; res += '-------------------------------------------------\n\t'; res += '-- Equation (self = f):\n' - j = _sage_const_0 ; - while ((j < self.getOrder()) and (self.getOperator().getCoefficient(j) == _sage_const_0 )): - j += _sage_const_1 ; + j = 0 ; + while ((j < self.getOrder()) and (self.getOperator().getCoefficient(j) == 0 )): + j += 1 ; res += self.__get_line__(self.getOperator().getCoefficient(j), j, detail, True); - for i in range(j+_sage_const_1 ,self.getOrder()+_sage_const_1 ): - if(self.getOperator().getCoefficient(i) != _sage_const_0 ): + for i in range(j+1 ,self.getOrder()+1 ): + if(self.getOperator().getCoefficient(i) != 0 ): res += '\n'; res += self.__get_line__(self.getOperator().getCoefficient(i), i, detail); res += '\n\t\t = 0 \n\t'; res += '-- Initial values:\n\t' - res += format(self.getInitialValueList(self.equation.get_jp_fo()+_sage_const_1 )); + res += format(self.getInitialValueList(self.equation.get_jp_fo()+1 )); res += '\n-------------------------------------------------'; @@ -2029,12 +2045,12 @@ class DDFunction (IntegralDomainElement): if(isinstance(element, DDFunction) and (not(element.__name is None))): crit = element.__critical_numbers__(); else: - ind = string.find(")")+_sage_const_1 ; + ind = string.find(")")+1 ; crit = string[:ind]; string = string[ind:]; ## Printing critical numbers if detail is required - if(detail and len(crit) > _sage_const_0 ): + if(detail and len(crit) > 0 ): res += crit + "\t"; if(len(res.expandtabs()) < len("\t\t".expandtabs())): res += "\t"; @@ -2042,21 +2058,21 @@ class DDFunction (IntegralDomainElement): res += "\t\t"; ## Adding the arithmetic symbol - if(not(first) and string[_sage_const_0 ] != '-'): + if(not(first) and string[0 ] != '-'): res += '+ '; - elif(string[_sage_const_0 ] == '-'): + elif(string[0 ] == '-'): res += '- '; - if(string[_sage_const_1 ] == '('): - string = string[_sage_const_1 :-_sage_const_1 ]; + if(string[1 ] == '('): + string = string[1 :-1 ]; else: - string = string[_sage_const_1 :]; + string = string[1 :]; else: res += ' '; ## Adding the function (with its derivative if needed) - if(der > _sage_const_1 ): + if(der > 1 ): res += "D^%d(f) " %der; - elif(der == _sage_const_1 ): + elif(der == 1 ): res += " D(f) "; else: res += " f "; @@ -2072,7 +2088,7 @@ class DDFunction (IntegralDomainElement): Representation method for DDFunctions. It prints basic information of the DDFunction. ''' if(self.is_constant): - return str(self.getInitialValue(_sage_const_0 )); + return str(self.getInitialValue(0 )); if(self.__name is None): return "(%s:%s:%s)DD-Function in (%s)" %(self.getOrder(),self.equation.get_jp_fo(),self.size(),self.parent()); else: @@ -2098,7 +2114,7 @@ class DDFunction (IntegralDomainElement): ##################################################### class DDRingFunctor (ConstructionFunctor): def __init__(self, depth, base_field): - self.rank = _sage_const_11 ; + self.rank = 11 ; self.__depth = depth; self.__base_field = base_field; super(DDRingFunctor, self).__init__(IntegralDomains(),IntegralDomains()); @@ -2122,7 +2138,7 @@ class DDRingFunctor (ConstructionFunctor): class ParametrizedDDRingFunctor (DDRingFunctor): def __init__(self, depth, base_field, var = set([])): - self.rank = _sage_const_12 ; + self.rank = 12 ; self.__vars = var; super(ParametrizedDDRingFunctor, self).__init__(depth, base_field); @@ -2156,11 +2172,11 @@ class DDSimpleMorphism (sage.categories.map.Map): def _call_(self, p): if(isinstance(self.codomain(), DDRing)): try: - return self.codomain().element([self.codomain().base()(coeff) for coeff in p.equation.getCoefficients()], p.getInitialValueList(p.equation.get_jp_fo()+_sage_const_1 )); + return self.codomain().element([self.codomain().base()(coeff) for coeff in p.equation.getCoefficients()], p.getInitialValueList(p.equation.get_jp_fo()+1 )); except: raise ValueError("Impossible the coercion of element \n%s\n into %s" %(p, self.codomain())) if(p.is_constant): - return self.codomain()(p.getInitialValue(_sage_const_0 )); + return self.codomain()(p.getInitialValue(0 )); raise ValueError("Impossible perform this coercion: non-constant element"); @@ -2171,22 +2187,116 @@ class DDSimpleMorphism (sage.categories.map.Map): ##################################### DDRing.Element = DDFunction; +##################################### +### Extra public methods +##################################### def zero_extraction(el): try: R = el.parent(); if(isinstance(R, DDRing)): return el.zero_extraction; elif(x in R.gens()): - n = _sage_const_0 ; - val = el(**{'x':_sage_const_0 }); - while(el != R.zero() and val == _sage_const_0 ): - n += _sage_const_1 ; + n = 0 ; + val = el(**{'x':0 }); + while(el != R.zero() and val == 0 ): + n += 1 ; el = R(el/x); - val = el(**{'x':_sage_const_0 }); + val = el(**{'x':0 }); return (n, el); except AttributeError as e: - return (_sage_const_0 ,el); + return (0 ,el); + +def polynomial_computation(poly, functions): + ''' + Method that compute a polynomial operation over a set of DDFunctions. + + INPUT: + - poly: a multivariate polynomial representing the operation we want to perform. + - functions: the set of functions related with the polynomial. + OUTPUT: + A DDFunction in the corresponding DDRing. + + WARNING: + We assume that for any pair of functions f,g in functions, there are not + constants c,d such that f = cg+d. Otherwise the result will still + be correct, but will be computed slower. + ''' + ### Checking the argument 'poly' + if(not (is_Polynomial(poly) or is_MPolynomial(poly))): + raise TypeError("First argument require to be a polynomial.\n\t- Got: %s (%s)" %(poly, type(poly))); + parent_poly = poly.parent(); gens = parent_poly.gens(); + + ### Checking the argument 'functions' + tfunc= type(functions); + if(not ((tfunc is list) or (tfunc is set) or (tfunc is tuple))): + functions = [functions]; + + if(len(functions) < parent_poly.ngens()): + raise TypeError("Not enough functions given for the polynomial ring of the polynomial."); + + ### Coarcing all the functions to a common parent + parent = functions[0].parent(); + for i in range(1,len(functions)): + parent = pushout(parent, functions[i].parent()); + + ### Base case: non-DDRing + if(not isinstance(parent, DDRing)): + return poly(**{str(gens[i]) : functions[i] for i in range(len(gens))}); + + functions = [parent(f) for f in functions]; + + ### Grouping the elements that are derivatives of each other + ### Using the assumption on the input functions, we have that for each + ### chain of derivatives, we only need to check if the function is + ### the derivative of the last or if the first element is the derivative + ### of the function. + chains = []; + for f in functions: + inPlace = False; + for i in range(len(chains)): + if(f.derivative() == chains[i][0]): + chains[i] = [f] + chains[i]; + inPlace = True; + break; + elif(chains[i][-1].derivative() == f): + chains[i] = chains[i] + [f]; + inPlace = True; + break; + if(not inPlace): + chains += [[f]]; + dics = {c[0] : [gens[functions.index(f)] for f in c] for c in chains}; + + ## We build the new operator + newOperator = None; + 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); + if(newOperator is None): + newOperator = _get_equation_poly(poly, dics); + + ### Getting the needed initial values for the solution + 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): + newInit = newInit[:newInit.index(None)]; + + ## Computing the new name + newName = None; + if(all(f.has_name() for f in functions)): + newName = DinamicString(_m_replace(str(poly), {str(gens[i]) : "_%d" %(i+1) for i in range(len(gens))}), [f._DDFunction__name for f in functions]); + + ## Building the element + result = parent.element(newOperator, newInit, check_init=False, name=newName); + result.change_built("polynomial", (poly, {str(gens[i]) : functions[i]})); + + return result; ################################################################################################### ### PRIVAME MODULE METHODS @@ -2210,24 +2320,39 @@ def _m_indices(string, *seps): def _m_replace(string, to_replace, escape=("(",")")): ## Checking if ther are scape characters if(not escape is None): - pos_init = string.find(escape[_sage_const_0 ]); - if(pos_init >= _sage_const_0 ): + pos_init = string.find(escape[0 ]); + if(pos_init >= 0 ): ## Escape initial character found, skipping outside the escape pos_end = len(string); - pos_end = string.rfind(escape[_sage_const_1 ]); + pos_end = string.rfind(escape[1 ]); if(pos_end <= pos_init): pos_end = len(string); - return string[:pos_init+_sage_const_1 ] + _m_replace(string[pos_init+_sage_const_1 :pos_end], to_replace, None) + string[pos_end:]; + return string[:pos_init+1 ] + _m_replace(string[pos_init+1 :pos_end], to_replace, None) + string[pos_end:]; ## No escape characters: doing normal computation all_index = _m_indices(string, *to_replace.keys()); - res = ""; current_index = _sage_const_0 ; + res = ""; current_index = 0 ; for el in all_index: - res += string[current_index:el[_sage_const_0 ]] + to_replace[el[_sage_const_1 ]]; - current_index = el[_sage_const_0 ]+len(el[_sage_const_1 ]); + res += string[current_index:el[0 ]] + to_replace[el[1 ]]; + current_index = el[0 ]+len(el[1 ]); res += string[current_index:] return res; + +def _get_equation_poly(poly, dic): + raise DevelopementError("_get_equation_poly"); + +def _get_initial_poly(poly, dic, i): + ## 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))); + + ## More than one variable + if(len(poly.variables()) > 1): + raise DevelopementError("_get_initial_poly"); + + raise DevelopementError("_get_initial_poly"); ################################################################################################### ### COMMAND CONVERSION OF DD_FUNCTIONS @@ -2246,10 +2371,10 @@ def command(e): def _command_list(elements): res = "["; - for i in range(len(elements)-_sage_const_1 ): + for i in range(len(elements)-1 ): res += command(elements[i]) + ", "; - if(len(elements) > _sage_const_0 ): - res += command(elements[-_sage_const_1 ]); + if(len(elements) > 0 ): + res += command(elements[-1 ]); res += "]"; return res; @@ -2267,7 +2392,6 @@ except ImportError: warnings.warn("Package ore_algebra was not found. It is hightly recommended to get this package built by M. Kauers and M. Mezzaroba (http://marc.mezzarobba.net/code/ore_algebra-analytic/ or http://www.kauers.de/software.html)", NoOreAlgebraWarning); DFinite = DDRing(PolynomialRing(QQ,x)); DDFinite = DDRing(DFinite); -#CFinite = DDRing(QQ, default_operator=w_OreOperator); DFiniteP = ParametrizedDDRing(DFinite, [var('P')]); def __get_recurrence(f): @@ -2282,29 +2406,29 @@ def __get_recurrence(f): num_of_bck = bound - op.getOrder(); m = None; - for i in range(num_of_bck, _sage_const_0 , -_sage_const_1 ): - if(op.backward(i) != _sage_const_0 ): + for i in range(num_of_bck, 0 , -1 ): + if(op.backward(i) != 0 ): m = -i; break; if(m is None): - for i in range(op.getOrder()+_sage_const_1 ): - if(op.forward(i) != _sage_const_0 ): + for i in range(op.getOrder()+1 ): + if(op.forward(i) != 0 ): m = i; break; - n = op._Operator__polynomialRing.gens()[_sage_const_0 ]; + n = op._Operator__polynomialRing.gens()[0 ]; - rec = [op.get_recursion_polynomial(i)(n=n+m) for i in range(m,op.forward_order + _sage_const_1 )]; + rec = [op.get_recursion_polynomial(i)(n=n+m) for i in range(m,op.forward_order + 1 )]; rec_gcd = gcd(rec); rec = [el/rec_gcd for el in rec]; return (rec, f.getSequenceList(op.forward_order-m)); -DFinite.get_recurrence = __get_recurrence; +DFinite._DDRing__get_recurrence = __get_recurrence; ################################################################################################### ### PACKAGE ENVIRONMENT VARIABLES ################################################################################################### -__all__ = ["DDRing", "DFinite", "DDFinite", "command", "zero_extraction" , "ParametrizedDDRing", "DFiniteP"]; +__all__ = ["DDRing", "DFinite", "DDFinite", "command", "zero_extraction", "polynomial_computation", "ParametrizedDDRing", "DFiniteP"]; diff --git a/ajpastor/operator/oreOperator.py b/ajpastor/operator/oreOperator.py index f65051d..16cc5f4 100644 --- a/ajpastor/operator/oreOperator.py +++ b/ajpastor/operator/oreOperator.py @@ -207,4 +207,13 @@ class w_OreOperator(Operator): ####################################################### + ####################################################### + ### OTHER WRAPPED METHODS + ####################################################### + def annihilator_of_associate(self, M): + M = w_OreOperator(self.base(), M).operator; + return w_OreOperator(self.base(), self.operator.annihilator_of_associate(M)); + + def annihilator_of_polynomial(self, poly): + return w_OreOperator(self.base(), self.operator.annihilator_of_polynomial(poly)); -- 2.1.4