Add more tests and examples.

Fix the orthogonality test.

diff --git a/mjo/orthogonal_polynomials.py b/mjo/orthogonal_polynomials.py
index edea32bf7a0788b7574f3dd9c7f7298e395fff1b..808a558eb925521148306cb7315fbd9df07bd32a 100644 (file)
--- a/mjo/orthogonal_polynomials.py
+++ b/mjo/orthogonal_polynomials.py
@@ -1,5 +1,4 @@
 from sage.all import *
-from sage.symbolic.expression import is_Expression
 
 def legendre_p(n, x, a = -1, b = 1):
     """
@@ -29,6 +28,36 @@ def legendre_p(n, x, a = -1, b = 1):
     returned. Otherwise, the value of the ``n``th polynomial at ``x``
     will be returned.
 
+    EXAMPLES:
+
+    Create the standard Legendre polynomials in `x`::
+
+        sage: legendre_p(0,x)
+        1
+        sage: legendre_p(1,x)
+        x
+
+    Reuse the variable from a polynomial ring::
+        sage: P.<t> = QQ[]
+        sage: legendre_p(2,t).simplify_rational()
+        3/2*t^2 - 1/2
+
+    If ``x`` is a real number, the result should be as well::
+
+        sage: legendre_p(3, 1.1)
+        1.67750000000000
+
+    Similarly for complex numbers::
+
+        sage: legendre_p(3, 1 + I)
+        7/2*I - 13/2
+
+    Even matrices work::
+
+        sage: legendre_p(3, MatrixSpace(ZZ, 2)([1, 2, -4, 7]))
+        [-179  242]
+        [-484  547]
+
     TESTS:
 
     We should agree with Maxima for all `n`::
@@ -55,10 +84,11 @@ def legendre_p(n, x, a = -1, b = 1):
         1
 
     Two different polynomials should be orthogonal with respect to the
-    inner product over [a,b]::
+    inner product over `[a,b]`. Note that this test can fail if QQ is
+    replaced with RR, since integrate() can return NaN::
 
-        sage: a = RR.random_element()
-        sage: b = RR.random_element()
+        sage: a = QQ.random_element()
+        sage: b = QQ.random_element()
         sage: j = ZZ.random_element(20)
         sage: k = j + 1
         sage: Pj = legendre_p(j, x, a, b)
@@ -81,6 +111,13 @@ def legendre_p(n, x, a = -1, b = 1):
         sage: bool(legendre_p(3, x, 0, 1) == p3)
         True
 
+    The zeroth polynomial is an element of the ring that we're working
+    in::
+
+        sage: legendre_p(0, MatrixSpace(ZZ, 2)([1, 2, -4, 7]))
+        [1 0]
+        [0 1]
+
     """
     if not n in ZZ:
         raise TypeError('n must be a natural number')
@@ -91,8 +128,15 @@ def legendre_p(n, x, a = -1, b = 1):
     if not (a in RR and b in RR):
         raise TypeError('both `a` and `b` must be real numbers')
 
-    a = RR(a)
-    b = RR(b)
+    if n == 0:
+        # Easy base case, save time. Attempt to return a value in the
+        # same field/ring as `x`.
+        return x.parent()(1)
+
+    # Even though we know a,b are real we use the symbolic ring. This
+    # lets us return pretty expressions where possible.
+    a = SR(a)
+    b = SR(b)
     n = ZZ(n) # Ensure that 1/(2**n) is not integer division.
     dn = 1/(2**n)
 
@@ -111,9 +155,4 @@ def legendre_p(n, x, a = -1, b = 1):
     # From Abramowitz & Stegun, (22.3.2) with alpha = beta = 0.
     P = dn * sum([ c(m)*g(m) for m in range(0,n+1) ])
 
-    # If `x` is a symbolic expression, we want to return a symbolic
-    # expression (even if that expression is e.g. `1`).
-    if is_Expression(x):
-        P = SR(P)
-
     return P