;;; Gaussian Quadrature, arbitrary precision.
;;; See http://www.cs.berkeley.edu/~fateman/generic/quad-ga.lisp or quad-fast.lisp 
;;; for basic ideas and comments.
;;; hacked to run in Maxima or Macsyma, but consequently slower than using mpfr

(eval-when (load compile eval)
  (if (string= "MAXIMA" (package-name *package*))
      (in-package :maxima) ;; for sourceforge maxima
    (in-package :climax))) ;; for commercial macsyma

;; I suspect that binding fpprec is not the right thing to do in commercial Macsyma..

;;;  Richard Fateman  Feb 18, 2007

;; legendre_pd returns legendre_p(k,x) and its derivative.
;;Note that d/dx P[n](x) =  (1/(x^2-1))*n*(x*P[n](x)-P[n-1](x)) except for x=+-1
;; This program is used with Newton iteration to refine roots of P[n](x).
;; useful only from lisp because it uses value construct to return.

(defun legendre_pd (k x) ;; use for double float 
  (assert (and (typep k 'fixnum)(<= 0 k)
	       (typep x 'double-float)))
  (case k 
    (0 (values 1 0))
    (1 (values x 1))
    (otherwise
     (let ((t0 1)
	   (t1 x)
	   (ans 0))
       (declare (fixnum k))
       (loop for i from 2 to k  do
	     (setf ans     (/ (- (* (1- (* 2 i)) x t1) (* (1- i) t0)) i)
		   t0 t1
		   t1 ans))
       (values t1
		      ;; (1/(x^2-1))*k*(x*P[k](x)-P[k-1](x))
		      ;; except if abs(x)=1 then use
		      (if (= 1 (abs x))
			  ;; 1/2*k*(k+1)*x^(k+1)
			  (/ (* k (1+ k) (if (oddp k) 1 x)) 2)
			(/ (* k (- (* x t1)t0))  (1-  (* x x)))
			
))))))

;;bigfloat version. Returns a macsyma list

(defun $legendre_pdbf (k x)
  (assert (and (typep k 'fixnum)(<= 0 k)))
  (setf x ($bfloat x))		;make sure it's a bigfloat
  (case k 
    (0 (list '(mlist) 1 0))
    (1 (list '(mlist) x 1))
    (otherwise
     (let ((t0 ($bfloat 1))
	   (t1 ($bfloat x))
	   (ans ($bfloat 0))
	   (i ($bfloat 1)))
       (declare (fixnum k))
       (loop for ii from 2 to k  do
	     (setf i (add i 1))
	     (setf ans 
		   (div (sub (mul (sub (mul 2 i) 1) x t1) (mul (sub i 1) t0)) i))
	     (setf t0 t1)
	     (setf t1 ans) )
       (list '(mlist)  t1
		;; (1/(x^2-1))*k*(x*P[k](x)-P[k-1](x))
		;; except if abs(x)=1 then use
		(if (equal (cdr bigfloatone) (fpabs (cdr x))) ;abs of the bf insides
		    ;; 1/2*k*(k+1)*x^(k+1)
			  (div (mul k (add 1 k) (if (oddp k) 1 x)) 2)
		  (div (mul k (sub (mul x t1)t0))  (sub  (mul x x) 1))))))))

;;The Gaussian quad weights are w[j]:=  -2/ ( (n+1)* P'[n](x[j])*P[n+1](x[j]))
;; this program computes w[k](x)

(defun $legendre_pdwt (k x &aux (kf ($bfloat k)))
  ;; compute -2 / (  k * P'[k-1](x)*P[k](x))
  (assert (and (typep k 'fixnum)(< 1 k))); k>=2
  (setf x ($bfloat x))
  (div -2 
     (mul kf
	(let 
	    ((t0 ($bfloat 1))
	     (t1 x)
	     (ans 0)
	     (i ($bfloat 0)))
	  (declare (fixnum k))
	  (loop for ii from 2 to (1- k)  do
		(setf i ($bfloat ii))
		(setf ans   
		      (div (sub (mul (sub (mul 2 i) 1) x t1) (mul (sub i 1) t0)) i)
		      t0 t1
		      t1 ans))
	  (mul				;deriv of p[k-1]
	   (div (mul (sub kf 1) (sub (mul x t1)t0)) (sub (mul x x) 1))
	   ;; legendre_p[k](x)
	    (div (sub (mul (sub (mul 2 kf) 1) x t1) (mul (sub kf 1) t0)) kf) )))))

(defun $ab_and_wts (n &optional(precision fpprec)) ;compute abscissae and weights upto bigfloat global precision
  (declare (special fpprec $float2bf ))
  (let ((a 0)
	(v 0.0d0)
	(np1 (+ n 1))
	(nph (/ 1.0d0(+ n 0.5d0)))
	(halfn (ash n -1))
	(myprec 38)
	(globalprec fpprec)
	($float2bf t) ; no warnings
	(abscissae  (make-array (ceiling n 2)))
	(weights   (make-array (ceiling n 2))))
    (declare (double-float v nph))
    (loop for i from 0 to (1- halfn) do
	  (setf v (cos(* pi (- (1+ i) 0.25) nph)));estimate of zero
	  (loop for i from 1 to 4  do	; 3 iterations is almost enough.
		(multiple-value-bind (val deriv)(legendre_pd n v)
		  (declare (double-float val deriv))
		  ;;new_guess= guess-f(x)/f'(x)
		  (setf v (- v (/ val deriv)))))
	  (setf  a ($bfloat v))
	  ;; compute more accurate abscissa in bigfloat
	  ;;  we build up gradually to goal precision.
	  (setf myprec 38)
	  (loop while (< myprec precision) do
		(setf fpprec  (min globalprec (setf myprec 
						       (* 2 myprec)  )))
		(let* ((valderiv($legendre_pdbf n a))
		       (val (cadr valderiv))
		       (deriv (caddr valderiv)))
		  ;;new_guess= guess-f(x)/f'(x)
		  (setf a  (sub a (div val deriv))) ))
	  ;; enough precision.
	  (setf (aref abscissae i) a)
	  (setf (aref weights i) ($legendre_pdwt np1 a)))
    (cond ((oddp n)
	   (setf (aref abscissae  halfn) 0)
	   (setf (aref weights  halfn)	; fill in middle element by subtraction 2-2*(sum of others)
		 (let ((sum  0))
		   (loop for i from 0 to (1- halfn)  do (setf sum (add sum (aref weights i ))))
		   (setf sum (sub 2 (add sum sum)))))))
    ;;(values abscissae weights) ;probably what we would prefer, for computing
    ;; to make it accessible to macsyma, return 2 lists.
    (list '(mlist) 
	  (cons '(mlist) (coerce abscissae 'list))
	  (cons '(mlist) (coerce weights 'list)))    
    ))

#|

A Maxima program

First define some function, e.g.  g(x):=exp(x).
Using gaussian integration to integrate g from -1 to 1 in the
current bigfloat precision, using n points, do gaussunit(g,n).

To see if your answer is correct, use n then 2n and 4n points, and/or
more and more digits. Compare answers. This will provide heuristic confirmation.


( gq1(%%g,n,aw):=
block([sum:0],
 map(lambda([%%a,%%w],sum:sum+%%w*(%%g(%%a)+%%g(-%%a))),aw[1],aw[2]),
 sum+ (if oddp(n) then -%%g(0)*aw[2][length(aw[1])] else 0 )),

 ab_and_wts[n,prec]:=   
/* memoize the abscissae and weights for Legendre zeros and Gauss quad */
  block([fpprec:prec],
         ab_and_wts(n)),

/* integrate function g(x) from x=-1 to x=1 using n points and fpprec */	 
gaussunit(%ggg,n):= gq1(%ggg,n,ab_and_wts[n,fpprec]),

/* integrate function g(x) from x=lo to x=hi using n points and fpprec */	 

gaussab(%hh,lo,hi,n):= 
  block([a:(hi-lo)/2, b:(hi+lo)/2], local(%zz),
     define (%zz(x),%hh(a*x+b)),
    a* gaussunit(%zz,n)),

/* integrate function g(x) from x=0 to x=inf using n points and fpprec */	 
gauss0inf(%gg,n):= 
 block([],local(%zz),
	       define(%zz(x),      
			 block([d:(1-x)], %gg((1+x)/d)/d^2)),
 2* gaussunit(%zz,n)),
/* integrate function g(x) from x=minf to x=inf using n points and fpprec */	 
gaussminfinf(%gg,n):= 
 block([],local(%zz),
      define(%zz(x),      
        block([d:(1-x),r:(1+x)/(1-x)], (%gg(r)+%gg(-r))/d^2)),
     2* gaussunit(%zz,n))

)$

			 |#