diff --git a/demos/demo1.py b/demos/demo1.py
index 1f7a3db..171f42a 100644
--- a/demos/demo1.py
+++ b/demos/demo1.py
@@ -26,24 +26,24 @@
 degree = 4
 
 # The first way of doing it is by directly supplying the weight function.
-wf = lambda(x): 1. / math.sqrt(2. * math.pi) * np.exp(-x ** 2 / 2.)
+wf = lambda x: 1. / math.sqrt(2. * math.pi) * np.exp(-x ** 2 / 2.)
 # Construct it:
 p = orthpol.OrthogonalPolynomial(degree,
                                 left=-np.inf, right=np.inf, # Domain
                                 wf=wf)
 # An orthogonal polynomial is though of as a function.
 # Here is how to get the number of inputs and outputs of that function
-print 'Number of inputs:', p.num_input
-print 'Number of outputs:', p.num_output
+print(('Number of inputs:', p.num_input))
+print(('Number of outputs:', p.num_output))
 # Test if the polynomials are normalized (i.e., their norm is 1.):
-print 'Is normalized:', p.is_normalized
+print(('Is normalized:', p.is_normalized))
 # Get the degree of the polynomial:
-print 'Polynomial degree:', p.degree
+print(('Polynomial degree:', p.degree))
 # Get the alpha-beta recursion coefficients:
-print 'Alpha:', p.alpha
-print 'Beta:', p.beta
+print(('Alpha:', p.alpha))
+print(('Beta:', p.beta))
 # The following should print a description of the polynomial
-print str(p)
+print((str(p)))
 # Now you can evaluate the polynomial at any points you want:
 X = np.linspace(-2., 2., 100)
 # Here is the actual evaluation
@@ -55,7 +55,7 @@
 plt.xlabel('$x$', fontsize=16)
 plt.ylabel('$p_i(x)$', fontsize=16)
 plt.legend(['$p_{%d}(x)$' % i for i in range(p.num_output)], loc='best')
-print 'Close the window to continue...'
+print('Close the window to continue...')
 plt.show()
 # You may also compute the derivatives of the polynomials:
 dphi = p.d(X)
@@ -65,5 +65,5 @@
 plt.xlabel('$x$', fontsize=16)
 plt.ylabel(r'$\frac{dp_i(x)}{dx}$', fontsize=16)
 plt.legend([r'$\frac{p_{%d}(x)}{dx}$' % i for i in range(p.num_output)], loc='best')
-print 'Close the window to end demo...'
+print('Close the window to end demo...')
 plt.show()
diff --git a/demos/demo10.py b/demos/demo10.py
index a316445..9fc3113 100644
--- a/demos/demo10.py
+++ b/demos/demo10.py
@@ -37,17 +37,17 @@
 p = orthpol.OrthogonalPolynomial(degree, rv=rv)
 # An orthogonal polynomial is though of as a function.
 # Here is how to get the number of inputs and outputs of that function
-print 'Number of inputs:', p.num_input
-print 'Number of outputs:', p.num_output
+print('Number of inputs:', p.num_input)
+print('Number of outputs:', p.num_output)
 # Test if the polynomials are normalized (i.e., their norm is 1.):
-print 'Is normalized:', p.is_normalized
+print('Is normalized:', p.is_normalized)
 # Get the degree of the polynomial:
-print 'Polynomial degree:', p.degree
+print('Polynomial degree:', p.degree)
 # Get the alpha-beta recursion coefficients:
-print 'Alpha:', p.alpha
-print 'Beta:', p.beta
+print('Alpha:', p.alpha)
+print('Beta:', p.beta)
 # The following should print a description of the polynomial
-print str(p)
+print(str(p))
 # Now you can evaluate the polynomial at any points you want:
 X = np.linspace(lower, upper, 100)
 # Here is the actual evaluation
@@ -59,7 +59,7 @@
 plt.xlabel('$x$', fontsize=16)
 plt.ylabel('$p_i(x)$', fontsize=16)
 plt.legend(['$p_{%d}(x)$' % i for i in range(p.num_output)], loc='best')
-print 'Close the window to continue...'
+print('Close the window to continue...')
 plt.show()
 # You may also compute the derivatives of the polynomials:
 dphi = p.d(X)
@@ -69,5 +69,5 @@
 plt.xlabel('$x$', fontsize=16)
 plt.ylabel(r'$\frac{dp_i(x)}{dx}$', fontsize=16)
 plt.legend([r'$\frac{p_{%d}(x)}{dx}$' % i for i in range(p.num_output)], loc='best')
-print 'Close the window to end demo...'
+print('Close the window to end demo...')
 plt.show()
diff --git a/demos/demo11.py b/demos/demo11.py
index 7d0c7c1..804e5cf 100644
--- a/demos/demo11.py
+++ b/demos/demo11.py
@@ -25,15 +25,15 @@
 # Then construct the product basis
 p = orthpol.ProductBasis(degree=degree, rvs=rvs)
 # Print info about the polynomials
-print str(p)
+print(str(p))
 # Evaluate the polynomials at some points
 X = np.hstack([rvs[0].rvs(size=(100, 1)), rvs[1].rvs(size=(100, 1))])
 # Look at the shape of X, it should be 100x2:
-print 'X shape:', X.shape
+print('X shape:', X.shape)
 # Evaluate the polynomials at X
 phi = p(X)
 # Look at the shape of phi, it should be 100xp.num_output
-print 'phi shape:', phi.shape
+print('phi shape:', phi.shape)
 # Take a look at the phi's also
-print 'phi:'
-print phi
+print('phi:')
+print(phi)
diff --git a/demos/demo2.py b/demos/demo2.py
index 6cae5ab..f7ec731 100644
--- a/demos/demo2.py
+++ b/demos/demo2.py
@@ -26,24 +26,24 @@
 degree = 4
 
 # The first way of doing it is by directly supplying the weight function.
-wf = lambda(x): np.exp(-x)
+wf = lambda x: np.exp(-x)
 # Construct it:
 p = orthpol.OrthogonalPolynomial(degree,
                                 left=0, right=np.inf, # Domain
                                 wf=wf)
 # An orthogonal polynomial is though of as a function.
 # Here is how to get the number of inputs and outputs of that function
-print 'Number of inputs:', p.num_input
-print 'Number of outputs:', p.num_output
+print('Number of inputs:', p.num_input)
+print('Number of outputs:', p.num_output)
 # Test if the polynomials are normalized (i.e., their norm is 1.):
-print 'Is normalized:', p.is_normalized
+print('Is normalized:', p.is_normalized)
 # Get the degree of the polynomial:
-print 'Polynomial degree:', p.degree
+print('Polynomial degree:', p.degree)
 # Get the alpha-beta recursion coefficients:
-print 'Alpha:', p.alpha
-print 'Beta:', p.beta
+print('Alpha:', p.alpha)
+print('Beta:', p.beta)
 # The following should print a description of the polynomial
-print str(p)
+print(str(p))
 # Now you can evaluate the polynomial at any points you want:
 X = np.linspace(0., 2., 100)
 # Here is the actual evaluation
@@ -55,7 +55,7 @@
 plt.xlabel('$x$', fontsize=16)
 plt.ylabel('$p_i(x)$', fontsize=16)
 plt.legend(['$p_{%d}(x)$' % i for i in range(p.num_output)], loc='best')
-print 'Close the window to continue...'
+print('Close the window to continue...')
 plt.show()
 # You may also compute the derivatives of the polynomials:
 dphi = p.d(X)
@@ -65,5 +65,5 @@
 plt.xlabel('$x$', fontsize=16)
 plt.ylabel(r'$\frac{dp_i(x)}{dx}$', fontsize=16)
 plt.legend([r'$\frac{p_{%d}(x)}{dx}$' % i for i in range(p.num_output)], loc='best')
-print 'Close the window to end demo...'
+print('Close the window to end demo...')
 plt.show()
diff --git a/demos/demo3.py b/demos/demo3.py
index 46209a9..a59d3a8 100644
--- a/demos/demo3.py
+++ b/demos/demo3.py
@@ -26,24 +26,24 @@
 degree = 4
 
 # The first way of doing it is by directly supplying the weight function.
-wf = lambda(x): 1. / np.sqrt(1. - x)
+wf = lambda x: 1. / np.sqrt(1. - x)
 # Construct it:
 p = orthpol.OrthogonalPolynomial(degree,
                                 left=-1., right=1., # Domain
                                 wf=wf)
 # An orthogonal polynomial is though of as a function.
 # Here is how to get the number of inputs and outputs of that function
-print 'Number of inputs:', p.num_input
-print 'Number of outputs:', p.num_output
+print('Number of inputs:', p.num_input)
+print('Number of outputs:', p.num_output)
 # Test if the polynomials are normalized (i.e., their norm is 1.):
-print 'Is normalized:', p.is_normalized
+print('Is normalized:', p.is_normalized)
 # Get the degree of the polynomial:
-print 'Polynomial degree:', p.degree
+print('Polynomial degree:', p.degree)
 # Get the alpha-beta recursion coefficients:
-print 'Alpha:', p.alpha
-print 'Beta:', p.beta
+print('Alpha:', p.alpha)
+print('Beta:', p.beta)
 # The following should print a description of the polynomial
-print str(p)
+print(str(p))
 # Now you can evaluate the polynomial at any points you want:
 X = np.linspace(-1., 1., 100)
 # Here is the actual evaluation
@@ -55,7 +55,7 @@
 plt.xlabel('$x$', fontsize=16)
 plt.ylabel('$p_i(x)$', fontsize=16)
 plt.legend(['$p_{%d}(x)$' % i for i in range(p.num_output)], loc='best')
-print 'Close the window to continue...'
+print('Close the window to continue...')
 plt.show()
 # You may also compute the derivatives of the polynomials:
 dphi = p.d(X)
@@ -65,5 +65,5 @@
 plt.xlabel('$x$', fontsize=16)
 plt.ylabel(r'$\frac{dp_i(x)}{dx}$', fontsize=16)
 plt.legend([r'$\frac{p_{%d}(x)}{dx}$' % i for i in range(p.num_output)], loc='best')
-print 'Close the window to end demo...'
+print('Close the window to end demo...')
 plt.show()
diff --git a/demos/demo4.py b/demos/demo4.py
index 51859de..996a9ce 100644
--- a/demos/demo4.py
+++ b/demos/demo4.py
@@ -28,24 +28,24 @@
 # Pick the alpha and beta of the Gegenbauer polynomials
 alpha = .5
 # The first way of doing it is by directly supplying the weight function.
-wf = lambda(x): (1. - x ** 2.) ** (alpha - 0.5)
+wf = lambda x: (1. - x ** 2.) ** (alpha - 0.5)
 # Construct it:
 p = orthpol.OrthogonalPolynomial(degree,
                                 left=-1., right=1., # Domain
                                 wf=wf)
 # An orthogonal polynomial is though of as a function.
 # Here is how to get the number of inputs and outputs of that function
-print 'Number of inputs:', p.num_input
-print 'Number of outputs:', p.num_output
+print('Number of inputs:', p.num_input)
+print('Number of outputs:', p.num_output)
 # Test if the polynomials are normalized (i.e., their norm is 1.):
-print 'Is normalized:', p.is_normalized
+print('Is normalized:', p.is_normalized)
 # Get the degree of the polynomial:
-print 'Polynomial degree:', p.degree
+print('Polynomial degree:', p.degree)
 # Get the alpha-beta recursion coefficients:
-print 'Alpha:', p.alpha
-print 'Beta:', p.beta
+print('Alpha:', p.alpha)
+print('Beta:', p.beta)
 # The following should print a description of the polynomial
-print str(p)
+print(str(p))
 # Now you can evaluate the polynomial at any points you want:
 X = np.linspace(-1., 1., 100)
 # Here is the actual evaluation
@@ -57,7 +57,7 @@
 plt.xlabel('$x$', fontsize=16)
 plt.ylabel('$p_i(x)$', fontsize=16)
 plt.legend(['$p_{%d}(x)$' % i for i in range(p.num_output)], loc='best')
-print 'Close the window to continue...'
+print('Close the window to continue...')
 plt.show()
 # You may also compute the derivatives of the polynomials:
 dphi = p.d(X)
@@ -67,5 +67,5 @@
 plt.xlabel('$x$', fontsize=16)
 plt.ylabel(r'$\frac{dp_i(x)}{dx}$', fontsize=16)
 plt.legend([r'$\frac{p_{%d}(x)}{dx}$' % i for i in range(p.num_output)], loc='best')
-print 'Close the window to end demo...'
+print('Close the window to end demo...')
 plt.show()
diff --git a/demos/demo5.py b/demos/demo5.py
index a377f7b..792e592 100644
--- a/demos/demo5.py
+++ b/demos/demo5.py
@@ -29,24 +29,24 @@
 alpha = 2.
 beta = 5.
 # The first way of doing it is by directly supplying the weight function.
-wf = lambda(x): (1. - x) ** alpha * (1 + x) ** beta
+wf = lambda x: (1. - x) ** alpha * (1 + x) ** beta
 # Construct it:
 p = orthpol.OrthogonalPolynomial(degree,
                                 left=-1., right=1., # Domain
                                 wf=wf)
 # An orthogonal polynomial is though of as a function.
 # Here is how to get the number of inputs and outputs of that function
-print 'Number of inputs:', p.num_input
-print 'Number of outputs:', p.num_output
+print('Number of inputs:', p.num_input)
+print('Number of outputs:', p.num_output)
 # Test if the polynomials are normalized (i.e., their norm is 1.):
-print 'Is normalized:', p.is_normalized
+print('Is normalized:', p.is_normalized)
 # Get the degree of the polynomial:
-print 'Polynomial degree:', p.degree
+print('Polynomial degree:', p.degree)
 # Get the alpha-beta recursion coefficients:
-print 'Alpha:', p.alpha
-print 'Beta:', p.beta
+print('Alpha:', p.alpha)
+print('Beta:', p.beta)
 # The following should print a description of the polynomial
-print str(p)
+print(str(p))
 # Now you can evaluate the polynomial at any points you want:
 X = np.linspace(-1., 1., 100)
 # Here is the actual evaluation
@@ -58,7 +58,7 @@
 plt.xlabel('$x$', fontsize=16)
 plt.ylabel('$p_i(x)$', fontsize=16)
 plt.legend(['$p_{%d}(x)$' % i for i in range(p.num_output)], loc='best')
-print 'Close the window to continue...'
+print('Close the window to continue...')
 plt.show()
 # You may also compute the derivatives of the polynomials:
 dphi = p.d(X)
@@ -68,5 +68,5 @@
 plt.xlabel('$x$', fontsize=16)
 plt.ylabel(r'$\frac{dp_i(x)}{dx}$', fontsize=16)
 plt.legend([r'$\frac{p_{%d}(x)}{dx}$' % i for i in range(p.num_output)], loc='best')
-print 'Close the window to end demo...'
+print('Close the window to end demo...')
 plt.show()
diff --git a/demos/demo6.py b/demos/demo6.py
index 992967e..78f8025 100644
--- a/demos/demo6.py
+++ b/demos/demo6.py
@@ -26,24 +26,24 @@
 degree = 4
 
 # The first way of doing it is by directly supplying the weight function
-wf = lambda(x): 1.
+wf = lambda x: 1.
 # Construct it:
 p = orthpol.OrthogonalPolynomial(degree,
                                 left=0., right=1., # Domain
                                 wf=wf)
 # An orthogonal polynomial is though of as a function.
 # Here is how to get the number of inputs and outputs of that function
-print 'Number of inputs:', p.num_input
-print 'Number of outputs:', p.num_output
+print('Number of inputs:', p.num_input)
+print('Number of outputs:', p.num_output)
 # Test if the polynomials are normalized (i.e., their norm is 1.):
-print 'Is normalized:', p.is_normalized
+print('Is normalized:', p.is_normalized)
 # Get the degree of the polynomial:
-print 'Polynomial degree:', p.degree
+print('Polynomial degree:', p.degree)
 # Get the alpha-beta recursion coefficients:
-print 'Alpha:', p.alpha
-print 'Beta:', p.beta
+print('Alpha:', p.alpha)
+print('Beta:', p.beta)
 # The following should print a description of the polynomial
-print str(p)
+print(str(p))
 # Now you can evaluate the polynomial at any points you want:
 X = np.linspace(0., 1., 100)
 # Here is the actual evaluation
@@ -55,7 +55,7 @@
 plt.xlabel('$x$', fontsize=16)
 plt.ylabel('$p_i(x)$', fontsize=16)
 plt.legend(['$p_{%d}(x)$' % i for i in range(p.num_output)], loc='best')
-print 'Close the window to continue...'
+print('Close the window to continue...')
 plt.show()
 # You may also compute the derivatives of the polynomials:
 dphi = p.d(X)
@@ -65,5 +65,5 @@
 plt.xlabel('$x$', fontsize=16)
 plt.ylabel(r'$\frac{dp_i(x)}{dx}$', fontsize=16)
 plt.legend([r'$\frac{p_{%d}(x)}{dx}$' % i for i in range(p.num_output)], loc='best')
-print 'Close the window to end demo...'
+print('Close the window to end demo...')
 plt.show()
diff --git a/demos/demo7.py b/demos/demo7.py
index 57b300b..32af6f7 100644
--- a/demos/demo7.py
+++ b/demos/demo7.py
@@ -33,17 +33,17 @@
 p = orthpol.OrthogonalPolynomial(degree, rv=rv)
 # An orthogonal polynomial is though of as a function.
 # Here is how to get the number of inputs and outputs of that function
-print 'Number of inputs:', p.num_input
-print 'Number of outputs:', p.num_output
+print('Number of inputs:', p.num_input)
+print('Number of outputs:', p.num_output)
 # Test if the polynomials are normalized (i.e., their norm is 1.):
-print 'Is normalized:', p.is_normalized
+print('Is normalized:', p.is_normalized)
 # Get the degree of the polynomial:
-print 'Polynomial degree:', p.degree
+print('Polynomial degree:', p.degree)
 # Get the alpha-beta recursion coefficients:
-print 'Alpha:', p.alpha
-print 'Beta:', p.beta
+print('Alpha:', p.alpha)
+print('Beta:', p.beta)
 # The following should print a description of the polynomial
-print str(p)
+print(str(p))
 # Now you can evaluate the polynomial at any points you want:
 X = np.linspace(0., 1., 100)
 # Here is the actual evaluation
@@ -55,7 +55,7 @@
 plt.xlabel('$x$', fontsize=16)
 plt.ylabel('$p_i(x)$', fontsize=16)
 plt.legend(['$p_{%d}(x)$' % i for i in range(p.num_output)], loc='best')
-print 'Close the window to continue...'
+print('Close the window to continue...')
 plt.show()
 # You may also compute the derivatives of the polynomials:
 dphi = p.d(X)
@@ -65,5 +65,5 @@
 plt.xlabel('$x$', fontsize=16)
 plt.ylabel(r'$\frac{dp_i(x)}{dx}$', fontsize=16)
 plt.legend([r'$\frac{p_{%d}(x)}{dx}$' % i for i in range(p.num_output)], loc='best')
-print 'Close the window to end demo...'
+print('Close the window to end demo...')
 plt.show()
diff --git a/demos/demo8.py b/demos/demo8.py
index d8113a3..f2c0324 100644
--- a/demos/demo8.py
+++ b/demos/demo8.py
@@ -33,17 +33,17 @@
 p = orthpol.OrthogonalPolynomial(degree, rv=rv)
 # An orthogonal polynomial is though of as a function.
 # Here is how to get the number of inputs and outputs of that function
-print 'Number of inputs:', p.num_input
-print 'Number of outputs:', p.num_output
+print('Number of inputs:', p.num_input)
+print('Number of outputs:', p.num_output)
 # Test if the polynomials are normalized (i.e., their norm is 1.):
-print 'Is normalized:', p.is_normalized
+print('Is normalized:', p.is_normalized)
 # Get the degree of the polynomial:
-print 'Polynomial degree:', p.degree
+print('Polynomial degree:', p.degree)
 # Get the alpha-beta recursion coefficients:
-print 'Alpha:', p.alpha
-print 'Beta:', p.beta
+print('Alpha:', p.alpha)
+print('Beta:', p.beta)
 # The following should print a description of the polynomial
-print str(p)
+print(str(p))
 # Now you can evaluate the polynomial at any points you want:
 X = np.linspace(-2., 2., 100)
 # Here is the actual evaluation
@@ -55,7 +55,7 @@
 plt.xlabel('$x$', fontsize=16)
 plt.ylabel('$p_i(x)$', fontsize=16)
 plt.legend(['$p_{%d}(x)$' % i for i in range(p.num_output)], loc='best')
-print 'Close the window to continue...'
+print('Close the window to continue...')
 plt.show()
 # You may also compute the derivatives of the polynomials:
 dphi = p.d(X)
@@ -65,5 +65,5 @@
 plt.xlabel('$x$', fontsize=16)
 plt.ylabel(r'$\frac{dp_i(x)}{dx}$', fontsize=16)
 plt.legend([r'$\frac{p_{%d}(x)}{dx}$' % i for i in range(p.num_output)], loc='best')
-print 'Close the window to end demo...'
+print('Close the window to end demo...')
 plt.show()
diff --git a/demos/demo9.py b/demos/demo9.py
index 98a1ac5..5e5410a 100644
--- a/demos/demo9.py
+++ b/demos/demo9.py
@@ -33,17 +33,17 @@
 p = orthpol.OrthogonalPolynomial(degree, rv=rv)
 # An orthogonal polynomial is though of as a function.
 # Here is how to get the number of inputs and outputs of that function
-print 'Number of inputs:', p.num_input
-print 'Number of outputs:', p.num_output
+print('Number of inputs:', p.num_input)
+print('Number of outputs:', p.num_output)
 # Test if the polynomials are normalized (i.e., their norm is 1.):
-print 'Is normalized:', p.is_normalized
+print('Is normalized:', p.is_normalized)
 # Get the degree of the polynomial:
-print 'Polynomial degree:', p.degree
+print('Polynomial degree:', p.degree)
 # Get the alpha-beta recursion coefficients:
-print 'Alpha:', p.alpha
-print 'Beta:', p.beta
+print('Alpha:', p.alpha)
+print('Beta:', p.beta)
 # The following should print a description of the polynomial
-print str(p)
+print(str(p))
 # Now you can evaluate the polynomial at any points you want:
 X = np.linspace(-9., -1., 100)
 # Here is the actual evaluation
@@ -55,7 +55,7 @@
 plt.xlabel('$x$', fontsize=16)
 plt.ylabel('$p_i(x)$', fontsize=16)
 plt.legend(['$p_{%d}(x)$' % i for i in range(p.num_output)], loc='best')
-print 'Close the window to continue...'
+print('Close the window to continue...')
 plt.show()
 # You may also compute the derivatives of the polynomials:
 dphi = p.d(X)
@@ -65,5 +65,5 @@
 plt.xlabel('$x$', fontsize=16)
 plt.ylabel(r'$\frac{dp_i(x)}{dx}$', fontsize=16)
 plt.legend([r'$\frac{p_{%d}(x)}{dx}$' % i for i in range(p.num_output)], loc='best')
-print 'Close the window to end demo...'
+print('Close the window to end demo...')
 plt.show()
diff --git a/orthpol/__init__.py b/orthpol/__init__.py
index 49d7b55..e01371f 100644
--- a/orthpol/__init__.py
+++ b/orthpol/__init__.py
@@ -11,6 +11,6 @@
 __all__ = ['OrthogonalPolynomial', 'ProductBasis', 'QuadratureRule']
 
 
-import _orthpol
+import orthpol._orthpol as _orthopol
 from ._quadrature_rule import *
 from ._orthogonal_polynomial import *
diff --git a/orthpol/_lancz.py b/orthpol/_lancz.py
index 7b2c49b..c201033 100644
--- a/orthpol/_lancz.py
+++ b/orthpol/_lancz.py
@@ -12,7 +12,7 @@
 
 
 import numpy as np
-import _orthpol as orthpol
+import orthpol._orthpol as orthpol
 
 
 def lancz(x, w, n):
@@ -33,7 +33,7 @@ def lancz(x, w, n):
 if __name__ == '__main__':
     x = np.linspace(0., 1, 100)
     w = np.ones(100)
-    print x, w
+    print(x, w)
     alpha, beta = lancz(x, w, 10)
-    print alpha
-    print beta
+    print(alpha)
+    print(beta)
diff --git a/orthpol/_orthogonal_polynomial.py b/orthpol/_orthogonal_polynomial.py
index 4a776a6..9678f21 100644
--- a/orthpol/_orthogonal_polynomial.py
+++ b/orthpol/_orthogonal_polynomial.py
@@ -1,363 +1,365 @@
-"""
-Describes an orthogonal polynomial.
-
-Author:
-    Ilias Bilionis
-
-Date:
-    7/25/2013
-"""
-
-
-__all__ = ['OrthogonalPolynomial', 'ProductBasis']
-
-
-import numpy as np
-import math
-import itertools
-from ._quadrature_rule import *
-from ._lancz import *
-import _orthpol as orthpol
-
-
-class OrthogonalPolynomial(object):
-
-    """1D Orthogonal Polynomial via recursive relation.
-
-    A polynomial is of course a function.
-    """
-
-    # Recurrence coefficient alpha
-    _alpha = None
-
-    # Recurrence coefficient beta
-    _beta = None
-
-    # Recurrence coefficient gamma
-    _gamma = None
-
-    # Is the polynomial normalized
-    _is_normalized = None
-
-    # The number of inputs
-    _num_input = None
-
-    # The number of outputs
-    _num_output = None
-
-    @property
-    def degree(self):
-        """Return the degree of the polynomial."""
-        return self.alpha.shape[0] - 1
-
-    @property
-    def alpha(self):
-        return self._alpha
-
-    @property
-    def beta(self):
-        return self._beta
-
-    @property
-    def gamma(self):
-        return self._gamma
-
-    @property
-    def is_normalized(self):
-        return self._is_normalized
-
-    @property
-    def num_input(self):
-        return self._num_input
-
-    @property
-    def num_output(self):
-        return self._num_output
-
-    def __init__(self, degree, rv=None, left=-1, right=1, wf=lambda(x): 1.,
-                 ncap=50, quad=None,
-                 name='Orthogonal Polynomial'):
-        """Construct the polynomial.
-
-        Keyword Arguments:
-            rv      ---     If not None, then it is assumed to be a
-                            RandomVariable object which used to define
-                            the support interval and the pdf.
-            degree  ---     The degree of the polynomial.
-            left    ---     The left end of the interval.
-            right   ---     The right end of the interval.
-            wf      ---     The weight function. The default is the identity.
-            ncap    ---     The number of quadrature points.
-            quad    ---     A quadrature rule you might want to use in
-                            case you are not satisfied with the default
-                            one.
-            name    ---     A name for the polynomial.
-        """
-        self.__name__ = name
-        if rv is not None:
-            left, right = rv.interval(1)
-            wf = rv.pdf
-        if quad is None:
-            quad = QuadratureRule(left=left, right=right, wf=wf, ncap=ncap)
-        self._alpha, self._beta = lancz(quad.x, quad.w, degree + 1)
-        self._gamma = np.ones(self.degree + 1, dtype='float64')
-        self.normalize()
-        self._num_input = 1
-        self._num_output = self.degree + 1
-
-    def __call__(self, x):
-        """Evaluate the function at x."""
-        return orthpol.poly_eval_all(x, self.alpha, self.beta, self.gamma)
-
-    def d(self, x):
-        return orthpol.poly_deval_all(x, self.alpha, self.beta, self.gamma)
-
-    def _eval(self, x):
-        """Evaluate the polynomial basis at x."""
-        return orthpol.poly_eval(x, self.alpha, self.beta, self.gamma)
-
-    def _d_eval(self, x):
-        """Evaluate the derivative of the polynomial.
-
-        Arguments:
-            x   ---     The input point(s).
-        """
-        return orthpol.poly_deval(x, self.alpha, self.beta, self.gamma)
-
-    def _evaluate_square_norms(self):
-        """Evaluate the square norms of the polynomials."""
-        s_norm = np.zeros(self.degree + 1)
-        s_norm[0] = self.beta[0] / (self.gamma[0] ** 2)
-        for i in range(1, self.degree + 1):
-            s_norm[i] = (self.beta[i] / self.gamma[i]) * s_norm[i - 1]
-        return s_norm
-
-    def normalize(self):
-        """Normalize the polynomials."""
-        self._beta, self._gamma = orthpol.poly_normalize(self.beta, self.gamma)
-        self._is_normalized = True
-
-    def __str__(self):
-        """Return a string representation of the object."""
-        s = self.__name__ + '\n'
-        s += ' alpha: ' + str(self.alpha) + '\n'
-        s += ' beta: ' + str(self.beta) + '\n'
-        s += ' gamma: ' + str(self.gamma) + '\n'
-        s += ' normalized: ' + str(self.is_normalized)
-        return s
-
-
-class ProductBasis(object):
-
-    """A multi-input orthogonal polynomial basis."""
-
-    # A container of polynomials
-    _polynomials = None
-
-    # The total order of the basis
-    _degree = None
-
-    # An array of basis terms
-    _terms = None
-
-    # The number of terms up to each order
-    _num_terms = None
-
-    # The number of inputs
-    _num_input = None
-
-    # The number of outputs
-    _num_output = None
-
-    @property
-    def polynomials(self):
-        return self._polynomials
-
-    @property
-    def degree(self):
-        return self._degree
-
-    @property
-    def terms(self):
-        return self._terms
-
-    @property
-    def num_terms(self):
-        return self._num_terms
-
-    @property
-    def num_input(self):
-        return self._num_input
-
-    @property
-    def num_output(self):
-        return self._num_output
-
-    def __init__(self, rvs=None, degree=1, polynomials=None, ncap=50,
-                 quad=None, name='Product basis'):
-        """Initialize the object.
-
-        Keyword Argument
-            rvs             ---     If not None, then it is assumed to
-                                    be a list of random variables.
-            degree          ---     The total degree of the basis. Each
-                                    one of the polynomials will have this
-                                    degree.
-            polynomials     ---     We only look at this if rv is None.
-                                    A collection of 1D orthogonal
-                                    polynomials.
-            ncap    ---     The number of quadrature points.
-            quad    ---     A quadrature rule you might want to use in
-                            case you are not satisfied with the default
-                            one.
-            name            ---     A name for the basis.
-        """
-        self.__name__ = name
-        assert isinstance(degree, int)
-        assert degree >= 0
-        if rvs is not None:
-            assert isinstance(rvs, list) or isinstance(rvs, tuple)
-            polynomials = [OrthogonalPolynomial(degree, rv=r, ncap=ncap,
-                                                quad=quad) for r in rvs]
-        assert (isinstance(polynomials, tuple) or
-                isinstance(polynomials, list))
-        for p in polynomials:
-            assert isinstance(p, OrthogonalPolynomial)
-        self._polynomials = polynomials
-        # Find the total order of the basis
-        self._degree = max([p.degree for p in polynomials])
-        # The number of inputs
-        self._num_input = len(polynomials)
-        # Compute the basis terms
-        self._num_output = self._compute_basis_terms()
-
-    def _compute_basis_terms(self):
-        """Compute the basis terms.
-
-        The following is taken from Stokhos.
-
-        The approach here for ordering the terms is inductive on the total
-        order p.  We get the terms of total order p from the terms of total
-        order p-1 by incrementing the orders of the first dimension by 1.
-        We then increment the orders of the second dimension by 1 for all of the
-        terms whose first dimension order is 0.  We then repeat for the third
-        dimension whose first and second dimension orders are 0, and so on.
-        How this is done is most easily illustrated by an example of dimension 3:
-
-        Order  terms   cnt  Order  terms   cnt
-        0    0 0 0          4    4 0 0  15 5 1
-                                 3 1 0
-        1    1 0 0  3 2 1        3 0 1
-             0 1 0               2 2 0
-             0 0 1               2 1 1
-                                 2 0 2
-        2    2 0 0  6 3 1        1 3 0
-             1 1 0               1 2 1
-             1 0 1               1 1 2
-             0 2 0               1 0 3
-             0 1 1               0 4 0
-             0 0 2               0 3 1
-                                 0 2 2
-        3    3 0 0  10 4 1       0 1 3
-             2 1 0               0 0 4
-             2 0 1
-             1 2 0
-             1 1 1
-             1 0 2
-             0 3 0
-             0 2 1
-             0 1 2
-             0 0 3
-        """
-        # Number of inputs
-        num_dim = len(self.polynomials)
-
-        # Temporary array of terms grouped in terms of same order
-        terms_order = [[] for i in range(self.degree + 1)]
-
-        # Store number of terms up to each order
-        self._num_terms = np.zeros(self.degree + 2, dtype='i')
-
-        # Set order zero
-        terms_order[0] = ([np.zeros(num_dim, dtype='i')])
-        self.num_terms[0] = 1
-
-        # The array cnt stores the number of terms we need to
-        # increment for each dimension.
-        cnt = np.zeros(num_dim, dtype='i')
-        for j, p in itertools.izip(range(num_dim), self.polynomials):
-            if p.degree >= 1:
-                cnt[j] = 1
-
-        cnt_next = np.zeros(num_dim, dtype='i')
-        term = np.zeros(num_dim, dtype='i')
-
-        # Number of basis functions
-        num_basis = 1
-
-        # Loop over orders
-        for k in range(1, self.degree + 1):
-            self.num_terms[k] = self.num_terms[k - 1]
-            # Stores the inde of the term we are copying
-            prev = 0
-            # Loop over dimensions
-            for j, p in itertools.izip(range(num_dim), self.polynomials):
-                # Increment orders of cnt[j] terms for dimension j
-                for i in range(cnt[j]):
-                    if terms_order[k - 1][prev + i][j] < p.degree:
-                        term = terms_order[k - 1][prev + i].copy()
-                        term[j] += 1
-                        terms_order[k].append(term)
-                        num_basis += 1
-                        self.num_terms[k] += 1
-                        for l in range(j + 1):
-                            cnt_next[l] += 1
-                if j < num_dim - 1:
-                    prev += cnt[j] - cnt[j + 1]
-            cnt[:] = cnt_next
-            cnt_next[:] = 0
-        self.num_terms[self.degree + 1] = num_basis
-        # Copy into final terms array
-        self._terms = []
-        for k in range(self.degree + 1):
-            num_k = len(terms_order[k])
-            for j in range(num_k):
-                self._terms.append(terms_order[k][j])
-        return num_basis
-
-    def __str__(self):
-        """Return a string representation of the object."""
-        s = self.__name__ + '\n'
-        s += ' sz = ' + str(self.num_output) + '\n'
-        for i in range(self.num_output):
-            s += ' ' + str(i) + ': '
-            for j in range(self.num_input):
-                s += str(self.terms[i][j]) + ' '
-            s += '\n'
-        s +=  ' num_terms = '
-        for i in range(self.degree + 1):
-            s += str(self.num_terms[i]) + ' '
-        return s
-
-    def __call__(self, x):
-        num_pt = x.shape[0]
-        basis_eval_tmp = [self.polynomials[j](x[:, j])
-                          for j in range(self.num_input)]
-        phi = np.ndarray((num_pt, self.num_output))
-        for k in range(self.num_output):
-            phi[:, k] = 1.
-            for j in range(self.num_input):
-                phi[:, k] *= basis_eval_tmp[j][:, self.terms[k][j]]
-        return phi
-
-    def _eval(self, x, hyp):
-        """Evaluate the polynomials at x."""
-        phi = np.ndarray(self.num_output)
-        basis_eval_tmp = [[] for j in range(self.num_input)]
-        for j in range(self.num_input):
-            basis_eval_tmp[j] = self.polynomials[j](x[j]).flatten()
-        for k in range(self.num_output):
-            phi[k] = 1.
-            for j in range(self.num_input):
-                phi[k] *= basis_eval_tmp[j][self.terms[k][j]]
-        return phi
+"""
+Describes an orthogonal polynomial.
+
+Author:
+    Ilias Bilionis
+
+Date:
+    7/25/2013
+"""
+
+
+__all__ = ['OrthogonalPolynomial', 'ProductBasis']
+
+
+import numpy as np
+import math
+import itertools
+from ._quadrature_rule import *
+from ._lancz import *
+import orthpol._orthpol as orthpol
+
+
+class OrthogonalPolynomial(object):
+
+    """1D Orthogonal Polynomial via recursive relation.
+
+    A polynomial is of course a function.
+    """
+
+    # Recurrence coefficient alpha
+    _alpha = None
+
+    # Recurrence coefficient beta
+    _beta = None
+
+    # Recurrence coefficient gamma
+    _gamma = None
+
+    # Is the polynomial normalized
+    _is_normalized = None
+
+    # The number of inputs
+    _num_input = None
+
+    # The number of outputs
+    _num_output = None
+
+    @property
+    def degree(self):
+        """Return the degree of the polynomial."""
+        return self.alpha.shape[0] - 1
+
+    @property
+    def alpha(self):
+        return self._alpha
+
+    @property
+    def beta(self):
+        return self._beta
+
+    @property
+    def gamma(self):
+        return self._gamma
+
+    @property
+    def is_normalized(self):
+        return self._is_normalized
+
+    @property
+    def num_input(self):
+        return self._num_input
+
+    @property
+    def num_output(self):
+        return self._num_output
+
+    def __init__(self, degree, rv=None, left=-1, right=1, wf=None,
+                 ncap=50, quad=None,
+                 name='Orthogonal Polynomial'):
+        """Construct the polynomial.
+
+        Keyword Arguments:
+            rv      ---     If not None, then it is assumed to be a
+                            RandomVariable object which used to define
+                            the support interval and the pdf.
+            degree  ---     The degree of the polynomial.
+            left    ---     The left end of the interval.
+            right   ---     The right end of the interval.
+            wf      ---     The weight function. The default is the identity.
+            ncap    ---     The number of quadrature points.
+            quad    ---     A quadrature rule you might want to use in
+                            case you are not satisfied with the default
+                            one.
+            name    ---     A name for the polynomial.
+        """
+        self.__name__ = name
+        if wf is None:
+            wf = lambda x: 1.0
+        if rv is not None:
+            left, right = rv.interval(1)
+            wf = rv.pdf
+        if quad is None:
+            quad = QuadratureRule(left=left, right=right, wf=wf, ncap=ncap)
+        self._alpha, self._beta = lancz(quad.x, quad.w, degree + 1)
+        self._gamma = np.ones(self.degree + 1, dtype='float64')
+        self.normalize()
+        self._num_input = 1
+        self._num_output = self.degree + 1
+
+    def __call__(self, x):
+        """Evaluate the function at x."""
+        return orthpol.poly_eval_all(x, self.alpha, self.beta, self.gamma)
+
+    def d(self, x):
+        return orthpol.poly_deval_all(x, self.alpha, self.beta, self.gamma)
+
+    def _eval(self, x):
+        """Evaluate the polynomial basis at x."""
+        return orthpol.poly_eval(x, self.alpha, self.beta, self.gamma)
+
+    def _d_eval(self, x):
+        """Evaluate the derivative of the polynomial.
+
+        Arguments:
+            x   ---     The input point(s).
+        """
+        return orthpol.poly_deval(x, self.alpha, self.beta, self.gamma)
+
+    def _evaluate_square_norms(self):
+        """Evaluate the square norms of the polynomials."""
+        s_norm = np.zeros(self.degree + 1)
+        s_norm[0] = self.beta[0] / (self.gamma[0] ** 2)
+        for i in range(1, self.degree + 1):
+            s_norm[i] = (self.beta[i] / self.gamma[i]) * s_norm[i - 1]
+        return s_norm
+
+    def normalize(self):
+        """Normalize the polynomials."""
+        self._beta, self._gamma = orthpol.poly_normalize(self.beta, self.gamma)
+        self._is_normalized = True
+
+    def __str__(self):
+        """Return a string representation of the object."""
+        s = self.__name__ + '\n'
+        s += ' alpha: ' + str(self.alpha) + '\n'
+        s += ' beta: ' + str(self.beta) + '\n'
+        s += ' gamma: ' + str(self.gamma) + '\n'
+        s += ' normalized: ' + str(self.is_normalized)
+        return s
+
+
+class ProductBasis(object):
+
+    """A multi-input orthogonal polynomial basis."""
+
+    # A container of polynomials
+    _polynomials = None
+
+    # The total order of the basis
+    _degree = None
+
+    # An array of basis terms
+    _terms = None
+
+    # The number of terms up to each order
+    _num_terms = None
+
+    # The number of inputs
+    _num_input = None
+
+    # The number of outputs
+    _num_output = None
+
+    @property
+    def polynomials(self):
+        return self._polynomials
+
+    @property
+    def degree(self):
+        return self._degree
+
+    @property
+    def terms(self):
+        return self._terms
+
+    @property
+    def num_terms(self):
+        return self._num_terms
+
+    @property
+    def num_input(self):
+        return self._num_input
+
+    @property
+    def num_output(self):
+        return self._num_output
+
+    def __init__(self, rvs=None, degree=1, polynomials=None, ncap=50,
+                 quad=None, name='Product basis'):
+        """Initialize the object.
+
+        Keyword Argument
+            rvs             ---     If not None, then it is assumed to
+                                    be a list of random variables.
+            degree          ---     The total degree of the basis. Each
+                                    one of the polynomials will have this
+                                    degree.
+            polynomials     ---     We only look at this if rv is None.
+                                    A collection of 1D orthogonal
+                                    polynomials.
+            ncap    ---     The number of quadrature points.
+            quad    ---     A quadrature rule you might want to use in
+                            case you are not satisfied with the default
+                            one.
+            name            ---     A name for the basis.
+        """
+        self.__name__ = name
+        assert isinstance(degree, int)
+        assert degree >= 0
+        if rvs is not None:
+            assert isinstance(rvs, list) or isinstance(rvs, tuple)
+            polynomials = [OrthogonalPolynomial(degree, rv=r, ncap=ncap,
+                                                quad=quad) for r in rvs]
+        assert (isinstance(polynomials, tuple) or
+                isinstance(polynomials, list))
+        for p in polynomials:
+            assert isinstance(p, OrthogonalPolynomial)
+        self._polynomials = polynomials
+        # Find the total order of the basis
+        self._degree = max([p.degree for p in polynomials])
+        # The number of inputs
+        self._num_input = len(polynomials)
+        # Compute the basis terms
+        self._num_output = self._compute_basis_terms()
+
+    def _compute_basis_terms(self):
+        """Compute the basis terms.
+
+        The following is taken from Stokhos.
+
+        The approach here for ordering the terms is inductive on the total
+        order p.  We get the terms of total order p from the terms of total
+        order p-1 by incrementing the orders of the first dimension by 1.
+        We then increment the orders of the second dimension by 1 for all of the
+        terms whose first dimension order is 0.  We then repeat for the third
+        dimension whose first and second dimension orders are 0, and so on.
+        How this is done is most easily illustrated by an example of dimension 3:
+
+        Order  terms   cnt  Order  terms   cnt
+        0    0 0 0          4    4 0 0  15 5 1
+                                 3 1 0
+        1    1 0 0  3 2 1        3 0 1
+             0 1 0               2 2 0
+             0 0 1               2 1 1
+                                 2 0 2
+        2    2 0 0  6 3 1        1 3 0
+             1 1 0               1 2 1
+             1 0 1               1 1 2
+             0 2 0               1 0 3
+             0 1 1               0 4 0
+             0 0 2               0 3 1
+                                 0 2 2
+        3    3 0 0  10 4 1       0 1 3
+             2 1 0               0 0 4
+             2 0 1
+             1 2 0
+             1 1 1
+             1 0 2
+             0 3 0
+             0 2 1
+             0 1 2
+             0 0 3
+        """
+        # Number of inputs
+        num_dim = len(self.polynomials)
+
+        # Temporary array of terms grouped in terms of same order
+        terms_order = [[] for i in range(self.degree + 1)]
+
+        # Store number of terms up to each order
+        self._num_terms = np.zeros(self.degree + 2, dtype='i')
+
+        # Set order zero
+        terms_order[0] = ([np.zeros(num_dim, dtype='i')])
+        self.num_terms[0] = 1
+
+        # The array cnt stores the number of terms we need to
+        # increment for each dimension.
+        cnt = np.zeros(num_dim, dtype='i')
+        for j, p in zip(list(range(num_dim)), self.polynomials):
+            if p.degree >= 1:
+                cnt[j] = 1
+
+        cnt_next = np.zeros(num_dim, dtype='i')
+        term = np.zeros(num_dim, dtype='i')
+
+        # Number of basis functions
+        num_basis = 1
+
+        # Loop over orders
+        for k in range(1, self.degree + 1):
+            self.num_terms[k] = self.num_terms[k - 1]
+            # Stores the inde of the term we are copying
+            prev = 0
+            # Loop over dimensions
+            for j, p in zip(list(range(num_dim)), self.polynomials):
+                # Increment orders of cnt[j] terms for dimension j
+                for i in range(cnt[j]):
+                    if terms_order[k - 1][prev + i][j] < p.degree:
+                        term = terms_order[k - 1][prev + i].copy()
+                        term[j] += 1
+                        terms_order[k].append(term)
+                        num_basis += 1
+                        self.num_terms[k] += 1
+                        for l in range(j + 1):
+                            cnt_next[l] += 1
+                if j < num_dim - 1:
+                    prev += cnt[j] - cnt[j + 1]
+            cnt[:] = cnt_next
+            cnt_next[:] = 0
+        self.num_terms[self.degree + 1] = num_basis
+        # Copy into final terms array
+        self._terms = []
+        for k in range(self.degree + 1):
+            num_k = len(terms_order[k])
+            for j in range(num_k):
+                self._terms.append(terms_order[k][j])
+        return num_basis
+
+    def __str__(self):
+        """Return a string representation of the object."""
+        s = self.__name__ + '\n'
+        s += ' sz = ' + str(self.num_output) + '\n'
+        for i in range(self.num_output):
+            s += ' ' + str(i) + ': '
+            for j in range(self.num_input):
+                s += str(self.terms[i][j]) + ' '
+            s += '\n'
+        s +=  ' num_terms = '
+        for i in range(self.degree + 1):
+            s += str(self.num_terms[i]) + ' '
+        return s
+
+    def __call__(self, x):
+        num_pt = x.shape[0]
+        basis_eval_tmp = [self.polynomials[j](x[:, j])
+                          for j in range(self.num_input)]
+        phi = np.ndarray((num_pt, self.num_output))
+        for k in range(self.num_output):
+            phi[:, k] = 1.
+            for j in range(self.num_input):
+                phi[:, k] *= basis_eval_tmp[j][:, self.terms[k][j]]
+        return phi
+
+    def _eval(self, x, hyp):
+        """Evaluate the polynomials at x."""
+        phi = np.ndarray(self.num_output)
+        basis_eval_tmp = [[] for j in range(self.num_input)]
+        for j in range(self.num_input):
+            basis_eval_tmp[j] = self.polynomials[j](x[j]).flatten()
+        for k in range(self.num_output):
+            phi[k] = 1.
+            for j in range(self.num_input):
+                phi[k] *= basis_eval_tmp[j][self.terms[k][j]]
+        return phi
diff --git a/orthpol/_quadrature_rule.py b/orthpol/_quadrature_rule.py
index edac13f..dbbcf17 100644
--- a/orthpol/_quadrature_rule.py
+++ b/orthpol/_quadrature_rule.py
@@ -14,7 +14,7 @@
 
 import numpy as np
 import math
-import _orthpol as orthpol
+import orthpol._orthpol as orthpol
 
 
 def symtr(t):
@@ -77,7 +77,7 @@ def w(self):
     def num_quad(self):
         return self._x.shape[0]
 
-    def __init__(self, left=-1, right=1, wf=lambda(x): 1., ncap=500,
+    def __init__(self, left=-1, right=1, wf=lambda x: 1., ncap=500,
                  name='Quadrature Rule'):
         """Construct a quadrature rule.
 
@@ -90,7 +90,7 @@ def __init__(self, left=-1, right=1, wf=lambda(x): 1., ncap=500,
         """
         x, w = fejer(ncap)
         if wf is None:
-            wf = lambda(x): np.ones(x.shape)
+            wf = lambda x: np.ones(x.shape)
         if math.isinf(left) and math.isinf(right):
             phi, dphi = symtr(x)
             self._x = phi
diff --git a/unittests/test.py b/unittests/test.py
index b2c704b..590caca 100644
--- a/unittests/test.py
+++ b/unittests/test.py
@@ -25,7 +25,7 @@ def test_normalization(self):
         p = orthpol.ProductBasis(rvs, degree = 5)
         U = st.uniform.rvs(loc = -1, scale = 2., size = (10000,2))
         P = p(U)
-        print np.sum(P**2, axis = 0) / 10000
+        print(np.sum(P**2, axis = 0) / 10000)
 
 if __name__ == '__main__':
     unittest.main()