diff --git a/svgpathtools/bezier.py b/svgpathtools/bezier.py
index 1919e68..7c82834 100644
--- a/svgpathtools/bezier.py
+++ b/svgpathtools/bezier.py
@@ -83,12 +83,13 @@ def bezier2polynomial(p, numpy_ordering=True, return_poly1d=False):
         coeffs = p
     else:
         # https://en.wikipedia.org/wiki/Bezier_curve#Polynomial_form
-        n = len(p) + 1
+        n = len(p) - 1
         coeffs = [fac(n)//fac(n-j) * sum(
-            (-1)**(i+j) * p[i] / (fac(i) * fac(j-i)) for i in xrange(j+1))
+            (-1)**(i+j) * p[i] / (fac(i) * fac(j-i)) for i in range(j+1))
             for j in range(n+1)]
-    if not numpy_ordering:
         coeffs.reverse()
+    if not numpy_ordering:
+        coeffs = coeffs[::-1]  # can't use .reverse() as might be tuple
     if return_poly1d:
         return poly1d(coeffs)
     return coeffs
diff --git a/test/test_bezier.py b/test/test_bezier.py
index 87517f2..4c82e93 100644
--- a/test/test_bezier.py
+++ b/test/test_bezier.py
@@ -1,21 +1,58 @@
 from __future__ import division, absolute_import, print_function
-# import unittest
-# try:
-#     from ..bezier import *
-# except ValueError:
-#     from svgpathtools.bezier import *
-#
-#
-# class Test_bezier(unittest.TestCase):
-#     def test_bernstein(self):
-#         self.fail()
-#
-#     def test_bezier_point(self):
-#         self.fail()
-#
-#     def test_split_bezier(self):
-#         self.fail()
-#
-#
-# if __name__ == '__main__':
-#     unittest.main()
+import numpy as np
+import unittest
+from svgpathtools.bezier import *
+from svgpathtools.path import bpoints2bezier
+
+
+class HigherOrderBezier:
+    def __init__(self, bpoints):
+        self.bpts = bpoints
+
+    def bpoints(self):
+        return self.bpts
+
+    def point(self, t):
+        return bezier_point(self.bpoints(), t)
+
+    def __repr__(self):
+        return str(self.bpts)
+
+
+def random_polynomial(degree):
+    return np.poly1d(np.random.rand(degree + 1))
+
+
+def random_bezier(degree):
+    if degree <= 3:
+        return bpoints2bezier(polynomial2bezier(np.random.rand(degree + 1)))
+    else:
+        return HigherOrderBezier(np.random.rand(degree + 1))
+
+
+class TestBezier2Polynomial(unittest.TestCase):
+    def test_bezier2polynomial(self):
+        tvals = np.linspace(0, 1, 10)
+        for d in range(1, 10):
+            b = random_bezier(d)
+            p = np.poly1d(bezier2polynomial(b.bpoints()))
+            for t in tvals:
+                msg = ("degree {}\nt = {}\nb(t) = {}\n = {}\np(t) = \n{}\n = {}"
+                       "".format(d, t, b, b.point(t), p, p(t)))
+                self.assertAlmostEqual(b.point(t), p(t), msg=msg)
+
+
+class TestPolynomial2Bezier(unittest.TestCase):
+    def test_polynomial2bezier(self):
+        tvals = np.linspace(0, 1, 10)
+        for d in range(1, 3):
+            p = random_polynomial(d)
+            b = HigherOrderBezier(polynomial2bezier(p))
+            for t in tvals:
+                msg = ("degree {}\nt = {}\nb(t) = {}\n = {}\np(t) = \n{}\n = {}"
+                       "".format(d, t, b, b.point(t), p, p(t)))
+                self.assertAlmostEqual(b.point(t), p(t), msg=msg)
+
+
+if __name__ == '__main__':
+    unittest.main()