; $Id: powerseries.scm,v 1.9 2008/03/08 03:52:39 uwe Exp $

; A small package for manipulating power series implemented as streams.
; There is nothing here that forces the use of floating-point, so if user
; inputs are integers or rationals, the results will also be rational; this
; is nice for getting out exact results.

; The stream of natural numbers (0 1 2 3 ...),
; useful as a building block for stuff below

(define ps-nat-numbers
  (letrec ((next
	    (lambda (n)
	      (cons n (delay (next (+ n 1)))))))
    (next 0)))

; Some functions for power-series specifically

; negation, addition, and subtraction are trivial

(define (ps-neg ps) (stream-scale -1 ps))
(define (ps-add ps . pss) (apply stream-map (append (list + ps) pss)))
(define (ps-sub ps . pss) (apply stream-map (append (list - ps) pss)))

; multiplication and reciprocal are a little less trivial

(define (ps-mul ps1 ps2)
  (let ((h1 (car ps1))
	(h2 (car ps2))
	(t1 (cdr ps1))
	(t2 (cdr ps2)))
    (cons (* h1 h2) (delay (ps-add (stream-scale h1 (force t2))
				   (stream-scale h2 (force t1))
				   (cons 0 (ps-mul (force t1) (force t2))))))))

; It is an error to take the reciprocal of a power series without a constant
; term because such a resulting series is not representable as a power
; series: it contains negative powers of x, and these cause infinities at x
; = 0. It would be a Laurent series, and we don't handle these.

(define (ps-recip ps)
  (when (zero? (car ps))
	(raise "ps-recip error: power series has no constant term!"))
  (letrec* ((r (cons 1 (delay (ps-neg (ps-mul r (fcdr ps)))))))
	   (stream-scale (/ (car ps)) r)))

; division is trivial when expressed in terms of reciprocal

(define (ps-div ps1 ps2) (ps-mul ps1 (ps-recip ps2)))

; Term-by-term differentiation and integration

(define (ps-derivative ps) (fcdr (stream-map * ps-nat-numbers ps)))

(define (ps-integral ps ctrm)
  (cons ctrm (stream-map / ps (fcdr ps-nat-numbers))))

; Semi-evaluate a power-series by multiplying each term by a specific x
; raised to the Nth power; the sum of all of these is the actual value of
; the expression, which may be approximated by the ps-sums routine below.

(define (ps-eval ps x)
  (stream-map * ps (stream-map (lambda (n) (expt x n)) ps-nat-numbers)))

; The stream of partial sums of the input

(define (ps-sums ps)
  (letrec ((accfn (lambda (acc str)
		    (set! acc (+ acc (car str)))
		    (cons acc (delay (accfn acc (fcdr str)))))))
    (accfn 0 ps)))

; TODO: Aitken acceleration formula (but is x_k here the kth term, or
; the kth partial sum??? check that!)
; new estimate = x_k - (x_{k+1} - x_k)^2/(x_{k+2} - 2*x_{k+1} + x_k)

; And some actual power series

; geometric series 1/(1 - x)

(define ps-geom (stream-map (lambda (n) 1) ps-nat-numbers))

(define ps-exp (stream-map (lambda (n) (/ (factorial n))) ps-nat-numbers))

; ln(1+x) -- valid only for -1 < x <= 1
; note that this converges *very* slowly for |x| near 1

(define ps-logxp1
  (cons 0 (stream-map
	   (lambda (n) (/ (expt -1 (+ n 1)) n)) (fcdr ps-nat-numbers))))

(define ps-sin (stream-map
		(lambda (n)
		  (if (even? n)
		      0
		      (/ (expt -1 (quotient (- n 1) 2))
			 (factorial n))))
		ps-nat-numbers))

(define ps-cos (stream-map
		(lambda (n)
		  (if (odd? n)
		      0
		      (/ (expt -1 (quotient n 2))
			 (factorial n))))
		ps-nat-numbers))

(define ps-tan (ps-div ps-sin ps-cos))

(define ps-atan (stream-map
		 (lambda (n)
		   (if (even? n)
		       0
		       (/ (expt -1 (quotient (- n 1) 2)) n)))
		 ps-nat-numbers))

(define ps-sinh (stream-map
		 (lambda (n)
		   (if (even? n)
		       0
		       (/ (factorial n)))) ps-nat-numbers))

(define ps-cosh (stream-map
		 (lambda (n)
		   (if (odd? n)
		       0
		       (/ (factorial n)))) ps-nat-numbers))

(define ps-tanh (ps-div ps-sinh ps-cosh))

(define ps-atanh (stream-map
		 (lambda (n)
		   (if (even? n)
		       0
		       (/ n))) ps-nat-numbers))

; exp(-x^2)

(define ps-gaussian (stream-map
		     (lambda (n)
		       (if (odd? n)
			   0
			   (begin (set! n (quotient n 2))
				  (/ (expt -1 n)
				     (factorial n))))) ps-nat-numbers))

; This is actually not quite Erf(x): there is a scale factor of 1/sqrt(pi)
; missing. I'm leaving that out so that stuff doesn't get forced to
; floating-point, since I don't have an infinite-precision rational version
; of (sqrt pi).

(define ps-erf (stream-scale 2 (ps-integral ps-gaussian 0)))

; Bessel functions J_n for integer n >= 0:
;
;                                (-1)^m x^(2*m)
; J_n = x^n * sum_m=0^infinity -------------------
;                              2^(2*m+n)*m!*(m+n)!

(define (ps-bessel-j n)
  (when (negative? n)
	(raise "ps-bessel-j can't handle negative n!"))
  (letrec ((term-fn (lambda (m2)
		      (if (odd? m2)
			  0
			  (let ((m (quotient m2 2)))
			    (/ (expt -1 m)
			       (* (expt 2 (+ m2 n))
				  (factorial m)
				  (factorial (+ m n))))))))
	   (cons-fn (lambda (m obj)
		      (if (zero? m)
			  obj
			  (cons 0 (cons-fn (- m 1) obj))))))
    (cons-fn n (stream-map term-fn ps-nat-numbers))))

; Bessel functions I_n for integer n >= 0:
;
;                                    x^(2*m)
; I_n = x^n * sum_m=0^infinity -------------------
;                              2^(2*m+n)*m!*(m+n)!

(define (ps-bessel-i n)
  (when (negative? n)
	(raise "ps-bessel-i can't handle negative n!"))
  (letrec ((term-fn (lambda (m2)
		      (if (odd? m2)
			  0
			  (let ((m (quotient m2 2)))
			    (/ (* (expt 2 (+ m2 n))
				  (factorial m)
				  (factorial (+ m n))))))))
	   (cons-fn (lambda (m obj)
		      (if (zero? m)
			  obj
			  (cons 0 (cons-fn (- m 1) obj))))))
    (cons-fn n (stream-map term-fn ps-nat-numbers))))

; Lambert W function: this satisfies the implicit equation W*exp(W) = x.
; This series converges for |x| < exp(-1). The function is multi-valued
; for -exp(-1) < x < 0; this series converges to the value closest to 0.

(define ps-lambert-w
  (cons 0 (stream-map
	   (lambda (n) (/ (expt (- n) (- n 1)) (factorial n)))
	   (fcdr ps-nat-numbers))))

; A couple of small utility functions to more easily show streams and
; tabulate values

(define (ps-show n p strm)
  (if (zero? p)
      (for-each (lambda (val)
		  (write-string (number->string val 10) #\linefeed))
		(stream-head strm n))
      (for-each (lambda (val)
		  (write-string (number->string val 10 p) #\linefeed))
		(stream-head strm n))))

(define (ps-table fn nterms lo hi step)
  (do ((x lo (+ x step)))
      ((> x hi) #t)
    (write-string (number->string x 10 8) #\tab
		  (number->string (last (stream-head
					 (ps-sums (ps-eval fn x)) nterms))
				  10 -8) #\linefeed)))