;;; Gaussian Quadrature, arbitrary precision.
;;; See quad-ga.lisp for basic ideas and comments
;;; This file is hacked to run faster and/or use less temp storage.
;;; Richard Fateman Jan 14, 2007.
(eval-when (compile load)
(declaim (optimize (speed 3) (safety 0) (space 0) (compilation-speed 0)))
(load "mpfr.fasl"))
(in-package :mpfr)
;; 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).
(defun legendre_pd (k x) ;; use for float or generic version
(assert (and (typep k 'fixnum)(<= 0 k)))
(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
(defun legendre_pdbf (k x)(legendre_pd k x))
(defun legendre_pdbf (k x)
(assert (and (typep k 'fixnum)(<= 0 k)))
(setf x (mpfr::into x)) ;make sure it's a bigfloat
(case k
(0 (values 1 0))
(1 (values x 1))
(otherwise
(let ((t0 (mpfr::into 1))
(t1 x); (mpfr::into x))
(ans (mpfr::into 0))
(i (mpfr::into 1)))
(declare (fixnum k))
(loop for ii from 2 to k do
(dsetv i (+ i 1)) ;; i is bigfloat version of ii
;; make a new copy for ans, so t1 and ans are different!
(setf ans (with-temps(/ (- (* (1- (* 2 i)) x t1) (* (1- i) t0)) i)))
(setf t0 t1)
(setf 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))))))))))
;; Starting from a point that approximates a root of a Legendre polynomial,
;; improve it.
(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)))))
;; the bigfloat version
;; 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)
(if (= deriv 0) (if (= val 0) val (error "newton iteration failed at ~s"guess))
;;new_guess= guess-f(x)/f'(x)
;;(- guess (/ val deriv)) ; example of open-coding mpfr calls.
(progn
(mpfr_div (gmpfr-f val)(gmpfr-f val)(gmpfr-f deriv) 0) ;;val:=val/deriv
(mpfr_sub (gmpfr-f guess)(gmpfr-f guess)(gmpfr-f val) 0)
;; this places the answer where the previous guess was.
;; we return it too, just for debugging.
guess))))
;;The weights are w[j]:= -2/ ( (n+1)* P'[n](x[j])*P[n+1](x))
;; this program computes 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 (mpfr::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 (- kf 1) (if (evenp k) 1 x)) 2)
(with-temps(/ (* (- kf 1) (- (* x t1)t0)) (- (* x x) 1))))
;; legendre_p[k](x)
(with-temps (/ (- (* (1- (* 2 kf)) x t1) (* (1- kf) t0)) kf))
)))))
#|
:pa :mpfr
(mpfr::set-prec 200) ; bits of fraction
(int_gs_l2bf #'exp 40)
;; this one is twice as fast:
(int_gs_l4bf #'exp 40)
;; running the integration formula repeatedly is very cheap, after weights etc computed.
;; compare to correct answer
(time (int_gs_l5bf #'exp 40))
(time (int_gs_l5bf #'exp 40)) ;; much faster 2nd time
;; the answer to this integration is computed by (right)
(defun right ()(let ((k (exp (mpfr::into 1)))) (- k (/ 1 k)))))
|#
;; To save allocation time and space for the array of weights etc, we do this:
;; 0. initialize s to 0. For i= 0 to n do lines 1-5:
;; 1. compute one abscissa x=x[i]
;; 2. compute one weight w=w[i]
;; 3. compute one value v=f(x[i])
;; 4. set s:= s+w*v.
;; 5. reusing x,w,v storage, repeat until n.
(defun int_gs_l2bf(fun n)
(let* ((np1 (+ n 1))
(v nil)
(fv (mpfr::into 0))
(myprec 38) ; might be 53, but maybe not all correct
(sum (mpfr::into 0))
(goal-precision (mpfr::get-prec))
(nph (/ 1.0d0(+ n 0.5d0)) ))
;; we step from point 1 to point n, summing along the way
(loop for j from 1 to n do
(setf v (cos(* pi (- j 0.25d0) nph))) ;initial approx. to a root of p[n]
(loop for k from 1 to 4 do ;; 4 iterations in double-float, from init refines it
(setf v (improve_lg_root n v )))
;; do this refinement some number of times more in bigfloat
(setf fv (mpfr::into v))
(setf myprec 38)
(loop while (< myprec goal-precision) do
(mpfr::set-prec (min goal-precision (setf myprec (* 2 myprec))))
(dsetv fv (improve_lg_rootbf n fv)))
;; now that fv is accurate enough put it into weighted sum
(dsetv sum (with-temps(+ sum (* (legendre_pdwt np1 fv) (funcall fun fv))))))
sum))
(defun ab_and_wts (n &optional (precision (get-prec))) ;mpfr precision. compute abscissae and weights.
(let ((a 0)
(v 0.0d0)
(np1 (+ n 1))
(nph (/ 1.0d0(+ n 0.5d0)))
(halfn (ash n -1))
(myprec 38)
(globalprec (get-prec))
(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 (- (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 (- v (/ val deriv)))))
(setf a (into v))
;; compute more accurate abscissa.
;; we build up gradually to goal precision.
(setf myprec 38)
(loop while (< myprec precision) do
(mpfr::set-prec (min globalprec (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
;; (format t "~%approx=~s" a)
(setf a (- a (/ val deriv))) ))
;; enough precision.
(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) ; fill in middle element by subtraction 2-2*(sum of others)
(let ((sum (into 0)))
(loop for i from 0 to (1- halfn) do (dsetv sum (+ sum (aref weights i ))))
(dsetv sum (- 2 (+ sum sum)))))))
(mpfr::set-prec globalprec)
(values abscissae weights)))
;; This is a version of quadrature computing weights etc as we go, like l2, but using symmetry.
(defun int_gs_l4bf(fun n)
(let* ((np1 (+ n 1))
(v nil)
(fv #.(mpfr::into 0))
(sum (into 0))
(myprec 38)
(nph (/ 1.0d0(+ n 0.5d0)) )
(goal-precision (mpfr::get-prec))
(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))
(setf myprec 38)
(loop while (< myprec goal-precision) do
(mpfr::set-prec (min goal-precision (setf myprec (* 2 myprec))))
(dsetv fv (improve_lg_rootbf n fv)))
;; now that v is accurate enough put it into weighted sum
(dsetv sum (with-temps(+ sum (* (legendre_pdwt np1 fv)
(+ (funcall fun fv)
(funcall fun (- fv))))))) )
(if (oddp n)
(+ sum (* (funcall fun 0)(legendre_pdwt np1 0)))
sum)))
;;;; memoizing stuff...
;; here's where we store precomputed wts and abscissae
(defparameter *legendabs* (make-hash-table))
(defparameter *legendwts* (make-hash-table))
(defparameter *legendprec* (make-hash-table))
;; call (clear-leg-has) to deallocate storage from hash tables for legendre polys
;; in case you think you need it for something else.
(defun clear-leg-hash()
(clrhash *legendabs*)
(clrhash *legendwts*)
(clrhash *legendprec*))
(defun int_gs_l5bf(fun n)
;; memoize the abscissae and weights
(let ((prec (or (gethash n *legendprec*) 0))
;;did we compute these weights previously
(weights nil)
(abscissae nil))
(cond ((cl::>= prec (get-prec))
;; Is enough precision already computed for this number of terms?
(setf abscissae (gethash n *legendabs*));use previous version
(setf weights (gethash n *legendwts*)))
(t
;; not found or insufficient precision, so compute and remember them now.
(multiple-value-setq (abscissae weights)(ab_and_wts n (get-prec)))
(setf (gethash n *legendabs*) abscissae)
(setf (gethash n *legendwts*) weights)
(setf (gethash n *legendprec*) (get-prec))))
;; 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)))
(defun gaussunit(fun n) (int_gs_l5bf fun n)) ;; integrate from -1 to 1
;;; alternative calling functions, based on versions in macsyma
;;; see quad-maxima.lisp comments.
;;;gaussab(%gg,lo,hi,n):=
;;; block([a:(hi-lo)/2, b:(hi+lo)/2],
;;; a* gq1(lambda([x],%gg(a*x+b)),n,ab_and_wts[n,fpprec])),
(defun gaussab(fun lo hi n) ;; integrate from lo to hi, using n points.
(let ((a (into (/ (- hi lo) 2)))
(b (into (/ (+ hi lo) 2)) ))
(* a (gaussunit #'(lambda(x)(funcall fun (with-temps (+ b (* a (into x)))))) n))))
;;;gauss0inf(%gg,n):=
;;; 2* gq1(lambda([t],block([d:(1-t)],
;;; %gg((1+t)/d)/d^2)),n,ab_and_wts[n,fpprec]),
(defun gauss0inf(fun n) ;; integrate from 0 to infinity
(* 2 (gaussunit
#'(lambda(z)
(let ((d (into (- 1 z))))
(/ (funcall fun
(with-temps
(/ (+ 1 z) d)))(* d d))))
n)))
;;;gaussminfinf(%gg,n):=
;;; 2* gq1(lambda([t],block([d: (1-t),
;;; r: (1+t)/(1-t)],
;;; (%gg(r)+%gg(-r))/d^2)),n,ab_and_wts[n,fpprec])
(defun gaussminfinf(fun n) ;; integrate from -infinity to infinity
(* 2 (gaussunit
#'(lambda(z)
(let* ((d (into (- 1 z)))
(r (with-temps (/ (+ 1 z) d ))))
(/ (+ (funcall fun r)
(funcall fun (- r))) (* d d))))
n)))
;; we could try playing with this:
;; accurate summation in a multiple-precision setting...
#+ignore
(defun acsum(a)
;;accurate summation of a[0]+ ...+a[last], all mpfr numbers
;; in an array a. Assume we are free to sort the array.
(sort a #'> :key #'abs)
(let ((p(set-prec (+ (get-prec) 48))))
(prog1 (reduce #'ga::two-arg-+ a)
(set-prec p))))
(defun acsum(a)
;;accurate summation of a[0]+ ...+a[last], all mpfr numbers
;; in an array a. Assume we are free to sort the array.
;;(sort a #'> :key #'abs)
(let ((p(set-prec (cl::+ (get-prec) 48)))); hack: just boost precision a bunch
(prog1 (reduce #'ga::two-arg-+ a)
(set-prec p))))
;;(setf r (vector (into 1)(into (expt 10 1000)) (into -1)(into (- (expt 10 1000)))))
;;; slower, more memory, acsum with sort has bug, different but maybe not much better when it counts.
(defun int_gs_l5bfa(fun n) ;like l5 but with accurate summation, stores terms
;; memoize the abscissae and weights
(let ((prec (or (gethash n *legendprec*) 0))
;;did we compute these weights previously
(weights nil)
(abscissae nil))
(cond ((cl::>= prec (get-prec))
;; Is enough precision already computed for this number of terms?
(setf abscissae (gethash n *legendabs*));use previous version
(setf weights (gethash n *legendwts*)))
(t
;; not found or insufficient precision, so compute and remember them now.
(multiple-value-setq (abscissae weights)(ab_and_wts n (get-prec)))
(setf (gethash n *legendabs*) abscissae)
(setf (gethash n *legendwts*) weights)
(setf (gethash n *legendprec*) (get-prec))))
;; compute the sum.
(let ((halfn (ash n -1))
(summands nil))
(loop for i from 0 to (1- halfn) do
(push (* (aref weights i)
(funcall fun (aref abscissae i))) summands)
(push (* (aref weights i)
(funcall fun (with-temps (-(aref abscissae i))))) summands))
(if (oddp n)
(push (* (aref weights halfn)(funcall fun (into 0))) summands))
(acsum summands)
)))
#| some testing in Mathematica regarding Legendre zeros
;; One zero of LegendreP[10,x] is about
;;0.973906528517171720077964012084452053428269946692382119231212066696595203234636159625723564956268556258233042518774211215022168601434477779920540958726
;; correct to the last rounded digit.
This program give about
;;0.97390652851717172007796401208425031143498830232617797549445377752271945185287852121950994663924446112734455480619806549061
;;................................^wrong
|#