Add more positive operator examples from the paper.

[sage.d.git] / mjo / orthogonal_polynomials.py

diff --git a/mjo/orthogonal_polynomials.py b/mjo/orthogonal_polynomials.py
index 808a558eb925521148306cb7315fbd9df07bd32a..e8e9210f46437755dc99333ce3159f9b7abd5a71 100644 (file)
--- a/mjo/orthogonal_polynomials.py
+++ b/mjo/orthogonal_polynomials.py
@@ -39,7 +39,7 @@ def legendre_p(n, x, a = -1, b = 1):
 
     Reuse the variable from a polynomial ring::
         sage: P.<t> = QQ[]
-        sage: legendre_p(2,t).simplify_rational()
+        sage: legendre_p(2,t)
         3/2*t^2 - 1/2
 
     If ``x`` is a real number, the result should be as well::
@@ -58,6 +58,31 @@ def legendre_p(n, x, a = -1, b = 1):
         [-179  242]
         [-484  547]
 
+    And finite field elements::
+
+        sage: legendre_P(3, GF(11)(5))
+        8
+
+    Solve a simple least squares problem over `[-\pi, \pi]`::
+
+        sage: a = -pi
+        sage: b = pi
+        sage: def inner_product(v1, v2):
+        ...       return integrate(v1*v2, x, a, b)
+        ...
+        sage: def norm(v):
+        ...       return sqrt(inner_product(v,v))
+        ...
+        sage: def project(basis, v):
+        ...       return sum([ inner_product(v, b)*b/norm(b)**2
+        ...                    for b in basis])
+        ...
+        sage: f = sin(x)
+        sage: legendre_basis = [ legendre_p(k, x, a, b) for k in range(0,4) ]
+        sage: proj = project(legendre_basis, f)
+        sage: proj.simplify_trig()
+        5/2*(7*(pi^2 - 15)*x^3 - 3*(pi^4 - 21*pi^2)*x)/pi^6
+
     TESTS:
 
     We should agree with Maxima for all `n`::
@@ -133,17 +158,30 @@ def legendre_p(n, x, a = -1, b = 1):
         # 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)
+    # Even though we know a,b are real we use a different ring. We
+    # prefer ZZ so that we can support division of finite field
+    # elements by (a-b). Eventually this should be supported for QQ as
+    # well, although it does not work at the moment. The preference of
+    # SR over RR is to return something attractive when e.g. a=pi.
+    if a in ZZ:
+        a = ZZ(a)
+    else:
+        a = SR(a)
+
+    if b in ZZ:
+        b = ZZ(b)
+    else:
+        b = SR(b)
+
+    # Ensure that (2**n) is an element of ZZ. This is used later --
+    # we can divide finite field elements by integers but we can't
+    # multiply them by rationals at the moment.
+    n = ZZ(n)
 
     def phi(t):
-        # This is an affine map from [a,b] into [-1,1] and so preserves
-        # orthogonality.
-        return (2 / (b-a))*t + 1 - (2*b)/(b-a)
+        # This is an affine map from [a,b] into [-1,1] and so
+        # preserves orthogonality.
+        return (a + b - 2*t)/(a - b)
 
     def c(m):
         return binomial(n,m)*binomial(n, n-m)
@@ -153,6 +191,7 @@ def legendre_p(n, x, a = -1, b = 1):
         return ( ((phi(x) - 1)**(n-m)) * (phi(x) + 1)**m )
 
     # From Abramowitz & Stegun, (22.3.2) with alpha = beta = 0.
-    P = dn * sum([ c(m)*g(m) for m in range(0,n+1) ])
+    # Also massaged to support finite field elements.
+    P = sum([ c(m)*g(m) for m in range(0,n+1) ])/(2**n)
 
     return P