(in-package :ga)
;;look in quad-ga for documentation

(defun legendre_pd (k x)
  (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)))))))))

;; Here is a Newton iteration to converge to a root of
;; legendre_p[n](x) starting from a particular guess.
;; The answer's precision will depend on the generic arithmetic
;; in  (- guess (/ val deriv))

(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 weights are w[j]:=  -2/ ( (n+1)* P'[n](x[j])*P[n+1](x))
;; this program computes them. 
(defun legendre_pdwt (k x)		; compute -2 / (  k * P'[k-1](x)*P[k](x))
  (assert (and (typep k 'fixnum)(< 1 k)))	; k>=2
  (/ -2 
     (* k
	(let 
	    ((t0 1)
	     (t1 x)
	     (ans 0) )
	  (declare (fixnum k))
	  (loop for i from 2 to (1- k)  do
		(setf ans     (/ (- (* (1- (* 2 i)) x t1) (* (1- i) t0)) i)
		      t0 t1
		      t1 ans))
	  (*				;deriv of p[k-1]
	   (let((k (1- k))) ;use formula above..
	     (if (= 1 (abs x))
		 (/ (* k (1+ k) (if (oddp k) 1 x)) 2)
	       (/ (* k (- (* x t1)t0)) (1- (* x x)))))
	   ;; legendre_p[k](x)
	   (/ (- (* (1- (* 2 k)) x t1) (* (1- k) t0)) k))))     ))


;; put almost all the pieces together:
;; Gaussian integration of the function fun using
;; n legendre points between -1 and 1, using whatever
;; precision is default for generic arithmetic.
(defun int_gs_l1(fun n) 

  (let* ((abscissae (improvedroots (setroots n) n))   ;; the points to evaluate fun at.
	 (np1 (+ n 1))
	 (weights (map 'array  #'(lambda(r)(legendre_pdwt np1 r))
		       abscissae)))
    ;; compute the sum.

    (loop for i from 0 to (1- n) sum
	  (* (funcall fun (aref abscissae i))(aref weights i)))
    
    ))



(defun int_gs_l1bf(fun n)  ;; bigfloats

  (let* ((np1 (+ n 1))
	 (v nil)
	 (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.
	  (loop for k from 1 to 4 do ;; 4 iterations in double-float, from init
		(setf v (improve_lg_root n v )))
	  ;;  do this refinement some number of times..
	  (setf v (mpfr::into v))
	  (setf myprec 38) 
	  (loop while (< myprec goal-precision) do
		(mpfr::set-prec (min goal-precision (setf myprec (* 2 myprec))))
		(setf v (improve_lg_root n v)))
	  ;; now that v is accurate enough put it into weighted sum
	  (setf sum (+ sum (* (legendre_pdwt np1 v) (funcall fun v)))))
    sum))


#|
:pa :ga
(mfpr::set-prec 200) ; bits of fraction
(int_gs_l1bf #'exp 10)			; 10 points
(int_gs_l1bf #'exp 20)			; 20 points
(int_gs_l1bf #'exp 40)

;; compare to correct answer 
(let ((k (exp (mpfr::into 1)))) (- k (/ 1 k)))
|#