recent changes to orthopoly and kron_delta, was: Two Patches for Maxima (on SF)

I would still find it useful to have a function

orthopoly_leadcoef(f,args)

similar to

orthopoly_weight(f,args)

returning the highest power of the respective polys.
as a substitute for an obviously not working

ratcoef(legendre_p(n,x),x,n);

Andre


http://www.ma.utexas.edu/pipermail/maxima/2008/012524.html

On 08/22/2012 08:50 PM, Robert Dodier wrote:
> I've made some changes recently (committed & pushed) to the orthopoly
> package and kron_delta. Details are shown in the attached log. In
> summary:
>
>   * fixed bugs in definition of kron_delta Lisp function
>   * make comments in orthopoly agree with code
>   * return simplifiable kron_delta expressions from orthopoly
>   * handle assoc_legendre_p in orthopoly_weight
>   * express legendre_{p,q} recurrences in terms of n + 1, n, and n - 1
>
> About the last one, I suppose it is to some degree a matter of taste.
> However, now the legendre_p & legendre_q recurrences are analogous to
> other recurrences known to orthopoly_recur, and agree with A&S 22.7.
>
> batch(test_orthopoly) reports no errors, nor does
> batch(rtest_pochhammer, test).
>
> Thanks to Kris Katterjohn for suggesting most of these changes.
>
> best,
>
> Robert Dodier
>
> PS.
> commit e5a8abaa23caf38de6ed397f8d6a79f5be46ae65
> Author: robert_dodier <robert_dodier at users.sourceforge.net>
> Date:   Wed Aug 22 12:15:05 2012 -0600
>
>      Express recurrences for legendre_p and legendre_q in terms of n + 1, n, and n - 1,
>      as other recurrences are expressed in orthopoly, and as shown in A&S 22.7. Thanks to Kris Katterjohn:
>      http://sourceforge.net/tracker/?func=detail&aid=3553317&group_id=4933&atid=304933
>
> diff --git a/doc/info/orthopoly.texi b/doc/info/orthopoly.texi
> index 12911e5..61cca4a 100644
> --- a/doc/info/orthopoly.texi
> +++ b/doc/info/orthopoly.texi
> @@ -824,10 +824,10 @@ respect to the polynomial degree.
>   @c ===end===
>   @example
>   (%i1) orthopoly_recur (legendre_p, [n, x]);
> -                (2 n - 1) P     (x) x + (1 - n) P     (x)
> -                           n - 1                 n - 2
> -(%o1)   P (x) = -----------------------------------------
> -         n                          n
> +                    (2 n + 1) P (x) x - n P     (x)
> +                               n           n - 1
> +(%o1)   P     (x) = -------------------------------
> +         n + 1                   n + 1
>   @end example
>   
>   The second argument to @code{orthopoly_recur} must be a list with the
> diff --git a/share/orthopoly/orthopoly.lisp b/share/orthopoly/orthopoly.lisp
> index d7c643a..9a851b3 100644
> --- a/share/orthopoly/orthopoly.lisp
> +++ b/share/orthopoly/orthopoly.lisp
> @@ -1260,6 +1260,8 @@ Maxima code for evaluating orthogonal polynomials listed in Chapter 22 of Abramo
>       
>   ;; For recursion relations, see A & S 22.7 page 782.
>   
> +;; legendre_p(n+1,x) = ((2*n+1)*legendre_p(n,x)*x-n*legendre_p(n-1,x))/(n+1)
> +
>   ;;  jacobi_p(n+1,a,b,x) = (((2*n+a+b+1)*(a^2-b^2) +
>   ;;    x*pochhammer(2*n+a+b,3)) * jacobi_p(n,a,b,x) -
>   ;;    2*(n+a)*(n+b)*(2*n+a+b+2)*jacobi_p(n-1,a,b,x))/(2*(n+1)*(n+a+b+1)*(2*n+a+b))
> @@ -1350,16 +1352,16 @@ Maxima code for evaluating orthogonal polynomials listed in Chapter 22 of Abramo
>   	 (let* ((n (nth 1 arg))
>   	       (x (nth 2 arg))
>   	       (z (if (eq fn '$legendre_q)
> -		      `((mtimes) -1 ((%kron_delta) ,n 1)) 0)))
> +		      `((mtimes) -1 ((%kron_delta) ,n 0)) 0)))
>   	   (simplify
> -	     `((mequal) ((,fn) ,n ,x)
> +	     `((mequal) ((,fn) ((mplus) 1 ,n) ,x)
>   	       ((mplus)
> -		((mtimes) ((mexpt) ,n -1)
> +		((mtimes) ((mexpt) ((mplus) 1 ,n) -1)
>   		 ((mplus)
> -		  ((mtimes) ((mplus) 1 ((mtimes) -1 ,n))
> -		   ((,fn) ((mplus) -2 ,n) ,x))
> -		  ((mtimes) ((mplus) -1 ((mtimes) 2 ,n))
> -		   ((,fn) ((mplus) -1 ,n) ,x) ,x)))
> +		  ((mtimes) ((mtimes) -1 ,n)
> +		   ((,fn) ((mplus) -1 ,n) ,x))
> +		  ((mtimes) ((mplus) 1 ((mtimes) 2 ,n))
> +		   ((,fn) ,n ,x) ,x)))
>   		,z)))))
>   
>   	((member fn '($assoc_legendre_p $assoc_legendre_q) :test 'eq)
> commit ff248ac1ee83a69eea03a762e365f1635bea30c8
> Author: robert_dodier <robert_dodier at users.sourceforge.net>
> Date:   Wed Aug 22 11:42:00 2012 -0600
>
>      Change recurrence relations for ultraspherical and Hermite polynomials to
>      agree with output of orthopoly_recur and A&S 22.7.
>
> diff --git a/share/orthopoly/orthopoly.lisp b/share/orthopoly/orthopoly.lisp
> index 7c63aea..d7c643a 100644
> --- a/share/orthopoly/orthopoly.lisp
> +++ b/share/orthopoly/orthopoly.lisp
> @@ -1265,7 +1265,7 @@ Maxima code for evaluating orthogonal polynomials listed in Chapter 22 of Abramo
>   ;;    2*(n+a)*(n+b)*(2*n+a+b+2)*jacobi_p(n-1,a,b,x))/(2*(n+1)*(n+a+b+1)*(2*n+a+b))
>   
>   ;; ultraspherical(n+1,a,x) = (2*(n+a)*x * ultraspherical(n,a,x) -
> -;;    (n+2*a+1)*ultraspherical(n-1,a,x))/(n+1)
> +;;    (n+2*a-1)*ultraspherical(n-1,a,x))/(n+1)
>   
>   ;; chebyshev_t(n+1,x) = 2*x*chebyshev_t(n,x) -chebyshev_t(n-1,x)
>   
> @@ -1276,7 +1276,7 @@ Maxima code for evaluating orthogonal polynomials listed in Chapter 22 of Abramo
>   ;; gen_laguerre(n+1,a,x) = (((2*n+a+1) - x)*gen_laguerre(n,a,x)
>   ;;  -(n+a)*gen_laguerre(n-1,a,x))/(n+1)
>   
> -;; hermite(n+1,x) = 2*x*hermite(n,x) +2*n*hermite(n-1,x)
> +;; hermite(n+1,x) = 2*x*hermite(n,x) -2*n*hermite(n-1,x)
>   
>   ;; See G & R 8.733.2; A & S 22.7.11 might be wrong -- or maybe I need
>   ;; reading glasses.
> commit 4743ec4737e9dba1b9cbededdf86edc8fa18c2be
> Author: robert_dodier <robert_dodier at users.sourceforge.net>
> Date:   Wed Aug 22 11:18:50 2012 -0600
>
>      Apply patch for assoc_legendre_p in orthopoly_weight. Thanks to Kris Katterjohn:
>      http://sourceforge.net/tracker/?func=detail&aid=3553910&group_id=4933&atid=304933
>
> diff --git a/share/orthopoly/orthopoly.lisp b/share/orthopoly/orthopoly.lisp
> index 84607e6..7c63aea 100644
> --- a/share/orthopoly/orthopoly.lisp
> +++ b/share/orthopoly/orthopoly.lisp
> @@ -1484,6 +1484,13 @@ variable ~:M" arg (car (last arg))))
>   	 (check-arg-length fn 2 (- (length arg) 1))
>   	 `((mlist) 1 -1 1))
>   
> +	; This is for a fixed order.  There is also an orthogonality
> +	; condition for fixed degree with weight function 1/(1-x^2).
> +	; See A & S 8.14.11 and 8.14.12.
> +	((eq fn '$assoc_legendre_p)
> +	 (check-arg-length fn 3 (- (length arg) 1))
> +	 `((mlist) 1 -1 1))
> +
>   	((eq fn '$laguerre)
>   	 (check-arg-length fn 2 (- (length arg) 1))
>   	 (let ((x (nth 2 arg)))
> commit 8bd381f1261024d64a3f4a84906a21699cbd3df0
> Author: robert_dodier <robert_dodier at users.sourceforge.net>
> Date:   Sat Aug 18 15:01:56 2012 -0600
>
>      In package orthopoly, return expressions containing %KRON_DELTA instead of $KRON_DELTA
>      since simplification property is attached to %KRON_DELTA. With this change, batch(test_orthopoly)
>      reports no errors.
>
> diff --git a/share/orthopoly/orthopoly.lisp b/share/orthopoly/orthopoly.lisp
> index a2ff70b..84607e6 100644
> --- a/share/orthopoly/orthopoly.lisp
> +++ b/share/orthopoly/orthopoly.lisp
> @@ -581,7 +581,7 @@ Maxima code for evaluating orthogonal polynomials listed in Chapter 22 of Abramo
>   	 '((n x)
>   	   ((unk) "$first" "$legendre_p")
>   	   ((mplus)
> -	    ((mtimes) -1 (($kron_delta) 0 n)
> +	    ((mtimes) -1 ((%kron_delta) 0 n)
>   	     ((mexpt) ((mplus) -1 ((mexpt) x 2)) -1))
>   	    ((mtimes)
>   	     ((mplus)
> @@ -997,7 +997,7 @@ Maxima code for evaluating orthogonal polynomials listed in Chapter 22 of Abramo
>   (defun $spherical_bessel_j (n x)
>     (cond ((and (eq '$zero (csign ($ratdisrep x)))
>   	      (or (integerp n) ($featurep n '$integer)))
> -	 `(($kron_delta) 0 ,n))
> +	 `((%kron_delta) 0 ,n))
>   
>   	((and (use-float x) (integerp n))
>   	 (let ((d 1) (xr) (xi) (z))
> @@ -1350,7 +1350,7 @@ Maxima code for evaluating orthogonal polynomials listed in Chapter 22 of Abramo
>   	 (let* ((n (nth 1 arg))
>   	       (x (nth 2 arg))
>   	       (z (if (eq fn '$legendre_q)
> -		      `((mtimes) -1 (($kron_delta) ,n 1)) 0)))
> +		      `((mtimes) -1 ((%kron_delta) ,n 1)) 0)))
>   	   (simplify
>   	     `((mequal) ((,fn) ,n ,x)
>   	       ((mplus)
> commit a82e9cbc10664d48a0eb012f4feb0e9029fa561a
> Author: robert_dodier <robert_dodier at users.sourceforge.net>
> Date:   Sat Aug 18 14:01:59 2012 -0600
>
>      Correct implementation of Lisp function kron_delta.
>      (1) In (DEFUN $KRON_DELTA ...), allow for multiple arguments
>      (2)  ... and punt directly to simplifier (no need for detour through TAKE)
>      (3)  ... and simplify with operator %KRON_DELTA (as $KRON_DELTA has no simplification property).
>      (4) Fix misspelling (knon_delta --> kron_delta) in error message.
>
> diff --git a/src/nset.lisp b/src/nset.lisp
> index e6e273b..875152e 100644
> --- a/src/nset.lisp
> +++ b/src/nset.lisp
> @@ -863,7 +863,7 @@
>   (setf (get '%kron_delta 'verb) '$kron_delta)
>   (setf (get '$kron_delta 'alias) '%kron_delta)
>   (setf (get '%kron_delta 'reversealias) '$kron_delta)
> -(defun $kron_delta (x) (take '($kron_delta) x))
> +(defun $kron_delta (&rest x) (simplifya `((%kron_delta) , at x) t))
>   (setf (get '%kron_delta 'real-valued) t) ;; conjugate(kron_delta(xxx)) --> kron_delta(xxx)
>   (setf (get '%kron_delta 'integer-valued) t) ;; featurep(kron_delta(xxx), integer) --> true
>   
> @@ -873,7 +873,7 @@
>     (declare (ignore y))
>   
>     (setq l (cdr l)) ;; remove (($kron_delta simp)
> -  (if (and l (null (cdr l))) (wna-err '$knon_delta)) ;; wrong number of arguments error for exactly one argument
> +  (if (and l (null (cdr l))) (wna-err '$kron_delta)) ;; wrong number of arguments error for exactly one argument
>   
>     ;; Checking both mnqp and meqp is convenient, but unnecessary. This code misses simplifications that
>     ;; involve three or more arguments. Example: kron_delta(a,b,a+b+1,a-b+5) could (but doesn't) simplify
>
> _______________________________________________
> Maxima mailing list
> Maxima at math.utexas.edu
> http://www.math.utexas.edu/mailman/listinfo/maxima
>