Repository: sympy/sympy

Issue:
Should the output of arbitrary_point be changed for Plane objects?
According to the docs, [arbitrary_point](http://docs.sympy.org/latest/_modules/sympy/geometry/plane.html#Plane.arbitrary_point) is supposed to parameterize a circle of radius 1 in the plane.  Maybe this what its output does, but it's hard to see how anyone could use it (scroll to the right to see the entire formula). 
```
>>> from sympy import *
>>> Plane(Point(1, 2, 3), Point(2, 3, 4), Point(0, 0, 0)).arbitrary_point()
Point3D(-(-2*sin(_t) + cos(_t) - 3)*(-2*sqrt(-84*sin(2*_t) - 216*sin(_t) - 171*cos(2*_t) - 1044*cos(_t) - 921) + 24*sin(_t) + 12*cos(_t) + 60)/(6*(4*sin(2*_t) + 24*sin(_t) + 3*cos(2*_t) + 36*cos(_t) + 49)) + (-2*sqrt(-84*sin(2*_t) - 216*sin(_t) - 171*cos(2*_t) - 1044*cos(_t) - 921) + 24*sin(_t) + 12*cos(_t) + 60)*cos(_t)/(4*sin(2*_t) + 24*sin(_t) + 3*cos(2*_t) + 36*cos(_t) + 49) + (-2*sqrt(-84*sin(2*_t) - 216*sin(_t) - 171*cos(2*_t) - 1044*cos(_t) - 921) + 24*sin(_t) + 12*cos(_t) + 60)/(4*sin(2*_t) + 24*sin(_t) + 3*cos(2*_t) + 36*cos(_t) + 49), (sin(_t) + 2)*(-2*sqrt(-84*sin(2*_t) - 216*sin(_t) - 171*cos(2*_t) - 1044*cos(_t) - 921) + 24*sin(_t) + 12*cos(_t) + 60)/(4*sin(2*_t) + 24*sin(_t) + 3*cos(2*_t) + 36*cos(_t) + 49) + (-2*sin(_t) + cos(_t) - 3)*(-2*sqrt(-84*sin(2*_t) - 216*sin(_t) - 171*cos(2*_t) - 1044*cos(_t) - 921) + 24*sin(_t) + 12*cos(_t) + 60)/(3*(4*sin(2*_t) + 24*sin(_t) + 3*cos(2*_t) + 36*cos(_t) + 49)), -(-2*sin(_t) + cos(_t) - 3)*(-2*sqrt(-84*sin(2*_t) - 216*sin(_t) - 171*cos(2*_t) - 1044*cos(_t) - 921) + 24*sin(_t) + 12*cos(_t) + 60)/(6*(4*sin(2*_t) + 24*sin(_t) + 3*cos(2*_t) + 36*cos(_t) + 49)))
```
Normally, parametric equation of a circle involves sin(t) and cos(t) in a linear way, not under square roots and mixed with sin(2t) and cos(2t). Simplifying didn't get make much difference.

My suggestion: instead of trying to fix the method, change it entirely so that it returns a point with two parameters s_, t_ where all three coordinates are linear functions of s_ and t_. This means changing a public method, but from what I can see, this method could hardly be used in practice. 

Also, #13789 proposes to implement a parameter finding routine, given a point on the entity. This was attempted in #13790 but since plane's arbitrary_point does not run through the entire plane, the parameter finding can't possibly work correctly. (It also won't work because `solve` will hang on such complex equations.)  

Discussion:

Patch:
diff --git a/sympy/geometry/plane.py b/sympy/geometry/plane.py
index 0f92f08efe..e3c2be83c9 100644
--- a/sympy/geometry/plane.py
+++ b/sympy/geometry/plane.py
@@ -10,7 +10,7 @@
 from sympy import simplify
 from sympy.core import Dummy, Rational, S, Symbol
 from sympy.core.compatibility import is_sequence
-from sympy.functions.elementary.trigonometric import acos, asin, sqrt
+from sympy.functions.elementary.trigonometric import cos, sin, acos, asin, sqrt
 from sympy.matrices import Matrix
 from sympy.polys.polytools import cancel
 from sympy.solvers import solve, linsolve
@@ -67,6 +67,8 @@ def __new__(cls, p1, a=None, b=None, **kwargs):
                 raise ValueError(filldedent('''
                     Either provide 3 3D points or a point with a
                     normal vector expressed as a sequence of length 3'''))
+            if all(coord.is_zero for coord in normal_vector):
+                raise ValueError('Normal vector cannot be zero vector')
         return GeometryEntity.__new__(cls, p1, normal_vector, **kwargs)
 
     def __contains__(self, o):
@@ -133,20 +135,22 @@ def angle_between(self, o):
             return acos(c/(d*e))
 
 
-    def arbitrary_point(self, t=None):
-        """ Returns an arbitrary point on the Plane; varying `t` from 0 to 2*pi
-        will move the point in a circle of radius 1 about p1 of the Plane.
+    def arbitrary_point(self, u=None, v=None):
+        """ Returns an arbitrary point on the Plane. If given two parameters,
+        the point ranges over the entire plane. If given one or no parameters,
+        returns a point with one parameter, varying which from 0 to 2*pi will
+        move the point in a circle of radius 1 about p1 of the Plane.
 
         Examples
         ========
 
         >>> from sympy.geometry.plane import Plane
-        >>> from sympy.abc import t
-        >>> p = Plane((0, 0, 0), (0, 0, 1), (0, 1, 0))
+        >>> from sympy.abc import u, v, t
+        >>> p = Plane((3, 1, 4), (2, 0, 1), (-1, -2, -2))
+        >>> p.arbitrary_point(u, v)
+        Point3D(-6*u - 3*v + 3, -3*u + 6*v + 1, 45*v + 4)
         >>> p.arbitrary_point(t)
-        Point3D(0, cos(t), sin(t))
-        >>> _.distance(p.p1).simplify()
-        1
+        Point3D(-sqrt(230)*sin(t)/230 - 2*sqrt(5)*cos(t)/5 + 3, sqrt(230)*sin(t)/115 - sqrt(5)*cos(t)/5 + 1, 3*sqrt(230)*sin(t)/46 + 4)
 
         Returns
         =======
@@ -154,21 +158,28 @@ def arbitrary_point(self, t=None):
         Point3D
 
         """
-        from sympy import cos, sin
-        t = t or Dummy('t')
+        circle = v is None
+        u = u or Dummy('u', real=True)
+        v = v or Dummy('v', real=True)
         x, y, z = self.normal_vector
         a, b, c = self.p1.args
-        if x == y == 0:
-            return Point3D(a + cos(t), b + sin(t), c)
-        elif x == z == 0:
-            return Point3D(a + cos(t), b, c + sin(t))
-        elif y == z == 0:
-            return Point3D(a, b + cos(t), c + sin(t))
-        m = Dummy()
-        p = self.projection(Point3D(self.p1.x + cos(t), self.p1.y + sin(t), 0)*m)
-
-        # TODO: Replace solve with solveset, when this line is tested
-        return p.xreplace({m: solve(p.distance(self.p1) - 1, m)[0]})
+        # x1, y1, z1 is a nonzero vector parallel to the plane
+        if x.is_zero and y.is_zero:
+            x1, y1, z1 = S.One, S.Zero, S.Zero
+        else:
+            x1, y1, z1 = -y, x, S.Zero
+        # x2, y2, z2 is also parallel to the plane, and orthogonal to x1, y1, z1
+        x2, y2, z2 = tuple(Matrix((x, y, z)).cross(Matrix((x1, y1, z1))))
+        if circle:
+            x1, y1, z1 = (w/sqrt(x1**2 + y1**2 + z1**2) for w in (x1, y1, z1))
+            x2, y2, z2 = (w/sqrt(x2**2 + y2**2 + z2**2) for w in (x2, y2, z2))
+            p = Point3D(a + x1*cos(u) + x2*sin(u), \
+                        b + y1*cos(u) + y2*sin(u), \
+                        c + z1*cos(u) + z2*sin(u))
+        else:
+            p = Point3D(a + x1*u + x2*v, b + y1*u + y2*v, c + z1*u + z2*v)
+        return p
+
 
     @staticmethod
     def are_concurrent(*planes):
@@ -219,6 +230,7 @@ def are_concurrent(*planes):
                     return False
             return True
 
+
     def distance(self, o):
         """Distance between the plane and another geometric entity.
 
@@ -280,6 +292,7 @@ def distance(self, o):
             f = sqrt(sum([i**2 for i in self.normal_vector]))
             return abs(e / f)
 
+
     def equals(self, o):
         """
         Returns True if self and o are the same mathematical entities.
@@ -305,6 +318,7 @@ def equals(self, o):
         else:
             return False
 
+
     def equation(self, x=None, y=None, z=None):
         """The equation of the Plane.
 
@@ -326,6 +340,7 @@ def equation(self, x=None, y=None, z=None):
         c = self.normal_vector
         return (sum(i*j for i, j in zip(b, c)))
 
+
     def intersection(self, o):
         """ The intersection with other geometrical entity.
 
@@ -406,6 +421,7 @@ def intersection(self, o):
                 for i in (x, y, z): result = result.subs(i, 0)
                 return [Line3D(Point3D(result), direction_ratio=c)]
 
+
     def is_coplanar(self, o):
         """ Returns True if `o` is coplanar with self, else False.
 
@@ -431,6 +447,7 @@ def is_coplanar(self, o):
         elif isinstance(o, GeometryEntity):  # XXX should only be handling 2D objects now
             return all(i == 0 for i in self.normal_vector[:2])
 
+
     def is_parallel(self, l):
         """Is the given geometric entity parallel to the plane?
 
@@ -471,6 +488,7 @@ def is_parallel(self, l):
             else:
                 return False
 
+
     def is_perpendicular(self, l):
         """is the given geometric entity perpendicualar to the given plane?
 
@@ -774,5 +792,6 @@ def random_point(self, seed=None):
             rng = random.Random(seed)
         else:
             rng = random
-        t = Dummy('t')
-        return self.arbitrary_point(t).subs(t, Rational(rng.random()))
+        u, v = Dummy('u'), Dummy('v')
+        params = {u: Rational(rng.random()), v: Rational(rng.random())}
+        return self.arbitrary_point(u, v).subs(params)


Test patch:
diff --git a/sympy/geometry/tests/test_plane.py b/sympy/geometry/tests/test_plane.py
index f57aac6180..4ce39aacc4 100644
--- a/sympy/geometry/tests/test_plane.py
+++ b/sympy/geometry/tests/test_plane.py
@@ -1,6 +1,6 @@
 from __future__ import division
 
-from sympy import Dummy, S, Symbol, pi, sqrt, asin
+from sympy import Dummy, S, Symbol, symbols, pi, sqrt, asin, sin, cos
 from sympy.geometry import Line, Point, Ray, Segment, Point3D, Line3D, Ray3D, Segment3D, Plane
 from sympy.geometry.util import are_coplanar
 from sympy.utilities.pytest import raises, slow
@@ -8,9 +8,7 @@
 
 @slow
 def test_plane():
-    x = Symbol('x', real=True)
-    y = Symbol('y', real=True)
-    z = Symbol('z', real=True)
+    x, y, z, u, v = symbols('x y z u v', real=True)
     p1 = Point3D(0, 0, 0)
     p2 = Point3D(1, 1, 1)
     p3 = Point3D(1, 2, 3)
@@ -80,6 +78,16 @@ def test_plane():
     assert pl6.is_perpendicular(pl7)
     assert pl6.is_perpendicular(l1) is False
 
+    assert pl6.distance(pl6.arbitrary_point(u, v)) == 0
+    assert pl7.distance(pl7.arbitrary_point(u, v)) == 0
+    assert pl6.distance(pl6.arbitrary_point(t)) == 0
+    assert pl7.distance(pl7.arbitrary_point(t)) == 0
+    assert pl6.p1.distance(pl6.arbitrary_point(t)).simplify() == 1
+    assert pl7.p1.distance(pl7.arbitrary_point(t)).simplify() == 1
+    assert pl3.arbitrary_point(t) == Point3D(-sqrt(30)*sin(t)/30 + \
+        2*sqrt(5)*cos(t)/5, sqrt(30)*sin(t)/15 + sqrt(5)*cos(t)/5, sqrt(30)*sin(t)/6)
+    assert pl3.arbitrary_point(u, v) == Point3D(2*u - v, u + 2*v, 5*v)
+
     assert pl7.distance(Point3D(1, 3, 5)) == 5*sqrt(6)/6
     assert pl6.distance(Point3D(0, 0, 0)) == 4*sqrt(3)
     assert pl6.distance(pl6.p1) == 0
@@ -106,6 +114,7 @@ def test_plane():
     assert Plane.are_concurrent(pl3, pl4, pl5) is False
     assert Plane.are_concurrent(pl6) is False
     raises(ValueError, lambda: Plane.are_concurrent(Point3D(0, 0, 0)))
+    raises(ValueError, lambda: Plane((1, 2, 3), normal_vector=(0, 0, 0)))
 
     assert pl3.parallel_plane(Point3D(1, 2, 5)) == Plane(Point3D(1, 2, 5), \
                                                       normal_vector=(1, -2, 1))