Why does SymPy calculate wrong intersections of planes?

Question

I have the strange problem, that the intersection of planes in SymPy works with simple examples, but fails for one with more complicated coordinates. I post both a simple example that works, and the one that fails. As the Povray images show, I have three planes that go through vertices of a polyhedron and are perpendicular to the line through the respective vertex and the center. I would like to calculate the point where these planes intersect, but SymPy gives wrong results for the lines in which pairs of planes intersect. In the images the correct intersections can be seen as short lines (created with CSG intersection). The long lines parallel to them are the ones calculated by SymPy.

Am I doing something wrong, or is this a bug in SymPy?

More images are here: http://paste.watchduck.net/1712/sympy_planes/
Does anyone know how to put many images on a page, without being kept from posting the question? ("Your post appears to contain code that is not properly formatted as code.")

Works

Code:

from sympy import Point3D, Plane


pointR = Point3D(1/2, 0, 1/2)
pointG = Point3D(1, 0, 0)

planeR = Plane(pointR, pointR)
planeG = Plane(pointG, pointG)

print('\n######## Intersection of the planes:')
lineRG = planeR.intersection(planeG)[0]  # yellow
print(lineRG)

print('\n######## Intersection of plane and contained line returns the line:')
lineRG_again = planeR.intersection(lineRG)[0]
print(lineRG_again.equals(lineRG))

Output:

######## Intersection of the planes:
Line3D(Point3D(1, 0, 0), Point3D(1, 1/2, 0))

######## Intersection of plane and contained line returns the line:
True

Fails

Code:

from sympy import sqrt, Point3D, Plane

pointR = Point3D(-1, 1 + sqrt(2), -2*sqrt(2) - 1)
pointG = Point3D(-sqrt(2) - 1, 1, -2*sqrt(2) - 1)
pointB = Point3D(-1, -sqrt(2) - 1, -2*sqrt(2) - 1)

planeR = Plane(pointR, pointR)
planeG = Plane(pointG, pointG)
planeB = Plane(pointB, pointB)

print('\n######## Intersections of the planes:')

lineRG = planeR.intersection(planeG)[0]  # yellow
lineRB = planeR.intersection(planeB)[0]  # magenta
lineGB = planeG.intersection(planeB)[0]  # cyan

print(lineRG)
print(lineRB)
print(lineGB)

print('\n######## Lines RG (yellow) and GB (cyan) intersect:')
print(lineRG.intersection(lineGB))
print('\n######## But line RB (magenta) is skew to both of them:')
print(lineRB.intersection(lineRG))
print(lineRB.intersection(lineGB))

print('\n######## Intersection of plane and contained line fails:')
lineRG_again = planeR.intersection(lineRG)

Output:

######## Intersections of the planes:
Line3D(Point3D(-1, 1, 0), Point3D((1 + sqrt(2))*(-2*sqrt(2) - 1) + 2*sqrt(2), -2*sqrt(2) + (-2*sqrt(2) - 1)*(-sqrt(2) - 1), -1 - (1 + sqrt(2))*(-sqrt(2) - 1)))
Line3D(Point3D(-1, 0, 0), Point3D((1 + sqrt(2))*(-2*sqrt(2) - 1) - (-2*sqrt(2) - 1)*(-sqrt(2) - 1) - 1, 0, 2 + 2*sqrt(2)))
Line3D(Point3D(-1, 1, 0), Point3D(-(-2*sqrt(2) - 1)*(-sqrt(2) - 1) - 2*sqrt(2) - 2, -(-2*sqrt(2) - 1)*(-sqrt(2) - 1) + 2 + 2*sqrt(2), 1 + (-sqrt(2) - 1)**2))

######## Lines RG (yellow) and GB (cyan) intersect:
[Point3D(-1, 1, 0)]

######## But line RB (magenta) is skew to both of them:
[]
[]

######## Intersection of plane and contained line fails:
Traceback (most recent call last):
  File "planes2.py", line 47, in <module>
    lineRG_again = planeR.intersection(lineRG)
  File "/home/tilman/Code/024_polyhedron/env/lib/python3.5/site-packages/sympy/geometry/plane.py", line 390, in intersection
    p = a.subs(t, c[0])
  File "/home/tilman/Code/024_polyhedron/env/lib/python3.5/site-packages/sympy/core/basic.py", line 916, in subs
    rv = rv._subs(old, new, **kwargs)
  File "/home/tilman/Code/024_polyhedron/env/lib/python3.5/site-packages/sympy/core/cache.py", line 93, in wrapper
    retval = cfunc(*args, **kwargs)
  File "/home/tilman/Code/024_polyhedron/env/lib/python3.5/site-packages/sympy/core/basic.py", line 1030, in _subs
    rv = fallback(self, old, new)
  File "/home/tilman/Code/024_polyhedron/env/lib/python3.5/site-packages/sympy/core/basic.py", line 1007, in fallback
    rv = self.func(*args)
  File "/home/tilman/Code/024_polyhedron/env/lib/python3.5/site-packages/sympy/geometry/point.py", line 1104, in __new__
    args = Point(*args, **kwargs)
  File "/home/tilman/Code/024_polyhedron/env/lib/python3.5/site-packages/sympy/geometry/point.py", line 159, in __new__
    raise ValueError('Imaginary coordinates are not permitted.')
ValueError: Imaginary coordinates are not permitted.

Images:

Edit: Works with SymPy 1.1.2

After installing the development version of SymPy (pip install git+https://github.com/sympy/sympy.git) I get the correct results:

######## Intersections of pairs of planes:
Line3D(Point3D(-7 + sqrt(2)/2, -sqrt(2)/2 + 7, 0), Point3D((1 + sqrt(2))*(-2*sqrt(2) - 1) - 6 + 5*sqrt(2)/2, -5*sqrt(2)/2 + 6 + (-2*sqrt(2) - 1)*(-sqrt(2) - 1), -1 - (1 + sqrt(2))*(-sqrt(2) - 1)))
Line3D(Point3D(-13 - 6*sqrt(2), 0, 0), Point3D(-13 + (1 + sqrt(2))*(-2*sqrt(2) - 1) - (-2*sqrt(2) - 1)*(-sqrt(2) - 1) - 6*sqrt(2), 0, 2 + 2*sqrt(2)))
Line3D(Point3D(-13/2 - 3*sqrt(2), -7*sqrt(2)/2 + 1/2, 0), Point3D(-(-2*sqrt(2) - 1)*(-sqrt(2) - 1) - 15/2 - 5*sqrt(2), -(-2*sqrt(2) - 1)*(-sqrt(2) - 1) - 3*sqrt(2)/2 + 3/2, 1 + (-sqrt(2) - 1)**2))

######## Intersection of all planes:
[Point3D(0, 0, -20*sqrt(2)/7 - 11/7)]