Opened 5 years ago
Closed 5 years ago
#20330 closed defect (fixed)
hyperbolic_geodesic midpoint bugfix
Reported by: | jhonrubia6 | Owned by: | |
---|---|---|---|
Priority: | major | Milestone: | sage-7.3 |
Component: | geometry | Keywords: | hyperbolic geometry, geodesic |
Cc: | Merged in: | ||
Authors: | Javier Honrubia González | Reviewers: | Vincent Delecroix |
Report Upstream: | N/A | Work issues: | |
Branch: | 2ba1c96 (Commits, GitHub, GitLab) | Commit: | 2ba1c96e25d7e128adc832bc4c4e4c51e8835b98 |
Dependencies: | Stopgaps: |
Description (last modified by )
midpoint() fails to end for some geodesic with finite length example:
sage: g=HyperbolicPlane().UHP().get_geodesic(2*I,2*I+1) sage: g.midpoint()' ... RuntimeError: maximum recursion depth exceeded in __instancecheck__
Change History (51)
comment:1 Changed 5 years ago by
- Component changed from PLEASE CHANGE to geometry
- Keywords hyperbolic geometry geodesic added
- Type changed from PLEASE CHANGE to defect
comment:2 Changed 5 years ago by
- Description modified (diff)
comment:3 Changed 5 years ago by
- Description modified (diff)
comment:4 Changed 5 years ago by
comment:5 Changed 5 years ago by
It seems that the matrix is sufficiently complicated that it just runs up against the Python function call limit. If I run simplify_full
on the matrix first, then it works (and I can even run in again to simplify it further). So this would be one way to hack around this.
Although we might benefit from having a special case for inverting 2x2 matrices using the explicit form:
[A B]^-1 = __1__ [ D -B] [C D] AD-BC [-C A]
This would give another way around this (it should be faster and in some ways, more robust) as this code only deals with 2x2 matrices.
comment:6 Changed 5 years ago by
I had tried simplify_full
(although I hadn't thought of applying it multiple times) but I am somehow disappointed with the solution even that it gets since the result is still too complicated (in the motivating example should have Re(z)=1/2 instead of the monster that it shows.
I am trying to understand the code to see how to get a simpler expression, if possible without resorting to simplify_full
the result.
The transformation S it uses to carry the geodesic to the imaginary axis it is not even recognized as an isometry by the model, fails the test self._model.isometry_in_model(S)
(although the 'full_simplified' expression successes)
Obviously works flawlessly with
h=HyperbolicPlane().UHP().get_geodesic(2.0*I,2.0*I+1.0) h.midpoint()
comment:7 Changed 5 years ago by
- Branch set to u/jhonrubia6/hyperbolic_geodesic_midpoint_bugfix
comment:8 Changed 5 years ago by
- Commit set to 89500eb900f245d091afcc4461f4d8f19f6db5b9
- Status changed from new to needs_review
The matrix of the isometry that sends the geodesic to the imaginary axis is a symbolic array, when this matrix is numeric the algorithm evaluates efficiently but when is purely symbolic it is necessary a double full_simplify() operation so that the posterior operation of inverting the matrix does not get a runtime error. After profiling the code, I do not see (surprisingly) any speed advantage in coding explicitly the inverse matrix, so I decided to left that unchanged.
New commits:
89500eb | If the S matrix is not numeric, apply full_simplify() in order to avoid the runtime error.
|
comment:9 Changed 5 years ago by
comment:10 follow-up: ↓ 12 Changed 5 years ago by
Hello,
Why in all this code the type of number is not preserved
sage: g = HyperbolicPlane().UHP().get_geodesic(CC(0,2),CC(1,2)) sage: type(g.endpoints()[0].coordinates()) <type 'sage.rings.complex_number.ComplexNumber'> sage: type(g.midpoint().coordinates()) <type 'sage.symbolic.expression.Expression'>
If I am using floating point number it is likely that I want to do floating point computations.
comment:11 follow-up: ↓ 13 Changed 5 years ago by
Related question: in your code your assume that the coordinates are symbolic. Why can you make such assumption (since I am able to construct geodesic with any kind of coordinates)?
comment:12 in reply to: ↑ 10 Changed 5 years ago by
- Status changed from needs_review to needs_work
Replying to vdelecroix:
Hello,
Why in all this code the type of number is not preserved
sage: g = HyperbolicPlane().UHP().get_geodesic(CC(0,2),CC(1,2)) sage: type(g.endpoints()[0].coordinates()) <type 'sage.rings.complex_number.ComplexNumber'> sage: type(g.midpoint().coordinates()) <type 'sage.symbolic.expression.Expression'>If I am using floating point number it is likely that I want to do floating point computations.
Since I was trying to fix a bug, I did not question the way the module was implemented in the beginning (I assumed it was discussed, and reached a consensus).
Taking a closer look, the effect you mention is caused by _crossratio_matrix()
as in
sage: (p1, p2, p3) = [UHP.random_point().coordinates() for k in range(3)] sage: A = HyperbolicGeodesicUHP._crossratio_matrix(p1, p2, p3) sage: type(A) <type 'sage.matrix.matrix_symbolic_dense.Matrix_symbolic_dense'>
even
sage: B = HyperbolicGeodesicUHP._crossratio_matrix(CC(1,1), CC(3,2), CC(2,2)) sage: type(B) <type 'sage.matrix.matrix_generic_dense.Matrix_generic_dense'>
One thing we can try is to check the type of p1,p2,p3 and construct the matrix in CDF if we have isinstance(p1,numbers.Complex)
using
return matrix(CDF,[[p1 - p2, (p1 - p2)*(-p0)],[p1 - p0, (p1 - p0)*(-p2)]])
What do you think of this approach? Any pitfalls ahead I am not aware of?
comment:13 in reply to: ↑ 11 Changed 5 years ago by
Replying to vdelecroix:
Related question: in your code your assume that the coordinates are symbolic. Why can you make such assumption (since I am able to construct geodesic with any kind of coordinates)?
Well, it seems that all the code is working on symbolics as in
sage: UHP = HyperbolicPlane().UHP() sage: P = UHP.random_point().coordinates() sage: type(P) <type 'sage.symbolic.expression.Expression'>
comment:14 follow-up: ↓ 16 Changed 5 years ago by
On a related topic. If I check with number.Complex I get unwanted results as in
sage: import numbers sage: from sage.rings.complex_number import ComplexNumber sage: Q = UHP.get_point(CC(sqrt(2)*(1+I))) sage: Q.coordinates() 1.41421356237309 + 1.41421356237309*I sage: isinstance(Q.coordinates(),numbers.Complex) False sage: isinstance(Q.coordinates(),ComplexNumber) True
comment:15 follow-up: ↓ 17 Changed 5 years ago by
I correct myself the crossratio matrix is
<type 'sage.matrix.matrix_generic_dense.Matrix_generic_dense'>
it gets converted to symbolic when multiplied by
B = matrix([[1, 0], [0,I]])
so if we change it to
matrix([[1r,0r],[0r,CC(0,1)]])
we mantain the type to
<type 'sage.matrix.matrix_generic_dense.Matrix_generic_dense'>
It is what you are requesting or it you require to get a
<type 'sage.matrix.matrix_complex_double_dense.Matrix_complex_double_dense'>
In that case, will we have to construct
return matrix(CDF,[[p1 - p2, (p1 - p2)*(-p0)],[p1 - p0, (p1 - p0)*(-p2)]])
as mentioned before when the points p1,p2,p3 are instances of ComplexNumber and leaving as it is in the rest of cases?
comment:16 in reply to: ↑ 14 ; follow-up: ↓ 18 Changed 5 years ago by
Replying to jhonrubia6:
On a related topic. If I check with number.Complex I get unwanted results as in
sage: import numbers sage: from sage.rings.complex_number import ComplexNumber sage: Q = UHP.get_point(CC(sqrt(2)*(1+I))) sage: Q.coordinates() 1.41421356237309 + 1.41421356237309*I sage: isinstance(Q.coordinates(),numbers.Complex) False sage: isinstance(Q.coordinates(),ComplexNumber) True
Which version of Sage are you using? On sage-7.2.beta6 complex numbers are perfectly recognized by Python numbers.Complex
:
sage: import numbers sage: isinstance(CC(0,1), numbers.Complex) True sage: isinstance(CDF(0,1), numbers.Complex) True
comment:17 in reply to: ↑ 15 ; follow-up: ↓ 19 Changed 5 years ago by
Replying to jhonrubia6:
I correct myself the crossratio matrix is
<type 'sage.matrix.matrix_generic_dense.Matrix_generic_dense'>
it gets converted to symbolic when multiplied byB = matrix([[1, 0], [0,I]])
so if we change it tomatrix([[1r,0r],[0r,CC(0,1)]])
we mantain the type to<type 'sage.matrix.matrix_generic_dense.Matrix_generic_dense'>
It is what you are requesting or it you require to get a<type 'sage.matrix.matrix_complex_double_dense.Matrix_complex_double_dense'>
In that case, will we have to construct
return matrix(CDF,[[p1 - p2, (p1 - p2)*(-p0)],[p1 - p0, (p1 - p0)*(-p2)]])as mentioned before when the points p1,p2,p3 are instances of ComplexNumber and leaving as it is in the rest of cases?
It would be bad to use CC(0,1)
by default for the complex imaginary unit. The only viable solution is to use the base ring of the coordinates
def get_I(a): from sage.structure.element import Element if isinstance(a, complex): return complex("j") elif isinstance(a, Expression): return SR("I") elif isinstance(a, Element): return a.parent().gen() else: raise ValueError("not a complex number")
It seems to work in most cases
sage: get_I(SR(1)) I sage: parent(_) Symbolic Ring sage: get_I(QQbar(1+I)) I sage: parent(_) Algebraic Field sage: get_I(CDF.an_element()) 1.0*I sage: parent(_) Complex Double Field sage: get_I(ComplexField(256).an_element()) 1.000000000000000000000000000000000000000000000000000000000000000000000000000*I sage: parent(_) Complex Field with 256 bits of precision
comment:18 in reply to: ↑ 16 Changed 5 years ago by
Replying to vdelecroix:
Replying to jhonrubia6:
On a related topic. If I check with number.Complex I get unwanted results as in
sage: import numbers sage: from sage.rings.complex_number import ComplexNumber sage: Q = UHP.get_point(CC(sqrt(2)*(1+I))) sage: Q.coordinates() 1.41421356237309 + 1.41421356237309*I sage: isinstance(Q.coordinates(),numbers.Complex) False sage: isinstance(Q.coordinates(),ComplexNumber) TrueWhich version of Sage are you using? On sage-7.2.beta6 complex numbers are perfectly recognized by Python
numbers.Complex
:sage: import numbers sage: isinstance(CC(0,1), numbers.Complex) True sage: isinstance(CDF(0,1), numbers.Complex) True
I am not at my installation right now (I'll check later). I tried the previous code in cloud.sagemath.com
sage: import numbers sage: isinstance(CC(0,1), numbers.Complex) False sage: isinstance(CDF(0,1), numbers.Complex) False
But then, maybe cloud.sagemath.com is not in the appropiate version. As I said I'll check later
comment:19 in reply to: ↑ 17 Changed 5 years ago by
Replying to vdelecroix:
Replying to jhonrubia6:
I correct myself the crossratio matrix is
<type 'sage.matrix.matrix_generic_dense.Matrix_generic_dense'>
it gets converted to symbolic when multiplied byB = matrix([[1, 0], [0,I]])
so if we change it tomatrix([[1r,0r],[0r,CC(0,1)]])
we mantain the type to<type 'sage.matrix.matrix_generic_dense.Matrix_generic_dense'>
It is what you are requesting or it you require to get a<type 'sage.matrix.matrix_complex_double_dense.Matrix_complex_double_dense'>
In that case, will we have to construct
return matrix(CDF,[[p1 - p2, (p1 - p2)*(-p0)],[p1 - p0, (p1 - p0)*(-p2)]])as mentioned before when the points p1,p2,p3 are instances of ComplexNumber and leaving as it is in the rest of cases?
It would be bad to use
CC(0,1)
by default for the complex imaginary unit. The only viable solution is to use the base ring of the coordinatesdef get_I(a): from sage.structure.element import Element if isinstance(a, complex): return complex("j") elif isinstance(a, Expression): return SR("I") elif isinstance(a, Element): return a.parent().gen() else: raise ValueError("not a complex number")It seems to work in most cases
sage: get_I(SR(1)) I sage: parent(_) Symbolic Ring sage: get_I(QQbar(1+I)) I sage: parent(_) Algebraic Field sage: get_I(CDF.an_element()) 1.0*I sage: parent(_) Complex Double Field sage: get_I(ComplexField(256).an_element()) 1.000000000000000000000000000000000000000000000000000000000000000000000000000*I sage: parent(_) Complex Field with 256 bits of precision
Understood. I appreciate your insight, I have not the experience in Sage programming to foresee these subleties. I'll introduce your code in the next patch
comment:20 Changed 5 years ago by
- Commit changed from 89500eb900f245d091afcc4461f4d8f19f6db5b9 to afdf46e9237f644889eb993a6236c31e91c6cf4d
Branch pushed to git repo; I updated commit sha1. New commits:
afdf46e | Modified the way it calculates the B matrix that transforms (0,1,inf)->(0,I,inf) trying to preserve the type of the calculations
|
comment:21 Changed 5 years ago by
- Reviewers set to Vincent Delecroix
- The following syntax is wrong
TESTS:: bla bla. sage: one_test()
It should beTESTS: bla bla.:: sage: one_test()
- You should get rid of trailing whitespaces, that is line that only contains spaces.
- The
.is_numeric()
method only applies to symbolic expressions. There is something weird with your codeif S[0,0].is_numeric()
since you might also want to manipulate floating point objects.
comment:22 Changed 5 years ago by
- Commit changed from afdf46e9237f644889eb993a6236c31e91c6cf4d to 663a76cb67fb1598301ce842dccd8341ab07218b
Branch pushed to git repo; I updated commit sha1. New commits:
663a76c | Removed trailing blank spaces. Fixed the TEST syntax. Modified the matrix element check for symbolic matrix check
|
comment:23 Changed 5 years ago by
I've also modified the doctoring in _to_std_geod()
method since it randomly failed (depending on the random UHP points, for example:
sage: UHP = HyperbolicPlane().UHP() sage: p1=0.869685163129173 + 0.897487924526729*I sage: p2=5.71917016919078 + 7.09666543326670*I sage: g = UHP.get_geodesic(p1, p2) sage: p= g.midpoint().coordinates() sage: A = g._to_std_geod(p);A # Send midpoint to I. sage: A = UHP.get_isometry(A) sage: [s, e]= g.complete().endpoints();s;e;p sage: bool(abs(A(s).coordinates()) < 10**-9) True sage: bool(abs(A(g.midpoint()).coordinates() - I) < 10**-9) True sage: bool(A(e).coordinates() == infinity) False sage: bool(abs(A(e).coordinates())>10**9) True
comment:24 Changed 5 years ago by
- Status changed from needs_work to needs_review
comment:25 Changed 5 years ago by
- the syntax is for refering to trac tickets is
:trac:`number`
- I don't understand your code. Instead of copy/paste you would better do
if isinstance(S, Matrix_symbolic_dense): S = S.simplify_full() S_1 = S.inverse() ... etc ...
comment:26 Changed 5 years ago by
You might want to add some doctest similar to the following to check that floating points remain floating points in any method.
sage: g = UHP.get_geodesic(CC(0,1), CC(2,2)) sage: UHP.dist(g.start(), g.end()) 1.45057451382258 sage: parent(_) Real Field with 53 bits of precision sage: g.midpoint() Point in UHP 0.666666666666666 + 1.69967317119759*I sage: parent(g.midpoint().coordinates()) Complex Field with 53 bits of precision sage: gc = g.complete() Geodesic in UHP from -0.265564437074637 to 3.76556443707464 sage: parent(gc.start().coordinates()) Real Field with 53 bits of precision sage: parent(gc._to_std_geod(g.start().coordinates())) Full MatrixSpace of 2 by 2 dense matrices over Complex Field with 53 bits of precision
comment:27 Changed 5 years ago by
Ok. I will try, but would it not be better to get its own trac ticket for that?
comment:28 Changed 5 years ago by
- Commit changed from 663a76cb67fb1598301ce842dccd8341ab07218b to 62159c386ed9ea95663b50606ae770f3e1afdbee
Branch pushed to git repo; I updated commit sha1. New commits:
62159c3 | Code simplified as noted by reviewer. docstring syntax error corrected. New test regarding floating points added in methods.
|
comment:29 Changed 5 years ago by
- Status changed from needs_review to needs_work
Thanks for incorporing some (unrelated) changes for floating point testing.
There are still some issues with documentation
- The code in the tests/examples must be indented
TESTS: This is a test:: sage: 1+1 2
orTESTS:: sage: 1+1 2
(the same holds for EXAMPLES)
- before any test there should be two colons (you forgot some at one place)
comment:30 Changed 5 years ago by
- Commit changed from 62159c386ed9ea95663b50606ae770f3e1afdbee to 929c7c73d9968adb10414dedf879060eb235728f
Branch pushed to git repo; I updated commit sha1. New commits:
929c7c7 | Revised rst syntax. Tried to adhere to pep8 recommendations. Some long lines reaming though
|
comment:31 Changed 5 years ago by
- Status changed from needs_work to needs_review
I tried to fix long lines too (as noted by flake8) , but there are some I do not find a way to reduce, so I left them.
comment:32 Changed 5 years ago by
- Status changed from needs_review to needs_work
Still some documentation issues
:meth:`midpoint()`
should be:meth`midpoint`
.
- before any tests there must be two colons. There is at least one missing occurrence after
... floating points in :meth:`dist()`.
comment:33 Changed 5 years ago by
- Commit changed from 929c7c73d9968adb10414dedf879060eb235728f to dd193839547305ed790c523e905c3e3f63ce563a
Branch pushed to git repo; I updated commit sha1. New commits:
dd19383 | Removed trailing blanks, added double colons before sage code, and revised indentation.
|
comment:34 Changed 5 years ago by
- Status changed from needs_work to needs_review
I think I got them all this time. Used pep8 checkers and a rest online checker.
comment:35 Changed 5 years ago by
Why did you do these kind of changes
-#*********************************************************************** +# ***********************************************************************
comment:36 Changed 5 years ago by
The line #********
produced an error
E265 block comment should start with '# '
in the PEP8 check.
The pep8 recommendation says
"Each line of a block comment starts with a # and a single space (unless it is indented text inside the comment).
Paragraphs inside a block comment are separated by a line containing a single # ."
comment:37 Changed 5 years ago by
As you noticed these are recommendations. Everywhere else in Sage, there is no space after the comment (especially in banners). It has to be consistent within Sage source code before being pep8 complient.
comment:38 Changed 5 years ago by
- Commit changed from dd193839547305ed790c523e905c3e3f63ce563a to 8fea9c2fd8182498e593fbd2e03e151ebefdbea0
Branch pushed to git repo; I updated commit sha1. New commits:
8fea9c2 | Removed blank spaces after #
|
comment:39 Changed 5 years ago by
Ok. I reverted those changes. What do you think about adding pictures illustrating the geodesics?. If so, I think we have to rename show
methods by plot
for sphinx_plot
to work.
comment:40 Changed 5 years ago by
It would be nice be open a new ticket for that purpose.
comment:41 Changed 5 years ago by
Ok. I've opened #20530 for that.
comment:42 Changed 5 years ago by
The method get_B
would be simpler and safer as follows
from sage.structure.element import Element from sage.symbolic.expression import Expression if isinstance(a, (int,float,complex)): # Python number a = CDF(a) if isinstance(a, Expression): # symbolic P = SR zero = SR.zero() one = SR.one() I = SR("I") elif isinstance(a, Element): # Sage number P = a.parent() zero = P.zero() one = P.one() I = P.gen() if I.is_one() or (I*I).is_one() or not (-I*I).is_one(): raise ValueError("invalid number") else: raise ValueError("not a complex number") return matrix(P, 2, [one, zero, zero, -I])
comment:43 Changed 5 years ago by
- Commit changed from 8fea9c2fd8182498e593fbd2e03e151ebefdbea0 to 2ba1c96e25d7e128adc832bc4c4e4c51e8835b98
Branch pushed to git repo; I updated commit sha1. New commits:
2ba1c96 | Repleaced the code for _get_B()
|
comment:44 Changed 5 years ago by
Ok. It seems nicer coding. I think you have contributed more code than me to this ticket Maybe we should exchange roles, put you as author and me as reviewer? ;-)
comment:45 follow-up: ↓ 46 Changed 5 years ago by
Once you've established which parent P to work over, you need a way to find I
in there. Just as getting 0, 1, -1 into P is to construct them in the simplest parent (ZZ in this case) and converting them into P via P(0),P(1),P(-1)
, constructing I
in P should be a similar procedure. There is one solution that works most of the time via cyclotomic fields: P(CyclotomicField(4).0)
In a way, if that fails for P
then P
is probably not a good parent for this kind of computations.
There are some numerical noise problems with CC(CyclotomicField(4),0)
which mightbe worth fixing.
Currently, GF(5)(CyclotomicField(4).0)
doesn't work.
Code along these lines could be a lot more robust for new parents than the special-casing approach taken now. I'm not quite sure we have suitable infrastructure for it to use it now, but perhaps we should (the main issue is of course that the choice of I isn't unique, whereas the integers convert uniquely into any integral domain)
comment:46 in reply to: ↑ 45 ; follow-up: ↓ 48 Changed 5 years ago by
Replying to nbruin:
Once you've established which parent P to work over, you need a way to find
I
in there. Just as getting 0, 1, -1 into P is to construct them in the simplest parent (ZZ in this case) and converting them into P viaP(0),P(1),P(-1)
, constructingI
in P should be a similar procedure. There is one solution that works most of the time via cyclotomic fields:P(CyclotomicField(4).0)
In a way, if that fails forP
thenP
is probably not a good parent for this kind of computations.There are some numerical noise problems with
CC(CyclotomicField(4),0)
which mightbe worth fixing.Currently,
GF(5)(CyclotomicField(4).0)
doesn't work.Code along these lines could be a lot more robust for new parents than the special-casing approach taken now. I'm not quite sure we have suitable infrastructure for it to use it now, but perhaps we should (the main issue is of course that the choice of I isn't unique, whereas the integers convert uniquely into any integral domain)
I am out of my field here, but if I understand you right, are you proposing to replace the elif block proposed by Vincent by something based on the CyclotomicField?? Why is it more robust? aside of the problems you mention of CC(CyclotomicField(4).0)
which may be should be accomplished on a different ticket
comment:47 Changed 5 years ago by
One way to check the sign would be to check the embedding in CC
as CC(custom_I) == CC.gen()
.
And an other way to find I
would be (-P.one()).sqrt()
... but currently there is no specification of this sqrt
operation and sometimes you get a symbolic element where you do not want to.
sage: (-QQ.one()).sqrt().parent() Symbolic Ring sage: (-RR.one()).sqrt().parent() Complex Field with 53 bits of precision sage: (-AA.one()).sqrt().parent() Algebraic Field sage: K = QuadraticField(-1) sage: (-K.one()).sqrt().parent() Number Field in a with defining polynomial x^2 + 1 sage: K = QuadraticField(2) sage: (-K.one()).sqrt().parent() Symbolic Ring
I guess it is reasonable to have a temporary function in sage/rings/
def imaginary_unit(ring): r""" Return the imaginary unit in ``ring`` or in its algebraic closure if `-1` has no square root. If ``ring`` has a defined embedding in some complex field, then the imaginary unit is normalized so that its image in this complex field is its canonical imaginary unit. """ ...
We could make it better later on and possibly remove if we find some other way. At the same time I would advocate that .sqrt()
should return an element in the same ring or, if it is not possible, in its algebraic closure.
comment:48 in reply to: ↑ 46 Changed 5 years ago by
Replying to jhonrubia6:
I am out of my field here, but if I understand you right, are you proposing to replace the elif block proposed by Vincent by something based on the CyclotomicField?? Why is it more robust? aside of the problems you mention of
CC(CyclotomicField(4).0)
which may be should be accomplished on a different ticket
There is no particular reason for P.gen(0) to be a fourth root of unity, even if P contains a fourth root of unity. Therefore, the proposed code can easily fail if it doesn't really have to. Hence, it would be nice to express intent a little more explicitly here. For CyclotomicField(4)
it's clear its generator is a fourth root of unity, so that's a candidate for a "universal" place to find such an element. I'm not so sure it's the best place, but it's at least a place to hook into the general conversion infrastructure.
comment:49 Changed 5 years ago by
The code
if isinstance(a, (int,float,complex)): # Python number a = CDF(a) if isinstance(a, Expression): # symbolic P = SR zero = SR.zero() one = SR.one() I = SR("I") elif isinstance(a, Element): # Sage number P = a.parent() zero = P(0) # Get 0 in Parent ring one = P(1) # Get 1 in Parent ring CF = CyclotomicField(4) # CF.gen() is guaranteed to be I I = P(CF.gen()) # Get I in Parent ring # In case there is some numerical 'noise' # e.g. P=CC ; P(CyclotomicField) == 6.123233995736766e-17 + 1.0*I I = 0.5*(I- I.conjugate()) if I.is_one() or (I*I).is_one() or not (-I*I).is_one(): raise ValueError("invalid number") else: raise ValueError("not a complex number") return matrix(P, 2, [one, zero, zero, -I])
fails when a = QQbar(1+I)
resulting in
Exception raised: Traceback (most recent call last): ... TypeError: Illegal initializer for algebraic number
The offending code being the matrix(P, 2 [one, zero, zero, -I])
sage: P = QQbar(1+I).parent() sage: zero = P(0); one = P(1) sage: matrix(P, 2, [one, zero, zero, -P(CyclotomicField(4).gen())]) Exception raised: Traceback (most recent call last): ... TypeError: Illegal initializer for algebraic number sage: matrix(P,2,[one,zero,zero,-P.gen()]) [ 1 0] [ 0 -I] sage: type(P(CyclotomicField(4).gen()));P(CyclotomicField(4).gen()) <class 'sage.rings.qqbar.AlgebraicNumber'> I sage: type(P.gen());P.gen() <class 'sage.rings.qqbar.AlgebraicNumber'> I
Do we have to open a ticket for that and leave this stopped or we leave the code as it is in the latest commit make it to the release and wait until the problems with CyclotomicField? (numerical noise, and this particular matrix generation) get solved to open a new ticket on this module? IMHO , the commit as it is solves the bug of the ticket (not introducing any new one) while both of your ideas go into the direction of a improvement of the code.
comment:50 Changed 5 years ago by
- Milestone changed from sage-7.2 to sage-7.3
- Status changed from needs_review to positive_review
You are right. We should think about improvements for later tickets! This one already does her part of the job.
comment:51 Changed 5 years ago by
- Branch changed from u/jhonrubia6/hyperbolic_geodesic_midpoint_bugfix to 2ba1c96e25d7e128adc832bc4c4e4c51e8835b98
- Resolution set to fixed
- Status changed from positive_review to closed
The problem seems to be triggered when calculating the inverse of the _to_std_geod(start) matrix