Opened 3 years ago

Closed 3 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) Commit: 2ba1c96e25d7e128adc832bc4c4e4c51e8835b98
Dependencies: Stopgaps:

Description (last modified by jhonrubia6)

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 3 years ago by jhonrubia6

  • Component changed from PLEASE CHANGE to geometry
  • Keywords hyperbolic geometry geodesic added
  • Type changed from PLEASE CHANGE to defect

comment:2 Changed 3 years ago by jhonrubia6

  • Description modified (diff)

comment:3 Changed 3 years ago by jhonrubia6

  • Description modified (diff)

comment:4 Changed 3 years ago by jhonrubia6

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

comment:5 Changed 3 years ago by tscrim

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 3 years ago by jhonrubia6

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 3 years ago by jhonrubia6 (previous) (diff)

comment:7 Changed 3 years ago by jhonrubia6

  • Branch set to u/jhonrubia6/hyperbolic_geodesic_midpoint_bugfix

comment:8 Changed 3 years ago by jhonrubia6

  • 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:

89500ebIf the S matrix is not numeric, apply full_simplify() in order to avoid the runtime error.

comment:9 Changed 3 years ago by jhonrubia6

  • Authors set to Javier Honrubia González

comment:10 follow-up: Changed 3 years ago by 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.

comment:11 follow-up: Changed 3 years ago by 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)?

comment:12 in reply to: ↑ 10 Changed 3 years ago by jhonrubia6

  • 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?

Last edited 3 years ago by jhonrubia6 (previous) (diff)

comment:13 in reply to: ↑ 11 Changed 3 years ago by jhonrubia6

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: Changed 3 years ago by 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

comment:15 follow-up: Changed 3 years ago by jhonrubia6

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: Changed 3 years ago by 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)
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: Changed 3 years ago by 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 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 3 years ago by jhonrubia6

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)
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 3 years ago by jhonrubia6

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 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 these subleties. I'll introduce your code in the next patch

Last edited 3 years ago by jhonrubia6 (previous) (diff)

comment:20 Changed 3 years ago by git

  • Commit changed from 89500eb900f245d091afcc4461f4d8f19f6db5b9 to afdf46e9237f644889eb993a6236c31e91c6cf4d

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

afdf46eModified 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 3 years ago by vdelecroix

  • Reviewers set to 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 3 years ago by git

  • Commit changed from afdf46e9237f644889eb993a6236c31e91c6cf4d to 663a76cb67fb1598301ce842dccd8341ab07218b

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

663a76cRemoved trailing blank spaces. Fixed the TEST syntax. Modified the matrix element check for symbolic matrix check

comment:23 Changed 3 years ago by jhonrubia6

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 3 years ago by jhonrubia6

  • Status changed from needs_work to needs_review

comment:25 Changed 3 years ago by vdelecroix

  • 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 3 years ago by vdelecroix

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 3 years ago by jhonrubia6

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

comment:28 Changed 3 years ago by git

  • Commit changed from 663a76cb67fb1598301ce842dccd8341ab07218b to 62159c386ed9ea95663b50606ae770f3e1afdbee

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

62159c3Code simplified as noted by reviewer. docstring syntax error corrected. New test regarding floating points added in methods.

comment:29 Changed 3 years ago by vdelecroix

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

  • Commit changed from 62159c386ed9ea95663b50606ae770f3e1afdbee to 929c7c73d9968adb10414dedf879060eb235728f

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

929c7c7Revised rst syntax. Tried to adhere to pep8 recommendations. Some long lines reaming though

comment:31 Changed 3 years ago by jhonrubia6

  • 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.

Last edited 3 years ago by jhonrubia6 (previous) (diff)

comment:32 Changed 3 years ago by vdelecroix

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

  • Commit changed from 929c7c73d9968adb10414dedf879060eb235728f to dd193839547305ed790c523e905c3e3f63ce563a

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

dd19383Removed trailing blanks, added double colons before sage code, and revised indentation.

comment:34 Changed 3 years ago by jhonrubia6

  • Status changed from needs_work to needs_review

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

Last edited 3 years ago by jhonrubia6 (previous) (diff)

comment:35 Changed 3 years ago by vdelecroix

Why did you do these kind of changes

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

comment:36 Changed 3 years ago by jhonrubia6

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 3 years ago by vdelecroix

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

  • Commit changed from dd193839547305ed790c523e905c3e3f63ce563a to 8fea9c2fd8182498e593fbd2e03e151ebefdbea0

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

8fea9c2Removed blank spaces after #

comment:39 Changed 3 years ago by jhonrubia6

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 3 years ago by vdelecroix

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

comment:41 Changed 3 years ago by jhonrubia6

Ok. I've opened #20530 for that.

comment:42 Changed 3 years ago by vdelecroix

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

  • Commit changed from 8fea9c2fd8182498e593fbd2e03e151ebefdbea0 to 2ba1c96e25d7e128adc832bc4c4e4c51e8835b98

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

2ba1c96Repleaced the code for _get_B()

comment:44 Changed 3 years ago by jhonrubia6

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 3 years ago by jhonrubia6 (previous) (diff)

comment:45 follow-up: Changed 3 years ago by 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 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: Changed 3 years ago by jhonrubia6

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 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 3 years ago by vdelecroix

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 3 years ago by nbruin

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 3 years ago by jhonrubia6

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 3 years ago by vdelecroix

  • 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 3 years ago by vbraun

  • Branch changed from u/jhonrubia6/hyperbolic_geodesic_midpoint_bugfix to 2ba1c96e25d7e128adc832bc4c4e4c51e8835b98
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.