From: Antonio-JP Date: Fri, 22 Mar 2019 12:10:53 +0000 (+0100) Subject: Middle point implementing the complete DAlgebraic X-Git-Url: http://git.risc.jku.at/gitweb/?a=commitdiff_plain;h=f933e5b715c19a58bc4c8a2a2e5054b67e1115d5;p=ajpastor%2Fdiff_defined_functions.git Middle point implementing the complete DAlgebraic Fixed some errors in lazyDDRing.py --- diff --git a/ajpastor/dd_functions/ddExamples.py b/ajpastor/dd_functions/ddExamples.py index eca1dfb..4dbe537 100644 --- a/ajpastor/dd_functions/ddExamples.py +++ b/ajpastor/dd_functions/ddExamples.py @@ -5,6 +5,7 @@ from sage.rings.polynomial.polynomial_ring import is_PolynomialRing; from sage.rings.polynomial.multi_polynomial_ring import is_MPolynomialRing; from ajpastor.dd_functions.ddFunction import *; +from ajpastor.dd_functions.lazyDDRing import LazyDDRing; from ajpastor.misc.dinamic_string import *; @@ -1862,24 +1863,33 @@ def DAlgebraic(polynomial, init=[], dR=None): ############################################### parent = polynomial.parent(); if(not (isPolynomial(parent) or isMPolynomial(parent))): - raise TypeError("The minimal polynomial is NOT a polynomial"); + raise TypeError("DAlgebraic: the input polynomial is NOT a polynomial"); base_ring = None; F = None; poly_ring = parent; if(isMPolynomial(parent)): - base_ring = PolynomialRing(parent.base(),parent.gens()[:-1]); - poly_ring = PolynomialRing(base_ring.fraction_field(), parent.gens()[-1]); + ## Only valid for 2 variables + if(parent.ngens() > 2): + raise TypeError("DAlgebraic: the input can not be a multivariate polynomial with more than 2 variables"); + base_ring = PolynomialRing(parent.base(),parent.gens()[0]); + F = base_ring.fraction_field(); else: if(isinstance(parent.base().construction()[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()[-1]); - - F = poly_ring.base(); + + if(is_DDRing(base_ring)): + F = LazyDDRing(base_ring); + elif(not parent.base().is_field()): + F = parent.base().fraction_field(); + else: + F = base_ring; + + poly_ring = PolynomialRing(F, parent.gens()[-1]); y = poly_ring.gens()[-1]; + ## At this point we have the following ## - F is a field ## - y is a variable @@ -1888,7 +1898,7 @@ def DAlgebraic(polynomial, init=[], dR=None): if(polynomial.degree() == 1): return -polynomial[0]/polynomial[1]; elif(polynomial.degree() <= 0): - raise TypeError("Constant polynomial given for algebraic function: IMPOSSIBLE!!"); + raise TypeError("DAlgebraic: constant polynomial given for algebraic function: IMPOSSIBLE!!"); ################################################# ## Building and checking the destiny ring @@ -1898,7 +1908,7 @@ def DAlgebraic(polynomial, init=[], dR=None): destiny_ring = DDRing(base_ring); else: destiny_ring = dR; - coercion = destiny_ring._coerce_map_from_(base_ring); + coercion = destiny_ring.base()._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)); @@ -1909,34 +1919,61 @@ def DAlgebraic(polynomial, init=[], dR=None): ################################################## ## Computing the derivative dy = polynomial.derivative(y); - - ## Getting its gcd with the polynomial - g,r,s = polynomial.xgcd(dy); - if((g != 1) or (not(g in poly_ring.base()))): - raise ValueError("No irreducible polynomial given"); ## Computing the coefficient-wise derivation of polynomial - mon = poly_ring(1); - ky = poly_ring(0); - for i in range(polynomial.degree()+1): - ky += (polynomial[i].derivative())*mon; - mon *= y; - - ## Getting the polynomial expression of y', y'',..., y^{(deg(polynomial))} - rows = [[0]*polynomial.degree()]; - mon = poly_ring(1); - for i in range(polynomial.degree()-1): - rows += [((-(i+1)*mon*s*ky)%polynomial).coefficients(False)]; - mon *= y; - + ky = sum(polynomial[i].derivative()*y**i for i in range(polynomial.degree()+1)); + + ### WARNING From this point until said otherwise, in the case base_ring is DDRing + ### we change the poylnomials to a finite setting because when this was code, + ### computing remainders with InfinitePolynomial coefficients broke the program. + ### + ### For further information about it, please check the track: (pending) + ### The code that should work once that is solved is the following: + if(is_DDRing(base_ring)): + # Changing to a finite setting + def get_sub_vars(*polys): + if(len(polys) > 1): + return get_sub_vars(polys); + return list(set(sum([sum([list(el[i].variables()) for i in range(el.degree()+1)],[]) for el in polys[0]],[]))); + + sub_vars = get_sub_vars([polynomial, dy, ky]); + poly_ring_fin = PolynomialRing(PolynomialRing(base_ring.base_ring(), sub_vars).fraction_field(), poly_ring.gens_dict().items()[0][0]); + polynomial = poly_ring_fin(str(polynomial)) + dy = poly_ring_fin(str(dy)); + ky = poly_ring_fin(str(ky)); + y = poly_ring_fin.gens()[0]; + + ## Getting its gcd with the polynomial + g,r,s = polynomial.xgcd(dy); + if((g != 1) and (g.degree() != 0)): + raise ValueError("DAlgebraic: no irreducible polynomial (%s) given" %polynomial); + + ## Getting the polynomial expression of y' + yp = (-ky*s)%polynomial; + ## Building the derivation matrix of <1,y,y^2,...> + rows = [[0]*polynomial.degree()]; + for i in range(1,polynomial.degree()): + # (y^i)' = i*y^(i-1)*y' + current = ((i*yp*y**(i-1))%polynomial).coefficients(False); + current += [0 for i in range(len(current),polynomial.degree())]; + rows += [current]; + M = Matrix(F, rows).transpose(); ## Creating the vector representing y y_vector = vector(F, [0,1] + [0]*(polynomial.degree()-2 )); - ## Building ans solving the system + + if(is_DDRing(base_ring)): raise RuntimeError("DEBUG Stop"); + + ## Building and solving the system to_solve = move(M, y_vector, lambda p : p.derivative(), M.ncols()+1); + + if(is_DDRing(base_ring)): raise RuntimeError("DEBUG Stop"); + v = to_solve.right_kernel_matrix()[0]; + if(is_DDRing(base_ring)): raise RuntimeError("DEBUG Stop"); + ## Cleaning denominators cleaning = lcm(el.denominator() for el in v); diff --git a/ajpastor/dd_functions/lazyDDRing.py b/ajpastor/dd_functions/lazyDDRing.py index f14e0eb..f8b71d4 100644 --- a/ajpastor/dd_functions/lazyDDRing.py +++ b/ajpastor/dd_functions/lazyDDRing.py @@ -69,7 +69,11 @@ class _LazyDDFunction(IntegralDomainElement): return self.__poly; def simplify(self): - self.__poly = self.parent().simplify(self.poly()); + if(self.__poly != 0): + # First we remove common factors in numerator and denominator + n_pol = self.__poly.parent()(prod(self.__poly.factor())); + # We then call the simplify method from parent + self.__poly = self.parent().simplify(n_pol); return self; def variables(self): @@ -269,7 +273,6 @@ class LazyDDRing (UniqueRepresentation, ConversionSystem, IntegralDomain): self.__poly_field = self.__poly_ring.fraction_field(); self.__gen = self.__poly_ring.gens()[0]; self.__used = 0; - self.stop = False; ## Initializing the map of variables self.__map_of_vars = {}; # Map where each variable is associated with a 'base' object @@ -449,8 +452,19 @@ class LazyDDRing (UniqueRepresentation, ConversionSystem, IntegralDomain): elif(times > 1): return self.derivative(self.derivative(el, times-1)); - ## Splitting the derivative into the monomials + ## COmputing the derivative of the polynomial associated poly = el.poly(); + ## Rational function case + try: + if(poly.denominator() != 1): + n = self(poly.numerator()); dn = self.derivative(n); + d = self(poly.denominator()); dd = self.derivative(d); + + return (dn*d - dd*n)/d**2; + except AttributeError: + pass; + + ## Splitting the derivative into the monomials if(poly == 0): return self.zero(); coeffs = poly.coefficients(); @@ -494,8 +508,9 @@ class LazyDDRing (UniqueRepresentation, ConversionSystem, IntegralDomain): ################################################################################################ ### Other Ring methods ################################################################################################ - def is_field(self): + def is_field(self, proof=False): return self.base().is_field(); + #return True; def is_finite(self): return self.base().is_finite(); @@ -505,6 +520,38 @@ class LazyDDRing (UniqueRepresentation, ConversionSystem, IntegralDomain): def is_noetherian(self): return self.base().is_noetherian(); + + def _xgcd_univariate_polynomial(self, a, b): + ''' + Return an extended gcd of ``a`` and ``b``. + + INPUT: + + - ``a``, ``b`` -- two univariate polynomials defined over ``self`` + + OUTPUT: + + A tuple ``(t, r, s)`` of polynomials such that ``t`` is the + greatest common divisor (monic or zero) of ``a`` and ``b``, + and ``r``, ``s`` satisfy ``t = r*a + s*b``. + ''' + SY = a.parent(); + ## Casting the polynomials a and b to polynomials with rational functions coefficients + # Getting the variables + coeffs = [self.to_poly(el) for el in a.coefficients() + b.coefficients()]; + gens = list(set(sum([[str(g) for g in poly.variables()] for poly in coeffs], []))); + # Building the rational functions + F = PolynomialRing(self.base_ring(), gens).fraction_field(); + # Building the polynomial ring + y = a.parent().gens_dict().items()[0][0]; + R = PolynomialRing(F, y); + a = R(str(a)); b = R(str(b)); + + ## Computing the xgcd in that extended field using the standard algorithm + t,r,s = F._xgcd_univariate_polynomial(a,b); + + ## Returning the result casted to Lazy Elements + return (SY(t.coefficients(False)), SY(r.coefficients(False)), SY(s.coefficients(False))); ################################################################################################ ### Coercion methods diff --git a/ajpastor/lazy/conversion.py b/ajpastor/lazy/conversion.py index 5a68890..74cff07 100644 --- a/ajpastor/lazy/conversion.py +++ b/ajpastor/lazy/conversion.py @@ -214,7 +214,7 @@ class ConversionSystem(object): try: return self._simplify(self.poly_ring()(element)); except AttributeError: - return self.poly_ring()(element); + return element; elif(element in self.poly_field()): element = self.poly_field()(element); n = self.simplify(element.numerator()); diff --git a/ajpastor/misc/sequence_manipulation.py b/ajpastor/misc/sequence_manipulation.py index 64eae8e..72b3ba3 100644 --- a/ajpastor/misc/sequence_manipulation.py +++ b/ajpastor/misc/sequence_manipulation.py @@ -1,6 +1,7 @@ from sage.functions.other import factorial; from sage.combinat.combinat import bell_polynomial; from sage.arith.misc import falling_factorial; +from sage.misc.cachefunc import cached_function ################################################################################ ################################################################################ @@ -24,6 +25,7 @@ def shift(f): def composition(f,g): if(g(0) != 0): raise ValueError("Impossible compose with g(0) != 0"); + @cached_function def _composition(n): if(n == 0): return f(0); @@ -59,6 +61,7 @@ def inv_lagrangian(f): if(f(1) == 0): raise NotImplementedError("Case with order higher than 1 not implemented"); f = ogf_egf(f); + @cached_function def _inverse_egf(n): if(n == 0): return 0; diff --git a/releases/diff_defined_functions.zip b/releases/diff_defined_functions.zip index cbf4346..f23ab03 100644 Binary files a/releases/diff_defined_functions.zip and b/releases/diff_defined_functions.zip differ diff --git a/releases/diff_defined_functions__0.6.zip b/releases/diff_defined_functions__0.6.zip index cbf4346..f23ab03 100644 Binary files a/releases/diff_defined_functions__0.6.zip and b/releases/diff_defined_functions__0.6.zip differ diff --git a/releases/old/diff_defined_functions__0.6__19.03.22_13:10:53.zip b/releases/old/diff_defined_functions__0.6__19.03.22_13:10:53.zip new file mode 100644 index 0000000..f23ab03 Binary files /dev/null and b/releases/old/diff_defined_functions__0.6__19.03.22_13:10:53.zip differ