haskeem-0.6.10: powerseries.scm
; $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)))