;; for Maxima, evaluate a power series at a point, fast.
;; Restriction: The point substituted for the one and only (MAIN) VAR in the series.
;; There are integer exponents. Rrational numerical coefficients or floats.
;; sample call:
;; epsp(taylor(sin(x),x,0,5),  1/2)

(defun check-for-epsp (p)
   (and (consp p)		; it's not a symbol or number
	(eq (caar p) 'mrat)		; it's a taylor series or a rat (CRE)
	(eq (cadr p) 'ps)		; it's a taylor series and not a constant
	;;(null(cdr (third (car p))))	; we would like to check that the body 
					;has only one variable. This won't do it .
	; varlist of taylor(sin(x),x,0,4)  has 2 elements, {sin(x), x}
		))

(defun $epsp(p r)

  ;; should do checking here that we have a power series with numeric coeffs etc
  
  (if (not (check-for-epsp p))  p ;; Don't change it if it doesn't pass muster.
	   
    (let* ((m (nthcdr 4 p)) ;; extract the "polynomial-ish" part of the data structure.
	   (lim (car(cadddr p))) ;; the highest exponent to retain, as a pair.
	   (limnum (/ (car lim)(cdr lim)))
	  (ans 0))  ;; see file hayat.lisp for truncated taylor series background
    (if (and (consp r)(eq (caar r) 'rat))
	(setf r (/ (cadr r)(caddr r)))) ; convert maxima rational to CL rational
    (map nil #'(lambda(q)
		 (let* 
		     ((coeff (cdr q))
		      (pow (car q))
		      (power (/ (car pow)(cdr pow)))) ; is (cdr power) always 1 for us?
		   (if (> power limnum) nil  ;; beyond trunc level
		 (incf ans (* (expt r power)
			      (/ (car coeff)(cdr coeff)))))))
	 m)
    (if(or(floatp ans)(= 1 (denominator ans)))	   ans
	  (list '(rat)(numerator ans)(denominator ans)))))) ; maxima rational ugh


(defun $epspf(p r)			; r is a double-float


  (declare (double-float r)(optimize (speed 3)(safety 0)))
  (if (not (check-for-epsp p))  p ;; Don't change it if it doesn't pass muster.
    (let ((m (nthcdr 4 p))
	  (ans 0.0d0))
      (declare (double-float ans))
      (map nil #'(lambda(q)
		   (let 
		       ((coeff (cdr q))
			(power (car q)))
		     (incf ans (* (expt r  (car power) ;; assume denom 1 (/ (car power)(cdr power))
					)
				  (/ (car coeff)(cdr coeff))))))
	   m)
      ans)))

;;; compare to using subst and ev ... 
;;; e.g. z:taylor(1/(x+1),x,0,1000)
;;; epsp(z,0.5d0) takes 0.009 seconds
;;; epspf(z,0.5d0) takes something less (this is GCL)
;;;subst(0.5d0,x,z) takes 0.012 seconds  
;; ev(z,x=0.5d0)  takes about 0.009  also!

;;; A better way might be to use Horner's rule.
;;; We have to figure the gap between exponents, e.g. for sin(x) series...
;;; In the program below, we try a kind of horner's rule.

;;; epsph(z,0.5d0) takes 0.0008 seconds. That's 10X faster than our other
;;; program or ev.
;;; but it is not as general.


(defun $epsph(p r) ; r is a double-float
  ;; should do checking here that we have a power series with numeric coeffs etc
  ;; we could do better if we knew the power series coefs were floats too
  ;; and if we used a cleverer optimized lisp. This is GCL.
  (declare (double-float r)(optimize (speed 3)(safety 0)))
  
  
  (if (not (check-for-epsp p))  p ;; Don't change it if it doesn't pass muster.
    
    (let* ((lim (car(cadddr p))) ;; the highest exponent to retain, as a pair.
	   (limnum (/ (car lim)(cdr lim))))
  (epspf2 (nthcdr 4 p) (* 1.0d0 r) 1.0 0 limnum))))   ;; r^n, n

(defun epspf2 (m r oldr^n oldn limnum)
  (declare (double-float r oldr^n)(fixnum oldn) (optimize (speed 3)(safety 0)))
  (cond ((null m) 0.0d0)
	(t  (let* ((q (car m))	; the lead term
		   (coeff (cdr q))	; its coefficient
		   (n (caar q)) ;; its power, assuming denom of power is 1
		   (increm (- n oldn))
		   (r^n (* oldr^n (if (= increm 1 r) r  (expt r increm)))))
	      (declare (double-float r^n)(fixnum increm n))
	      (if (> n limnum) 0
		(+ (* (/ (car coeff)(cdr coeff))
		      r^n)
		   (epspf2 (cdr m) r r^n n limnum)))))))