Opened 7 years ago

Closed 7 years ago

# hyperbolic_geodesic midpoint bugfix

Reported by: Owned by: Javier Honrubia González major sage-7.3 geometry hyperbolic geometry, geodesic Javier Honrubia González Vincent Delecroix N/A 2ba1c96 2ba1c96e25d7e128adc832bc4c4e4c51e8835b98

### Description (last modified by Javier Honrubia González)

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__

### comment:1 Changed 7 years ago by Javier Honrubia González

Component: PLEASE CHANGE → geometry hyperbolic geometry geodesic added PLEASE CHANGE → defect

### comment:2 Changed 7 years ago by Javier Honrubia González

Description: modified (diff)

### comment:3 Changed 7 years ago by Javier Honrubia González

Description: modified (diff)

### comment:4 Changed 7 years ago by Javier Honrubia González

The problem seems to be triggered when calculating the inverse of the _to_std_geod(start) matrix

### comment:5 Changed 7 years ago by Travis Scrimshaw

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 7 years ago by Javier Honrubia González

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 double 'full_simplified' expression successes) Obviously works flawlessly with

h=HyperbolicPlane().UHP().get_geodesic(2.0*I,2.0*I+1.0)
h.midpoint()
Last edited 7 years ago by Javier Honrubia González (previous) (diff)

### comment:7 Changed 7 years ago by Javier Honrubia González

Branch: → u/jhonrubia6/hyperbolic_geodesic_midpoint_bugfix

### comment:8 Changed 7 years ago by Javier Honrubia González

Commit: → 89500eb900f245d091afcc4461f4d8f19f6db5b9 new → 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 7 years ago by Javier Honrubia González

Authors: → Javier Honrubia González

### comment:10 follow-up:  12 Changed 7 years ago by Vincent Delecroix

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 7 years ago by Vincent Delecroix

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 7 years ago by Javier Honrubia González

Status: needs_review → needs_work

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?

Last edited 7 years ago by Javier Honrubia González (previous) (diff)

### comment:13 in reply to:  11 Changed 7 years ago by Javier Honrubia González

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 7 years ago by Javier Honrubia González

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 7 years ago by Javier Honrubia González

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 7 years ago by Vincent Delecroix

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 7 years ago by Vincent Delecroix

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?

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 7 years ago by Javier Honrubia González

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

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 7 years ago by Javier Honrubia González

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?

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

Understood. I appreciate your insight, I have not the experience in Sage programming to foresee this subleties. I'll introduce your code in the next patch

Version 1, edited 7 years ago by Javier Honrubia González (previous) (next) (diff)

### comment:20 Changed 7 years ago by git

Commit: 89500eb900f245d091afcc4461f4d8f19f6db5b9 → 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 7 years ago by Vincent Delecroix

Reviewers: → Vincent Delecroix
1. The following syntax is wrong
TESTS::

bla bla.

sage: one_test()

It should be
TESTS:

bla bla.::

sage: one_test()

1. You should get rid of trailing whitespaces, that is line that only contains spaces.
1. The .is_numeric() method only applies to symbolic expressions. There is something weird with your code if S[0,0].is_numeric() since you might also want to manipulate floating point objects.

### comment:22 Changed 7 years ago by git

Commit: afdf46e9237f644889eb993a6236c31e91c6cf4d → 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 7 years ago by Javier Honrubia González

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 7 years ago by Javier Honrubia González

Status: needs_work → needs_review

### comment:25 Changed 7 years ago by Vincent Delecroix

• 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 7 years ago by Vincent Delecroix

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 7 years ago by Javier Honrubia González

Ok. I will try, but would it not be better to get its own trac ticket for that?

### comment:28 Changed 7 years ago by git

Commit: 663a76cb67fb1598301ce842dccd8341ab07218b → 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 7 years ago by Vincent Delecroix

Status: needs_review → 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

or
TESTS::

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 7 years ago by git

Commit: 62159c386ed9ea95663b50606ae770f3e1afdbee → 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 7 years ago by Javier Honrubia González

Status: needs_work → 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.

Last edited 7 years ago by Javier Honrubia González (previous) (diff)

### comment:32 Changed 7 years ago by Vincent Delecroix

Status: needs_review → 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 7 years ago by git

Commit: 929c7c73d9968adb10414dedf879060eb235728f → 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 7 years ago by Javier Honrubia González

Status: needs_work → needs_review

I think I got them all this time. Used pep8 checkers and a rest online checker.

Last edited 7 years ago by Javier Honrubia González (previous) (diff)

### comment:35 Changed 7 years ago by Vincent Delecroix

Why did you do these kind of changes

-#***********************************************************************
+# ***********************************************************************

### comment:36 Changed 7 years ago by Javier Honrubia González

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 7 years ago by Vincent Delecroix

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 7 years ago by git

Commit: dd193839547305ed790c523e905c3e3f63ce563a → 8fea9c2fd8182498e593fbd2e03e151ebefdbea0

Branch pushed to git repo; I updated commit sha1. New commits:

 ​8fea9c2 Removed blank spaces after #

### comment:39 Changed 7 years ago by Javier Honrubia González

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 7 years ago by Vincent Delecroix

It would be nice be open a new ticket for that purpose.

### comment:41 Changed 7 years ago by Javier Honrubia González

Ok. I've opened #20530 for that.

### comment:42 Changed 7 years ago by Vincent Delecroix

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 7 years ago by git

Commit: 8fea9c2fd8182498e593fbd2e03e151ebefdbea0 → 2ba1c96e25d7e128adc832bc4c4e4c51e8835b98

Branch pushed to git repo; I updated commit sha1. New commits:

 ​2ba1c96 Repleaced the code for _get_B()

### comment:44 Changed 7 years ago by Javier Honrubia González

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? ;-)

Last edited 7 years ago by Javier Honrubia González (previous) (diff)

### comment:45 follow-up:  46 Changed 7 years ago by Nils Bruin

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 7 years ago by Javier Honrubia González

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)

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 7 years ago by Vincent Delecroix

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 7 years ago by Nils Bruin

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 7 years ago by Javier Honrubia González

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 7 years ago by Vincent Delecroix

Milestone: sage-7.2 → sage-7.3 needs_review → 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 7 years ago by Volker Braun

Branch: u/jhonrubia6/hyperbolic_geodesic_midpoint_bugfix → 2ba1c96e25d7e128adc832bc4c4e4c51e8835b98 → fixed positive_review → closed
Note: See TracTickets for help on using tickets.