Formal power series: Difference between revisions

=={{header|EchoLisp}}==

We implement infinite formal power series (FPS) using '''streams'''. No operator overloading in EchoLisp, so we provide the operators '''s-add, s-mul''' ,.. which implement the needed operations. '''poly->stream''' converts a finite polynomial into an infinite FPS, and '''s-value''' gives the value of a FPS at x.

<lang lisp>

(require 'math)

;; converts a finite polynomial (a_0 a_1 .. a_n) to an infinite serie (a_0 ..a_n 0 0 0 ...)

(define (poly->stream list)

(make-stream (lambda(n) (cons (if (< n (length list)) (list-ref list n) 0) (1+ n))) 0))

;; c = a + b , c_n = a_n + b_n

(define (s-add a b)

(make-stream (lambda (n) (cons (+ (stream-ref a n) (stream-ref b n)) (1+ n))) 0))

;; c = a * b , c_n = ∑ (0 ..n) a_i * b_n-i

(define (s-mul-coeff n a b) (sigma (lambda(i) (* (stream-ref a i)(stream-ref b (- n i)))) 0 n))

(define (s-mul a b)

(make-stream (lambda(n) (cons (s-mul-coeff n a b) (1+ n))) 0))

;; b = 1/a ; b_0 = 1/a_0, b_n = - ∑ (1..n) a_i * b_n-i / a_0

(define (s-inv-coeff n a b)

(if (zero? n) (/ (stream-ref a 0))

(- (/ (sigma (lambda(i) (* (stream-ref a i)(stream-ref b (- n i)))) 1 n)

(stream-ref a 0)))))

;; note the self keyword which refers to b = (s-inv a)

(define (s-inv a)

(make-stream (lambda(n) (cons (s-inv-coeff n a self ) (1+ n))) 0))

;; b = (s-k-add k a) = k + a_0, a_1, a_2, ...

(define (s-k-add k a)

(make-stream (lambda(n) (cons

(if(zero? n) (+ k (stream-ref a 0)) (stream-ref a n)) (1+ n))) 0))

;; b = (s-neg a) = -a_0,-a_1, ....

(define (s-neg a)

(make-stream (lambda(n) (cons (- (stream-ref a n)) (1+ n))) 0))

;; b = (s-int a) = ∫ a ; b_0 = 0 by convention, b_n = a_n-1/n

(define (s-int a)

(make-stream (lambda(n) (cons (if (zero? n) 0 (/ (stream-ref a (1- n)) n)) (1+ n))) 0))

;; value of power serie at x, n terms

(define (s-value a x (n 20))

(poly x (take a n)))

;; stream-cons allows mutual delayed references

;; sin = ∫ cos

(define sin-x (stream-cons 0 (stream-rest (s-int cos-x))))

;; cos = 1 - ∫ sin

(define cos-x (stream-cons 1 (stream-rest (s-k-add 1 (s-neg (s-int sin-x))))))

</lang>

{{out}}

<lang lisp>

(take cos-x 16)

→ (1 0 -1/2 0 1/24 0 -1/720 0 1/40320 0 -1/3628800 0 1/479001600 0 -1.1470745597729725e-11 0)

(take sin-x 16)

→ (0 1 0 -1/6 0 1/120 0 -1/5040 0 1/362880 0 -1/39916800 0 1.6059043836821613e-10 0 -7.647163731819816e-13)

;; compute (cos PI)

(s-value cos-x PI)

→ -1.0000000035290808

;; check that 1 / (1 - x) = 1 + x + x^1 + x^2 + ...

(define fps-1 (poly->stream '( 1 -1)))

(take fps-1 13)

→ (1 -1 0 0 0 0 0 0 0 0 0 0 0)

(define inv-fps-1 (s-inv fps-1))

(take inv-fps-1 13)

→ (1 1 1 1 1 1 1 1 1 1 1 1 1)

(s-value inv-fps-1 0.5) ;; check that 1 / (1 - 0.5) = 2

→ 1.9999980926513672

(s-value inv-fps-1 0.5 100) ;; 100 terms

→ 2

</lang>