title: Why does SymPy calculate wrong intersections of planes?
body: <p>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 href="https://en.wikipedia.org/wiki/Truncated_cuboctahedron" rel="noreferrer">a polyhedron</a> 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 <a href="https://en.wikipedia.org/wiki/Constructive_solid_geometry" rel="noreferrer">CSG</a> intersection). The long lines parallel to them are the ones calculated by SymPy.</p>
  
  <p>Am I doing something wrong, or is this a bug in SymPy?</p>
  
  <p>More images are here: <a href="http://paste.watchduck.net/1712/sympy_planes/" rel="noreferrer">http://paste.watchduck.net/1712/sympy_planes/</a> <br>
  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.")</p>
  
  <h1>Works</h1>
  
  <p>Code:</p>
  
  <pre><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))
  </code></pre>
  
  <p>Output:</p>
  
  <pre><code>######## Intersection of the planes:
  Line3D(Point3D(1, 0, 0), Point3D(1, 1/2, 0))
  
  ######## Intersection of plane and contained line returns the line:
  True
  </code></pre>
  
  <h1>Fails</h1>
  
  <p>Code:</p>
  
  <pre><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)
  </code></pre>
  
  <p>Output:</p>
  
  <pre><code>######## 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 &lt;module&gt;
      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.
  </code></pre>
  
  <p>Images:</p>
  
  <p><a href="https://i.stack.imgur.com/jbwCE.png" rel="noreferrer"><img src="https://i.stack.imgur.com/jbwCE.png" alt="planes"></a>
  <a href="https://i.stack.imgur.com/2ozkT.png" rel="noreferrer"><img src="https://i.stack.imgur.com/2ozkT.png" alt="wrong intersections"></a></p>
  
  <h1>Edit: Works with SymPy 1.1.2</h1>
  
  <p>After installing the development version of SymPy (<code>pip install git+https://github.com/sympy/sympy.git</code>) I get the correct results:</p>
  
  <pre><code>######## 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)]
  </code></pre>
  
  <p><a href="https://i.stack.imgur.com/RUHtS.png" rel="noreferrer"><img src="https://i.stack.imgur.com/RUHtS.png" alt="enter image description here"></a></p>