;;; Gaussian Quadrature, quad-double precision
;;; See quad-ga.lisp for basic ideas and comments
;;; quad-oqd.lisp for older versions.
;;; Jan 10, 2007.
;;; by Richard Fateman
;; 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).
(in-package :qd)
(eval-when (compile load)
(declaim (optimize (speed 3) (safety 0) (space 0) (compilation-speed 0)))
(load "qd.dll")
)
(defun legendre_pd (k x)
(declare(double-float x) (fixnum k) (optimize (speed 3)(safety 0)))
(assert (and (typep k 'fixnum)(<= 0 k)))
(case k
(0 (values 1.0d0 0.0d0))
(1 (values x 1.0d0))
(otherwise
(let ((t0 1.0d0)
(t1 x)
(ans 0.0d0))
(declare(double-float t0 t1 ans x))
(declare (fixnum k))
(loop for i from 2 to k do
(setf ans (cl::/ (cl::- (cl::* (cl::1- (* 2.0d0 i)) x t1) (cl::* (cl::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)
(cl::/ (cl::* k (cl::- (cl::* x t1)t0)) (cl::1- (cl::* x x)))))))))
;;qdfloat version
(defun legendre_pdbf (k x)
(assert (and (typep k 'fixnum)(<= 0 k)))
(case k
(0 (values 1 0))
(1 (values (into x) 1))
(otherwise
(let ((t0 (into 1))
(t1 (into x))
(ans (into 0))
(i (into 1)))
(declare (fixnum k))
(loop for ii from 2 to k do
(dsetv i (+ i 1))
(setf ans (with-temps(/(- (* (- (* 2 i)1) x t1) (* (- i 1) t0)) i)))
(dsetv t0 t1)
(dsetv 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)
(with-temps (/ (* (into k) (- (* x t1)t0)) (1- (* x x))))
))))))
;; the float version
(defun improve_lg_root(n guess)
(multiple-value-bind (val deriv)(legendre_pd n guess)
(if (= deriv 0) (if (= val 0)
val
(error "newton iteration failed at ~s"guess))
;;new_guess= guess-f(x)/f'(x)
(- guess (/ val deriv)))))
;; Here is a Newton iteration to converge to a root of legendre_p, bigfloat.
(defun improve_lg_rootbf(n guess)
(multiple-value-bind (val deriv)(legendre_pdbf n guess)
;;new_guess= guess-f(x)/f'(x)
;; dsetv alters array in place
(dsetv guess (- guess (/ val deriv)))) )
;;The weights are w[j]:= -2/ ( (n+1)* P'[n](x[j])*P[n+1](x))
;; this program computes one of them
(defun legendre_pdwt (k x &aux (kf (into k))) ; compute -2 / ( k * P'[k-1](x)*P[k](x))
(assert (and (typep k 'fixnum)(< 1 k))) ; k>=2
(setf x (into x))
(/ -2
(* kf
(let
((t0 (into 1))
(t1 x)
(ans 0)
(i (into 0)))
(declare (fixnum k))
(loop for ii from 2 to (1- k) do
(setf i (into ii))
(setf ans
;; make a new number here
(/ (with-temps (- (* (1- (* 2 i)) x t1) (* (1- i) t0))) i)
t0 t1
t1 ans))
(* ;deriv of p[k-1]
(if (= 1 (abs x)) ;not going to happen..
(/ (* kf (1- kf) (if (evenp k) 1 x)) 2)
(with-temps(/ (* (- kf 1) (- (* x t1)t0)) (- (* x x) 1))))
;; legendre_p[k](x)
(with-temps (/ (- (* (- (* 2 kf) 1) x t1) (* (- kf 1) t0)) kf))))) ))
#|
:pa :qd
(int_gs_l4bf #'exp 40)
;; the next 2 lines are about as expensive as the previous one
(multiple-value-setq (*abs* *vals*) (ab_and_wts 40)) ;; compute constants
(int_gs_l3bf #'exp *abs* *vals*) ;; run the integration formula
;; but running the integration formula repeatedly is very cheap, after weights etc computed.
;; compare to correct answer
(int_gs_l5bf #'exp 40) ; combines both ideas: remembers the weights, once they are computed.
|#
(defun right ()(let ((k (exp (into 1)))) (- k (/ 1 k))))
;; Yet another approach, assuming you can pre-compute arrays of
;; abscissae and weights. Then you can compute the function values and
;; add everything together.
(defun int_gs_l3bf(fun abscissae weights) ;; bigfloats
;; compute the sum.
(let ((sum (into 0))
(n (1-(length abscissae))))
(loop for i from 0 to n do (setf sum
(+ sum(*(funcall fun (aref abscissae i))
(aref weights i)))))
sum))
#| ;; version 4
;; to save allocation time and space we
;; 0. initialize s to 0. For i= 0 to n/2 do lines 1-5: ;; assume n is even.
;; 1. compute one abscissa x=x[i]
;; 2. compute one weight w=w[i]
;; 3. compute two values v=f(x)+f(-x))
;; 4. set s:= s+w*v.
;; 5. reusing x,w,v storage, repeat until n.
;; depending on whether n was even or odd, there may be one more point, at f(0)
|#
(defun int_gs_l4bf(fun n)
(let* ((np1 (+ n 1))
(v nil)
(fv #.(qd::into 0))
(sum (into 0))
(nph (/ 1.0d0(+ n 0.5d0)) )
(halfn (ash n -1)))
;; we step from point 1 to point n, summing along the way
(loop for j from 1 to halfn do
(setf v (cos(* pi (- j 0.25d0) nph))) ;initial approx.
(loop for k from 1 to 4 do ;; 4 iterations in double-float, from init
(setf v (improve_lg_root n v )))
(setf fv (into v))
;; 4 iterations in qd
(loop for j from 1 to 4 do
(dsetv fv (improve_lg_rootbf n fv)))
;; now that v is accurate enough put it into weighted sum
(dsetv sum (+ sum (* (legendre_pdwt np1 fv)
(+ (funcall fun fv)
(funcall fun (- fv))))))
)
(if (oddp n)
(+ sum (* (funcall fun 0)(legendre_pdwt np1 0)))
sum)))
;; yet other posible improvements. Only compute half the data and
;; generate the symmetric parts when needed. Also, for the middle
;; weight of an odd number of terms, realize that the sum of all the
;; weights will be 2. So if we compute S, the sum weights up to n/2,
;; the remaining weight will be 2*(1 - S).
;; Also we can perhaps skip the full weight computation n is even,
;; since one can be deduced by subtraction.
;; look for precomputed guys, else compute new ones
(defparameter *legendabs* (make-hash-table))
(defparameter *legendwts* (make-hash-table))
;; double-up on weights and abscissae
(defun int_gs_l5bf(fun n) ; memoize the abscissae and weights
(let ((abscissae (gethash n *legendabs*))
(weights nil))
(cond (abscissae (setf weights (gethash n *legendwts* )))
(t;; compute them now.
(multiple-value-setq (abscissae weights)(ab_and_wts n))
(setf (gethash n *legendabs*) abscissae)
(setf (gethash n *legendwts*) weights)))
;; compute the sum.
(let ((sum (into 0))
(halfn (ash n -1)))
(loop for i from 0 to (1- halfn) do (dsetv sum
(+ sum (*(+ (funcall fun (aref abscissae i))
(funcall fun (with-temps (-(aref abscissae i)))))
(aref weights i)))))
(if (oddp n)
(dsetv sum (+ sum (* (aref weights halfn)(funcall fun (into 0))))))
sum)))
;; pre-compute abscissae and weights
(defun ab_and_wts (n );precision is somewhat less than 215, qd. Stores only half the values (maybe +1)
(let ((a 0)
(v 0.0d0)
(np1 (+ n 1))
(nph (/ 1.0d0(+ n 0.5d0)))
(halfn (ash n -1))
(myprec 38)
(abscissae (make-array (ceiling n 2)))
(weights (make-array (ceiling n 2))))
(loop for i from 0 to (1- halfn) do
(setf v (cos(* pi (- (cl::1+ i) 0.25) nph)))
(loop for i from 1 to 4 do ; 3 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 (cl::- v (cl::/ val deriv)))))
(setf a (into v))
;; compute more accurate abscissa.
;; we build up gradually to goal precision.
(setf myprec 38)
;;takes a while
(loop while (< myprec 208) do
(setf myprec (* 2 myprec) )
(multiple-value-bind (val deriv)(legendre_pdbf n a)
;;new_guess= guess-f(x)/f'(x)
;; dsetv alters a in place
(dsetv a (- a (/ val deriv))) ))
(setf (aref abscissae i) a)
(setf (aref weights i) (legendre_pdwt np1 a)) )
(cond ((oddp n)
(setf (aref abscissae halfn) (into 0))
(setf (aref weights halfn)
(legendre_pdwt np1 (aref abscissae halfn)))))
(values abscissae weights)))
(defun clear-leg-hash() ; call to deallocate storage from hash tables for legendre polys
(clrhash *legendabs*)
(clrhash *legendwts*))