-# 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
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
R = p.parent();
return derivative(p,R(x));
except Exception:
- return _sage_const_0 ;
+ return 0 ;
#####################################################
### Class for DD-Rings
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;
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;
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):
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;
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] = {};
## 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));
# 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:
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()):
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
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
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):
- 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)):
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'])));
##########################################
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];
##########################################
## 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;
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:
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:
return ();
def ngens(self):
- return _sage_const_0 ;
+ return 0 ;
def construction(self):
return (DDRingFunctor(self.depth(), self.base_field), self.__original);
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();
'''
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.
'''
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):
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);
@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);
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)):
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;
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);
# 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);
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);
#####################################
### 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'
## 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
# 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;
## 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
## 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);
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:
@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;
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]);
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 = [];
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 = [];
else:
for coeff in self.equation.getCoefficients():
try:
- if(coeff != _sage_const_0 ):
+ if(coeff != 0 ):
to_sum += [coeff.degree()]
except Exception:
pass;
def change_name(self, new_name):
self.__name = new_name;
+
+ def has_name(self):
+ return not(self.__name is None);
#####################################
### Arithmetic methods
'''
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();
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):
'''
### 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]);
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
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);
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
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);
'''
### 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;
- ``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]);
@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;
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;
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()):
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
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);
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`.
### 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)]);
## 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
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
######################################
## 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
######################################
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):
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):
## 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;
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
## 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
### 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);
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));
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
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)):
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];
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);
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);
'''
Missing method for DDFuntions.
'''
- return _sage_const_0 ;
+ return 0 ;
def __str__(self, detail=True):
'''
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-------------------------------------------------';
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";
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 ";
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:
#####################################################
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());
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);
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");
#####################################
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
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
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;
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):
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"];
projects / ajpastor / diff_defined_functions.git / commitdiff