Ask Your Question
3

Sagemath not evaluating complicated expression

asked 2019-05-29 13:19:14 -0600

cka gravatar image

updated 2019-05-31 15:08:39 -0600

tmonteil gravatar image

I'm working with Sagemath for an in-depth project for the first time. I'm experiencing what I perceive to be some very strange behavior and I'm totally stumped on resolving it. I have a class that is encoding some geometric objects.

It works perfectly fine for relatively simple inputs. The strange thing is that when self.center and/or C get sufficiently complicated, I get an error:

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
<ipython-input-11-4e5b3134f826> in <module>()
      1 c1 = circle(I+Integer(1),Integer(4))
      2 m1 = c1.mobToRealLine()
----> 3 c2 = m1.imageCircle(c1)
      4 print(c2)

<ipython-input-8-5610dc2742b9> in imageCircle(self, circle)
     70         print("new A, B, C = {}, {}, {}".format(newA, newB, newC))
     71 
---> 72         return Circle(newA, newB, newC)
     73 
     74 def circle(center, radius):

<ipython-input-8-5610dc2742b9> in __init__(self, A, B, C)
      9         if A != Integer(0):
     10             self.center = -conjugate(B)
---> 11             self.radius = abs(sqrt(abs(self.center)**Integer(2) - C))
     12         elif A == Integer(0):
     13             self.center = None

/Applications/SageMath/local/lib/python2.7/site-packages/sage/functions/other.pyc in sqrt(x, *args, **kwds)
    869             return sqrt(x)
    870         try:
--> 871             return x.sqrt(*args, **kwds)
    872         # The following includes TypeError to catch cases where sqrt
    873         # is called with a "prec" keyword, for example, but the sqrt

/Applications/SageMath/local/lib/python2.7/site-packages/sage/symbolic/expression.pyx in sage.symbolic.expression.Expression.sqrt (build/cythonized/sage/symbolic/expression.cpp:44675)()
   8075         """
   8076         return new_Expression_from_GEx(self._parent,
-> 8077                 g_hold2_wrapper(g_power_construct, self._gobj, g_ex1_2, hold))
   8078 
   8079     def sin(self, hold=False):

RuntimeError: unsupported type in numeric::lcm

The error occurs when it's trying to initialize a the class with

A, B, C, = 1/2*abs(sqrt(16) - I - 1)^2/sqrt(16) + (1/2*I + 1/2)*(sqrt(16) + I - 1)/sqrt(16) - (1/2*I - 1/2)*(sqrt(16) - I - 1)/sqrt(16) - 7/sqrt(16),1/2*I*(sqrt(16) - I + 1)*(sqrt(16) - I - 1)/sqrt(16) + (1/2*I - 1/2)*(sqrt(16) - I + 1)/sqrt(16) - (1/2*I + 1/2)*(sqrt(16) - I - 1)/sqrt(16) + 7*I/sqrt(16),1/2*abs(sqrt(16) + I + 1)^2/sqrt(16) + (1/2*I - 1/2)*(sqrt(16) + I + 1)/sqrt(16) - (1/2*I + 1/2)*(sqrt(16) - I + 1)/sqrt(16) - 7/sqrt(16)

If I (by hand) initialize the class with these values of A, B, and C, it works just fine. If, on the other hand, I use the "imagecircle" method that I've written (which is implemented as follows. prints are there for debugging), it fails, despite the fact that the A, B, and C given by the imagecircle method are exactly as above. Here is code that if run, gives the error:

class Circle:
    """Unified line and circle class. If A=0,
        then this represents a line."""

    def __init__(self, A, B, C):
        self.A = A
        self.B = B
        self.C = C
        if A != 0:
            self.center = -conjugate(B)
            self.radius = abs(sqrt(abs(self.center)^2 - C))
        elif A == 0:
            self.center = None
            self.radius = None
    def isCircle(self):
        """Am I a geometric circle in the complex plane?"""
        return self.A != 0

    def isLine(self):
        """Am I a geometric line in the complex plane?"""
        return self.A == 0

    def mobToRealLine(self):
        if self.isCircle():
            """return the Mobius that takes self to the real line. It takes z0+r to 0, and z0 +- Ir to +-1"""
            z0 = self.center
            r = self.radius
            a = 1/sqrt(r) * (1/2 - I/2)
            b = -1/sqrt(r) * (1/2 - I/2) * (r+z0)
            c = 1/sqrt(r) * (1/2 + I/2)
            d = 1/sqrt(r) * (1/2 + I/2) * (r-z0)
            return Mobius(a,b,c,d)
        if self.isLine() and imag(self.B) != 0: #This takes iC/2Im(B), 1/2 +- I(ReB + ImC)/(2ImB) to 0, +-1
            B, C = self.B, self.C
            a = (1-I) / sqrt( conjugate(B)/imag(B) )
            b = -( (1/2+I/2)*C*sqrt( conjugate(B)/imag(B) ) ) / ( conjugate(B) )
            c = 0
            d = (1/2+I/2)*sqrt( conjugate(B)/imag(B) )
            return Mobius(a,b,c,d)
        else:#Vertical line case x=x0
            B, C = self.B, self.C
            x0 = -C/(2*B)
            a = -sqrt(2)/2 + I*sqrt(2)/2
            b = (sqrt(2)/2 - I*sqrt(2)/2)*x0
            c = 0
            d = -sqrt(2)/2 - I*sqrt(2)/2
            return Mobius(a,b,c,d)



class Mobius:
    """Mobius transformation class"""
    def __init__(self, a, b, c, d):
        self.a   = a
        self.b   = b
        self.c   = c
        self.d   = d
        self.matrix = Matrix([[a,b],[c,d]])

    def imageCircle(self, circle):#This seems broken for a vertical line...
        """Return the image of circle under self."""
        A, B, C = circle.A, circle.B, circle.C
        cA, cB, cC = conjugate(A), conjugate(B), conjugate(C)
        a, b, c, d = self.a, self.b, self.c, self.d
        ca, cb, cc, cd = conjugate(a), conjugate(b), conjugate(c), conjugate(d)

        newA = A*abs(d)^2 - cc*d*B - c*cd*cB + C*abs(c)^2
        newB = -A*cb*d + ca*d*B + cb*c*cB - ca*c*C
        newC = A*abs(b)^2 - ca*b*B - a*cb*cB + abs(a)^2*C
        print("new A, B, C = {}, {}, {}".format(newA, newB, newC))

        return Circle(newA, newB, newC)

def circle(center, radius):
    """returns Circle(A,B,C) with appropriate A, B, C (for simplicity of use)"""
    A = 1
    B = -conjugate(center)
    C = abs(center)^2 - radius^2
    return Circle(A,B,C)

c1 = circle(I+1,4)
m1 = c1.mobToRealLine()
c2 = m1.imageCircle(c1)
print(c2)
edit retag flag offensive close merge delete

Comments

1

Please provide the exact code one can use to reproduce the error.

Try using algebraic numbers instead of the symbolic ring.

slelievre gravatar imageslelievre ( 2019-05-30 01:16:20 -0600 )edit

Ok, I updated my post with code to reproduce the error. If using algebraic numbers is the key, can you lead me to how to get that working here? Thanks in advance for your help.

cka gravatar imagecka ( 2019-05-30 08:51:16 -0600 )edit
1

There seems to be some weird "holding" going on, e.g. sqrt(abs(1+I)^2 + 14) evaluates to sqrt(16) instead of 4.

rburing gravatar imagerburing ( 2019-05-30 09:56:08 -0600 )edit

Yes, exactly! I don't understand why this is happening.

cka gravatar imagecka ( 2019-05-30 09:58:17 -0600 )edit
2

I reported this (sub)issue as trac ticket #27897.

rburing gravatar imagerburing ( 2019-05-30 10:19:34 -0600 )edit

2 answers

Sort by ยป oldest newest most voted
2

answered 2019-05-31 15:08:16 -0600

tmonteil gravatar image

This seems indeed to be a bug, thanks for reporting, and thanks @rburing for the bug report, it is now trac ticket 27897.

edit flag offensive delete link more
0

answered 2019-06-02 16:04:57 -0600

slelievre gravatar image

The code in the question seems to be designed to work with exact algebraic numbers.

In that case, it would seem natural to work with Sage's algebraic numbers, using QQbar (the field of algebraic numbers) and AA (the field of real algebraic numbers).

Below is an adaptation of the code in the question to do that (additionally, this works around the problem of algebraic numbers too deeply nested in Sage's symbolic ring).

The exactify method on elements of QQbar and AA obtained from computations can simplify their internal representation (and thus make them display nicer).

class Circle:
    """Unified line and circle class. If A is zero,
        then this represents a line."""

    def __init__(self, A, B, C):
        self.A = AA(A)
        self.B = QQbar(B)
        self.C = AA(C)
        if A != 0:
            self.center = -conjugate(B)
            self.radius = abs(sqrt(abs(self.center)^2 - C))
            self.center.exactify()
            self.radius.exactify()
        elif A == 0:
            self.center = None
            self.radius = None

    def isCircle(self):
        """Am I a geometric circle in the complex plane?"""
        return self.A != 0

    def isLine(self):
        """Am I a geometric line in the complex plane?"""
        return self.A == 0

    def mobToRealLine(self):
        I = QQbar.gen()
        if self.isCircle():
            """Return the Mobius that takes self to the real line. It takes z0 + r to 0, and z0 +- Ir to +-1"""
            z0 = self.center
            r = self.radius
            a = 1/sqrt(r) * (1/2 - I/2)
            b = -1/sqrt(r) * (1/2 - I/2) * (r+z0)
            c = 1/sqrt(r) * (1/2 + I/2)
            d = 1/sqrt(r) * (1/2 + I/2) * (r-z0)
            a.exactify()
            b.exactify()
            c.exactify()
            d.exactify()
            return Mobius(a, b, c, d)
        if self.isLine() and imag(self.B) != 0:
            # This takes iC/2Im(B), 1/2 +- I(ReB + ImC)/(2ImB) to 0, +-1
            B, C = self.B, self.C
            a = (1-I)/sqrt(conjugate(B)/imag(B))
            b = -(1/2+I/2)*C*sqrt(conjugate(B)/imag(B))/conjugate(B)
            c = QQbar.zero()
            d = (1/2+I/2)*sqrt(conjugate(B)/imag(B))
            a.exactify()
            b.exactify()
            c.exactify()
            d.exactify()
            return Mobius(a, b, c, d)
        else:
            # Vertical line x = x0
            B, C = self.B, self.C
            x0 = -C/(2*B)
            a = -sqrt(2)/2 + I*sqrt(2)/2
            b = (sqrt(2)/2 - I*sqrt(2)/2)*x0
            c = QQbar.zero()
            d = -sqrt(2)/2 - I*sqrt(2)/2
            a.exactify()
            b.exactify()
            c.exactify()
            d.exactify()
            return Mobius(a, b, c, d)

class Mobius:
    """Mobius transformation class"""
    def __init__(self, a, b, c, d):
        self.a = QQbar(a)
        self.b = QQbar(b)
        self.c = QQbar(c)
        self.d = QQbar(d)
        self.a.exactify()
        self.b.exactify()
        self.c.exactify()
        self.d.exactify()
        self.matrix = Matrix(QQbar, [[self.a, self.b], [self.c, self.d]])

    def imageCircle(self, circle):  # Seems broken for a vertical line
        """Return the image of circle under self."""
        A, B, C = circle.A, circle.B, circle.C
        cA, cB, cC = conjugate(A), conjugate(B), conjugate(C)
        a, b, c, d = self.a, self.b, self.c, self.d
        ca, cb, cc, cd = conjugate(a), conjugate(b), conjugate(c), conjugate(d)
        newA = A*abs(d)^2 - cc*d*B - c*cd*cB + C*abs(c)^2
        newB = -A*cb*d + ca*d*B + cb*c*cB - ca*c*C
        newC = A*abs(b)^2 - ca*b*B - a*cb*cB + abs(a)^2*C
        newA.exactify()
        newB.exactify()
        newC.exactify()
        print("new A, B, C = {}, {}, {}".format(newA, newB, newC))
        return Circle(newA, newB, newC)

def circle(center, radius):
    """returns Circle(A,B,C) with appropriate A, B, C (for simplicity of use)"""
    A = AA.one()
    B = -conjugate(QQbar(center))
    C = abs(QQbar(center))^2 - AA(radius)^2
    B.exactify()
    C.exactify()
    return Circle(A, B, C)

With the above definitions, we get:

sage: c1 = circle(I + 1, 4)
sage: m1 = c1.mobToRealLine()
sage: c2 = m1.imageCircle(c1)
new A, B, C = 0, 4*I, 0
edit flag offensive delete link more

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower

Stats

Asked: 2019-05-29 10:06:30 -0600

Seen: 124 times

Last updated: Jun 02