Opened 8 years ago
Closed 8 years ago
#10771 closed enhancement (fixed)
gcd and lcm for fraction fields
Reported by: | SimonKing | Owned by: | AlexGhitza |
---|---|---|---|
Priority: | major | Milestone: | sage-4.7.2 |
Component: | basic arithmetic | Keywords: | gcd lcm fraction fields |
Cc: | burcin | Merged in: | sage-4.7.2.alpha0 |
Authors: | Simon King | Reviewers: | Marco Streng, Mariah Lenox |
Report Upstream: | N/A | Work issues: | |
Branch: | Commit: | ||
Dependencies: | Stopgaps: |
Description (last modified by )
At sage-devel, the question was raised whether it really is a good idea that the gcd in the rational field should return either 0
or 1
.
Since any non-zero element of QQ
qualifies as gcd of two non-zero rationals, it should be possible to define gcd and lcm, so that gcd(x,y)*lcm(x,y)==x*y
holds for any rational numbers x,y, and so that gcd(QQ(m),QQ(n))==gcd(m,n)
and lcm(QQ(m),QQ(n))==lcm(m,n)
for any two integers m,n.
Moreover, it should be possible to provide gcd/lcm for any fraction field of a PID
: Note that currently gcd raises a type error for elements of Frac(QQ['x'])
.
The aim is to implement gcd and lcm as ElementMethods
of the category QuotientFields()
.
Approach
Let R be an integral domain, assume that it provides gcd and lcm, and let F be its fraction field. Since R has gcd, we can assume that x.numerator()
and x.denominator()
are relatively prime for any element x of F.
Then, define
gcd(x,y) = gcd(x.numerator(),y.numerator())/lcm(x.denominator(),y.denominator()) lcm(x,y) = lcm(x.numerator(),y.numerator())/gcd(x.denominator(),y.denominator())
Benefits
If that approach is mathematically sober, we obtain the following equalities up to units in R:
gcd(x,y)*lcm(x,y)==x*y
, for any x,y in F, provided that the equality holds for any x,y in R.gcd(F(x),F(y))==gcd(x,y)
andlcm(F(x),F(y))==lcm(x,y)
for any x,y in R.
Attachments (1)
Change History (39)
comment:1 Changed 8 years ago by
- Description modified (diff)
comment:2 Changed 8 years ago by
- Status changed from new to needs_review
- Type changed from PLEASE CHANGE to defect
comment:3 follow-up: ↓ 6 Changed 8 years ago by
- Reviewers set to Marco Streng
A couple of comments:
- I don't see the point of keeping a special function
gcd_rational(self, other, **kwds):
that returns a GCD in the set {0,1}. Why only for QQ? Why at all? Also, the difference betweengcd
andgcd_rational
is not explained by the word "rational". Perhapsgcd_zero_one
would be a more informative name.
- Everything generalizes from principal ideal domains to unique factorization domains (UFD's) (such as multivariate polynomial rings over unique factorization domains) as long as they have
gcd
andlcm
methods implemented. Why not write "unique factorization domain" in the documentation instead of "principal ideal domain"?
- Suppose I have a fraction field F of a ring of type R that does not have lcm and gcd methods, or these methods exist, but raise other kinds of errors, e.g. because the ring is not a UFD or the methods have not been implemented. Let a and b be elements of F. Then a and b have a gcd in F because F is a field, so I would expect a.gcd(b) to return something (anything basically). After applying your patch, if I do
a.gcd(b)
, it is very confusing to get anAttributeError: 'RElement' object has no attribute 'gcd'
: I'm not interested in gcd's of RElements, only of(Fraction)FieldElements
. You could put your entiregcd
andlcm
code betweentry:
andexcept (AttributeError, NotImplementedError, TypeError, ValueError):
to return the same 0 or 1 thatgcd_rational
would (which is a mathematically correct gcd in F after all).
- trac has a
t
too much in line 1610 of ring.pyx:the case of the rational field. However, since tract ticket #10771,
- you write "quotient field" in the documentation. You could write "fraction field" to avoid confusion with quotient rings, which may be fields. This makes it more clear that you refer to the mathematical counterpart of Sage's "
FractionField
"?
All tests passed!
with -long, well documented, looks good.
comment:4 Changed 8 years ago by
comment:5 Changed 8 years ago by
Trivial: typo "tract" for "trac" in patch.
comment:6 in reply to: ↑ 3 Changed 8 years ago by
Hi Marco,
Replying to mstreng:
A couple of comments:
- I don't see the point of keeping a special function
gcd_rational(self, other, **kwds):
that returns a GCD in the set {0,1}. Why only for QQ? Why at all? Also, the difference betweengcd
andgcd_rational
is not explained by the word "rational". Perhapsgcd_zero_one
would be a more informative name.
Perhaps I am too 'conservative': I would rather rename something than to delete it.
But I suppose it would be better to implement a gcd and an lcm as element methods of Fields()
-- see #9819. In that way, one could easily provide a correct gcd/lcm behaviour for all fields, and for fraction fields of PID one would still obtain a gcd/lcm behaviour that is not only correct but nice.
- Everything generalizes from principal ideal domains to unique factorization domains (UFD's) (such as multivariate polynomial rings over unique factorization domains) as long as they have
gcd
andlcm
methods implemented. Why not write "unique factorization domain" in the documentation instead of "principal ideal domain"?
Good question. I did it by accident. However, it matches the sad fact that the categories are a little askew here:
sage: PrincipalIdealDomains().is_subcategory(UniqueFactorizationDomains()) False # shouldn't that be "True"? sage: PrincipalIdealDomains().is_subcategory(GcdDomains()) True sage: UniqueFactorizationDomains().is_subcategory(GcdDomains()) True
Perhaps one should write "fraction field of an integral domain with gcd and lcm"? Because that's what is duck typed.
- Suppose I have a fraction field F of a ring of type R that does not have lcm and gcd methods, or these methods exist, but raise other kinds of errors, e.g. because the ring is not a UFD or the methods have not been implemented. Let a and b be elements of F. Then a and b have a gcd in F because F is a field, so I would expect a.gcd(b) to return something (anything basically). After applying your patch, if I do
a.gcd(b)
, it is very confusing to get anAttributeError: 'RElement' object has no attribute 'gcd'
: I'm not interested in gcd's of RElements, only of(Fraction)FieldElements
. You could put your entiregcd
andlcm
code betweentry:
andexcept (AttributeError, NotImplementedError, TypeError, ValueError):
to return the same 0 or 1 thatgcd_rational
would (which is a mathematically correct gcd in F after all).
If one adds a gcd and lcm for field elements (#9819) returning 0 or 1, then your suggestion certainly makes sense.
- trac has a
t
too much in line 1610 of ring.pyx:the case of the rational field. However, since tract ticket #10771,
- you write "quotient field" in the documentation. You could write "fraction field" to avoid confusion with quotient rings, which may be fields. This makes it more clear that you refer to the mathematical counterpart of Sage's "
FractionField
"?
Thanks! I am about to post a new patch, and I hope that Luis as author of #9819 does not mind if I add the case of arbitrary fields to my patch.
comment:7 follow-up: ↓ 11 Changed 8 years ago by
I just updated my patch, and I hope it addresses all your remarks. The typos are fixed, the old gcd of the rational field is now removed, and: I also added gcd/lcm for general fields. I am sorry that this makes #9819 a duplicate.
And I hope that I did not confuse gcd and lcm in the following setting:
sage: GF(2)(1).gcd(GF(2)(1)) 1 sage: GF(2)(1).gcd(GF(2)(0)) 1 sage: GF(2)(0).gcd(GF(2)(0)) 0 sage: GF(2)(1).lcm(GF(2)(0)) 0 sage: GF(2)(1).lcm(GF(2)(1)) 1
That's implemented as element methods for the category of Fields()
. Somehow, the category framework is cool, isn't it?
comment:8 follow-up: ↓ 9 Changed 8 years ago by
Oops, I forgot one aspect: If the gcd/lcm of the base ring raises an error, it should be caught. Implementing it now...
comment:9 in reply to: ↑ 8 Changed 8 years ago by
Replying to SimonKing:
Oops, I forgot one aspect: If the gcd/lcm of the base ring raises an error, it should be caught. Implementing it now...
Here it is:
sage: R.<q> = ZZ.extension(x^2+5) sage: gcd(q,q) 1 sage: gcd(q,0) 1 sage: gcd(R.zero(),0) 0 sage: lcm(q,q) 1 sage: lcm(q,0) 0
Correct me if I am wrong, but I think that now every complaint is addressed.
Best regards,
Simon
comment:10 follow-up: ↓ 22 Changed 8 years ago by
- Status changed from needs_review to needs_work
New complaints, sorry.
- I don't really like this:
sage: M = Matrix(GF(3),2,2,[1,2,0,1]); M [1 2] [0 1] sage: GF(2)(1).lcm(M) 1
I think you should first try to coerce other to self.parent() in gcd and lcm for Field elements.
- Line 64 of quotient_fields.py "
both GCD and LCM, it is possible to be a bit more specific
" uses capital letters GCD and LCM, while the methods that it (quasi-)refers to are lower case:AttributeError: 'sage.rings.rational.Rational' object has no attribute 'GCD'
- Next line: "define the GCD uniquely up to a unit in the base ring". Write "of" instead of "in" so that it is clear that you are not saying that the GCD is in the base ring.
- Your tests with
R.<q> = ZZ.extension(x^2+5)
is unrelated to fraction field: q is 1sage: R.<q> = ZZ.extension(x^2+5) sage: g = R.gens(); g [1, q] sage: q 1 sage: g[1] q sage: g[1] == q False
All your examples involving q get coerced to ZZ by the global GCD and LCM, so you should find other examples. This seems to work:sage: gcd(g[1],1) TypeError: unable to find gcd of q and 1 sage: gcd(g[1]/1,1) 1
(Correct behaviour both times)
Time to catch a bus, may have more comments later.
comment:11 in reply to: ↑ 7 ; follow-up: ↓ 12 Changed 8 years ago by
Replying to SimonKing:
I just updated my patch, and I hope it addresses all your remarks. The typos are fixed, the old gcd of the rational field is now removed, and: I also added gcd/lcm for general fields. I am sorry that this makes #9819 a duplicate.
Please.
And I hope that I did not confuse gcd and lcm in the following setting:
sage: GF(2)(1).gcd(GF(2)(1)) 1 sage: GF(2)(1).gcd(GF(2)(0)) 1 sage: GF(2)(0).gcd(GF(2)(0)) 0 sage: GF(2)(1).lcm(GF(2)(0)) 0 sage: GF(2)(1).lcm(GF(2)(1)) 1That's implemented as element methods for the category of
Fields()
. Somehow, the category framework is cool, isn't it?
I do not full understand this, but I am happy that works.
Some thoughts:
For QQ and the like, could it be that gcd and lcm should only take care of coercing to a good setting and then the real algorithm should be in _lcm, _gcd? Note that QQ already has ._gcd and ._lcm so these methods has to be taken into account.
For generic Fields, It should appear in the documentation that the methods return 0 if both arguments are zero and a non-zero element otherwise. The user should not suppose that the non-zero element is actually one, since this is changed by subclasses.
comment:12 in reply to: ↑ 11 Changed 8 years ago by
Hi Luis!
Replying to lftabera:
That's implemented as element methods for the category of
Fields()
. Somehow, the category framework is cool, isn't it?I do not full understand this, but I am happy that works.
Then let me explain:
Any category (in Sage) comes with a parent class and an "element class". That very often is an abstract thing, but you can inherit from it, and in particular you can implement methods for it.
If things are properly done, then a parent P is declared as an object in some category, and its class automatically(!) inherits from the parent class of its category:
sage: QQ.__class__ <class 'sage.rings.rational_field.RationalField_with_category'> sage: issubclass(QQ.__class__,QuotientFields().parent_class) True
If things are very properly done, then the parent has an attribute Element
, which is supposed to be a class. That class will (again automatically) be translated into a subclass of the element class of the parent's category, and is available as the attribute element_class
:
sage: m = matrix(GF(2),[[1]]) sage: G = MatrixGroup([m]) sage: G.Element <class 'sage.groups.matrix_gps.matrix_group_element.MatrixGroupElement'> sage: G.__class__ <class 'sage.groups.matrix_gps.matrix_group.MatrixGroup_gens_finite_field_with_category'> sage: G.element_class <class 'sage.groups.matrix_gps.matrix_group_element.MatrixGroup_gens_finite_field_with_category.element_class'> sage: issubclass(G.element_class,G.Element) True sage: issubclass(G.element_class,G.category().element_class) True sage: isinstance(G.random_element(),G.category().element_class) True
So, if things are very properly done, then a method defined as an "element method" of a category, will be available by looking up the method resolution order, since the elements of the parent are actual instances of the category's element class.
Sometimes things are done properly, but not very properly:
sage: QQ.category() Category of quotient fields sage: isinstance(QQ, QQ.category().parent_class) True sage: hasattr(QQ,'Element') False sage: isinstance(QQ.one(), QQ.category().element_class) False
Hence, the methods that I defined for the element class of the "category of quotient fields" is not available to elements of QQ
by simply looking up the method resolution order. But: There also is a __getattr__
method implemented for sage.structure.element.Element
, and this can access the element methods of the category of the parent of an element (uff!) even if the element is no instance of the "proper" element class.
This is why the new gcd works for the rationals, but it's slower than with the method resolution order, and this is why I also did not remove the existing lcm for rationals.
Some thoughts:
For QQ and the like, could it be that gcd and lcm should only take care of coercing to a good setting and then the real algorithm should be in _lcm, _gcd? Note that QQ already has ._gcd and ._lcm so these methods has to be taken into account.
That's a nice observation. By searching a little, I found that _gcd
is used in the method sage.structure.element.PrincipalIdealDomainElement.gcd
.
However, there are only very few classes that define a _gcd
method: sage.structure.element.EuclideanDomainElement
, sage.structure.element.FieldElement
, sage.rings.rational.Rational
, and then there are hits in ...polynomial_element_generic
and ...polynomial_modn_dense_ntl
.
In other words, I guess that the _gcd
method is an artefact of earlier attempts to reflect mathematics in the class hierarchy -- this should now be done in the category framework. _gcd
will only play a role if
- the class of an element derives from
PrincipalIdealDomainElement
, and
- the class of an element additionally defines
_gcd
.
It seems to me that there is precisely one class that directly inherits from PrincipalIdealDomainElement
, and there is precisely one class that directly inherits from that class:
sage: search_src("\(PrincipalIdealDomainElement") structure/element.pxd:83:cdef class EuclideanDomainElement(PrincipalIdealDomainElement): structure/element.pyx:2370:cdef class EuclideanDomainElement(PrincipalIdealDomainElement): sage: search_src("\(EuclideanDomainElement") structure/element.pyx:2362:PY_SET_TP_NEW(EuclideanDomainElement, Element) rings/integer.pxd:9:cdef class Integer(EuclideanDomainElement):
So, _gcd
is only used for the Integer
class. I guess it could easily be removed.
For generic Fields, It should appear in the documentation that the methods return 0 if both arguments are zero and a non-zero element otherwise. The user should not suppose that the non-zero element is actually one, since this is changed by subclasses.
Right. I'll take care of that, but now I have a bus to catch...
comment:13 follow-up: ↓ 14 Changed 8 years ago by
- Did you run a full doctest? I'm getting doctest failures that may have something to do with your patch, as they probably simplify fractions somewhere deep inside in some way. Or my installation is messed up again.
sage -t -long devel/sage/sage/symbolic/expression.pyx # 1 doctests failed sage -t -long devel/sage/sage/symbolic/integration/integral.py # 1 doctests failed sage -t -long devel/sage/sage/stats/basic_stats.py # 2 doctests failed
- The documentation of content for multivariate polynomials says "Returns the content of this polynomial.
Here, we define content as the gcd of the coefficients in the base ring."
Your changed doctest (1 becomes 2) is correct, but perhaps it would be more informative to
add
f.content().parent()
, and to give an example over another field. (Just a suggestion if you are editing the code anyway. Otherwise, don't bother.) More important:# Since trac ticket #10771, the gcd in QQ restricts to the # gcd in ZZ.
I don't think this comment is needed in this part of the code. If you want to include it, it would better fit in the doctest a few lines above it. And I think junk like the following should not be in the code at all:#,integer=self.parent() is ZZ)
comment:14 in reply to: ↑ 13 Changed 8 years ago by
Replying to mstreng:
- Did you run a full doctest? I'm getting doctest failures that may have something to do with your patch, as they probably simplify fractions somewhere deep inside in some way. Or my installation is messed up again.
sage -t -long devel/sage/sage/symbolic/expression.pyx # 1 doctests failed sage -t -long devel/sage/sage/symbolic/integration/integral.py # 1 doctests failed sage -t -long devel/sage/sage/stats/basic_stats.py # 2 doctests failed
Hum. Apparently I did not run these tests. The problem is that I currently work on various patches, frequently switch between them, and it may happen that I test one of them but forget the other.
Anyway, the failure in "expression.pyx" seems to be related with the patch. I'll try to understand what is happening.
- The documentation of content for multivariate polynomials says "Returns the content of this polynomial. Here, we define content as the gcd of the coefficients in the base ring." Your changed doctest (1 becomes 2) is correct, but perhaps it would be more informative to add
f.content().parent()
, and to give an example over another field. (Just a suggestion if you are editing the code anyway. Otherwise, don't bother.)
OK.
More important:
# Since trac ticket #10771, the gcd in QQ restricts to the # gcd in ZZ.I don't think this comment is needed in this part of the code. If you want to include it, it would better fit in the doctest a few lines above it.
I thought of it as a comment to people who know the original code. It could be that it is a method that uses gcd only internally, so that commenting in the code seems better than in the documentation.
And I think junk like the following should not be in the code at all:
#,integer=self.parent() is ZZ)
OK, that was an artefact of editing.
But I think most urgent are the doctest failures. It seems related to gcd for fields that are not fraction fields (such as RR
). That would also explain why all text passed with the original version of my patch, whereas the new version fails.
Cheers, Simon
comment:15 Changed 8 years ago by
Yep, I understand where the first error comes from.
The original gcd codes says:
if b is not None: if hasattr(a, "gcd"): return a.gcd(b, **kwargs) else: try: return ZZ(a).gcd(ZZ(b)) except TypeError: raise TypeError, "unable to find gcd of %s and %s"%(a,b)
Moreover, elements of RR
did not have a gcd method. Hence, it was tried to coerce them into ZZ, and so we had
sage: gcd(3.0,6.0) 3
But with my patch, all field elements (including elements of RR
) have a gcd method, so, a conversion to ZZ
is not attempted.
I am not sure yet how I would solve it: Change the code of `sage.symbolic.expression"? Change the gcd implemented for fields that are no fraction fields, such that a special case is made for inexact fields?
comment:16 Changed 8 years ago by
I suggest to do the following, which would generalize the current attempt to convert stuff into ZZ
:
For inexact fields, evaluation in the prime subfield is attempted:: sage: lcm(15.2,12.0); lcm(15.2,12.0).parent() 228 Rational Field sage: RR(76/5) 15.2000000000000 sage: lcm(76/5,12) 228 If this fails, we resort to the default we see above:: sage: lcm(6.0*CC.0,8*CC.0); lcm(6.0*CC.0,8*CC.0).parent() 1.00000000000000 Complex Field with 53 bits of precision
Similarly for gcd.
In that way, the first of the three doctest failures vanishes. But the other two remain problems.
comment:17 follow-up: ↓ 18 Changed 8 years ago by
The second failure can be triggered as follows:
sage: F(x) = 1/sqrt(2*pi*1^2)*exp(-1/(2*1^2)*(x-0)^2) sage: G(x) = 1/sqrt(2*pi*n(1)^2)*exp(-1/(2*n(1)^2)*(x-n(0))^2) sage: ((F(x)-G(x))^2).expand() /mnt/local/king/SAGE/sage-4.6.2.alpha4/local/bin/sage-sage: Zeile 300: 18366 Speicherzugriffsfehler sage-ipython "$@" -i
Segmentation fault while expanding. I could imagine that gcd==1
results in an infinite loop.
comment:18 in reply to: ↑ 17 Changed 8 years ago by
Replying to SimonKing:
Segmentation fault while expanding. I could imagine that
gcd==1
results in an infinite loop.
Ouch, it is worse: It already occurs before the expansion!
comment:19 Changed 8 years ago by
The third error boils down to:
Without patch
sage: x = -sqrt(2)-1/5*I sage: x*x 1/25*(5*sqrt(2) + I)^2
With patch
sage: x = -sqrt(2)-1/5*I sage: x*x 1/25*(-5*sqrt(2) - I)^2
The segfault in the second doctest problem also comes from multiplication of symbolic expressions. And the example above indicates that it is related with pulling common factors out of a list of expressions.
comment:20 Changed 8 years ago by
Just FYI, #8111 is possibly also related.
comment:21 follow-up: ↓ 23 Changed 8 years ago by
- Cc burcin added
- Status changed from needs_work to needs_review
With the new patch, all tests pass for me. The behaviour is now (mutatis mutandis for lcm):
- If
ZZ
is a subring of a field F that is not a fraction field thena.gcd(b)
first tests if a and b can be converted intoZZ
. If it is the case, their gcd inZZ
as element ofZZ
is returned. This is not nice, but it is compatible with the "unpatched" behaviour of Sage and is assumed in some other parts of Sage:# unpatched! sage: gcd(RR(1),RR(1)) 1 sage: _.parent() Integer Ring
- If
ZZ
is no subring of F or conversion fails thena.gcd(b)
returnsF(0)
orF(1)
. This is, e.g., the case forgcd(CC(2),CC(2))
, sinceCC(2)
can not be converted intoZZ
.
- a.gcd(b) only raises an error if
b
can not be converted intoa.parent()
.
- gcd and lcm in fraction fields of a unique factorization domain R restrict to R, and moreover they satisfy
gcd(a,b)*lcm(a,b)==a*b
up to multiplication with a unit of R.
There is a price to pay, and that's why I put Burcin on Cc:
The distribution of minus signs in symbolic expression changes, apparently since now gcd raises no errors in examples like the following:
sage: (I - 1/3*sqrt(2))^2 1/9*(-sqrt(2) + 3*I)^2 # Without patch, we would get # 1/9*(sqrt(2) - 3*I)^2
The question is whether a changed minus sign distribution is such a serious thing.
comment:22 in reply to: ↑ 10 Changed 8 years ago by
Hi Marco,
sorry, I forgot to explain how I addressed your complaints:
Replying to mstreng:
New complaints, sorry.
- I don't really like this:
sage: M = Matrix(GF(3),2,2,[1,2,0,1]); M [1 2] [0 1] sage: GF(2)(1).lcm(M) 1I think you should first try to coerce other to self.parent() in gcd and lcm for Field elements.
Yep, that was a mistake. We now have:
sage: M = Matrix(GF(3),2,2,[1,2,0,1]); M [1 2] [0 1] sage: GF(2)(1).lcm(M) --------------------------------------------------------------------------- ArithmeticError Traceback (most recent call last) ... ArithmeticError: The second argument can not be interpreted in the parent of the first argument. Can't compute the lcm sage: GF(2)(1).gcd(M) ERROR: An unexpected error occurred while tokenizing input The following traceback may be corrupted or invalid The error message is: ('EOF in multi-line statement', (4, 0)) --------------------------------------------------------------------------- ArithmeticError Traceback (most recent call last) ... ArithmeticError: The second argument can not be interpreted in the parent of the first argument. Can't compute the gcd
- Line 64 of quotient_fields.py "
both GCD and LCM, it is possible to be a bit more specific
" uses capital letters GCD and LCM, while the methods that it (quasi-)refers to are lower case:
I changed the doc string to lower case.
- Next line: "define the GCD uniquely up to a unit in the base ring". Write "of" instead of "in" so that it is clear that you are not saying that the GCD is in the base ring.
Done.
- Your tests with
R.<q> = ZZ.extension(x^2+5)
is unrelated to fraction field:
You are right. My intention was to start with R
, show that there is no gcd in R, and then expose what happens in Frac(R)
. The test is changed accordingly.
Best regards,
Simon
comment:23 in reply to: ↑ 21 ; follow-up: ↓ 24 Changed 8 years ago by
There is a price to pay, and that's why I put Burcin on Cc:
The distribution of minus signs in symbolic expression changes, apparently since now gcd raises no errors in examples like the following:
sage: (I - 1/3*sqrt(2))^2 1/9*(-sqrt(2) + 3*I)^2 # Without patch, we would get # 1/9*(sqrt(2) - 3*I)^2The question is whether a changed minus sign distribution is such a serious thing.
What happens with the already existing doctests? That would be what I would be most concerned about. Also, is the infinite loop now gone? IIRC you had fixed that somehow.
Really, in some sense the new version is 'better' because it keeps the minus where the user put it. Do (-I - 1/3*sqrt(2))^2
and (I + 1/3*sqrt(2))^2
do anything different from before?
comment:24 in reply to: ↑ 23 ; follow-up: ↓ 25 Changed 8 years ago by
Replying to kcrisman:
What happens with the already existing doctests? That would be what I would be most concerned about.
There are two tests (really in fact only one) in sage.stats
where I had to change the expected output. The old and new output are equivalent:
sage: std([I, sqrt(2), 3/5]) # result with patch sqrt(1/450*(-5*sqrt(2) + 10*I - 3)^2 + 1/450*(-5*sqrt(2) - 5*I + 6)^2 + 1/450*(10*sqrt(2) - 5*I - 3)^2) # Without patch, we had sqrt(1/450*(5*sqrt(2) + 5*I - 6)^2 + 1/450*(5*sqrt(2) - 10*I + 3)^2 + 1/450*(10*sqrt(2) - 5*I - 3)^2) # Testing equality: sage: bool(sqrt(1/450*(5*sqrt(2) + 5*I - 6)^2 + 1/450*(5*sqrt(2) - 10*I + 3)^2 + 1/450*(10*sqrt(2) - 5*I - 3)^2) == std([I, sqrt(2), 3/5])) True
There is essentially the same happening with variance([I, sqrt(2), 3/5])
.
Also, is the infinite loop now gone? IIRC you had fixed that somehow.
Yep. The problem was that in an old version of my patch I would not try conversion to ZZ
but conversion to the prime subfield. A real number with finitely many digits can of course always be converted into RR.prime_subfield()
, which is QQ
. So, I guess that the segfault came from rationals with very large numerator and denominator. Here is the proof that the segfault has gone:
sage: F(x) = 1/sqrt(2*pi*1^2)*exp(-1/(2*1^2)*(x-0)^2) sage: G(x) = 1/sqrt(2*pi*n(1)^2)*exp(-1/(2*n(1)^2)*(x-n(0))^2) sage: (F-G)**2 x |--> 1/4*(sqrt(2)*e^(-1/2*x^2)/sqrt(pi) - 1.41421356237309*e^(-0.500000000000000*x^2)/sqrt(pi))^2
which coincides with the answer of unpatched Sage.
Really, in some sense the new version is 'better' because it keeps the minus where the user put it. Do
(-I - 1/3*sqrt(2))^2
and(I + 1/3*sqrt(2))^2
do anything different from before?
They do.
sage: (-I - 1/3*sqrt(2))^2 1/9*(-sqrt(2) - 3*I)^2 # Was: 1/9*(sqrt(2) + 3*I)^2 sage: (I + 1/3*sqrt(2))^2 1/9*(sqrt(2) + 3*I)^2 # Was: 1/9*(sqrt(2) + 3*I)^2
Really I don't understand where that comes from. Namely, as far as I know, the sign of gcd or lcm did not change. The only difference in gcd/lcm examples involving I
and sqrt(2)
is the fact that without my patch an error would be raised.
comment:25 in reply to: ↑ 24 Changed 8 years ago by
Really, in some sense the new version is 'better' because it keeps the minus where the user put it. Do
(-I - 1/3*sqrt(2))^2
and(I + 1/3*sqrt(2))^2
do anything different from before?They do. sage: (-I - 1/3*sqrt(2))^{2 1/9*(-sqrt(2) - 3*I)}2 # Was: 1/9*(sqrt(2) + 3*I)^{2 }
Okay, that's what I figured would happen.
sage: (I + 1/3*sqrt(2))^{2 1/9*(sqrt(2) + 3*I)}2 # Was: 1/9*(sqrt(2) + 3*I)^{2 }
Okay, that is the same as before in any case.
Really I don't understand where that comes from. Namely, as far as I know, the sign of gcd or lcm did not change. The only difference in gcd/lcm examples involving
I
andsqrt(2)
is the fact that without my patch an error would be raised.
Yes,
sage: gcd(I + 1/3*sqrt(2),I + 1/3*sqrt(2)) --------------------------------------------------------------------------- RuntimeError Traceback (most recent call last)
Interesting. Like I said earlier, Sage now sends such symbolic things to Pynac/Ginac?, which however needs to use Sage again in order to do certain calculations like this (Mike Hansen explained it well in your thread on sage-devel). So apparently in the past it decided to factor out a -1 in certain cases based on some gcd - maybe just a default of 1 or something - whereas now it does not, due to the new gcd. What does
sage: gcd(I + 1/3*sqrt(2),I + 1/3*sqrt(2))
do after the patch? What happens with this?
sage: gcd(4,2+2*I*sqrt(3))
Notice that I am purposely taking two things NOT in a polynomial or extension thing, where it can be calculated whether it's a UFD. In the symbolic ring, who knows what's "right"?
comment:26 Changed 8 years ago by
I think I've figured out hg queues now, so I can actually try out the patch! [No guarantees, though.] Was fuzz-testing it and noticed the following:
sage: var("x") x sage: a = 1/3*x**0 sage: b = 2/5*x**0 sage: type(a) <type 'sage.symbolic.expression.Expression'> sage: gcd(a,b) 1/15 sage: lcm(a,b) --------------------------------------------------------------------------- TypeError Traceback (most recent call last) /Applications/sage_devel/devel/sage-main/<ipython console> in <module>() /Applications/sage_devel/local/lib/python2.6/site-packages/sage/rings/arith.pyc in lcm(a, b) 1562 return ZZ(a).lcm(ZZ(b)) 1563 except TypeError: -> 1564 raise TypeError, "unable to find lcm of %s and %s"%(a,b) 1565 return LCM(b) 1566 TypeError: unable to find lcm of 1/3 and 2/5
All bets are off in SR, but there is an Expression.gcd (no .lcm, strangely). Now that QQ reduces to ZZ, maybe switch the attempt to try ZZ to QQ, so that rationals Just Work(tm)?
This surprised me a little too:
sage: R.<x>=QQ[] sage: sage: a = 1/3*x**0 sage: b = 2/5*x**0 sage: type(a) <type 'sage.rings.polynomial.polynomial_rational_flint.Polynomial_rational_flint'> sage: gcd(a,b), lcm(a,b) (1, 1) sage: sage: a = (1*x**0)/(3*x**0) sage: b = (2*x**0)/(5*x**0) sage: type(a) <class 'sage.rings.fraction_field_element.FractionFieldElement_1poly_field'> sage: gcd(a,b), lcm(a,b) (1, 1)
given that the seemingly harder cases
sage: a = 1/3+x sage: b = 2/5+x sage: type(a) <type 'sage.rings.polynomial.polynomial_rational_flint.Polynomial_rational_flint'> sage: gcd(a,b), lcm(a,b), a*b-gcd(a,b)*lcm(a,b) (1, x^2 + 11/15*x + 2/15, 0) sage: sage: a = (1*x**0)/(3*x**0)+x sage: b = (2*x**0)/(5*x**0)+x sage: type(a) <class 'sage.rings.fraction_field_element.FractionFieldElement_1poly_field'> sage: gcd(a,b), lcm(a,b), a*b-gcd(a,b)*lcm(a,b) (1, x^2 + 11/15*x + 2/15, 0)
behave a little nicer. Is it possible to recover the rational results?
comment:27 follow-up: ↓ 28 Changed 8 years ago by
Ah, I just read the comments in #8111 (re: having to dive into flint). Well, at least changing ZZ to QQ in rings.arith.lcm would be straightforward.
comment:28 in reply to: ↑ 27 Changed 8 years ago by
Replying to dsm:
Ah, I just read the comments in #8111 (re: having to dive into flint). Well, at least changing ZZ to QQ in rings.arith.lcm would be straightforward.
It could be that it is not straight forward. If my analysis of the segfault occuring with a previous version of my patch is correct, changing ZZ
to QQ
could very well trigger the segfault again.
comment:29 follow-up: ↓ 31 Changed 8 years ago by
- Owner changed from AlexGhitza to dsm
Well, even after changing ZZ to QQ there everything seems to work as expected for me (no segfaults, a quick subset of tests all passed; I'll try testall long later and look for more crashers), and if it did explode it sort of feels like it'd be due to a bug somewhere else that it revealed.. but that can wait for another ticket later.
The real meat is already in your patch, and I'm already using it, so I'm very grateful for the work, and hope someone gives the okay soon! :^)
comment:30 Changed 8 years ago by
- Owner changed from dsm to AlexGhitza
comment:31 in reply to: ↑ 29 ; follow-up: ↓ 32 Changed 8 years ago by
Replying to dsm:
Well, even after changing ZZ to QQ there everything seems to work as expected for me (no segfaults, a quick subset of tests all passed; I'll try testall long later and look for more crashers), and if it did explode it sort of feels like it'd be due to a bug somewhere else that it revealed.. but that can wait for another ticket later.
What came out of the tests?
The real meat is already in your patch, and I'm already using it, so I'm very grateful for the work, and hope someone gives the okay soon!
:^)
You (after finishing the tests, of course)?? Or you could post a patch for changing ZZ
to QQ
in rings.arith.lcm
, and perhaps I can have a look on it? If I remember correctly, it is ok if you are reviewer for my patch and I am reviewer for your patch.
comment:32 in reply to: ↑ 31 Changed 8 years ago by
Replying to SimonKing:
What came out of the tests?
I seem to recall it passed just fine but I can't find the output now :-(
so I'm rerunning them (takes me about ~3h).
You (after finishing the tests, of course)??
Following was's recommendations (http://ask.sagemath.org/question/31/should-i-_really_-review-trac-tickets) I don't think I've got enough experience with the Sage codebase yet to set something to positive review; I'm still finding my way around.
comment:33 Changed 8 years ago by
- Reviewers changed from Marco Streng to Marco Streng, Mariah Lenox
- Status changed from needs_review to positive_review
Applied patch to sage-4.7.1.alpha2, did 'sage -b' followed by 'make testlong'. All tests passed. Positive review!
comment:34 Changed 8 years ago by
<Mr Burns> Excellent. </Burns> Many thanks to Simon for the work!
comment:35 Changed 8 years ago by
- Milestone changed from sage-4.7.1 to sage-4.7.2
- Type changed from defect to enhancement
comment:36 follow-up: ↓ 37 Changed 8 years ago by
- Status changed from positive_review to needs_work
Could you please change the commit message such that the lines are not so long? You should wrap long lines but make sure the first still makes sense by itself, that will appear in hg log
.
Changed 8 years ago by
Implement gcd and lcm for general fields, with the special case of the fraction field of a unique factorization domain
comment:37 in reply to: ↑ 36 Changed 8 years ago by
- Status changed from needs_work to positive_review
Replying to jdemeyer:
Could you please change the commit message such that the lines are not so long? You should wrap long lines but make sure the first still makes sense by itself, that will appear in
hg log
.
Admittedly the log line was "a bit" two long. It is now spread over three lines. I hope it is fine now.
comment:38 Changed 8 years ago by
- Merged in set to sage-4.7.2.alpha0
- Resolution set to fixed
- Status changed from positive_review to closed
Without the patch, we had:
With the patch, one has