1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
|
;;; "limit.scm" Scheme Implementation of one-side limit algorithm.
;Copyright 2005 Aubrey Jaffer
;
;Permission to copy this software, to modify it, to redistribute it,
;to distribute modified versions, and to use it for any purpose is
;granted, subject to the following restrictions and understandings.
;
;1. Any copy made of this software must include this copyright notice
;in full.
;
;2. I have made no warranty or representation that the operation of
;this software will be error-free, and I am under no obligation to
;provide any services, by way of maintenance, update, or otherwise.
;
;3. In conjunction with products arising from the use of this
;material, there shall be no use of my name in any advertising,
;promotional, or sales literature without prior written consent in
;each case.
;;@code{(require 'limit)}
(define (inv-root f1 f2 f3 prec)
(define f1^2 (* f1 f1))
(define f2^2 (* f2 f2))
(define f3^2 (expt f3 2))
(require 'root) ; SLIB
(newton:find-root (lambda (f0)
(+ (- (* (expt f0 2) f1))
(* f0 f1^2)
(* (- (* 2 (expt f0 2)) (* 3 f1^2)) f2)
(* (+ (- (* 2 f0)) (* 3 f1)) f2^2)
(* (- (+ (- (expt f0 2)) (* 2 f1^2)) f2^2)
f3)
(* (+ (- f0 (* 2 f1)) f2) f3^2)))
(lambda (f0)
(+ (- (+ (* -2 f0 f1) f1^2 (* 4 f0 f2))
(* 2 f2^2)
(* 2 f0 f3))
f3^2))
f1
prec))
(define (invintp f1 f2 f3)
(define f1^2 (* f1 f1))
(define f2^2 (* f2 f2))
(define f3^2 (expt f3 2))
(let ((c (+ (* -3 f1^2 f2)
(* 3 f1 f2^2)
(* (- (* 2 f1^2) f2^2) f3)
(* (- f2 (* 2 f1)) f3^2)))
(b (+ (- f1^2 (* 2 f2^2)) f3^2))
(a (- (* 2 f2) f1 f3)))
(define disc (- (* b b) (* 4 a c)))
;;(printf "discriminant: %g\n" disc)
(if (negative? (real-part disc))
(/ b -2 a)
(let ((sqrt-disc (sqrt disc)))
(define root+ (/ (- sqrt-disc b) 2 a))
(define root- (/ (+ sqrt-disc b) -2 a))
(if (< (magnitude (- root+ f1)) (magnitude (- root- f1)))
root+
root-)))))
(define (extrapolate-0 fs)
(define n (length fs))
(define (choose n k)
(do ((kdx 1 (+ 1 kdx))
(prd 1 (/ (* (- n kdx -1) prd) kdx)))
((> kdx k) prd)))
(do ((k 1 (+ 1 k))
(lst fs (cdr lst))
(L 0 (+ (* -1 (expt -1 k) (choose n k) (car lst)) L)))
((null? lst) L)))
(define (sequence->limit proc sequence)
(define lval (proc (car sequence)))
(if (finite? lval)
(let ((val (proc (cadr sequence))))
(define h_n*nsamps (* (length sequence) (magnitude (- val lval))))
(if (finite? val)
(let loop ((sequence (cddr sequence))
(fxs (list val lval))
(trend #f)
(ldelta (- val lval))
(jdx (+ -1 (length sequence))))
(cond ((null? sequence)
(case trend
((diverging) (and (real? val) (/ ldelta 0.0)))
((bounded) (invintp val lval (caddr fxs)))
(else (cond ((zero? ldelta) val)
((not (real? val)) #f)
(else (extrapolate-0 fxs))))))
(else
(set! lval val)
(set! val (proc (car sequence)))
;;(printf "f(%12g)=%12g; delta=%12g hyp=%12g j=%3d %s\n" (car sequence) val (- val lval) (/ h_n*nsamps jdx) jdx (or trend ""))
(if (finite? val)
(let ((delta (- val lval)))
(define h_j (/ h_n*nsamps jdx))
(cond ((case trend
((converging) (<= (magnitude delta) h_j))
((bounded) (<= (magnitude ldelta) (magnitude delta)))
((diverging) (>= (magnitude delta) h_j))
(else #f))
(loop (cdr sequence) (cons val fxs) trend delta (+ -1 jdx)))
(trend #f)
(else
(loop (cdr sequence) (cons val fxs)
(cond ((> (magnitude delta) h_j) 'diverging)
((< (magnitude ldelta) (magnitude delta)) 'bounded)
(else 'converging))
delta (+ -1 jdx)))))
(and (eq? trend 'diverging) val)))))
(and (real? val) val)))
(and (real? lval) lval)))
(define (limit proc x1 x2 . k)
(set! k (if (null? k) 8 (car k)))
(cond ((not (finite? x2)) (slib:error 'limit 'infinite 'x2 x2))
((not (finite? x1))
(or (positive? (* x1 x2)) (slib:error 'limit 'start 'mismatch x1 x2))
(limit (lambda (x) (proc (/ x))) 0.0 (/ x2) k))
((= x1 (+ x1 x2)) (slib:error 'limit 'null 'range x1 (+ x1 x2)))
(else (let ((dec (/ x2 k)))
(do ((x (+ x1 x2 0.0) (- x dec))
(cnt (+ -1 k) (+ -1 cnt))
(lst '() (cons x lst)))
((negative? cnt)
(sequence->limit proc (reverse lst))))))))
|
