##################################################################################
##################################################################################
###
+### Algebraic functions
+###
+##################################################################################
+##################################################################################
+def DAlgebraic(polynomial, init=[], dR=None):
+ '''
+ Method that transform an algebraic function to a DD-Function.
+
+ INPUT:
+ - polynomial: the minimal polynomial of the function we want to transform.
+ - init: the initial values that the function will have. Two options are
+ possible: a list is given, then we will use it directly to build the final
+ result, or a value is given an we will compute the others using the polynomial
+ equation.
+ - dR: the ring where we want to include the result. If None, an automatic
+ destiny ring will be computed.
+
+ OUTPUT:
+ - A DDFunction in a particuar DDRing.
+
+ WARNINGS:
+ - There is no control about the irreducibility of the polynomial.
+
+ ERRORS:
+ - If the function can not be represented in dR a TypeError will be raised
+ - If any error happens with the initial values, a ValueError will be raised
+ '''
+ ## Local imports
+ from sage.rings.polynomial.polynomial_ring import is_PolynomialRing as isPolynomial;
+ from sage.rings.polynomial.multi_polynomial_ring import is_MPolynomialRing as isMPolynomial;
+ from sage.categories.pushout import FractionField;
+ from ajpastor.misc.matrix import matrix_of_dMovement as move;
+
+ ###############################################
+ ## Dealing with the polynomial input
+ ###############################################
+ parent = polynomial.parent();
+ if(not (isPolynomial(parent) or isMPolynomial(parent))):
+ raise TypeError("The minimal polynomial is NOT a polynomial");
+
+ base_ring = None;
+ F = None;
+ poly_ring = parent;
+ if(isMPolynomial(parent)):
+ base_ring = PolynomialRing(parent.base(),parent.gens()[:-_sage_const_1 ]);
+ poly_ring = PolynomialRing(base_ring.fraction_field(), parent.gens()[-_sage_const_1 ]);
+ else:
+ if(isinstance(parent.base().construction()[_sage_const_0 ], FractionField)):
+ base_ring = parent.base().base();
+ else:
+ base_ring = parent.base();
+ if(not parent.base().is_field()):
+ poly_ring = PolynomialRing(parent.base().fraction_field(), parent.gens()[-_sage_const_1 ]);
+
+ F = poly_ring.base();
+ y = poly_ring.gens()[-_sage_const_1 ];
+ ## At this point we have the following
+ ## - F is a field
+ ## - y is a variable
+ ## - poly_ring == F[y]
+ polynomial = poly_ring(polynomial); ## Now the structure is univariate
+ if(polynomial.degree() == _sage_const_1 ):
+ return -polynomial[_sage_const_0 ]/polynomial[_sage_const_1 ];
+ elif(polynomial.degree() <= _sage_const_0 ):
+ raise TypeError("Constant polynomial given for algebraic function: IMPOSSIBLE!!");
+
+ #################################################
+ ## Building and checking the destiny ring
+ #################################################
+ destiny_ring = None;
+ if(dR is None):
+ destiny_ring = DDRing(base_ring);
+ else:
+ destiny_ring = dR;
+ coercion = destiny_ring._coerce_map_from_(base_ring);
+ if((coercion is None) or (coercion is False)):
+ raise TypeError("Incompatible polynomial with destiny ring:\n\t- Coefficients in: %s\n\t- Destiny Ring: %s" %(base_ring, destiny_ring));
+
+ dest_var = repr(destiny_ring.variables()[_sage_const_0 ]);
+
+ ##################################################
+ ## Computing the differential equation
+ ##################################################
+ ## Computing the derivative
+ dy = polynomial.derivative(y);
+
+ ## Getting its gcd with the polynomial
+ g,r,s = polynomial.xgcd(dy);
+ if((g != _sage_const_1 ) or (not(g in poly_ring.base()))):
+ raise ValueError("No irreducible polynomial given");
+
+ ## Computing the coefficient-wise derivation of polynomial
+ mon = poly_ring(_sage_const_1 );
+ ky = poly_ring(_sage_const_0 );
+ for i in range(polynomial.degree()+_sage_const_1 ):
+ ky += (polynomial[i].derivative())*mon;
+ mon *= y;
+
+ ## Getting the polynomial expression of y', y'',..., y^{(deg(polynomial))}
+ rows = [[_sage_const_0 ]*polynomial.degree()];
+ mon = poly_ring(_sage_const_1 );
+ for i in range(polynomial.degree()-_sage_const_1 ):
+ rows += [((-(i+_sage_const_1 )*mon*s*ky)%polynomial).coefficients(False)];
+ 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 ];
+
+ ## Cleaning denominators
+ cleaning = lcm(el.denominator() for el in v);
+
+ equation = destiny_ring.element([F.base()(el*cleaning) for el in v]).equation;
+
+ ##################################################
+ ## Getting the initial values
+ ##################################################
+ if(not (type(init) is list)):
+ ## We try to compute the new initial values
+ init = [init];
+ go_on = True;
+ for i in range(_sage_const_1 ,min(equation.get_jp_fo()+_sage_const_2 , to_solve.ncols())):
+ try:
+ init += [sum(to_solve[j,i](**{dest_var:_sage_const_0 })*init[_sage_const_0 ]**j for j in range(polynomial.degree()))];
+ except ZeroDivisionError:
+ go_on = False;
+ break;
+
+ if(go_on and (equation.get_jp_fo()+_sage_const_2 > to_solve.ncols())):
+ extra = move(M, vector(F,[el[_sage_const_0 ] for el in to_solve[:,-_sage_const_1 ]]), equation.get_jp_fo()+_sage_const_2 -to_solve.ncols());
+ init += [sum(extra[j,i](**{dest_var:_sage_const_0 })*init[_sage_const_0 ]**j for j in range(polynomial.degree())) for i in range(extra.ncols())];
+
+ ##################################################
+ ## Returning the DDFunction
+ ##################################################
+ return destiny_ring.element(equation, init);
+
+##################################################################################
+##################################################################################
+###
### Particular differential Equations
###
##################################################################################
# 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_4 = Integer(4); _sage_const_100 = Integer(100); _sage_const_12 = Integer(12); _sage_const_11 = Integer(11); _sage_const_1en10 = RealNumber('1e-10')
+_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')
# Python imports
import warnings;
### Static methods
#################################################
@staticmethod
+ def __algebraic_equivalent__(R,S):
+ '''
+ Method that checks if two NumberField extensions are equivalent in SAGE.
+
+ We say that R~S if they have the same polynomials defining their elements
+ and they have exactly the same name.
+
+ If R or S are not an algebraic extension then the method returns R==S
+ '''
+ 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))):
+ return R==S;
+
+ if(const_S[_sage_const_1 ] != const_R[_sage_const_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;
+
+ if(len(polys_R) != len(polys_S)):
+ return False;
+
+ for i in range(len(polys_R)):
+ try:
+ index = polys_S.index(polys_R[i]);
+ if(names_R[i] != names_S[index]):
+ return False;
+ except ValueError:
+ return False;
+ return True;
+
+ @staticmethod
def __get_gens__(parent):
+ from sage.categories.pushout import AlgebraicExtensionFunctor as algebraic;
+ from sage.categories.pushout import PolynomialFunctor as polynomial;
+ from sage.categories.pushout import MultiPolynomialFunctor as multi_polynomial;
+ from sage.categories.pushout import FractionField as fraction;
+
current = parent;
const = parent.construction();
gens = [el for el in parent.gens()];
+ 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)):
+ 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];
+ found = None;
+ for gen in current.gens():
+ if(current(name) == gen):
+ found = 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['other'] += list(current.gens());
+
current = const[_sage_const_1 ];
- gens += [el for el in current.gens()];
const = current.construction();
- while(_sage_const_1 in gens):
- gens.remove(_sage_const_1 );
+ ## Cleaning "1" from the result
+ ## Put everything as tuples
+ if(_sage_const_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 );
+ res['polynomial'] = tuple(set(res['polynomial']));
+ while(_sage_const_1 in res['other']):
+ res['other'].remove(_sage_const_1 );
+ res['other'] = tuple(set(res['other']));
- return tuple(set(gens));
+ return res, current;
#################################################
### Builder
- ### Setting new coercions
- current = self;
+ ### Setting new conversions
+ current = self.base();
+ morph = DDSimpleMorphism(self, current);
+ current.register_conversion(morph);
while(not(current.base() == current)):
current = current.base();
morph = DDSimpleMorphism(self, current);
current.register_conversion(morph);
-
-
-
+
#################################################
### Coercion methods
#################################################
def _coerce_map_from_(self, S):
+ '''
+ Method to get the coerce map from the SAGE structure 'S' (if possible).
+
+ To allow the agebraic numbers, we use the method '__get_gens__' to compare how the ring 'S' and the ring 'self' where built. If at some point we can not use the behaviour of the generators, we will rely on the usual _coerce_map_from_ with 'self.base()'.
+ '''
+ ## Checking the easy cases
coer = None;
if(isinstance(S, DDRing)):
coer = self.base()._coerce_map_from_(S.base());
if(not(coer is False) and not(coer is None)):
return True;
- return None;
+ #return None;
+ ## If not coercion is found, we check deeper using generators
+ gens_self, p_self = DDRing.__get_gens__(self);
+ try:
+ gens_S, pS = DDRing.__get_gens__(S);
+ except RuntimeError:
+ return None;
+
+ ## Using the 'other' generators
+ if(len(gens_S['other']) > _sage_const_0 ):
+ return self.base()._coerce_map_from_(S);
+
+ ## Comparing the algebraic construction
+ for i in range(len(gens_S['algebraic'])):
+ 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 ])):
+ found = True;
+ break;
+ if(not found):
+ return self.base()._coerce_map_from_(S);
+
+ ## Comparing the parametric space
+ if(any(str(el) not in [str(ob) for ob in gens_self['polynomial']] for el in gens_S['polynomial'])):
+ return self.base()._coerce_map_from_(S);
+
+ ## Comparing the depth of the structure S
+ if(isinstance(S, DDRing)):
+ return S.depth() <= self.depth();
+
+ #coer = None;
+ #if(isinstance(S, DDRing)):
+ # coer = self.base()._coerce_map_from_(S.base());
+ #elif(S == self.base()):
+ # coer = True;
+ #elif(isinstance(S, sage.symbolic.ring.SymbolicRing)):
+ # coer = True;
+ #else:
+ # coer = self.base()._coerce_map_from_(S);
+ #
+ #if(not(coer is False) and not(coer is None)):
+ # return True;
+ #return None;
def __is_polynomial(self, S):
from sage.rings.polynomial.polynomial_ring import is_PolynomialRing as isUniPolynomial;
return isUniPolynomial(S) or isMPolynomial(S);
def _pushout_(self, S):
- if(isinstance(S, DDRing)):
- from sage.categories.pushout import pushout;
- if(isinstance(S, ParametrizedDDRing) and isinstance(self, ParametrizedDDRing)):
- if(len(set(self.parameters(True)).intersection(S.parameters(True))) > _sage_const_0 ):
- import traceback;
- entry = traceback.extract_stack(limit=_sage_const_3 )[-_sage_const_3 ];
- warnings.warn_explicit("Merging rings with common parameters. \n\t* %s\n\t* %s\n--Are you sure they are the same?" %(repr(self),repr(S)), MergingRingWarning, *(list(entry[:_sage_const_3 ])));
- pushout_ddRing = pushout(self.base_ddRing(),S.base_ddRing());
- return ParametrizedDDRing(pushout_ddRing, set(self.parameters(True)).union(S.parameters(True)));
- elif(isinstance(S, ParametrizedDDRing)):
- return S._pushout_(self);
- elif(isinstance(self, ParametrizedDDRing)):
- pushout_ddRing = pushout(self.base_ddRing(),S);
- return ParametrizedDDRing(pushout_ddRing, self.parameters(True));
- else:
- pushout_base = pushout(self.original_ring(),S.original_ring());
- pushout_field = pushout(self.base_field, S.base_field);
- return DDRing(pushout_base,depth = max(self.depth(),S.depth()),base_field=pushout_field);
+ '''
+ This method compute the pushout of 'self' with the parent class 'S'. This method returns None if it is not possible to get the pushout or a DDRing with all the properties required to include 'self' and 'S'.
+
+ The method goes as follows:
+ 1 - Compute the simplest field over everything is based (let us call it 'F'). This will be 'QQ' or an algebraic extension of 'QQ' (see the class 'NumberField')
+ 2 - Compute the list of parameters involved. This will be the constant trascendent elements that will extend the ring 'F' as a rational function ring. Let 'B' be such extension.
+ 3 - Compute the list of variables involved. We usually expect just one variable but allow to have more of them. We will build 'R' the polynomial ring over 'F' with those variables.
+ 4 - Compute the depth required for include 'self' and 'S'. Then we build the corresponding DDRing of the appropriate depth. If we had parameters, we build the ParametrizedDDRing.
- const = S.construction();
- if((not(const is None)) and (isinstance(const[_sage_const_0 ], sage.categories.pushout.FractionField))):
- S = S.base();
- if(self.__is_polynomial(S)):
- all_params = set(str(el) for el in DDRing.__get_gens__(S));
- params = all_params - set([str(el) for el in self.variables(True)]);
- variables = all_params - params;
-
- base_ring = self;
- if(len(variables) > _sage_const_0 ):
- base_ring = DDRing(PolynomialRing(self.base_field, list(variables)), depth=self.depth());
-
- if(len(params) > _sage_const_0 ):
- return ParametrizedDDRing(base_ring, params);
- else:
- return base_ring;
+ INPUT:
+ - A parent structure of SAGE 'S'.
- return None;
+ OUTPUT:
+ - A DDRing or ParametrizedDDRing such that any element of 'self' and 'S' is in it.
+
+ WARNINGS:
+ - A warning will pop up if we merge rings with the same parameter names.
+
+ ERRORS:
+ - 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)):
+ return None;
+
+ ## We get a list of the generators of self and S with their types
+ gens_self, pself = DDRing.__get_gens__(self);
+ try:
+ gens_S, pS = DDRing.__get_gens__(S);
+ except RuntimeError:
+ return None;
+
+ if(len(gens_S['other']) > _sage_const_0 or len(gens_self['other']) > _sage_const_0 ):
+ raise TypeError("Impossible to compute a pushout: generators not recognized found.\n\t- %s" %(list(gens_S['other']) + list(gens_self['other'])));
+
+ ##########################################
+ ## Compute the basic field F
+ ##########################################
+ ## Computing the original field
+ F = None;
+ try:
+ F = pushout(pself, pS);
+ except:
+ pass;
+ if(F is None):
+ raise TypeError("Incompatible original structures:\n\t- %s\n\t- %s" %(pself, pS));
+
+ if(not F in IntegralDomains()):
+ raise TypeError("Pushout of the original structures is not an integral domain:\n\t- %s" %F);
+ if(not F.is_field()):
+ 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']}
+ 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 ];
+
+ sorted_names = sorted(polys);
+ sorted_polys = [polys[el] for el in sorted_names];
+
+ F = NumberField([PolynomialRing(F,x)(el) for el in sorted_polys], sorted_names);
+
+ ##########################################
+ ## Compute the list of parameters and variables
+ ##########################################
+ all_params = set(str(el) for el in (gens_S['polynomial']+gens_self['polynomial']));
+ params = all_params - set([str(el) for el in self.variables(True)]);
+ variables = all_params - params;
+
+ ##########################################
+ ## Compute the required depth
+ ##########################################
+ depth = self.depth();
+ if(isinstance(S, DDRing)):
+ depth = max(depth,S.depth());
+
+ ##########################################
+ ## Building the final DDRing
+ ##########################################
+ if(len(variables) > _sage_const_0 ):
+ F = PolynomialRing(F,[str(el) for el in variables]);
+ R = DDRing(F, depth = depth);
+ if(len(params) > _sage_const_0 ):
+ R = ParametrizedDDRing(R, params);
+
+ return R;
def _has_coerce_map_from(self, S):
coer = self._coerce_map_from_(S);
return self.element([coeff for coeff in X.equation.getCoefficients()], X.getInitialValueList(X.equation.get_jp_fo()+_sage_const_1 ), name=X._DDFunction__name);
else:
try:
- X = self.base()(X);
+ try:
+ X = self.base()(X);
+ except TypeError as e:
+ try:
+ X = self.base()(str(X));
+ except:
+ raise e;
el = self.element([-self.base_derivation(X), X]);
name = str(X);
try:
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);
except AttributeError:
- return self.element([_sage_const_1 ],[], self.base()(X), name=str(X));
+ try:
+ return self.element([_sage_const_1 ],[], self.base()(X), name=str(X));
+ except Exception:
+ return self(str(element));
except TypeError as e:
- #try:
- # if(isinstance(X.parent(),sage.symbolic.ring.SymbolicRing)):
- # #print "Casting a symbolic"
- # new_parent = ParametrizedDDRing(self,X.variables())
- # el = new_parent(X);
- # print new_parent;
- # print el;
- # return el
- #except:
raise TypeError("Cannot cast an element to a DD-Function of (%s):\n\tElement: %s (%s)\n\tReason: %s" %(self, X, type(X), e));
def gens(self):
def __contains__(self, X):
try:
- return (X.parent() is self) or (self._has_coerce_map_from(X.parent()));
- except AttributeError:
- try:
- self(X)
+ if((X.parent() is self) or (self._has_coerce_map_from(X.parent()))):
return True;
- except Exception:
- return False;
+ except AttributeError:
+ pass;
+ try:
+ self(X)
+ return True;
+ except Exception:
+ return False;
#################################################
### Magic python methods
def to_depth(self, dest):
return DDRing(self.original_ring(), depth = dest, base_field = self.base_field, invertibility = self.base_inversionCriteria, derivation = self.base_derivation, default_operator = self.default_operator);
- def is_field(self):
+ def is_field(self, **kwds):
return False;
- def is_finite(self):
+ def is_finite(self, **kwds):
return False;
- def is_noetherian(self):
+ def is_noetherian(self, **kwds):
return True;
#################################################
@staticmethod
def __classcall__(cls, *args, **kwds):
## In order to call the __classcall__ of DDRing we treat the arguments received
- if(len(args) != _sage_const_2 ):
- raise TypeError("Incompatible number of arguments for ParametrizedDDRing");
-
- base_ddRing = args[_sage_const_0 ]; parameters = args[_sage_const_1 ];
+ base_ddRing = args[_sage_const_0 ];
+ if(len(args) > _sage_const_1 ):
+ parameters = args[_sage_const_1 ];
+ else:
+ parameters = kwds.get('parameters',None);
+ names = kwds.get('names',None);
+
+ ## Using the "names" approach of SAGE
+ if(parameters is None and names is None):
+ raise ValueError("Invalid parameters input format: no parameters given");
+ elif(parameters is None):
+ parameters = names;
+ elif(not(names is None)):
+ raise ValueError("Invalid syntax creating a ParametrizedDDRing. Use one of the following syntaxes:\n\t- D = ParametrizedDDRing(R,['a','b'])\n\t- D.<a,b> = ParametrizedDDRing(R)");
## First we get the new variables
if ((not(type(parameters) == tuple)) and (not(type(parameters) == list)) and (not(type(parameters) == set))):
def gens(self):
return self.parameters(True);
+ def _first_ngens(self, n):
+ return self.parameters(False)[:n];
+
def ngens(self):
return len(self.parameters());
while(self.getInitialValue(n) == _sage_const_0 ):
n = n+_sage_const_1 ;
+ X = self.parent().variables()[_sage_const_0 ];
if(n == _sage_const_0 ):
return (_sage_const_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;
+ 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 )];
### Computing the new name
newName = None;
if(not(self.__name is None)):
- newName = DinamicString("(_1)/(_2^%d)" %n, [self.__name, "x"]);
+ newName = DinamicString("(_1)/(_2^%d)" %n, [self.__name, repr(X)]);
result = self.parent().element(newEquation, newInit, check_init=False, name=newName);
- result.change_built("polynomial", (PolynomialRing(self.parent().base_field,['x','x1']).fraction_field()('x1/(x^%d)' %n), {'x':self.parent()(x),'x1':self}));
+ result.change_built("polynomial", (PolynomialRing(self.parent().base_field,[repr(X),'x1']).fraction_field()('x1/(%s^%d)' %(repr(X),n)), {repr(X):self.parent()(X),'x1':self}));
return (n,result);
def divide(self, other):
return other;
elif(other.is_null):
return self;
+
+ X = self.parent().variables()[_sage_const_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[_sage_const_0 ],oe[_sage_const_0 ])*gcd(se[_sage_const_1 ].getInitialValue(_sage_const_0 ),oe[_sage_const_1 ].getInitialValue(_sage_const_0 )));
#####################################
### Differential methods
return self.parent().element(newOperator, newInit, check_init=False, name=newName);
- def compose(self, element):
+ def compose(self, other):
'''
- Method to get a DDFunction `g` that is the composition of `self` with `other`. It return the composition in the corresponding ring.
- '''
- ## Checking we fall in teh cases already implemented
- if((not (element in self.parent().original_ring())) and (self.getOrder() > _sage_const_3 )):
- raise NotImplementedError("Composition not yet implemented with elements not in %s for functions of order greater than 2 (got %d)" %(self.parent().original_ring(),self.getOrder()));
-
+ Method to compute the composition of 'self' with 'other' (if possible).
+
+ The method first compute the new ring where the composition will belong and then relies on the method 'compose_solution' of the Operator class.
+
+ Then, it computes the new initial values using the Faa di Bruno's formula.
+
+ INPUT:
+ - 'self': a DDFunction
+ - 'other': an element of SAGE
+
+ OUTPUT:
+ - A new DDFunction 'f' in the corresponding ring such f == self(other)
+
+ ERRORS:
+ - ValueError: if other(0) != 0
+ - TypeError: if we can not compute a destination ring
+ - Any error from Operator.compose_solution
+
+ WARNINGS:
+ - If the depth of the resulting DDFunction is greater than 3, a warning will pop-up
+ '''
+ ######################################
+ ## Initial checking and considerations
+ ######################################
+ ## Checking that the 'other' element is NOT in the Symbolic Ring
+ if(isinstance(other.parent(), sage.symbolic.ring.SymbolicRing)):
+ try:
+ other = self.parent().original_ring()(str(other));
+ except Exception as e:
+ raise TypeError("Impossible to perform a composition with a symbolic element. Try to cast (%s) to some field, polynomial ring or some DDRing" %(other));
+
## 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 ];
- if(self_var == element or self.is_constant):
+ if(self_var == other or self.is_constant):
return self;
-
- ## Otherwise, we check that the value at 0 is zero
- init = element(**{str(self_var):_sage_const_0 });
- if(init != _sage_const_0 ):
- raise ValueError("Impossible to compose a power series with other that has non-zero initial value.\n\t- Get: %s" %init);
-
- ## We compute the differential operator case by case
- newOperator = None;
+
+ ## Checking that 'other' at zero is zero
+ value = None;
+ try:
+ value = other(**{str(self_var):_sage_const_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 );
+ 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 })));
+
+ ######################################
+ ## Computing the destiny ring
+ ######################################
destiny_ring = None;
- if(not(element in self.parent().original_ring())):
- ## Particular cases: we have a formula
- ## First, we get the pushout of both parents
- ## This will be a DDRing always
- from sage.categories.pushout import pushout;
- mix_ring = pushout(self.parent(), element.parent());
-
- ## We compute the depth in the DDRing structure of 'element'
- element_depth = _sage_const_0 ;
- if(isinstance(element.parent(),DDRing)):
- element_depth = element.parent().depth();
-
- ## We build the destiny DDRing
- destiny_ring = mix_ring.to_depth(self.parent().depth() + element_depth);
- if(destiny_ring.depth() >= _sage_const_3 ):
- warnings.warn("Reached a too high depth in the recursion (%d). This may not finish." %destiny_ring.depth(), TooDeepWarning, stacklevel=_sage_const_2 );
- #element = destiny_ring.base()(element);
-
- if(element in destiny_ring.original_ring()):
- return destiny_ring(self).compose(element);
- elif(self.getOrder() == _sage_const_0 ):
- raise ValueError("Impossible Error: Non-constant element of order zero!!");
- elif(self.getOrder() == _sage_const_1 ):
- a = self.equation.getCoefficient(_sage_const_1 )(x=element);
- b = self.equation.getCoefficient(_sage_const_0 )(x=element);
- dg = destiny_ring.base_derivation(element);
- newOperator = destiny_ring.element([dg*b,a]).equation;
- elif(self.getOrder() == _sage_const_2 ):
- a = self.equation.getCoefficient(_sage_const_2 )(x=element);
- b = self.equation.getCoefficient(_sage_const_1 )(x=element);
- c = self.equation.getCoefficient(_sage_const_0 )(x=element);
- dg = destiny_ring.base_derivation(element);
- ddg = destiny_ring.base_derivation(dg);
- newOperator = destiny_ring.element([dg**_sage_const_3 *c,dg**_sage_const_2 *b-ddg*a,dg*a]).equation;
- elif(self.getOrder() == _sage_const_3 ):
- a = self.equation.getCoefficient(_sage_const_3 )(x=element);
- b = self.equation.getCoefficient(_sage_const_2 )(x=element);
- c = self.equation.getCoefficient(_sage_const_1 )(x=element);
- d = self.equation.getCoefficient(_sage_const_0 )(x=element);
- dg = destiny_ring.base_derivation(element);
- ddg = destiny_ring.base_derivation(dg);
- dddg = destiny_ring.base_derivation(ddg);
- newOperator = destiny_ring.element([dg**_sage_const_5 *d, dg**_sage_const_4 *c-ddg*dg**_sage_const_2 *b-dg*dddg*a-_sage_const_3 *ddg**_sage_const_2 *a,dg**_sage_const_3 *b-_sage_const_3 *ddg*dg*a,dg**_sage_const_2 *a]).equation;
+ ## First, compute the pushout fo the parents
+ from sage.categories.pushout import pushout;
+ sp = self.parent(); op = other.parent();
+ push = pushout(sp, op);
+
+ ## Second, compute the final depth of the DDRing
+ if(not isinstance(push, DDRing)):
+ raise TypeError("A non DDRing obtained for the composition: that is impossible -- review method _pushout_ of DDRing class");
+ if(isinstance(op, DDRing)):
+ destiny_ring = push.to_depth(op.depth()+sp.depth());
else:
- newOperator = self.equation.compose_solution(element);
- destiny_ring = self.parent();
-
-
- ### We compute the initial values required
- required = newOperator.get_jp_fo()+_sage_const_1 ;
-
- f = self.getInitialValueList(required);
+ 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 );
+
+ ######################################
+ ## Computing the new operator
+ ######################################
+ ## HERE IS THE MAIN RECURSION STEP
+ equation = destiny_ring.element([coeff(**{str(self_var) : other}) for coeff in self.equation.getCoefficients()]).equation; ## Equation with coefficients composed with 'other'
+ g = destiny_ring.base()(other); ## Casting the element 'other' to the base ring
+
+ new_equation = equation.compose_solution(g);
+
+ ######################################
+ ## Computing the new initial values
+ ######################################
+ required = new_equation.get_jp_fo()+_sage_const_1 ;
+ ## Getting as many initial values as we can and need
+ init_f = self.getInitialValueList(required);
+ init_g = None;
try:
- ## Case element is simpler than self
- g = [factorial(i)*self.parent().getSequenceElement(element,i) for i in range(required)];
- except TypeError:
- ## Case element has higher complexity than self
- g = [factorial(i)*element.parent().getSequenceElement(element,i) for i in range(required)];
+ 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)];
+ ## 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
- newInit = [f[_sage_const_0 ]]+[sum([f[j]*bell_polynomial(i,j)(*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 ,required)]; ## See Faa di Bruno's formula
- ### Computing the new name
- newName = None;
+ ######################################
+ ## Computing the new name
+ ######################################
+ new_name = None;
if(not(self.__name is None)):
- # if(isinstance(element, DDFunction) and not(element.__name is None)):
- # newName = self.__name.replace("exp", "####").replace("x", element.__name).replace("####", "exp");
- # else:
- # newName = self.__name.replace("exp", "####").replace("x", repr(element)).replace("####", "exp");
- if(isinstance(element, DDFunction) and (not (element.__name is None))):
- newName = din_m_replace(self.__name, {"x":element.__name}, True);
- elif(not isinstance(element, DDFunction)):
- newName = din_m_replace(self.__name, {"x":repr(element)}, True);
-
- return destiny_ring.element(newOperator, newInit, check_init=False, name=newName);
+ if(isinstance(other, DDFunction) and (not (other.__name is None))):
+ new_name = din_m_replace(self.__name, {"x":other.__name}, True);
+ elif(not isinstance(other, DDFunction)):
+ new_name = din_m_replace(self.__name, {"x":repr(other)}, True);
+ return destiny_ring.element(new_equation, new_init, name=new_name);
+
#####################################
### Equality methods
#####################################
newName = None;
if(not(self.__name is None)):
newName = DinamicString("(_1)^%d" %(other), self.__name);
- self.__pows[other].__name = newName;
+ if(isinstance(self.__pows[other], DDFunction)):
+ self.__pows[other].__name = newName;
else:
try:
inverse = self.inverse;
return self.__depth;
class ParametrizedDDRingFunctor (DDRingFunctor):
- def __init__(self, depth, base_field, var = {}):
+ def __init__(self, depth, base_field, var = set([])):
self.rank = _sage_const_12 ;
self.__vars = var;
super(ParametrizedDDRingFunctor, self).__init__(depth, base_field);
raise TypeError("The object [%s] is not an element of [%s]" %(x, self.domain()));
return x;
+ def _repr_(self):
+ return "ParametrizedDDRing(*,%s)" %(self.__vars);
+
def _apply_functor(self, x):
return ParametrizedDDRing(x, self.__vars);
super(DDSimpleMorphism, self).__init__(domain, codomain);
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 ));
+ 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 ));
projects / ajpastor / diff_defined_functions.git / commitdiff